An empirical model of heliospheric cosmic ray modulation on long-term time scale

Galactic Cosmic Rays (GCRs) entering the heliosphere are subject to modulation processes due to variable solar magnetic activity. Finding a relationship between cosmic-ray variations and the heliospheric parameters is important for reconstruction of solar activity in the past. Here, we develop a semi-empirical model describing the heliospheric modulation of GCRs in terms of heliospheric parameters such as the open solar magnetic flux, the tilt angle of the heliospheric current sheet and the polarity of the large scale solar magnetic field. Our model is fitted using annual data obtained for the period 1976–2013, which includes the very weak solar minimum during 2008–2010. The model shows a good agreement with the data, and therefore, can be used for reconstructions of the modulation potential at different levels of solar activity. The model’s validity is also tested using the cosmogenic radionuclides C and Be stored in terrestrial archives. The tilt angle used to fit the parameters in our semi-empirical modulation model is reconstructed by a mathematical model described here.


Introduction
A theory describing Galactic Cosmic Rays (GCRs) transport in the heliosphere (Parker 1965;Krymskij 1969) includes four main processes: diffusion of particles along and perpendicular to the heliospheric magnetic field (HMF), outward convection and adiabatic deceleration in the expanding solar wind (SW), gradient and curvature drifts in the large scale HMF and drift due to the tilted heliospheric current sheet (HCS).This interplay results in spatial and temporal variations of GCRs intensity and energy.
Despite essential progress in theoretical development of the heliosperic transport models (see a review by Potgieter 2013 and references therein), a semi-empirical approach (e.g.Cliver 1993;Belov 2000;Alanko-Huotari et al. 2006) is still useful in many practical studies.Fully developed numerical 3D timedependent models of cosmic rays heliospheric transport do exist (e.g.Jokipii & Thomas 1981;Kota & Jokipii 1983;Hattingh & Burger 1995;Potgieter et al. 2001;Ferreira & Potgieter 2004;Potgieter 2013).However, such models include parameters that cannot be directly measured, such as, e.g., the diffusion tensor or the spectrum of turbulence of the solar wind magnetic inhomogeneities.Therefore, it is difficult to apply realistic computations using the sophisticated 3D models which are often based on ad hoc parameterisations (e.g.Ferreira & Potgieter 2004).As an alternative way, a number of (semi-)empirical models are based on observations, aiming to associate cosmic-ray variations with various solar and heliospheric parameters.Some of these empirical models are limited either by using cosmic rays intensity at a fixed energy (Stozhkov et al. 2004), or time series obtained by a single neutron monitor (Sabbah & Rybansk 2006), or neglecting possible long-term solar variability (Belov et al. 2006).
In this work, we use an approach similar to that by Alanko-Huotari et al. (2006), improved by using revised relations and the extended database of the heliospheric parameters for the period 1951-2013.Earlier models were developed based on observations corresponding to the period of high levels of solar activity, the so-called Modern Grand maximum of solar activity (Solanki et al. 2004).Accordingly, they may not well represent a possible centennial trend and periods with low solar activity.We fit our model for the period including the latest very weak solar cycle, which is similar to low activity levels in the past, such as the Dalton Minimum (DM) ca.1800 or the Modern minimum ca. 1900(McCracken & Beer 2014).Thus, our model aims at reconstructing the GCRs behaviour during periods with different solar activity levels.
We reconstruct, by applying our model, the modulation potential since 1616, a parameter that describes the mean energy loss of GCR particles within the heliosphere due to solar modulation, and using that, we model the global production of radionuclides 14 C and 10 Be, adopting the method described in Kovaltsov et al. (2012) and Kovaltsov & Usoskin (2010), respectively.The modelled series are next compared with the records from terrestrial archives such as tree rings (Roth & Joos 2013) and ice cores (Berggren et al. 2009).

Reconstruction of the HCS tilt angle
The tilt of the HCS is an important factor affecting modulation of galactic cosmic rays in the heliosphere (e.g., Jokipii & Thomas 1981;Kota & Jokipii 1983).Therefore, it is important to know its variability in order to study large scale HMF dynamics and the heliospheric modulation of cosmic rays on long time scale including the centennial trends.Although continuous observations of the tilt angle are made by the Wilcox Solar Observatory (WSO) since 1976, for the purpose of centennial reconstructions of the CR modulation it is important to know the tilt angle variability on longer time scale.Here we develop a model to reconstruct the HCS tilt angle using the phase of the solar cycle.This model describes the cyclic behaviour of the tilt angle (Hoeksema 1991;Suess et al. 1993) and its dependence solely on the phase of the solar cycle.
Variability of the HCS tilt angle follows the 11-year solar cycle along with the sunspot number (Fig. 1).However, although the latter varies significantly from cycle to cycle (as reflected in different cycle magnitudes varying by a factor of 2), the tilt angle exhibits roughly the same variations over all cycles, in accord with the idea of a regular cyclic behaviour (Suess et al. 1993;Cliver & Ling 2001;Alanko-Huotari et al. 2007), depending only on the solar cycle phase and not on its strength.This makes it possible to reconstruct the tilt angle variability using only the solar cycle phase.Moreover, strong variations in sunspot number amplitude during the solar maxima indicate that the sunspot number cannot be a good direct proxy for the tilt angle.
Figure 2 shows a superposed plot of the tilt angle as a function of the cycle phase for solar cycles 21 through 24.By comparing the tilt angle values from cycle to cycle, one can notice that the tilt angle has a similar behaviour during the ascending phase with a fast and smooth increase, for all the cycles.During the descending phase, the tilt angle fluctuates essentially, and these variations differ from cycle to cycle.However, the decrease is gradual in all cycles supporting the idea of a cyclic behaviour.Cliver & Ling (2001) concluded, from a similar analysis, that the descending phase of the tilt angle cycle is more gradual for odd cycles, but that is not observed in Figure 2, and moreover the statistic is too low (only two odd cycles) to substantiate such a conclusion.Alanko-Huotari et al. (2007) developed an empirical relation describing the cyclic behaviour of the tilt angle for the period 1976-2005.With the WSO tilt angle calculations extending until present, we revisited the model and fitted new parameters.Here we used, as the HCS tilt angle, the maximum latitudinal extent of the HCS, as provided by the WSO (http://wso.stanford.edu/Tilts.html).WSO offers two different approaches to estimate these values, using their photospheric magnetic field observations: One is the so-called ''classic'' model, and the other is the newer radial model.The ''classic'' model applies line-of-sight boundary conditions to the photosphere and therefore, it requires polar field corrections.Such corrections are not necessary for the radial model, since it applies radial conditions.Also, the radial model computes results with a higher source surface radius, and accordingly, the current sheet appears flatter and has a lower maximum extent.Here we use the values defined by the radial model, similar to Alanko-Huotari et al. (2007), in view of its more robust definition (Wang & Sheeley 1992).Averages of the values produced by WSO over a calendar year are derived and considered as annual values.
The best-fit empirical model, describing the annual value of the HCS tilt angle based on the progress of the solar cycle, takes the form: where is the year in the cycle and N is the cycle length in years.The value X i represents the phase of a cycle.We note that this relation is applied only to annual values.The length of the solar cycle 24 was defined, assuming the minimum of the cycle 24 in 2021 (Uzal et al. 2012).For this model, the ascending and maximum phases are considered to be shorter than the descending phase, reflecting the asymmetry of the tilt angle cyclic shape (Hathaway 2010).
Although observational limitations restrict the maximum angle to be 70° (Suess et al. 1993), this does not affect the use of this model's reconstructions for the study of GCRs heliospheric modulation, since convective-diffusive processes are the main contributors to cosmic rays modulation for tilt angle values greater than 50° (Lopate & Simpson 1991;Potgieter & Le Roux 1992;Cliver et al. 2013).The cyclic behaviour of the model is consistent with the cyclic shape followed by the WSO estimated values (Fig. 2).In order to test how well the tilt angle values, derived by this empirical relation, agree with the ones obtained by WSO, the root mean square (RMS) error between the two sets was estimated.The model was applied for the years 1976-2013, and the corresponding R rms is 6.6°.We note that this is a significant improvement with respect to the earlier model (Alanko-Huotari et al. 2007) which yields R rms = 11.2°forannual values.The fact that the model with more data gives a better fit to the data confirms its convergence and validity.Figure 3 shows the annual variations of the (radial) tilt angle provided by WSO and the one calculated here (Eq.( 1)).The two quantities are well correlated (R ¼ 0:956 þ0:012 À0:017 ) with the best correlation found during the ascending phase (Fig. 4).Some indirect sparse tilt angle estimates, based on image analysis of solar eclipses from 1870 through 2002, are done by Pishkalo (2006).We also show them in the same figure for comparison.They all appear to be in good agreement (see Fig. 3), confirming the robustness of the model.

Reconstruction of the modulation potential
Observable heliospheric parameters involved in the modulation processes described in Parker's theory are the solar wind velocity and density, the HCS tilt angle, the polarity, the strength and the level of turbulence of the HMF.These parameters vary with solar activity and the solar cycle phase.
In our model, we consider the heliospheric parameters that satisfy the following two main criteria (Alanko-Huotari et al. 2006).Firstly, they should describe the global heliosphere and secondly, they should be recorded or calculated continuously for a long period.Parameters that satisfy these criteria are the open solar magnetic flux (F), the tilt angle of the HCS (a) and the HMF polarity (p).The solar wind velocity is not considered within these parameters since its variations at 1 AU, in the ecliptic plane, do not represent the global behaviour of the solar wind, especially at high heliolatitudes.The open solar magnetic flux is used in the model as a global parameter of the HMF, instead of the HMF strength, B, measured in the ecliptic plane, since the latter is local and poorly represents the global HMF distribution.The time variations of these parameters are shown in Figure 5.
Parker's transport equation can be reduced, under some simplifying yet realistic assumptions, to the so-called force field approximation (Gleeson & Axford 1968;Caballero-Lopez & Moraal 2004), where the modulation is described by a single parameter, the modulation potential, /, which parametrises the shape of the GCR energy spectrum (see formalism in Usoskin et al. 2005) and is expected to be inversely proportional to the diffusion coefficient, j, of the heliospheric transport of GCRs, to some power n.Taking into consideration that the diffusion coefficient, j, is supposed to be roughly inversely related to the strength of the interplanetary magnetic field (Wibberenz et al. 2001b), which in turn is directly related, but not necessarily linearly proportional, to the OSF, one can assume that the modulation potential is roughly proportional to the OSF, / / F n (Usoskin et al. 2002).In some studies n is taken to be equal to unity (e.g.Kta & Jokipii 2001), resulting in a linear relationship between the modulation potential and OSF.On the other hand, some studies suggest that the relation is not exactly linear or quasilinear (Mursula et al. 2003;Alanko-Huotari et al. 2006).Alanko-Huotari et al. (2006) studied the correlation between / and F over periods of flat or highly tilted HCS and showed that the dependence of / on F varies with the tilt angle.This is supported by other studies according to which n can differ from unity (e.g.Caballero-Lopez & Moraal 2004), having a dependence on the solar cycle phase (Wibberenz et al. 2001b;Ferreira & Potgieter 2004) and vary with the HCS tilt angle (Wibberenz et al. 2001a;Ferreira & Potgieter 2004).
The role of the HMF polarity in the GCR modulation within the heliosphere varies with the phase of solar cycle and subsequently with the flatness/waviness of the HCS.During periods of minimum solar activity, when the HCS is flat, particles drift along the equatorial plane and away from the Sun for positive polarity periods, while the opposite occurs for negative polarity periods.For the negative polarity periods the modulation becomes stronger with increasing tilt angle (Alanko-Huotari et al. 2007; Dorman 2006 and references therein).During positive polarity periods an increasing tilt angle leads to declining modulation.However, the effects of the HCS on the modulation potential are not significant when the solar cycle reaches its maximum.
Accordingly, we use a model relating the heliospheric parameters to the modulation potential in the following form: where the mid-term defines the dependence of / on the OSF, that is slightly modulated by the HCS tilt angle, and the last term parametrises the drift effect of the HCS.In this model / 0 , a 0 , n and b are free parameters.
In order to find the best-fit parameters that minimise the mean square differences between the observed and the modelled /, we performed a non-linear model fitting using the annual modulation potential values, calculated from ground based cosmic ray observations (Usoskin et al. 2005(Usoskin et al. , 2011)), the OSF reconstructed by Lockwood et al. (2013aLockwood et al. ( , 2013bLockwood et al. ( , 2014aLockwood et al. ( , 2014b) ) and the tilt angle as described in Section 2. The polarity used for the fitting is defined as p = 1 for positive polarity periods, p = À1 for negative polarity periods and p = 0 for the years when polarity reversal occurs.Based on the Wilcox Solar Observatory (WSO) polar field observations since 1976, the polarity reversal years are taken to be 1980, 1990, 2000, 2013.For the years prior to 1976, when no direct data are available, we assume that the polarity reverses at the year of sunspot maximum.The fitting was done for the period 1951-2013, and the best fitted parameters are / 0 = 1473.9MV, a 0 = 150°, n = 1.03 and b = 0.095.
The annual variation of the computed modulation potential is shown in Figure 6 (magenta dashed curve).In the same figure the annual reconstructed values of the modulation potential by ground based cosmic ray data are plotted (blue curve).The correlation coefficient between the two is R = 0.88 ± 0.03 (Fig. 7).The modelled curve appears to follow well the one based on neutron monitor observations.However, a discrepancy appears during the maximum of the solar cycle 22.This may be related to the high plasma flow pressure during years 1991 and 1992.Though solar wind streams are an important factor of cosmic rays modulation, they could not be evaluated in the past and therefore it is difficult to include them in the model, which is intended for long-term studies.
The discrepancy during the maximum of solar cycle 24 can be explained in terms of the polarity dependence of the modulation potential in Eq. ( 2), and the polar field reversal around that period.According to WSO, the polarity reversed in the end of 2013.Therefore, in our model p was set to zero on that year and negative prior.Both from a theoretical and observational perspective, for large tilt angles, i.e. around solar maxima, GCR particles of all energies experience larger modulation during negative polarity epochs comparing to positive polarity epochs.And as we expect, our model yields high modulation,  for p = À1.However, based on neutron monitor measurements the modulation potential was lower.This can be explained considering the unusual polarity reversal during this maximum.More precisely, the maximum of solar cycle 24 was characterised by a rather slow and intermittent reversal of the polar field, with north-south asymmetry, meaning that the Northern polar field reversed on November 2012 and the Southern on March 2014, with both polar fields being simultaneously positive for more than a year (Sun et al. 2015).Therefore, although the polarity is officially reversed from negative to positive at the end of 2013, it was positive for a rather long time since 2012, leading to the lower level of modulation recorded by neutron monitors.
Here, we explain the discrepancy over the last solar maximum considering the parameters involved in our model, and more specifically the polarity and the model's response to such a change in the input data.Considering that during periods of high solar activity, other time-dependent transient phenomena, such as CMEs and propagating interaction regions (Burlaga et al. 1993;Wang et al. 2006), have an essential contribution to GCRs modulation, it is worth mentioning that they might be a source of differences between the modelled and observed series in Figure 6.More sophisticated models are based on these observable parameters (e.g.Burlaga et al. 1993) and describe well the modulation during solar maxima over the past few solar cycles.Although we acknowledge their importance, such phenomena were not recorded in the past, over time scales for which we want to reconstruct the modulation potential and thus they are not considered by our model.

Centennial reconstructions and radionuclide production
Cosmic rays particles interact in the atmosphere with the nuclei of atmospheric gases, leading to the production of cosmogenic radionuclides, such as 14 C and 10 Be.Solar variability and the geomagnetic field affect the flux of cosmic rays entering the atmosphere and subsequently, the radionuclide production rate (Beer et al. 2012).After their production, the radionuclides, following different distribution and deposition processes, are eventually stored in terrestrial archives such as tree rings or ice cores.The radionuclide signal reflects mainly production (Berggren et al. 2009), making them a useful tool of solar variability reconstructions in the past (Usoskin 2013).However, contribution from local climate to 10 Be deposition may be  2009) (magenta curve) and modelled here (blue curve).Although the two curves show a very good agreement after 1900, they vary from each other significantly before that period.This is likely due to atmospheric distribution and deposition processes that affect the stored amount of the radionuclide in ice.E. Asvestari and I.G.Usoskin: Empirical model of cosmic ray modulation essential at interannual-decadal time scales (Pedro et al. 2006(Pedro et al. , 2011;;Usoskin et al. 2009).
We reconstructed the modulation potential since 1616, using the semi-empirical modulation model (Eq.( 2)), where for the OSF we used the one published by Lockwood et al. (2014a) and Lockwood & Owens (2014), the polarity p was used as defined in Section 3 and the tilt angle a used, is a series reconstructed using Eq.(1).To define the length of the solar cycles since 1616, and subsequently the cycle phase input in Eq. ( 1), we considered the sunspot number maxima and minima as listed in SIDC, version 1.0.The obtained modulation potential series was applied first to the global radiocarbon production rate model (Kovaltsov et al. 2012), and then to the 10 Be flux model (Kovaltsov & Usoskin 2010).In Figure 8, we compare the calculated global radiocarbon production rate with the reconstruction by Roth & Joos (2013) based on tree ring records.The two curves agree well prior to 1900.Due to fossil fuel burning (the so-called Suess effect), since the late 19th century, and atomic bomb testing after 1950, the 14 C records cannot be used for the modern epoch.
We also compared (Fig. 9) the calculated 10 Be flux, with the one reconstructed using the ice core records from Greenland NGRIP (Berggren et al. 2009).Our computations show similar long-term trends to the reconstruction by Berggren et al. (2009).It is worthy to note that the two curves show an excellent agreement with each other over the 20th century.Prior to this period, the two curves appear not to show similar solar cycle variations, which is likely due to regional climatic effects (Usoskin et al. 2009) leading to strong fluctuations of the 10 Be deposition (Pedro et al. 2006(Pedro et al. , 2011)).

Conclusions
In this paper we introduce an empirical model to reconstruct the HCS tilt angle annual variations and a semi-empirical model for centenial reconstructions of the modulation potential.The HCS tilt angle model depends only on the solar cycle phase and describes the cyclic behaviour of the HCS tilt angle, reflecting its asymmetric shape, with short length and fast ascending phase in contrast to the gradual descending phase.The reconstructed series agree well with the observations by WSO as well as with sparse estimates from image analysis of solar eclipses (Pishkalo 2006).Due to observational constraints there is an upper limit in the tilt angle maximum at 70°, but this has no effect on the application of the reconstructed tilt angle in studying the heliospheric modulation of GCRs.
The semi-empirical modulation model relates the GCR's modulation potential with those heliosperic parameters measured at 1 AU and describes the global heliospheric solar activity.These are the open solar magnetic flux, the HCS tilt angle and the solar magnetic field polarity.The model fitting is done for the period 1951-2013 and gives a good correlation between the modelled and the reconstructed by ground based measurements' modulation potential.However, there appears to be a small deviation during the solar maximum of cycle 22, which is likely related to high pressure of solar wind plasma flow velocity.Such events are hard to predict and subsequently to include in the model, but their effects may be significant.The second discrepancy occurring during the maximum of solar cycle 24 can be explained in terms of the unusually extended polar field reversal, with both northern and southern polar fields being simultaneously positive for over a year, leading to a higher flux of GCR particles at Earth.
To test the validity of our model, we reconstructed the modulation potential since the 1616, which was subsequently used to compute the global radiocarbon production rate and the 10 Be flux.The computations show similar long-term variations with the global radionuclides production records from terrestrial archives such as tree rings and ice cores which validate the approach.The period studied includes different levels of solar activity and therefore the semi-empirical model can be considered as a good approach for reconstructions in centennial and millennia scales.

Fig. 2 .Fig. 1 .
Fig. 2. Yearly averaged HCS tilt angle variations based on the model (Eq.(1)) and the WSO observations of the photospheric magnetic field in relation to the cycle phase, for solar cycles 21-24 as indicated in the legend.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3. Temporal variations of the HCS tilt angle based on the radial model by WSO (blue line) and reconstructed here using the empirical model (Eq.(1)) (magenta steplike line).Black stars represent reconstructed HCS tilt angles by an analysis of total solar eclipse images (courtesy of M.I.Pishkalo).

Fig. 6 .
Fig. 6.Temporal variations of the modulation potential reconstructed by ground based observations of cosmic rays (blue curve) and the modelled (magenta dashed curve) over the period 1951-2013.

Fig. 7 .
Fig. 7. Scatter plot of the data shown in Figure 6.The solid line represents the diagonal.The correlation coefficient, R, between the two parameters is approximately 0.88 ± 0.03.