Nonlinear solar cycle forecasting: theory and perspectives

In this paper we develop a modern approach to solar cycle forecasting, based on the mathematical theory of nonlinear dynamics. We start from the design of a static curve fitting model for the experimental yearly sunspot number series, over a time scale of 306 years, starting from year 1700 and we establish a least-squares optimal pulse shape of a solar cycle. The cycle-to-cycle evolution of the parameters of the cycle shape displays different patterns, such as a Gleissberg cycle and a strong anomaly in the cycle evolution during the Dalton minimum. In a second step, we extract a chaotic mapping for the successive values of one of the key model parameters – the rate of the exponential growthdecrease of the solar activity during the n-th cycle. We examine piece-wise linear techniques for the approximation of the derived mapping and we provide its probabilistic analysis: calculation of the invariant distribution and autocorrelation function. We find analytical relationships for the sunspot maxima and minima, as well as their occurrence times, as functions of chaotic values of the above parameter. Based on a Lyapunov spectrum analysis of the embedded mapping, we finally establish a horizon of predictability for the method, which allows us to give the most probable forecasting of the upcoming solar cycle 24, with an expected peak height of 93±21 occurring in 2011/2012.


Introduction
The sunspot number time series provides the longest existing record of solar activity, and is thus the best available data set for studying the long-term evolution of solar activity and, in particular, of the 11-year activity cycle.Here, we apply the mathematical concepts of the nonlinear dynamical systems, in order to derive dynamical properties of the solar Correspondence to: A. L. Baranovski (alexander.baranovski@gmail.com)cycle and to derive short-and long-term predictions of the sunspot number.
Many attempts have already been made in the past to predict future solar cycles.They can be grouped into two main categories (Hathaway et al., 1999): -Regression Techniques: standard auto-regression, curve-fitting and neural networks.
-Precursor Techniques: a combination of sunspot number indicators (tracers of the closed solar magnetic field) and geomagnetic indicators (proxies of the Sun's open dipolar magnetic field).
Regression techniques provide short-term extrapolations using observed indices of solar activity from the recent past to extrapolate into the near future.They suffer from a rather low performance during the transition from one cycle to the other, while precursor techniques can deliver an estimate of the amplitude for the next solar cycle one or two years before the previous cycle has ended.However, both approaches lack a physical basis.Next to these, magnetohydrodynamical models have been developed during the last decades to simulate the solar dynamo processes that generate the solar magnetic activity, but even now, they cannot reproduce the detailed evolution of the solar cycle (for a recent review of this effort, see Charbonneau, 2005).Therefore, a mathematical approach of the problem is still useful and productive, in order to progress in this field of research.However, so far, most of the published methods are either based on crude statistical techniques and purely empirical relations or they use classical mathematical tools adapted only to deterministic time series (e.g.Fourier analysis) or adopt a linear stochastic approach.The irregularity of the sunspot time series, both in period and amplitude, as well as the existence of grand minima, where the cyclicity vanishes, suggests a more complex evolution.Recently, this was confirmed by advanced models Those original numbers, derived from numerous visual observations, show not only the 11-year pulsation but also its irregular variations in amplitude and periodicity.
of the solar dynamo that include a meridional circulation of the weak surface fields (Charbonneau, 2005;Dikpati et al., 2006).Although they can only qualitatively reproduce the solar evolution, such physical models of a nonlinear dynamo show that the system is highly sensitive to the initial observational constraints and takes on a chaotic or intermittent behaviour.This is why we apply here two methods adapted to the analysis of such nonlinear dynamical systems with a limited number of degrees of freedom.The S plane surface and Poincare section techniques, as well as a nonlinear mapping analysis applied to a chaotic time series, can bring qualitative and quantitative insights into the underlying physical dynamics.Contrary to earlier studies, we make use here of the entire sunspot time series, instead of one or just a few past cycles.

Curve-fitting modelling
Figure 1 shows the annual sunspot number (SN) time series s(t) for the years 1700 through 2005.From the above SN time evolution, one can easily observe the cyclic behavior and the cycle-to-cycle variations in amplitude, shape and duration.This led Waldmeier (1935) to first propose the idea of representing the course of SN x(t) during each cycle, from minimum to minimum, by a standard curve belonging to a family of pre-defined standard profiles.Stewart and Panofsky (1938) proposed a more advanced function for representing the shape of a cycle, defined as: where t 0 is the starting time of a sunspot cycle and F , m, α are parameters that vary from cycle to cycle.In this representation, the ascending and descending phases of the cycle are described by a power and exponential law, respectively.By contrast, we first modify the pulse shape, as defined in Eq. ( 1), and detect the nonlinear dynamics of the leastsquares optimal shape parameters in their cycle-to-cycle evolution.This will allow us to achieve an optimal performance in a retroactive "forecast" of past cycles and to obtain a most probable prediction of the upcoming cycles.
We will assume that the time variations of SN can be described as a train of pulses using the following signal model: For any k-th pulse corresponding to sunspot cycle x k (t) occurring in the time interval τ k−1 ≤t<τ k between two successive SN minima, the model is characterized by 1. the deterministic pulse shape (Fig. 2) , 2. the pulse amplitude A k ; 3. the pulse duration T k =τ k −τ k−1 , where τ k is the starting time of the (k+1)-th pulse and All parameters A k , T k , ε k , as well as the damping factor α k vary from cycle to cycle.β=11 4 is a scale factor.
Now we are able to extract the values of the model parameters for each k-th pulse, thus obtaining a parametric description of the shape of each individual solar cycle: 1. sunspot number at the beginning of the cycle 2. time of extremum (maximum value of the SN) 3. maximum sunspot number (amplitude of the cycle) 4. sunspot number at the end of the cycle We can also formulate the estimation of the model parameters {A k , ε k , α k } by a least-squares fitting x(t) to the sunspot number time series s(t).Here, we assume that the lengths T k and starting times τ k−1 of the sunspot cycles are given a priori.A conventional approach for solving this nonlinear optimization problem is based on numerical methods, e.g. the maximum neighborhood algorithm of Marquardt.In our case, we are able to demonstrate an analytical solution of the following distance function minimization with respect to the 28 triples of the parameters {A k , ε k , α k } of k-th spot cycle, where k=−4, −3, . . ., 23 are given according to the Zürich cycle numbering system.Rewriting Eq. ( 8) in the form where and z k (t) = t−τ k−1 T k , we easily find a first equation which must be solved with respect to the amplitude A k .Introducing the new notations and also taking into account that Eq. ( 10) is then given by A second equation: after substitution of Eq. ( 11) leads to the solution where Finally, a third equation with parameters A k and ε k , respectively, expressed as Eqs.( 11) and ( 13), leads to a nonlinear equation for the last unknown parameter α k : The roots α * k , k=−4, −3, . . ., 23 can be found with high accuracy by the simplest numerical methods, e.g. the Newton method.Thus having α * k , we first obtain for any year t in the k-th sunspot cycle, i.e. τ k−1 ≤t<τ k .We then calculate the parameters A * k , ε * k according to Eqs. ( 7) and ( 9).By this way, we collect all 28 triples {A * k , ε * k , α k } and plot the model signal as shown in Fig. 4, where it is compared with the original sunspot number observations s(t).
The performance of the pulse-shape fitting model can be evaluated by a statistical analysis of residuals ξ (t)=s(t)−x(t).Figure 5 shows this ξ (t) function for the years 1700 to 2005.Its elementary statistical characteristics are: Table 1 lists an estimation of the probability density function (PDF) and of the autocorrelation function (ACF) of ξ (t).
The Goodness-of-Fit test detects a Gaussian distribution ξ (t).It is important to note that ξ (t) does not behave as white noise but rather as correlated or "coloured" noise, with an autocorrelation function of the form: The power spectral density of a noise with such an ACF can be calculated by using the Wiener-Chinchine theorem.This leads to: where ω 0 =π/6.This power distribution of the noise in the SN data series has thus a Lorentzian profile.

Nonlinear mapping technique
Now we rewrite the signal model for x(t) as: where First, from Eq. ( 5), we note that the length of the ascending phase T + k relates to the new parameter γ k by a simple law: Similarly, the µ-amplitude dependency versus γ , shown in Fig. 6, with points µ * k , γ * k , k=−4, −3, ..., 23, can be well approximated by the least-squares fitting function: From the least-squares optimal model Eq. ( 16) with a set of 28 triples {A * k , ε * k , α * k }, we also calculate the offset values η * k from Eq. ( 21).Finally, the neighboring points from a set γ * k , η * k , k=−4, ..., 23 are linked together by a broken line, as shown in Fig. 7.We thus finally obtain a piece-wise linear function which generates the approximated values of η k for arbitrary input values of the parameter γ in the interval A≡0.642≤γ n ≤B≡1.235.
By analogy, we can design a generator for the descending phases T − k (Fig. 8) for any arbitrary number k: based on the set γ * k , T − k , k=−4, ..., 23 , where the lengths T − k , k=−4, ..., 23 are easily obtained from the historical data series s(t).It leads to the following formula for the estimate of the length of any cycle k: A peculiarity was precisely detected by Sonett (1983) and Wilson (1988) in the solar cycle evolution at the beginning of the Dalton minimum.They explained it by a possible error or misplacement of SN minima in the Wolf data, respectively.Vitinsky, Kopecky and Kuklin (1986) and Serre and Nesme-Ribes (2000) link this anomaly of SN behavior to a phase catastrophe over the period 1790-1798 when the solar cycle evolution was not cyclic but roughly linear, with a greatly reduced phase evolution rate along the trajectory (Usoskin and Mursula, 2003).We also detect this phenomenon and in addition, we extend the time limits of its manifestation to the period {1775-1823}.

Nonlinear mapping for parameter γ and its properties
We now assume that γ * n is a one-dimensional chaotic process generated by a first-order map φ linking successive values γ * n−1 , γ * n .As all γ n are distributed within the interval A≡0.642≤γ n ≤B≡1.235(see Fig. 7), a map φ is used to transform (A,B) into itself.
The plot of all 27 pairs γ * n−1 , γ * n is given in Fig. 10.Points are numbered according to the Zürich cycle numbers, i.e. the pair γ * k−1 , γ * k corresponds to cycle number k. Linking couples of neighbouring points by a line, we obtain a piece-wise linear approximation of a map φ, as depicted in Fig. 11.
A map where coefficients a j , b j , j =1, 26 and a partition is positive.It indicates an average rate of the exponential divergence of the initially closed trajectories-sequences {γ n }.
The piece-wise linear map φ obeys a piece-wise constant probability density function p ϕ (x) (Baranovski and Daems, 1995), which can be found from the Perron-Frobenius equation One can show that φ is an ergodic map and using Birkhoff's ergodic theorem Eq. ( 27), this leads to and the autocorrelation function is given by It is important to note that Eq. ( 29) is in good agreement with direct estimations of the Lyapunov exponent λ max obtained from the monthly sunspot number time series (Mundt et al. 1992): This confirms that the extracted mapping φ precisely reflects a key dynamic property: a rate of exponential divergence in the sunspot time series.An estimation of the probability density and the autocorrelation functions for the chaotic trajectory {γ n , n=1, 10 000} are collected in Table 2.We believe that the small asymmetry in the shape of the PDF reflects the incompleteness of the map φ, due to the limited existing data record.In other words, we expect that new parameters {γ n } derived from the future observations of the annual sunspot numbers will be not lying on the lines of the current map φ and will define new nodes in φ, implying a splitting of some lines in φ.We think that mainly lines with the slanting slopes dϕ dγ <1 may be split into pairs of expanding lines with slopes dϕ dγ >1.As we see from Fig. 11, there are only four short lines with low slanting slopes and we believe that the main structure of the φ mapping is already included.This observation thus prompts us to already use the current mapping for attempting a forecast of future solar activity.
Currently, the Sun is approaching its next activity minimum, and most other predictions suggest that in five or sixyears time, the next solar maximum will be weak.Still, Dikpati et al. ( 2006) used a flux-transport dynamo model and concluded instead that cycle 24 will have a 30-50% higher peak than cycle 23.David Hathaway and Robert Wilson (NASA's Marshall Space Flight Center) also predict a maximum of the smoothed sunspot number of about 145±30 in 2010, while the following cycle should have a maximum of about 70±30 in 2023.Thus we agree only with their prediction for cycle 25 (see Fig. 16).

Auto-regression algorithm for γ -parameter
As we have seen earlier, γ is a central parameter in our nonlinear mapping technique.We now continue to investigate this parameter by using an auto-regression methodology.We thus consider {γ k , k=−4, ..., 23} as an auto-regressive process of order p, i.e.
where the auto-regression coefficients {a i , i=1, ..., p} can be easily calculated by a least-squares technique.
We establish that every j-th cycle has a different optimal value p, which minimizes the difference between    Making such plots for every cycle, starting from the second cycle (in this case p=1), we obtain the following empirical distribution for the orders: shape of the k-th cycle, we are able to plot all seven scenarios of evolution of this cycle.Figure 15 shows this for the 8 last cycles in comparison with the actual historical record s(t) (red curve).
In order to predict the value γ 24 for the upcoming cycle, we calculate the average and compare it with γ nonlinear 24 =ϕ γ * 23 =0.8354 found above, as shown in Fig. 16.
Figure 17 illustrates all possible scenarios for the upcoming cycle.
By including the average value for γ 24 as a 29th node in a updated design of the nonlinear mapping φ, we modify a bit the φ mapping constructed above.Then we can use the  17.Modelled average sunspot evolution during cycle 24: the 7 black curves and their average (green curve) result from the autoregression method, while the red curve is obtained using the nonlinear mapping technique.Comparatively, the latter curve suggests a lower maximum of about 90, which also peaks about one year later.new mapping for the generation of the upcoming cycles, as shown in Fig. 18.

Statistical analysis of the modelled annual sunspot number
Now we can study the statistical properties of the future evolution of the annual sunspot number modeled by the dynamical system approach (mapping φ) based on the continuoustime pulse shape representation of Eq. (3) for a cycle.A statistical analysis can be done analytically in terms of the invariant density of the mapping φ (see appropriate techniques in Baranovski andSchwarz, 1999 and2003) and theory of the semi-markovian processes (Nollau, 1980).Here, we just show in Table 3 the resulting plots of the estimates for the PDF and ACF of the annual sunspot number in a very long For such long-term extrapolations, due to the sensitivity of the solution to initial conditions, only qualitative information about the dynamlical evolution is meaningful beyond a few solar cycles.In this respect, this extrapolated sunspot time series contains long-term modulations of the cycle amplitude, in good accordance with observations from the last two centuries.Note also that this multi-century extrapolation does not display systematic secular trends or periods of Grand Minima (intermittent interruptions of the solar cycle, with constantly low activity, as observed in past centuries) in the upcoming 3 centuries.
time interval of 10 000 years.We note that the periodicity of the ACF is equal to 11 years and the estimates of the PDF and ACF are close to the ones of the observed daily sunspot number over the past observed period 1700-2005.

Conclusions
In this work, we first applied a nonlinear, least-squares fitting (NLF) to the annual sunspot time series R i (t) (or x(t), as used in text for mathematical convenience) for the years 1700<t<2005, which allowed us to characterize each cycle by a small set of parameters.We then applied two methods to the NLF output parameters, where each value represents one of the 23 available solar cycles: 1.A nonlinear mapping technique (NMT) has been constructed based on the output parameters 2. A linear regression technique (LRT) has been applied on the same NLF parameters 3. Finally, we synthesized an hybrid method based on the combination of both (LRT+NMT).
This allowed us to derive a prediction of the upcoming solar cycle (#24).Each method produces a different value for the activity maximum: -NMT method: maximum of 93±21 in 2011.79±0.5,However, those values span a limited range and are all similar to the current cycle amplitude (110).Although there is some discrepancy between the results, the range is restricted enough to exclude either a very strong cycle, similar to cycle #19 that peaked in 1959 (Rmax=200), or a very weak cycle.Finally, we modelled the future evolution of the sunspot number over several millenia.Over such a long term, the synthesized times series can only reproduce qualitatively the actual evolution of the sunspot number, given the limited number of past cycles available to build the nonlinear mapping and also given the intrinsically high sensitivity of the system to initial conditions.However, it allowed us to derive some global statistics of the time series, which contains many more cycles than in the available past observations.Those preliminary results illustrate quite well the potential offered by those mathematical methods as diagnostic tools for the physical processes underlying the observed solar activity.

Fig. 1 .
Fig. 1.Historical yearly sunspot number time series 1700-2005.Those original numbers, derived from numerous visual observations, show not only the 11-year pulsation but also its irregular variations in amplitude and periodicity.

Fig. 12 .
Fig. 12. Cycle-to-cycle evolution of the parameter γ n , modeled by the nonlinear mapping φ.The horizontal time range spans the period 1700-2800.

Fig. 13 .
Fig. 13.Lamerey-Königs diagram showing how predictions of the next few solar cycles can be iteratively derived from the mapping of Fig. 11, starting from the last observed cycle (#23).
j −k and γ * j .As an example, for the current cycle 23, we see from Fig.14that orders 8 and 9 give a very close auto-regression value γ 23 for γ * 23 =0.9707.

Fig. 15 .
Fig. 15.The 8 last cycles #16-23 plotted for the 7 most probable autoregression values in comparison with the actual historical record s(t) (red curve).

Fig. 18 .
Fig.18.Evolution of the sunspot cycles 24-56 modelled by the modified nonlinear mapping.For such long-term extrapolations, due to the sensitivity of the solution to initial conditions, only qualitative information about the dynamlical evolution is meaningful beyond a few solar cycles.In this respect, this extrapolated sunspot time series contains long-term modulations of the cycle amplitude, in good accordance with observations from the last two centuries.Note also that this multi-century extrapolation does not display systematic secular trends or periods of Grand Minima (intermittent interruptions of the solar cycle, with constantly low activity, as observed in past centuries) in the upcoming 3 centuries.

Table 1 .
Histogram and normalized ACF of residuals ξ (t) which shows Gaussian colored noise properties.

Table 1 .
Histogram and normalized acf of residuals ξ(t) which shows Gaussian colored noise

Table 1 .
Histogram and normalized acf of residuals ξ(t) which shows Gaussian colored noise 1 Fig.6.Scatterplot of the (γ n , µ n ) pairs for available values n=−4, . . ., 23.The points display a distinctive trend which allows one to derive a least-squares relation between those two parameters.

Table 2 .
Estimates of the PDF and ACF of the parameter γ .

Table 2 .
Estimates of the pdf and acf of the parameter γ

Table 3 .
Histogram and ACF of the modelled annual sunspot number time series over 10 000 years.The statistical distribution and the 11-year periodicity are quite similar to the past observed but much shorter time series.

Table 3 .
Histogram and acf of the modelled annual sunspot number time series over 10000 years.The statistical distribution and the 11-year periodicity are quite similar to the past observed but much shorter time series.