Modelling solar cycle length based on Poincaré maps for Lorenz-type equations

Two systems of Lorenz-type equations modelling solar magnetic activity are studied: Firstly a low order dynamic system in which the toroidal and poloidal fields are represented by xand y-coordinates respectively, and the hydrodynamical information is given by the z coordinate. Secondly a complex generalization of the three ordinary differential equations studied by Lorenz. By studying the Poincaré map we give numerical evidence that the flow has an attractor with fractal structure. The period is defined as the time needed for a point on a hyperplane to return to the hyperplane again. The periods are distributed in an interval. For large values of the Dynamo number there is a long tail toward long periods and other interesting comet-like features. These general relations found for periods can further be physically interpreted with improved helioseismic estimates of the parameters used by the dynamical systems. Solar Dynamic Observatory is expected to offer such improved measurements.


Introduction
Understanding and being able to predict the length of a solar cycle is of great importance.Not only per se, but also due to its relation to the solar magnetic activity.
Many indicators of solar magnetic activity may be used to describe the solar cycle length.In Lundstedt et al. (2005) and Fligge et al. (1999) wavelet studies of the sunspot number were carried out.The length can be estimated from the wavelet coefficient maximum (WCM).A length of about Correspondence to: H. Lundstedt (henrik@lund.irf.se)4300 days corresponds to an 11-year period.In Fig. 1 the length roughly shows an inverse relation to the activity.However sometimes a more complicated relation (Lundstedt et al., 2005) also appears.
The sunspot number is a coarse indicator of the variation of the toroidal magnetic field of the sun.During the Maunder minimum (1645-1715) very few sunspot appeared and no cyclicity is visible.The cycle length was also increasing just before the start of the MM (Frick et al., 2001).However, using 14 C as and 10 Be indicator, cyclicity and periods are found during the MM (Beer et al., 1998;Lundstedt et al., 2006;Knudsen et al., 2009).
What can solar dynamo models tell us about the solar activity cycle length, amplitude and processes behind them?Briefly, a modulation of the 11-year cycle dynamo and cycle length can be achieved by changing the three parameters, ω,α and the meridional circulation rate of the dynamo (Dikpati and Gilman, 2001).Dikpati and Charbonneau (1999) emphasize that the velocity of the meridional flow is a critical factor in determining the period of the dynamo cycle.A least-squares fit on their numerical data gave the following scaling law governing the dependence of the dynamo period on model parameters: where the time period T is measured in years and u 0 , s 0 and η T , i.e. the meridional flow speed, the source coefficient (strength of source term representing the surface generation of poloidal field due to the decay of tilted bipolar active regions) and the turbulent diffusivity, are all measured in cgs units.When they introduced a random variation of the meridional velocity, they found that a long cycle followed a short one, but the long-term average was fixed at 11 years by the long-term average velocity of about 17 m/s.
The strength of the polar field, at minimum, has been used as an indicator of the strength of the upcoming cycle (Svalgaard et al., 2005).It is based on the idea that the polar field Published by Copernicus Publications on behalf of the European Geosciences Union.(as proxy of the poloidal field) is submerged, and strengthened by the differential rotation into a toroidal magnetic field below the surface.When the toroidal field is strong enough it starts to move upward and finally emerges through the solar surface as a bipolar region and sunspots appear.Svalgaard et al. (2005) compared the strengths of the polar fields at sunspot minima and estimated the next cycle (Cycle 24) to be weak, as weak as around 1900.This estimate is in accordance to what we found in wavelet ampligram studies Lundstedt (2006b).
The most recent (2008) forecast by the Cycle 24 Prediction Panel, asked for by NASA and organized by NOAA, is that the sunspot maximum will reach 90 in May 2013.That will be the weakest since 1928.A comprehensive description of all the predictions on which the forecast is based is given in Pesnell (2008).In this article we try to evaluate what can be stated in general about a solar cycle length/period and about the next.We therefore do not wish to restrict our study by using limited observations of indicators of the solar magnetic activity, such as the sunspot number and to use statistical methods.Instead we start with two different physics-based dynamical systems (Tobias et al., 1995;Weiss et al., 1984) modelling the solar magnetic activity.On the other hand these two systems have also been shown to describe many of the features in the sunspot variation.These two systems are then studied mathematically, using topological methods, to examine what can be said in general about a period and the following period.

Topology
Topology is formally defined in terms of set operations.Let X be a set.A topology on X is a collection T of subsets of X, called the open sets.The set X together with a topology T is called a topological space.A subset A is dense in a space X if and only if A intersects every nonempty open set in X.

Dynamical systems
The basic goal of the theory of dynamical systems (Katok and Hasselblatt, 2006) is to understand the eventual or asymptotic behaviour of an interative process.If the process is a differential equation whose independent variable is time, then the theory attempts to predict the ultimate behaviour of solutions of the equation in either the distant future (t → +∞) or distant past (t → −∞).
A dynamical system is said to be chaotic (Devaney chaotic) (Devaney, 2003) if there exists at least one dense orbit and the set of periodic orbits is dense.Topological invariants of periodic orbits, such as the linking number (Lundstedt, 2009), can be used to identify the strange attractor and the stretching and squeezing mechanisms (Gilmore and Lefranc, 2002).

Solar cycles and dynamical systems
The regeneration of the solar cycle is given by, where B T is the toroidal, B P the poloidal field, the rotation and α is a pseudo-tensor requiring a turbulent model.This is by the so called ω and the α effects (Ossendrijver, 2003).
In α mean field models α has been suggested to be estimated either only by the kinetic helicity or by, where a = b √ 4πρ is the Alfvén speed based on the smallscale magnetic component and τ c the correlation time of the turbulent motion (Charbonneau, 2005).Tobias et al. (1995) describe the modulation and occurrence, of the solar cycle by using a low order dynamic system of Lorenz equation type (Eqs.6a, b, c).The toroidal (B T ) and poloidal fields (B P ) are represented by the x-and the ycoordinates, respectively.The hydrodynamical information is given by the z-coordinate.They suggest that as the parameters are varied, a single fixed point loses stability in a Hopf bifurcation to a periodic orbit, which then undergoes a secondary Hopf bifurcation creating a two-torus, and with the torus in turn breaking down to give chaotic motion.

Low dimension Lorenz equation
We solved this Lorenz system of equations for different values of the parameters (a,c,d,λ,ω,µ, )  10 in Tobias et al. (1995).λ represents the dynamo state and µ the hydrodynamic state.represents the rotation rate and is related to λ and µ as The ω and the α effects are not explicitly given, but introduced through .
The system of differential equations is given by where a, c, d, ω, λ and µ are parameters.

Existence of a strange attractor
Let M µ = {(x,y,z) : z = 0 and x 2 + y 2 > µ} be a plane with a hole of radius √ µ.The set M µ consists of the set of points (x,y,z) in the plane z = 0 such that ż < 0 according to Eq. (6c).Hence, the flow defined by Eq. ( 6) intersects the manifold M µ transversely.We may therefore consider the Poincaré map P : M µ → M µ of the flow with respect to the manifold M µ .Figure 4 shows a picture of the set M µ .Any point (x,y,0) on M µ flows downward and if it returns to M µ it flows up through the hole and then hits M µ from above in the point P (x,y,0), as indicated at the picture.
By studying the Poincaré map we will give numerical evidence that the flow (6) has an attractor with fractal structure.We will concentrate on three sets of parameters, but emphasise that our method can be used for any set of parameters.The three sets of parameters are shown in Table 1.
If U ⊂ M µ is a non empty and open set (open as a subset of M µ ), such that P (U ) ⊂ U , then there exists an attracting  closed set ⊂ U defined by This means that is the smallest closed set, such that if p is a point in U then d(P n (p), ) → 0 as n → ∞, where d denotes the distance.In other words, is the smallest set such that the orbit of all points in U is attracted to .Thus, in order to prove that there exists an attractor of the flow (6) it is sufficient to find a set U ⊂ M µ such that P (U ) ⊂ U .We will do this as follows.
Using the method in Tucker (2002), we use the classical Runge-Kutta method to integrate the flow of points on the circles C r 1 and C r 2 , where until they return to M µ .Note that C r 1 and C r 2 are the boundary of the set A(r 1 ,r 2 ).We thus calculate P (C r 1 ) and P (C r 2 ).
If P (C r 1 ) ⊂ A(r 1 ,r 2 ) and P (C r 2 ) ⊂ A(r 1 ,r 2 ), it follows by the continuity and invertibility of P that P (A(r 1 ,r 2 )) ⊂ A(r 1 ,r 2 ).This means that any point in A(r 1 ,r 2 ) returns to A(r 1 ,r 2 ) as indicated in Fig. 6.Since the numerical integration is associated with some errors, it is preferable to have numerical evidence that P (C r 1 ),P (C r 2 ) ⊂ A(r 1 ,r 2 ), with r 1 < r 1 < r 2 < r 2 and that the differences r 1 − r 1 and r 2 − r 2 are not too small.
After a few tries one easily finds values of r 1 and r 2 for each set of parameters, such that the numerical approximation P maps a set of 1024 equally spaced points on C r 1 and C r 2 into A(r 1 ,r 2 ).Table 2-Table 4 show, for different step sizes, the largest and the smallest possible values or r 1 and r 2 compatible with the calculated returns of 1024 points each on C r 1 and C r 2 .
Since in all three sets of parameters, we have good nummerical evidence that the Poincaré map P maps A(r 1 ,r 2 ) into itself, we feel confident that there is an attractor inside A(r 1 ,r 2 ) defined by Let us now investigate some of the properties of this attractor.Let us first observe that because of Eq. ( 8), the set P n (A(r 1 ,r 2 )) is an approximation of the attractor if n is sufficiently large.We will therefore visualise the set P n (A(r 1 ,r 2 )) for a different n.

Prediction of period length
Even though all points in A(r 1 ,r 2 ) return to A(r 1 ,r 2 ), different points may need different amounts of time to return.We will refer to this time as the period.If we have observed that a point has a certain period, what be said about the period of the next return.We pursue this problem in the following way.For any point (x,y,0) in the attractor we want to calculate the time needed for this point to return to P (x,y,0) and then calculate the time needed for P (x,y,0) to return to P 2 (x,y,0).We then want to plot the second period as a function of the first period in a diagram to see if there are any correlations.
We cannot know exactly what points are in , nor can we calculate the period for all infinitely many points in .We therefore take the finitely many points in P n (A(r 1 ,r 2 )) which we have found and are shown in Fig. 8-10, and we calculate the periods for these points.This was done, and the result is shown in Fig. 11-13.As is evident from these figures, there are some correlations between a period and the following period.This means that if one observes a period of some certain length, it is possible to make non-trivial predictions of the length of the following period, using only the observed period.For instance, if we observe a period of length 2.5 in the case of the third set of parameters, then we can predict that the next period will be of a length between 2.39 and 2.84.

Derivatives and Lyapunov exponents
Since the Poincaré map P maps points in the plane z = 0 onto themselves, the derivative d p P of P at a point p = (x,y,0) is a linear map that maps vectors parallel to the plane z = 0 onto vectors parallel to the plane z = 0. Hence the derivative can be represented as a 2 × 2 matrix.Similarly, if we consider any power P n of P , where n is a natural number, then we may consider the derivative d p (P n ) of P n at a point p.As with d p P we can represent d p (P n ) as a 2 × 2 matrix.One can show that for almost all points p ∈ A(r 1 ,r 2 ) the limit exists.The number λ + (p) is called the largest Lyapunov exponent of P at the point p.If λ + (p) is positiv, this means that in average, small neighbourhoods around p are stretched a factor ∼ e λ + (p)n in some direction under the action of P n .So if λ + (p) is positive, then there are arbitrary points close to p that move away from the orbit of p with exponential speed.A positive Lyapunov exponent is therefore related to a chaotic behaviour of the map P , and it is then impossible to make long time predictions since errors grow with exponential speed.We estimated λ + (p) for all points p that is in P n (A(r 1 ,r 2 )), which are our approximation of the attractor .Using the methods from Tucker (2002), the derivative of P k was numerically calculated for points in P n (A(r 1 ,r 2 )).By Eq. ( 9) we then obtain an approximation of the Lyapunov exponent by Table 5 shows the obtained approximations of the Lyapunov exponents.

Complex generalization of the Lorenz equation
Another equation system, describing the solar magnetic state, is considered: Weiss et al. (1984) developed a complex generalization of the Lorenz equation: The dynamo number D is given by:  Here α is explicitly mentioned.Kitiashvili and Kosovichev (2008) also include magnetic helicity and herewith α m .The dynamo number measures the strength of the two induction effects α 0 and 0 , relative to the diffusivity.We studied the Poincaré map of the flow returning to the hyperplane Re(A) = 0, for the parameters D = 2 and ν = 0.4, using the Runge-Kutta method.Since this flow is in 6 real dimensions (3 complex dimensions), the analysis of this Poincaré map is much more numerically involved than that of the system (6).Starting with a point (A,B, ) we calculated a big number of returns to the hyperplane Re(A) = 0.The returns were situated in a 5-dimensional box of size 4.1 × 3.4 × 4.0 × 3.0 × 6.0.This box was divided into a grid of n 5 boxes.The aim was to find a subset of these boxes such that any point in any of the boxes returns to some of the boxes.If such a subset is found, then we let U be the interior of the union of the boxes.Then, as in the previous case, there exists an attractor satifying Eq. ( 7).However we cannot numerically check the return of any point in the boxes, so instead we checked the returns of m 5 uniformly distributed points in each of the boxes.Some of these points did not return to the boxes, so instead we added new boxes containing these returns.For these new boxes, we again checked the returns, and this process was repeated several times.This process eventually stopped, and no more boxes were needed.Hence this is a numerical evidence for the existance of an attractor.These calculations were carried out for m = 4 and n = 256, and smaller values of these parameters.Each case resulted in that a trapping region was found.
To estimate the dimension of the attractor, we estimated the attractor with a big number of returns of one point and checked how many boxes from the above mentioned grid that were needed to cover the returns.This was done with different numbers of returns, and different grid sizes.
If the grid consists of n 5 boxes, then it is expected that asymptotically n s boxes are needed to cover the attractor.If we let N denote the number of boxes in the cover, by plotting logN against logn and fitting a straight line to the data, the slope of this line is an estimate of the dimension of the attractor.Table 6 shows the results from an estimate of the attractor using 10 8 points.There is a plot of this in Fig. 14.The slope of the line suggests a dimension of the attractor (as a subset of the hyperplane) of about 1.7.Hence this suggests that the dimension of the attractor of the flow is about 2.7.
We again studied the relationship between the period n and period n+1.In both cases ν < 0.4.In Fig. 15, D was chosen to be 2.0.Only two values for the following periods are expected.When D is sufficiently large and ν < 1, then the periodic solution becomes unstable, first by multiply periodic and finally by chaotic solutions.The periods are distributed in intervals A comet-like distribution of points appears for D = 5.0 and higher.Most of the points (period, following) appear in the head and are therefore similar.Figure 16 shows when D was set to 5.0.
When D is increased shorter periods exist and also interesting, pronounced comet-like tail structures.Figure 17 shows when D was set to 9.0.These tails illustrate how long periods are followed by short periods and vice versa.

Conclusions and future plans
The fundamental question we asked was, "What can be said in general about a solar cycle/period and the following one?"We used Poincaré maps to explore two dynamical systems: A low-dimension Lorenz system, constructed by Tobias et al. (1995), and a six-dimensional complex generalization (Weiss et al., 1984).The Poincaré maps gave numerical evidence that the flow has an attractor with fractal structure.The period was defined as the time needed for a point on a hyperplane to return the hyperplane again.The periods were distributed in an interval.For large values of the Dynamo number there is a long tail toward long periods and also other interesting features.There is a tendency the long periods are followed by short periods, and short periods are followed by long periods.
These results can be compared with what Dikpati and Charbonneau (1999) found.They emphasize that the velocity of the meridional flow is a critical factor in determining the period of the dynamo cycle.When introducing a random  variation of the meridional velocity, they also found that a long cycle followed a short.
In Fig. 18 the period of a sunspot cycle vs. the following is plotted.There might be similar discerned structures; long periods can be followed by a short or an average and a short can be followed by a long.However, the number of estimated solar cycle periods is very small.Some of the estimates are also uncertain.The estimates were produced by National Geophysical Data Center (NGDC) in Boulder, USA.
We expect upcoming helioseismic observations with Solar Dynamics Observatory (SDO) (Hoeksema and the HMI magnetic Team, 2008) will give us improved physical values for the dynamical systems.More advanced equation systems, including terms describing helicity, will also be explored.
The found interesting structures will herewith be further studied and neural networks will be used to learn these patterns for predictions (Lundstedt, 2005(Lundstedt, , 2008)).

Fig. 1 .
Fig. 1.Upper panel time of wavelet coefficient maximum representing the length vs. time in year.Lower panel: Wavelet coefficient maximum representing solar activity vs. time.

Fig. 2 .
Fig. 2. Upper panel shows sunspot number.Lower panel shows 14 C. MM stands for Maunder Minimum and DM for Dalton Minimum.

Fig. 3 .
Fig. 3. Trajectories of a dynamical system that are suggested to mimic present solar sunspot cycles (Third set of parameters).

Fig. 4 .
Fig. 4. The action of the map P is that it flows the point (x,y,0) until it hits M µ at the point P (x,y,0).

Fig. 5 .
Fig. 5. Picture of the set M µ and the set A(r 1 ,r 2 ).

Fig. 11 .
Fig. 11.Periods of the first parameter set.

1002H.
Lundstedt and T. Persson: Solar magnetic cycles and topology
, Table 1.The cases correspond to increasing values and Figs. 5, 7, and

Table 1 .
Three sets of parameters.

Table 2 .
Values of r 1 and r 2 for different step sizes for the first parameter set.

Table 3 .
Values of r 1 and r 2 for different step sizes for the second parameter set.

Table 4 .
Values of r 1 and r 2 for different step sizes for the third parameter set.

Table 5 .
Approximations of the largest Lyapunov exponent for the three sets of parameters.

Table 6 .
Estimating the dimension of the attractor of Eq. (11).