© ESO, 2013

1. Introduction

Solar activity has been observed for centuries (if not millennia) through the occurrence of sunspots. Seen as relatively darker circular regions on the solar disk, they have been accurately measured since the invention of telescope. However, exactly how and why they form has only begun to be understood in the last century. Sunspots are now understood to be regions where ribbons of intense magnetic field have risen up, e.g. due to magnetic buoyancy and breached the surface of the sun, which then eventually expand out into the solar atmosphere. The mechanism through which these ribbons of intense magnetic field occur can be explained by using a dynamo theory in the mean field electrodynamics. In this theory, solar magnetic fields are viewed as the sum of toroidal and poloidal fields where the toroidal fields are the components wrapping around the sun in the torodial direction, while the poloidal fields are the fields orthogonal to these toroidal magnetic fields.

As such, a conventional theory for the origin and evolution of the solar dynamo is often based on the two processes of the conversion of toroidal fields to poloidal fields and vice versa. Given poloidal fields, prominent latitudinal and/or radial differential rotation of the sun provides a simple mechanism to generate toroidal fields from poloidal fields, and this is so-called Ω effect. However, the conversion of poloidal fields from the toroidal fields is far more intricate and is believed to require the special property of small-scale turbulent convective motion in the solar interior such as helicity of turbulent flow u, e.g., given by (see, e.g. Parker 1957; Biskamp 2003). Here τ_C is the correlation of time of turbulent flow u. This is the so-called α effect from the pioneering work of Steenbeck et al. (1966). The combination of α and Ω thus constitutes a close loop for the generation of magnetic fields, providing us an intuitive explanation of how a large scale magnetic field can be generated from differential rotation and (helical) turbulent flows. As a basic mechanism for the observed 22 yr solar period, this so-called α-Ω dynamo has been popular and widely used in modelling solar dynamo.

Besides the regular 22 yr period, over the centuries, solar activity, a measure of magnetic field strength on the solar surface, has however exhibited irregularity on many different time scales. Roughly, depending on the 10 year average sunspot numbers (SN), the sun has three distinct behaviours: minima activity (<15 SN), normal activity (15 < SN < 50), and maxima activity (SN > 50) with the probability of 0.17, 0.75 and 0.08, respectively (see Usoskin 2007). Furthermore, the solar magnetism is known to enter periods of high or low activity apparently at random, with the typical example of Mauder Minimum. In order to explain these different time scales of solar magnetic activity, various attempts have been made to modify the conventional α-Ω dynamo model, for instance, by including nonlinearity (e.g. mode-coupling, saturation, etc.). In particular, a pioneering work by Weiss et al. (1984) proposed a simplfied nonlinear dynamic system where magnetic fields are shown to exhibit complex temporal behavior, quantitatively capturing observed irregularity in solar activity.

The effect of fluctuating α effect on solar activity has also been studied by previous authors (e.g. Hoyng et al. 1993). In fact, given that the α can arise from the total helicity, which is obtained as spatial, or temporal average, it is not unreasonable that alpha can fluctuate in time and/or space, even if slowly, over the solar cycle. The investigation of fluctuating alpha has also been stimulated by various recent works that seriously question the validity of the traditional α-Ω dynamo model because of significant quenching in alpha effect even by a very weak magnetic field due to the back-reaction of magnetic fields on turbulent flow (see, e.g. Cattaneo & Hughes 2006; Hughes & Proctor 2009), or because of rapidly oscillating alpha, with fluctuating part about three orders of magnitude greater than the constant mean part (see Hughes & Proctor 2009) of α. It has also been shown that the helicity alone in rotating convection is insufficient for dynamo action (see Cattaneo & Hughes 2006; Hughes & Proctor 2009). As such, the investigation of stochastic alpha effect is of great importance, and has received a considerable attention (e.g. Proctor 2007; Richardson & Proctor 2010; Mitra & Brandenburg 2012; Heinemann et al. 2011; Richardson & Proctor 2012). We also note that the limitation of mean field theory with constant transport coefficients has also been recognised in other systems (e.g. fusion plasmas) where the temporal and spatial variation of these transport coefficients has been invoked to improve the mean field theory (for instance see Hahm et al. 2004).

The purpose of this paper is to infer the statistical characteristics of the alpha effect by using observational data. To this end, we perform a systematic statistical comparison between irregularity in solar activity and stochasticity in α − Ω solar dynamo model. Specifically, we use the observational data of time history of solar activity inferred from C14 data in the past 11 000 years on a long time scale and direct observations of the sun spot numbers obtained in recent years 1795−1995 on a short time scale, and perform Monte Carlo simulations on these data to obtain probability distribution functions (PDFs) of the solar activity on both long and short time scales. To elucidate the key effect of stochastic dynamo process on solar activity, we adopt a simple nonlinear dynamical α-Ω model (see Weiss et al. 1984) and focus on the effect of temporal fluctuation in α effect. Compared with previous works on stochastic alpha dynamo, we undertake a systematic statistical study by varying the three key parameters of alpha: its magnitude, temporal correlation time and ratio of fluctuating to mean parts. For all these cases, the statistics of fluctuating alpha is assumed to be Gaussian with finite correlation time τ_α (i.e. coloured noise). We construct the PDFs of magnetic energy and then solar activity from our simulations which are compared with those from observational data, mentioned above. We show that the results from our dynamo model can match the PDFs of solar activity on time scales shorter than the 22 year solar cycle period.

The remainder of this paper is organized as follows. We present our stochastic dynamical model and numerical methods for solar dynamo in Sect. 2. Section 3 details the methods used to construct solar activity index from observational data and numerical simulations of our solar dynamo model on long and short time scales. From these solar activity index, we compute probability distribution functions (PDFs) of solar activity for observational data and our solar dynamo model in Sect. 4, and perform comparison of these PDFs in Sect. 5. Conclusions and discussions are presented in Sect. 6.

2. The stochastic dynamo model

To perform a systematic statistical study, we adopt a simplified non-linear α-Ω dynamo model represented by a dynamical system for the toroidal T and poloidal P magnetic fields (see Weiss et al. 1984). The equations are derived following Weiss et al. (1984) from the mean field dynamo equations in spherical polar co-ordinates for an axisymmetric magnetic field. Under the assumption that dynamo action occurs in a thin shell and the alpha effect is spatially homogeneous and isotropic, Cartesian co-ordinates are used to approximate an arbitrary element within this shell (x and y pointing east and north respectively, and z normal to surface). In this plane layer the velocity is assumed to vary with height in the y direction. Averaging vertically and considering the temporal evolution of a single (k = 1) fourier mode propagating in the x direction yields the following equations in non-dimensionalised form: (1)Specifically, this model describes the evolution of a single mode of toroidal and poloidal magnetic fields where the α effect acts to convert toroidal field T into poloidal field P, the Ω effect to converts the poloidal field into toroidal field and η is rate of the magnetic dissipation. The α is taken to have both mean and fluctuating parts as α(t) = α₀ + α′ = |α|(n₀ + n′(t)). Here, α₀ =  ⟨ α(t) ⟩  = |α|n₀ is the mean part while α′ = |α|n′ is the fluctuating part with ⟨ α′ ⟩  =  ⟨ n′(t) ⟩  = 0, where angular brackets ⟨  ⟩ represent the ensemble average. For fluctuating α′ = |α|n′, we assume it to have a normal distribution by taking n′ to be Gaussian and to have the following exponentially decaying temporal correlation

with correlation time τ_α. Note that the amplitude of the correlation function is chosen to make ⟨ (n₀ + n′)² ⟩  = 1. The ratio of the fluctuation to mean alpha is thus given by

In line with previous work (for example Richardson & Proctor 2010) we use the alpha effect quenching as whilst assuming that the Ω effect and turbulent dissipation of magnetic fields are unaffected by magnetic fields.

In dimensionless form, the equations are: (2)Here, η is the turbulent diffusivity of magnetic fields which will be taken to be η = 1 in our dimensionless unit with no loss of generality. In our dimensionless units, the typical oscillation timescale of the magnetic field is found to be 2π (see Fig. 1) for all parameter values. Therefore, by converting this dimensionless 2π to the 22 year solar cycle period, we obtain the magnetic dissipation time scale of 3.5 years of solar magnetic fields in physical units from η = 1. Note that the time resolution of our simulation is 4 days in this unit. To be consistent with current theories (see Dikpati & Charbonneau 1999), we choose ⟨ T² ⟩/⟨ P² ⟩  = 1000 by using | α| = 1000^−1/2.

For the accurate simulation of Gaussian noise with finite correlation time, we use the numerical technique by Honneycutt (the so-called Honneycutt SRK4 technique) (see Honeycutt 1992). Using this method, the ensemble average is taken over 50 different realizations obtained by simulating our models 50 times by using a cluster of 50 computers to determine PDFs accurately. Note that 50 different realizations for ensemble are obtained by using random initial conditions. We compute the statistical properties after our results reach a statistically stationary state.

Fig. 1 Frequency spectra of poloidal field, showing peaks at . Note that peaks continue to 3 × 10-5 under in red lines.

3. Activity index as a proxy for solar activity

For the comparison of our simulation results and observational data, we assume that sunspot number is linearly related to the magnetic energy (contained in the toroidal field) and introduce a normalized activity index (AI) by: (3)Data in Eq. (3) refers to (i) magnetic energy obtained from numerical simulation of the solar dynamo model and (ii) sunspot numbers obtained from both the C14 observation and recent group sunspot number data. Here, the angular brackets ⟨  ⟩ represent the mean value of the data. It is important to note that Eq. (3) is a linear mapping from which the original data can be recovered. The utilized observational data are the decadal averaged data of C14 infered fom tree rings since 10 000 B.C. (see Solanki et al. 2004) and monthly averaged data of group sunspot number from direct observations in 1605−1995 (see Hoyt & Schatten 1998). These are explained in details in the following two subsections.

3.1. C14 proxy data

For long time behaviour we use the solar activity reconstruction inferred from trapping of C14 in tree rings (see Solanki et al. 2004). The C14 record was obtained by carefully calibrating C14 concentrations in tree rings to the 350 year group sunspot number (henceforth referred to as the R_Gsunspot number) determined by Hoyt & Schatten (1998). Used in combination with models of the earth’s magnetic field, this permits us to determine the solar activity from the C14 concentration which extend back 11 000 years. The original C14 record is plotted in Fig. 2a. The C14 data is converted to AI by using Eq. (3), and noise in the data is then reduced by applying the Gleissburg 1-2-2-2-1 filter to increase the accuracy in determining extreme (maxima/minima) events (see Usoskin et al. 2007). Note that this filter has a slight effect on the mean solar activity. Negative sunspot numbers in this record represent data artefacts and are replaced with zero Solanki et al. (2004).

Fig. 2 Decadal and monthly sunspot activity taken from C14 and RG data respectively.

3.2. Sun sunspot number data

The group sunspot number is the best quality record of solar activity obtained between 1605 and 1995 (see Hoyt & Schatten 1998). This record has been collated from more than 60 000 observations by early to modern astronomers. Recent careful analysis of these observations has created a coherent index of solar activity known as the Group sunspot number. For our analysis, only the data between 1795 and 1995 are used since this period is the longest interval which has the continuous data without interruption. This data is shown in Fig. 2b, together with AI computed from Eq. (3). The red line in Fig. 2b represents the 10 year average value, which mimics the 10 year averaged C14 data. Note that 10 year average value in red line smooths out extreme events of both minima and maxima.

Fig. 3 Probability density function from C14 and RG sunspot number data. Solid lines show the mean PDFs while upper and lower dashed lines indicate the range of PDFs within one standard deviation.

3.3. Numerical simulation results

In order to compare our numerical simulation results of dynamo model (2) with solar activity reconstructions discussed above, we also convert our simulation results into AI as follows. First, we compute time series of magnetic energy B² = RE(T)² and smooth this time series on the two different time scales corresponding to the time intervals over which the two observational data were averaged. That is, we take the time average of magnetic energy from our simulation over 10 years and one month to be able to compare those with decadal averaged C14 data and monthly averaged R_G group sunspot number data, respectively. Second, the smoothing with the Gleissburg filter (see Solanki et al. 2004) is performed for both averaged time series. Finally, the resulting time series are then converted to AI by using Eq. (3).

4. Construction of PDFs

The original observational data are time traces of C14 and R_G sunspot activity with no direct information on the PDFs of solar activity. However, these data provide an estimate of the 1 standard deviation error bar associated with each data point. We perform Monte Carlo error estimation by using this 1 standard deviation error to obtain PDFs. This involves rather complicated technique and is thus detailed in Sect. 4.1. In comparison, the PDFs of magnetic energy from our dynamo model are numerically computed straightforwardly by considering 50 ensembles, as discussed in Sect. 4.2. Note that for all data and simulation results, 33 bins were used to obtain PDFs.

Here, we compare the two PDFs obtained from observational data. Figures 3a and b show the PDFs from the C14 and group sunspot numbers data, respectively. We have also analyzed the international sunspot number data in this paper and obtained almost identical PDFs to PDFs obtained previously from Group sunspot number. Figure 3b shows the same PDF of the monthly smoothed solar activity using the Group sunspot number and those using International sunspot number. Note that in this figure, the PDF of the International Sunspot number data are within the error range of the Group sunspot number. (Note that PDFs from numerical simulation are shown in Figs. 5 and 6) Dashed lines, solid lines, and crosses in Fig. 3 represent ± 1σ errors, mean PDF and the original data (i.e. original observational data with no error analysis), respectively. It is worthy noting that the mean PDFs are similar to original data besides that they are little smoother. Importantly, the original data typically lie within ± 1σ of the mean PDF. This thus confirms that our constructed PDFs of data are compatible with the original data, validating our statistical technique of computing PDFs from the time trace of the observational data by using Monte Carlo error estimate (see Sect. 4.1).

It is interesting that there are notable differences between the C14 data and Group sunspot data in the PDFs due to the different time interval over which the time average of these two data are taken (i.e. 10 yrs vs. 1 month), which can have a significant effect on extrema events, noted in Sect. 3.2. That is, on the C14 data, the average over 10 years combined with the Gleissburg filter smooths out extreme events, with the major effect of decreasing probability of minima and maxima events in comparison to the monthly averaged data. Similar effect was also manifested in the time trace shown in Fig. 2b. Furthermore, it has been demonstrated that the Gleissburg filter is very efficient at differentiating between prolonged minima/maxima events and it is precisely why it was used in our analysis (see Solanki et al. 2004), as noted previously in Sect. 3.1. It should also be noted that mean PDFs for both the C14 and sunspot data are strongly non Gaussian. Over plotted solid vertical lines in Fig. 3a represent the region of minima and maxima solar activity marked by SN < 15 and SN > 50, respectively.

4.1. Monte Carlo error estimation simulation for C14 and R_G sunspot number data

Both sets of observational data provide an estimate of the 1 standard deviation error bar associated with each data point. To construct probability density functions (PDFs), we first assume that the errors are independent random Gaussian noise with zero mean and standard deviation equal to the error provided at each data point. This enable us to create an ensemble of 100 000 different realizations of of AI for C14 data and R_G sunspot number data, which are then used to compute the PDFs. From these ensembles of 100 000 PDFs, we determine the mean PDF and the range of PDFs within one standard deviation of the ensembles in each bin. This process is the so-called Monte Carlo error estimation.

Fig. 4 AI of C14 data and dynamo model for τα = 4π, αR = 0.21.

4.2. Dynamo model results

It is much easier to compute PDFs from numerical simulations of our stochastic dynamo modl as we calculate them directly from the smoothed AI data for each of 50 ensembles and compute mean PDFs. The estimation of errors is obtained by computing one standard deviation in the PDFs of 50 ensembles in each bin. This allows us to reduce the error in numerically determined PDFs by running the simulation longer and making the error bars arbitrarily small.

5. Comparison of activity index and PDFs

5.1. C14 data and dynamo model

Before examining PDFs, we first compare the time trace of AI constructed from the magnetic energies obtained from the simulation of our dynamo model and the C14 data over 4000 yrs, which are shown in black and red lines in Fig. 4a. It is interesting to see that our dynamo model quantitatively reproduce the range of solar variability and periods of high and low activity. Furthermore, the time scale of maxima and minima events from our dynamo simulation appears to be consistent with the C14 reconstruction. Figure 4b shows the simulation results on a shorter time interval of 200 yrs where the black and red lines show the magnetic energy with and without 10 year average, respectively. Note again that the 10 yr averaged magnetic energy in red line deaccentuates extreme events. Note that also the oscillation shown in Fig. 4b is similar to that shown in Fig. 2b.

In order to infer the statistical properties of fluctuating α from observational data, we now compare PDFs from observational data and dynamo simulations. For comparison of PDFs, we vary the parameter values of numerical simulations within the range of τ_α ∈  [0.01,100], Ω|α| ∈  [0.001,100] (for the fixed value of ⟨ T² ⟩/⟨ P² ⟩  ~ 1000 and | α| = 1000^−1/2), and α_R ∈  [0.01,1000]. Note that the PDFs obtained by using parameter values outside these ranges have noticeable difference from the observational PDFs. We now make detailed comparisons of PDFs obtained by using these parameter values in order to further constrain these parameter values, which will enable us to infer the statistical properties of fluctuating α such as τ_α, Ω, and α_R.

Figure 5a−i shows the PDF of AI for Ω = 1000^1/2/3.5, τ_α = π,2π,4π, α_R = 0.21,0.24,0.29 (for fixed ⟨ T² ⟩/⟨ P² ⟩  ~ 1000 and | α| = 1000^−1/2). Different columns are for different values of the correlation time τ_α while different rows are for different values of α_R. We see that the agreement between the C14 data and our simulation results is good for a broad range of parameter values (see Figs. 5c, e, g). A close examination of Figs. 5c, e, g reveals a general tendency that a good agreement is achieved when τ_α increases (decreases) with decreasing (increasing) α_R. Specifically, for α_R = 0.24 and 0.29 with τ_α = π and 2π respectively, the agreement is particularly good. The case shown in Fig. 5g is on the boarder of overestimating low activity (see inset of Fig. 6g enlargement). This trend continues as τ_α further increases It is worth noting that the PDFs from dynamo simulations shown in Figs. 5c, e all lie very close to the one standard deviation errors of the C14 data for all values of the activity index.

However, increasing τ_α for fixed α_R or α_R for fixed τ_α leads to too fat tails, and/or too broadened PDFs (see Figs. 5f, h, i) to be compatible with C14 data; similarly, short τ_α or small α_R causes PDFs to become too narrow and Gaussian (see Figs. 5a, b, d), respectively, in disagreement with the PDFs from C14 data.

Fig. 5 Comparison of PDFs from C14 data and dynamo simulation for difference values of τα and αR. Red, solid and dashed line represent PDFs from dynamo simulation, mean, one standard deviation errors. Black, solid and dashed lines show the original C14 data and one standard deviation errors.

Fig. 6 Comparison of PDFs from C14 data and dynamo simulation for difference values of τα and αR. Red, solid and dashed line represent PDFs from dynamo simulation, mean, one standard deviations. Black, solid and dashed lines show the original C14 data and one standard deviation errors.

5.2. Group (R_G) sunspot number data and dynamo model

We now compare PDFs of the group sunspot data with those from our dynamo simulations for the same parameter values as used in the previous Sect. 5.1. Figures 6a, i show that the agreement in PDFs between dynamo simulation and R_G sunspot data are generally good for the same parameter values obtained for the best agreement between PDFs of C14 data and simulation results, as shown in Figs. 5c, e, g. For the best fit cases from Sect. 5.1, good agreement is seen in the τ_α = 4π in Fig. 6g. The cases of τ_α = 2π and τ_c = π overestimate the occurrence of high solar activity, as shown in Figs. 6c, e. Specifically the case τ_α = π case severely overestimates solar activity and is inconsistent with the monthly observed data.

5.3. Optimal parameter values of fluctuating alpha

The best agreement between the PDFs from the two observational data and those from our simulation is shown for parameter values of fluctuating α with τ_α = 2π (22 years) and α_R = 0.24. We remark that the Group sunspot number data represents a snapshot of solar activity as it covers only 18 solar cycles. As such, the PDF of activity taken from this data set is lacking the crucial information on the long term variation of solar activity. This could lead the to an underestimate of solar activity, explaining the overestimate of solar activity by our model. We note that outside the ranges considered for τ_α and α_R, features emerge in the PDF that are inconsistent with both C14 and R_G data. Specifically, significant overestimates of minima and underestimates of maxima are observed in comparison to both the C14 and R_G data. We thus conclude that the largest reasonable value of τ_α to be 44 years as the PDF shown in Fig. 5g lies at the boundary of overestimating the occurrence of minima events compared to the C14 data (see inset in Fig. 5g). On the other hand a lower limit seems to be about 22 years as the PDF in Fig. 6e slightly overestimates the occurrence of maxima events. In summary, we conclude that the α effect is constrained to be within the range τ_α ∈  [22,44] years with α_R ∈  [0.21,0.24].

6. Conclusion

We have inferred the statistical characteristics of the alpha effect by using observational data by performing a systematic statistical comparison between irregularity in solar activity and stochasticity in α − Ω solar dynamo model. Specifically, we used the observational data of time history of solar activity inferred from C14 data in the past 11 000 years on a long time scale and direct observations of the sun spot numbers obtained in recent years 1795−1995 on a short time scale, and performed Monte Carlo error estimate simulations on these data to obtain probability distribution functions (PDFs) of the solar activity on both long and short time scales. To elucidate the key effect of stochastic dynamo process on solar activity, we simulated a simple nonlinear dynamical α-Ω model (see Weiss et al. 1984) by assuming that α has both mean α₀ and fluctuating α′ parts. We have done a systematic statistical study by varying the three key parameters of alpha: its magnitude, temporal correlation time and ratio of fluctuating to mean parts. For all these cases, the statistics of fluctuating alpha are assumed to be Gaussian with finite correlation time τ_α (i.e. coloured noise). We construct the PDFs of magnetic energy and then solar activity from our simulations which are compared with those from observational data, mentioned above.

Through careful comparison of the PDFs averaged over different time scales, we have tuned our solar dynamo model to reproduce the distribution of observed activity. Specifically, we have shown that the magnitude of the fluctuating α is smaller than the mean part () with a correlation time between 22 and 44 years, with longer correlation times requiring a lower fluctuation ratio.

Despite the simplicity of our dynamo model, which was necessary for comprehensive statistical investigation, it has allowed us to elucidate key characteristics of the α effect. The fact that solar activity can be described by such a simple model is not only novel but also indicative of how the distribution of extreme event in complex system can be elucidated by reduced models.

We note that the α − Ω model lacks meridonal circulation modelling and it is not entirely clear as to how our results might apply to the flux transport model. It has however been recognised by both communities that understanding the effects of stochastically changing parameters are an important aspect of the dynamo problem. Therefore, the generic results of our results may also be valid in other stochastic dynamo models.

Finally, we have introduced and demonstrated novel statistical methodologies including the computation of activity index, the extraction of statistical properties by taking appropriate average using ensemble, different time intervals, etc. and computation of PDFs, especially by utilising Monte Carlos error estimation. These could stimulate further statistical study of solar dynamos and the comparison of observational data with theoretical prediction, by using more sophisticated models.

References

All Figures

	Fig. 1 Frequency spectra of poloidal field, showing peaks at . Note that peaks continue to 3 × 10-5 under in red lines.
In the text

	Fig. 2 Decadal and monthly sunspot activity taken from C14 and RG data respectively.
In the text

	Fig. 3 Probability density function from C14 and RG sunspot number data. Solid lines show the mean PDFs while upper and lower dashed lines indicate the range of PDFs within one standard deviation.
In the text

	Fig. 4 AI of C14 data and dynamo model for τα = 4π, αR = 0.21.
In the text

	Fig. 5 Comparison of PDFs from C14 data and dynamo simulation for difference values of τα and αR. Red, solid and dashed line represent PDFs from dynamo simulation, mean, one standard deviation errors. Black, solid and dashed lines show the original C14 data and one standard deviation errors.
In the text

	Fig. 6 Comparison of PDFs from C14 data and dynamo simulation for difference values of τα and αR. Red, solid and dashed line represent PDFs from dynamo simulation, mean, one standard deviations. Black, solid and dashed lines show the original C14 data and one standard deviation errors.
In the text