Propagation of Cosmic Rays in Heliosphere: the HelMod Model

TL;DR

The exttt{HelMod} model simulates cosmic ray propagation in the heliosphere, incorporating drift effects and matching observational data from multiple space missions during solar cycles 23 and 24.

Contribution

This paper presents the current version 3 of the exttt{HelMod} model, including detailed implementation of drift effects and validation against experimental data.

Findings

exttt{HelMod} reproduces measured cosmic ray spectra within experimental errors.

Model results agree with Ulysses measurements at high solar latitudes.

The model accurately simulates cosmic ray modulation during solar cycles 23 and 24.

Abstract

The heliospheric modulation model \helmod{} is a two dimensional treatment dealing with the helio-colatitude and radial distance from Sun and is employed to solve the transport-equation for the GCR propagation through the heliosphere down to Earth. This work presents the current version 3 of the \helmod{} model and reviews how main processes involved in GCR propagation were implemented.~The treatment includes the so-called particle drift effects --~e.g., those resulting, for instance, from the extension of the neutral current sheet inside the heliosphere and from the curvature and gradient of the IMF --, which affect the transport of particles entering the solar cavity as a function of their charge sign.~The \helmod{} model is capable to provide modulated spectra which well agree within the experimental errors with those measured by AMS-01, BESS, PAMELA and AMS-02 during the solar cycles 23 and 24.~Furthermore, the counting rate measured by Ulysses at +/- 80degree of solar latitude and 1 to 5\,AU was also found in agreement with that expected by \helmod{} code version 3.

Tables

Table 1. Table 1: Parameters c i subscript 𝑐 𝑖 c_{i} of the polynomial expression ( 6 ) as a function of solar polarity and phase. In the last column the rms value of the relative differences is shown.

	$c_0$	$c_1$	$c_2$	$c_3$	rms
A0 Ascending	0.0003059	-2.51e-6	1.284e-8	-2.838e-11	0.1097
A0 Descending	0.0002876	-3.715e-6	2.534e-8	-5.689e-11	0.1400
A$>$0 Ascending	0.0002262	-5.058e-7	–	–	0.1153
A$>$0 Descending	0.0002267	-7.118e-7	–	–	0.1607

Table 2. Table 2: Parameters p i subscript 𝑝 𝑖 p_{i} of the exponential expression ( 7 ) for the high solar activity. In the last column the rms value of the relative differences is shown.

	$p_0$	$p_1$	$p_2$	rms
MCMU	0.003753	-0.04791	0.0001365	0.100
OULU	0.001354	-0.10070	0.0007697	0.094

Table 3. Table 3: Parameters of the analytical fits to the proton and He LIS’s (from Table 5 of Boschini et al. ( 2017 ) .

	$a_0$	$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$b$	$c$	$d_1$	$d_2$	$e_1$	$e_2$	$f_1$	$f_2$	$g$
p	$94.1$	$-831$	0	$16700$	$-10200$	0	$10800$	$8590$	$-4230000$	$3190$	$274000$	$17.4$	$-39400$	$0.464$	0
He	$1.14$	0	$-118$	$578$	0	$-87$	$3120$	$-5530$	$3370$	$1.29$	$134000$	$88.5$	$-1170000$	$861$	$0.03$

Equations

$$\frac{\partial U}{\partial t}=-\nabla \cdot (U\vec{V})+\nabla \cdot \left[K^S\cdot \nabla U\right]+\frac{(\nabla \cdot {\vec{V}}_{sw})}{3}\frac{\partial }{\partial T}\left({\alpha }_{rel}TU\right),$$

$$\vec{V}={\vec{V}}_{sw}+{\vec{v}}_{drift},$$

$${\vec{v}}_{drift}=\nabla \cdot K^A$$

$${\alpha }_{rel}=\frac{T+2T_0}{T+T_0}.$$

$$\frac{\partial f}{\partial t}=-\nabla \cdot (f\vec{V})+\nabla \cdot \left[K^S\cdot \nabla f\right]+\frac{(\nabla \cdot {\vec{V}}_{sw})}{3p^2}\frac{\partial }{\partial p}\left(p^{3}f\right).$$

$$J=p^{2}f=\frac{\beta cU}{4\pi },$$

$$K_{||}=\frac{\beta }{3}{K}_0\left[\frac{P}{1\text{GV}}+g_{low}\right]\left(1+\frac{r}{\text{1AU}}\right),$$

$$P=\frac{pc}{|Z|e}$$

$$\frac{K_{\perp ,i}}{K_{||}}={\rho }_i,$$

$$0.022<{\rho }_i<0.083$$

$$K_0^{SSN}=c_0+c_1\mathrm{SSN}+c_2{\mathrm{SSN}}^2+c_3{\mathrm{SSN}}^3$$

$$K_0^{NMCR}=p_0\exp \left(p_1\mathrm{NMCR}+p_2{\mathrm{NMCR}}^2\right)$$

$${\vec{B}}_{Par}$$

$${\vec{B}}_{Pol}$$

$$\Gamma =\frac{\Omega (r-r_0)\sin \vartheta }{V_{sw}}.$$

$$\delta (\vartheta )=\frac{{\delta }_m}{\sin \vartheta },$$

$$J_{LIS}(P)\times P^{2.7}=$$

$$\displaystyle\left\{\begin{array}[]{ll}\sum_{i=0}^{5}a_{i}P^{i},&P\leq 1\ {\rm GV},\vskip 3.0pt plus 1.0pt minus 1.0pt\\ b+\frac{c}{P}+\frac{d_{1}}{d_{2}+P}+\frac{e_{1}}{e_{2}+P}+\frac{f_{1}}{f_{2}+P}+gP,&P>1\ {\rm GV},\end{array}\right.$$

$$\Delta r=$$

$$-\left[V_{SW}+v_{drift,r}\right]\Delta t$$

$$+(2K_{rr}{)}^{1/2}{\omega }_r\sqrt{\Delta t},$$

$$\Delta \theta =$$

$$+\frac{2K_{r\theta }}{r\sqrt{2K_{rr}}}{\omega }_r\sqrt{\Delta t}$$

$$+{\left[\frac{2K_{\theta \theta }}{r^2}-\frac{2K_{r\theta }^2}{r^2{K}_{rr}}\right]}^{1/2}{\omega }_{\theta }\sqrt{\Delta t},$$

$$\Delta T=$$

$$L=$$

$$J_{Earth}(T_i)=\frac{\beta (T_i)}{N}\sum _{k=1}^N\frac{J_{LIS}(T_{r,k})}{\beta (T_{r,k})}\cdot \exp (-L_{Int,k}).$$

$$J_{\mathrm{Earth}}(P_0)={\int }_0^{\infty }{J}_{\mathrm{LIS}}(P)G(P_0|P)dP.$$

$$\frac{\partial F}{\partial t}=\frac{1}{2}\sum _{i,j}[B{B}^{\top }{]}^{ij}\frac{{\partial }^{2}F}{\partial x_i\partial x_j}+\sum _i{A}_B^i\frac{\partial F}{\partial x_i}+LF.$$

$$d{x}_i=A_B^{i}dt+\sum _j{B}^{ij}d{\omega }_j.$$

$$\frac{\partial U}{\partial t}=$$

$$+\frac{1}{r\sin \theta }\frac{\partial K_{\varphi r}}{\partial \varphi }\frac{\partial U}{\partial r}-V_{sw}\frac{\partial U}{\partial r}-V_{dr}\frac{\partial U}{\partial r}$$

$$+\frac{1}{r^2}\frac{\partial r{K}_{r\theta }}{\partial r}\frac{\partial U}{\partial \theta }+\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }\left(\frac{K_{\theta \theta }\sin \theta }{r}\right)\frac{\partial U}{\partial \theta }$$

$$+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \varphi }\left(K_{\varphi \theta }\right)\frac{\partial U}{\partial \theta }-\frac{V_{d\theta }}{r}\frac{\partial U}{\partial \theta }$$

$$+\frac{1}{r^2}\frac{\partial }{\partial r}\left(\frac{r{K}_{r\varphi }}{\sin \theta }\right)\frac{\partial U}{\partial \varphi }+\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }\left(\frac{K_{\theta \varphi }}{r}\right)\frac{\partial U}{\partial \varphi }$$

$$+\frac{1}{r^2{\sin }^2\theta }\frac{\partial }{\partial \varphi }\left(K_{\varphi \varphi }\right)\frac{\partial U}{\partial \varphi }-\frac{V_{d\varphi }}{r\sin \theta }\frac{\partial U}{\partial \varphi }+\frac{2}{3}\frac{\alpha T{V}_{sw}}{r}\frac{\partial U}{\partial T}$$

$$+K_{rr}\frac{{\partial }^{2}U}{\partial r\partial r}+\frac{K_{r\theta }}{r}\frac{{\partial }^{2}U}{\partial r\partial \theta }+\frac{K_{r\varphi }}{r\sin \theta }\frac{{\partial }^{2}U}{\partial r\partial \varphi }$$

$$+\frac{K_{\theta r}}{r}\frac{{\partial }^{2}U}{\partial \theta \partial r}+\frac{K_{\theta \theta }}{r^2}\frac{{\partial }^{2}U}{\partial {\theta }^2}+\frac{K_{\theta \varphi }}{r^2\sin \theta }\frac{{\partial }^{2}U}{\partial \theta \partial \varphi }$$

Full text

††journal: Advances in Space Research

1 Introduction

Modulated omni-directional intensities of galactic cosmic rays (GCRs) were observed during different phases of solar activity using both balloon flights (for instance, see Boezio et al., 1999; Menn et al., 2000; Haino et al., 2004; Shikaze et al., 2007; Abe et al., 2008, 2016) and space-borne missions (e.g., see Alcaraz et al., 2000d, c, a, b; Aguilar et al., 2002, 2007; Adriani et al., 2009a, b, 2010, 2011, 2013, 2015, 2016; Aguilar et al., 2014, 2015b, 2015a, 2016), in particular during the latest solar cycles. The increased performance of on-board spectrometers was and is currently enabling to enhance the accuracy of the observed spectra. Thus, it was opening the way to a better understanding of processes related to the transport of GCRs through the Heliosphere – the so-called modulation effect – and, ultimately, to the capability a) to unveil local interstellar spectra (LIS) of GCR species (e.g., see Bisschoff and Potgieter, 2014, 2016; Della Torre et al., 2017a; Boschini et al., 2017, and also references therein), b) to investigate their generation, acceleration and diffusion process within the Milky Way (e.g., see Boella et al., 1998; Strong et al., 2007; Evoli et al., 2008; Putze et al., 2009), and, in turn, c) to possibly untangle features due to new physics – i.e., dark matter (e.g., see Bottino et al., 1998; Cirelli and Cline, 2010; Ibarra et al., 2010; Salati, 2011; Weniger, 2011, and references therein) – or additional astrophysical sources so far not taken into account (e.g., see Chang et al., 2008; Abdo et al., 2009; Adriani et al., 2009a; Cernuda, 2011; Mertsch and Sarkar, 2011; Della Torre et al., 2015; Rozza et al., 2015, and references therein).

Among space missions currently observing GCRs, AMS-02 – on-board of the International Space Station since May 2011 – is continuously providing data with unprecedented measurement accuracy. In fact, this spectrometer allowed one to determine the most accurate differential intensities of protons (Aguilar et al., 2015b), helium nuclei (Aguilar et al., 2015a), antiprotons (Aguilar et al., 2016), electrons and positrons (Aguilar et al., 2014). The high precision of these experimental data together with those from Ulysses spacecraft (e.g., see Simpson et al., 1992, 1996; Heber et al., 1996; Ferrando et al., 1996; De Simone et al., 2011; Gieseler and Heber, 2016) constitute a challenge for any modulation model of the inner part of heliosphere. Actually such a treatment has to reproduce the observed GCR spectra transported – during different solar activity phases – down to Earth and, also, outside the ecliptic plane at distances from Sun ranging from about 1 to 5 AU. In fact, the observations made using Ulysses spacecraft allowed one to determine both the latitudinal and radial dependence of GCR intensity. Furthermore, the data taken during Ulysses fast latitudinal scan exhibited a latitudinal dependence on i) the charge sign of the GCR species (i.e., protons and electrons, which are the dominant positively and negatively charged component, respectively), ii) solar activity and iii) polarity of the interplanetary magnetic field (IMF). It is worth to remark that these data may, in addition, allow a better understanding of space radiation environment close to Earth, thus extending our capability of predicting radiation hazards for astronauts and device damages in space missions (e.g., see Leroy and Rancoita, 2007; Golge et al., 2015, and Chapters 7 and 8 of Leroy and Rancoita (2016)).

In this work we present the version 3 of 2-D heliospheric modulation (HelMod) model – i.e., a two dimensional treatment dealing with the helio-colatitude and radial distance from Sun (Gervasi et al., 1999; Bobik et al., 2012, 2013b) – currently employed to solve the transport-equation for the GCR propagation through the heliosphere down to Earth. The relevant GCR propagation processes are described in Sects. 2 and 3 in order to illustrate how they are implemented in his latest version of HelMod model. Furthermore, details on Monte Carlo technique used to solve the stochastic integration are treated in Sects. 4, 5, and allow one a better understanding on HelMod capabilities to deal with solar modulation, within the inner part of the heliosphere. At present, the model only treats GCRs with energies $\gtrsim 0.5$ GeV/nucleon, thus modulation effects occurring in the outer heliosphere – i.e., beyond the termination shock (TS) (see, e.g., Langner et al., 2003; Langner and Potgieter, 2004; Bobik et al., 2008; Potgieter, 2008; Florinski and Pogorelov, 2009; Luo et al., 2013; Senanayake and Florinski, 2013) – are not accounted for.  It has to be pointed out that the HelMod model is capable of describing the current large set of observation data, which were collected during solar cycles 23 and 24 with the occurrence of two solar minimum. For this purpose, the model includes the so-called particle drift effects – e.g., those resulting, for instance, from the extension of the neutral current sheet inside the heliosphere and from the curvature and gradient of the IMF –, which affect the transport of particles entering the solar cavity as a function of their charge sign. These effects are particularly relevant when IMF exhibits a well-defined large-scale structure or this latter is still relevant. In fact, at the solar minimum and when the solar activity is not too far from such a condition, GCR modulated intensities exhibit a dependence on charge sign (e.g., see Garcia-Munoz et al., 1986; Clem et al., 1996, 2000; Boella et al., 2001). In fact, the IMF polarity reversal causes charge sign dependent modulation effects, for instance, those observed in particle over anti-particle intensities ratio at rigidities lower than about 10–20 GV (e.g., Adriani et al., 2016). These effects are treated in the Parker transport equation through the terms including the drift velocity. The analysis on Ulysses out-of-ecliptic observations (e.g., see Simpson, 1996; Simpson et al., 1996; Heber et al., 1996, 1998, 2008; Ferrando et al., 1996; De Simone et al., 2011; Gieseler and Heber, 2016) provided, so far, a unique point of view highlighting the presence of latitudinal gradients in the spatial distribution of GCRs, during period of low solar activity, i.e., when the combination of particle charge ($q$) and solar magnetic polarity ($A$) is positive ($qA>0$); while a more uniform distribution of GCRs in the inner part of the heliosphere occurs for $qA<0$.

As discussed in Bobik et al. (2011b, 2012), the model exhibits a smooth time dependence introduced by the parameters – related to solar activity and adopted within the model itself, as described in Sect. 2.3 –, which are averaged over time durations corresponding to Carrington rotations, i.e., a) solar wind speed ($V_{\rm sw}$, see Sect. 2.2), b) tilt angle ($\alpha_{\rm t}$, see Sect. 3) of the neutral current sheet, and c) diffusion parameter ($K_{0}$, see Sect. 2.3). Furthermore, it has to be remarked that the solar wind usually takes one year or even more to reach the border of heliosphere. The above parameters – usually determined at 1 AU – are transferred to describe the properties of any distant heliospheric sector, according to the time required by the solar wind to reach such a region from Sun (see discussion in Sect. 3).

In the present article, the LIS fluxes are those derived in Della Torre et al. (2017a); Boschini et al. (2017) by means of GALPROP v55. In that article leptons were not yet treated, i.e., no electron LIS is available yet using the new GALPROP version. Thus, the discussion on the modulated spectra obtained using HelMod model is mostly restricted to comparisons with experimental data regarding protons and helium nuclei. As discussed in Della Torre et al. (2017a); Boschini et al. (2017), the LIS, presented in Sect. 3, accommodate both the low energy interstellar CR spectra measured by Voyager 1 and the high energy observations by BESS, Pamela, AMS-01, and AMS-02 over solar cycles 23–24.

Finally, we have to remark that Engelbrecht and Burger (2013) exploited an ab-initio approach for a three dimensional steady state GCR modulation model, in which the effects of turbulence on both the diffusion and drift of these cosmic-rays are treated in a self-consistent description; Strauss et al. (2013) uses a hybrid modeling approach incorporating the plasma flow from a magnetohydrodynamics model with the particle transport; and, finally, Vos and Potgieter (2016, and reference therein) computed spatial gradients and absolute flux variations for GCR protons in the heliosphere for solar minimum. Although these models provide encouraging results, they still depends on parameters whose time evolution is not yet measurable or fully understood. So far, the found agreement among the modulated spectra from HelMod code and experimental data collected over a long period (e.g., see Bobik et al., 2012, 2013b; Della Torre et al., 2017a; Boschini et al., 2017) motivated the choice, in HelMod, to reduce the complexity of diffusion process using a unique time dependent variable, as described in Sect. 2.2.

2 Heliospheric Propagation of Cosmic Rays

2.1 Parker Equation

The cosmic rays propagation trough the heliosphere was treated by Parker (1965), who demonstrated that - in the framework of statistical physics - the random walk of cosmic ray particles is a Markoff process, describable by a Fokker–Planck equation (FPE).  In his original formulation, Parker’s transport-equation was expressed in terms of particle density for unit space and energy, i.e., $U(\vec{x},T)$ (e.g., see Jokipii and Parker, 1970; Fisk, 1971; Bobik et al., 2012, and also Sections 8.2–8.2.5 of Leroy and Rancoita, 2016):

[TABLE]

with $\vec{x}$ the 3D-spatial position in Cartesian coordinates,

[TABLE]

$\vec{V}_{sw}$ the solar wind (SW) velocity,

[TABLE]

the drift velocity (e.g., see Jokipii et al., 1977; Jokipii and Levy, 1977; Bobik et al., 2012, and references therein), $K^{A}$ and $K^{S}$ the antisymmetric and symmetric part of the diffusion tensor, respectively. In Eq. (1) $T$ is particle kinetic energy, $T_{0}$ is particle rest energy and, finally,

[TABLE]

We should note that some authors prefer to describe the modulation of GCRs re-expressing Eq. (1) in terms of the so-called omni-directional distribution function $f(\vec{x},p)$, where $p$ is particle momentum (e.g., see Jokipii and Kopriva, 1979; Yamada et al., 1998; Pei et al., 2010, and Section 8.2.4 of Leroy and Rancoita, 2016):

[TABLE]

Although Eqs. (1) and (3) describe the same process (see, e.g., Schlickeiser, 2002), they are commonly used as alternative formulation; therefore the SDE integration technique (described in Sect. 4) used in this work should be performed differently as described in Bobik et al. (2016).

Parker’s transport-equation describes i) the diffusion of GCRs by magnetic irregularities, ii) the so-called adiabatic-energy changes associated with expansions and compressions of cosmic radiation, iii) an effective convection resulting from the solar wind (SW, with velocity $\vec{V}_{{\rm sw}}$) convection effect and iv) the drift effects related to the drift velocity ($\vec{v}_{drift}$). In turn, the drift velocity is determined by the antisymmetric part of the diffusion tensor [see Eq. (2)] which accounts for gradient, curvature and current sheet drifts of particles in the IMF, i.e., it depends on the charge sign of particles.

As discussed in Sect. 2.2, the low scale irregularities are the main responsible of particle diffusion in the heliosphere. Nevertheless, large scale structures due to SW expansion contribute to convective and magnetic drift motion in Eq. (1). Although present in the original formulation of Parker’s equation, in early modulation models the drift velocity was neglected. In recent years, it was found that models accounting for charge dependent effect by means of drift velocity are, indeed, capable to describe the observations (e.g., see Della Torre et al., 2012; Bobik et al., 2013a; Maccione, 2013).

In HelMod model, we use the drift treatment originally developed by Potgieter and Moraal (1985) and refined using Parker’s magnetic field with polar correction described in Sect. 2.4. The full description and derivation can be found in Bobik et al. (2013b). Since during high activity periods the heliospheric magnetic field is far from being considered regular, we introduced a correction factor suppressing any drift velocity at solar maximum. For sake of completeness, we have to note that the presence of turbulence in the interplanetary medium should reduce the global effect of CR drift in the heliosphere (e.g., see discussion in Minnie et al. 2007) and this is usually taken into account introducing a drift suppression factor (e.g., see Strauss et al., 2011) that is effective at rigidities below 1 GV (for a general discussion about the energy range of drift effects in the solar modulation see, e.g., Nndanganeni and Potgieter, 2016).

As remarked by Jokipii et al. (1977), since $\nabla\cdot\vec{v}_{drift}=0$, one finds that drift velocity is added to that of solar wind resulting effective convection speed but do not contribute to the adiabatic-energy losses [third right-hand term of Eqs. (1, 3)]. Nevertheless, his contribution causes an interplay with the diffusion process resulting in different modulated spectral slope for period of same solar activity but opposite magnetic field polarity (see, e.g., Potgieter, 2013).

In HelMod model version 3, the description of the solar wind is that presented in Bobik et al. (2012), i.e., it is constant (from Earth up to the TS) and its value is that measured at Earth. During periods of low solar activity, the solar wind speed increases by almost a factor two from the ecliptic plane to the poles, thus subdividing the heliosphere into two regions with slow and fast solar wind (McComas et al., 2000); no latitudinal dependence is considered during high activity periods (Bobik et al., 2012).

We should note that second-order Fermi acceleration mechanisms are usually neglected in general treatment of the propagation of GCRs towards the inner part of heliosphere, while additional energy-loss processes – beside that adiabatic included in Eq. (1) – (like e.g. inverse-Compton) are negligible (e.g., see Bobik et al., 2011a).

Finally, the differential intensity of GCRs, $J$, is related to the particle density ($U$) and to the omni-directional distribution function ($f$) [see Eqs. (1, 3)] as:

[TABLE]

where $\beta$ is the particle velocity expressed in units of the speed of light, $c$.

2.2 The Diffusion Tensor

When Parker (1965) introduced his transport equation, he underlined how the diffusion process experienced by cosmic rays is mainly resulting from the magnetic irregularities of the IMF. Such magnetic-field irregularities, he remarked in that work, are transported rigidly in the SW and appear with dimensions comparable with the gyro-radius of about 0.1–10 GeV proton moving in an IMF of a few nT. The scattering of energetic particles on them causes a random walk, which is accounted for by a diffusion tensor having a parallel component ($K_{||}$) – with respect to the direction of large scale IMF – larger than that perpendicular [$K_{\perp,i}$, where the subscript $i$ refer to radial ($r$) or latitudinal ($\theta$) direction components]. For rigidity greater than 1 GV, the HelMod code version 3 implements a functional form that is linearly dependent to particle rigidity ($P$) and linearly proportional to solar distance ($r$) (e.g., see Boschini et al., 2017):

[TABLE]

where $K_{0}$ is the diffusion parameter – described in Sect. 2.3 – which dependence on time reflect the variability of interplanetary medium properties for the different phases of solar activity,

[TABLE]

is the particle rigidity expressed in GV with $p$ the particle momentum, $r$ is the heliocentric distance from Sun in AU and, finally, $g_{low}$ is a parameter (see a further discussion in this Section), which depends on the level of solar activity and allows the description of the flattening with rigidity below few GV. In the present model, the spatial dependence is proportional to the distance ($r$) from Sun; it is consistent with that used in Bobik et al. (2013b) for dealing with CR latitudinal gradients and no further latitudinal dependence appears to be needed (see also discussions in Jokipii and Kota, 1989; McDonald et al., 1997; Strauss et al., 2011).

It is important to remark how from the observation of fluctuations in the average magnetic field, Jokipii (1966) and Jokipii (1971) put the basis for the first successful description of diffusion tensor with the so-called *quasi-linear theory *(QLT) under the approximation of a weak turbulence. However, a complete parametrization for the components of diffusion tensor, from low to high rigidities, are still an open question. Nowadays, it is commonly accepted that (see, e.g., Shalchi, 2009; Engelbrecht and Burger, 2015, and reference therein): i) at higher rigidities the diffusion coefficient should have a quasi-linear dependence (e.g., see Gloeckler and Jokipii, 1966; Gleeson and Axford, 1968; Jokipii, 1966, 1971; Perko, 1987; Potgieter and Le Roux, 1994; Strauss et al., 2011), and ii) the diffusion coefficient should be more “flat” at lower rigidity (e.g., see Palmer, 1982; Bieber et al., 1994). In previous models of CR propagation in the heliosphere, at intermediate energies, to account for the above discussed conditions the parallel diffusion coefficient ($K_{||}$) was expressed with a sharp transition at $\sim 1$ GV between the two above mentioned regimes (e.g., see Perko, 1987; Alanko-Huotari et al., 2007; Strauss et al., 2011; Bobik et al., 2012). However, with increasing the experimental accuracy of collected data such a simple approach needs to be revised allowing a smoother transition between the two regimes. In fact, the present functional form of such a transition [Eq. (5)] is consistent with those presented in Burger and Hattingh (1998) for the same rigidity interval.

In the current model the perpendicular diffusion coefficient is taken to be proportional to $K_{||}$ following the ratio

[TABLE]

for both $r$ and $\theta$ $i$-coordinates(e.g., see Potgieter, 2000; Burger and Hattingh, 1998; Boschini et al., 2017, and references therein).  Palmer (1982) suggested that

[TABLE]

at Earth. The above description for $K_{\perp,i}$ is consistent at high rigidity with those from QLTs. Although this description was improved using more complicated approach (see, e.g., Engelbrecht and Burger, 2013), it remains one of the fundamental observational reference for transport theories (see discussion in Section 1.7.1 of Shalchi, 2009).  In the current version model, we use $\rho_{i}=0.06$ as discussed in Boschini et al. (2017). As discussed in Bobik et al. (2012), we used an enhanced $K_{\perp,\theta}$ in the polar regions; this enhancement is an implicit way of reducing drift effects by changing the CR intensity gradients significantly (Potgieter, 2013) in order to reproduce the amplitude and rigidity dependence of the latitudinal gradients of GCR differential intensities for protons (see Sect. 5.4 and, e.g., Vos and Potgieter, 2016).

Moving towards the solar maximum the rate of coronal mass ejections (CME) increases leading to a more chaotic structure of the solar magnetic field and to a stronger turbulence regime. In such conditions QLTs and other theories derived for a weak turbulence are no longer valid. In particular non linear effects are expected to be stronger for particles scattered at pitch angles around 90∘(Shalchi, 2009). As an example, simulations performed in strong turbulence condition showed a linear rigidity dependence of $K_{||}$ that extends to lower rigidity with respect to the value evaluated with QLTs (e.g., see Figure 3.5 and 6.5 of Shalchi, 2009). The effects of such a modification in the theories regard mainly rigidities lower than those considered in this work. Therefore, for sake of simplicity, the parametrization implemented in HelMod model was that expressed in Eq. (5), where $g_{low}$ decreases down to zero during high activity periods, while $g_{low}$ reaches its maximum value ($g_{low}=0.3$) for low solar activities111In current approach, $g_{low}$ was separately tuned for each set of observation. The resulting best values were then fitted to an empirical function, which, in turn, was tuned similarly to the other HelMod parameters, discussed in Sect. 3. A further discussion on $g_{low}$ values during the transition from low to high activity is available in Boschini et al. (2017).. For sake of completeness, as reported by, e.g., Guo and Florinski (2014), merged interaction region (MIR)222A MIR is the buildup of multiple interplanetary ejecta with enhanced solar wind speed, magnetic field, and plasma density can additionally reduce the diffusion and drift coefficients (le Roux and Fichtner, 1999), although not explicitly treated in this work, the time dependence of diffusion parameter includes such effect as it is linked to the real variation of GCR fluxes.

Finally, it has to be remarked that in the so-called force-field model (FFM) (e.g., see Gleeson and Axford 1968; Gleeson and Urch 1971, Section 2.1 of Bobik et al., 2012, also Section 8.2.4 of Leroy and Rancoita 2016 and references therein), the diffusion tensor is reduced to a scalar for spherically symmetric modulated number density of CR particles and steady-state modulation conditions. The FFM is an approximated way for treating solar modulation in many practical applications, but its intrinsic assumptions do not allow to account for relevant effects, like those related to the charge drift reported, for instance, from the observation of GCR modulation during periods with opposite field polarities of Sun (for instance, by Evenson and Meyer, 1984; Garcia-Munoz et al., 1986; Clem et al., 2000; Boella et al., 2001). Besides, evidences on how drift mechanisms can modify both the radial and (solar) latitude gradients were reported by Cummings et al. (1987, and reference therein); McKibben (1975, and reference therein); Simpson (1996, and reference therein); Heber et al. (1996, 2008, and reference therein); De Simone et al. (2011, and reference therein); Gieseler and Heber (2016, and reference therein). For a further discussion about limitations of the FFM to describe the modulation of cosmic rays, one can see, for instance, Caballero-Lopez and Moraal (2004). Nevertheless, it is worth to remark that the time variation of FFM modulation potential can be considered a good time variation proxy for the global behavior of heliosphere and, in turn, it can be exploited for determining the diffusion parameter (e.g., see Section 2.1 of Bobik et al., 2012, Section 8.2.4 of Leroy and Rancoita 2016 and Sect. 2.3).

2.3 The Diffusion Parameter

The diffusion parameter $K_{0}$ introduced in Eq. ((5)) (e.g., see Section 2.1 of Bobik et al. (2012)) is a scaling factor that defines the global behavior of particle flux modulation in the heliosphere and its dependence on time reflect the variability of interplanetary medium properties (like the actual solar magnetic field transported by SW and its turbulence) during the different phases of solar cycles (e.g., see Equation 4 in Manuel et al., 2014). $K_{0}$ was expressed in Section 2.1 of Bobik et al. 2012 and afterwards in Bobik et al. 2013b by means of a practical relationship with respect to the monthly smoothed sunspot numbers (SSN); in those papers, such a relationship was demonstrated to be adequate for the description of the dependence of the diffusion parameter on solar activity and polarity.

The current $K_{0}$ employed values are derived by means of the procedure discussed in Section 2.1 of Bobik et al. (2012) using the data from (Usoskin et al., 2011) and the SSN, updated to the most recent data series (World Data Center SILSO, 1964-2015; Clette et al., 2015); they are subdivided into four subsets, that is, ascending and descending phases for both negative and positive solar magnetic field polarities and are shown in Fig. 1. The updated practical relationship between $K_{0}$ in ${\rm AU^{2}GV^{-1}s^{-1}}$ and SSN values for $2.2\leq{\rm SSN}\leq 266.9$ for those periods (see Fig. 1) found is:

[TABLE]

with the parameters $c_{i}$ listed in Table 1. Furthermore, the root mean square (rms) values of the relative differences between the values obtained using Eq. (6) – following the procedure discussed in Bobik et al. (2012) – and current $K_{0}$ data are also reported in Table 1. It can be shown that the practical relationship (6) provides an overall agreement between calculated diffusion parameters, as function of SSN, and $K_{0}$ values.

This description is good enough to deal with CR modulation with HelMod code for periods of low solar activity. However, as soon as the high solar activity periods are considered – that is, in current approach, when tilt angles are $\gtrsim 50^{\circ}$ –, the rate and intensity of disturbances from Sun, such as CMEs and their short-term impacts, become increasingly larger, i.e., are resulting in a more chaotic structure with i) stronger perturbations of the magnetic field large-scale structure, which is hardly describable by means of a simple dipole approximation, and ii) local modulation effects which can be only related empirically to the actual SSN value. Therefore, a different or an additional solar activity indicator should be investigated. In fact, the monthly smoothed CR counting rates recorded by neutron monitors (NM) allows to better reproduce the short-time variation of $K_{0}$ that are needed to correctly calculate the CR modulation.  Among the available set of neutron monitor counting rate (NMCR) data series (Klein et al., 2009; Mavromichalaki et al., 2011), those from McMurdo NM (MCMU) in Antarctica (with effective vertical cutoff rigidity $P_{c}\simeq 0.3$ GV) and Oulu NM (Kananen et al., 1991) in Finland (with $P_{c}\simeq 0.8$ GV) were considered. These two mentioned NM stations have a low cutoff rigidity with a long enough period of data taking. The vertical geomagnetic cutoff rigidity is an adequate approximation of the rigidity lowest-limit of the primary spectrum to which the NM is sensitive333This is generally true for geographic position which have a geomagnetic cutoff rigidity below the atmospheric threshold, that for proton $\sim 1$ GV.. However, the maximum of NM sensitivity – i.e. the maximum of the response function (Clem and Dorman, 2000) – with respect to interactions of primary cosmic ray particles with the atmosphere, occurs in the rigidity interval 3–15 GV, depending on the NM site (e.g., see Artamonov et al., 2016, for a further discussion). Using monthly NMCR data, the new practical relationship for the high solar activity periods becomes:

[TABLE]

with the parameters $p_{i}$ listed in Tab. 2. For sake of completeness, we investigated also the usage of NMCR data for low activity periods. For such a purpose, a fit to a relationship similar to that provided by Eq. (7) was performed. It can be shown that the so-obtained diffusion parameters allow to derive modulated spectra similar to those in which $K_{0}$ is obtained by means of Eq. (6). Therefore, since the use of NMCR during low solar activity does not result in an appreciable difference, we kept the approach employing SSNs, as activity indicator for such periods.

It has to be noted that, although those two stations have different sensitivities and values of vertical rigidity cutoff, in a relative scale the relationship between NMCRs and the diffusion parameter is similar. This can be shown, for example, by computing the rms’s of the relative differences between values obtained using Eq. (7) and $K_{0}$ data. For the two NM stations considered the rms values found are shown in Tab. 2 for MCMU and OULU, respectively. The found results are compatible and, thus, only the relation (7) with MCMU data (shown in Fig. 2) was employed in this work, since it refers to the station with the lowest vertical cutoff.

2.4 The Interplanetary Magnetic Field

The GCR propagation in the heliosphere is affected by the outwards flowing SW with its embedded magnetic-field and magnetic-field irregularities. The Solar magnetic field is transported by the non-relativistic streaming particles of the SW, which carries the field into interplanetary space, producing the large-scale structure of the interplanetary (or heliospheric) magnetic field and the heliosphere geometry.  The so-generated and transported IMF is characterized by both the large scale structure (SW expansion from a rotating source) and low scale irregularities that change with time according to solar activity (e.g., due to variations of SW velocity or to plasma perturbations related to CMEs).

In the current version of HelMod we used IMF and drift as implemented in Bobik et al. (2013b). The heliosphere is divided into polar and equatorial regions, where different descriptions of IMF are applied. In the equatorial region, we used the Parker’s IMF ($\vec{B}_{Par}$) in the parametrization of Hattingh and Burger (1995), while in the polar regions we used a modified IMF ($\vec{B}_{Pol}$) that includes a latitudinal component, accounting for large scale fluctuations, dominant at high heliolatitudes, as suggested by Jokipii and Kota (1989):

[TABLE]

with

[TABLE]

In Eqs. (8, 9), $A$ is a coefficient that determines the IMF polarity and allows $|\vec{B}_{Par}|$ to be equal to $B_{\oplus}$, i.e., the value of the IMF at Earth’s orbit – as extracted from NASA/GSFC’s OMNI data set through OMNIWeb (King and Papitashvili, 2005); $\vec{e}_{r}$, $\vec{e}_{\theta}$ and $\vec{e}_{\varphi}$ are unit vector components in the radial, latitudinal and azimuthal directions, respectively; $\vartheta$ is the co-latitude (polar angle); $\vartheta^{\prime}$ is the polar angle determining the position of the heliospheric current sheet (HCS); $H$ is the Heaviside function: thus, the term, $[1-2H(\vartheta-\vartheta^{\prime})]$ allows $\vec{B}_{Par}$ to change sign in the two regions above and below the HCS (Jokipii and Thomas, 1981) of the heliosphere. Furthermore, $\Omega$ is the angular solar rotation speed and is assumed to be independent on the heliographic latitude and equal to the sidereal rotation at Sun’s equator. The separation between equatorial and polar regions were set to $\theta=30^{\circ}$ and $\theta=150^{\circ}$ according to Bobik et al. (2013b). Finally, in order to have a divergence free magnetic field, we require that the perturbation factor [$\delta(\vartheta)$] must be (Bobik et al., 2013b, and reference therein):

[TABLE]

where $\delta_{m}$ is the minimum perturbation factor of the field. The perturbation parameter is let to grow with decreasing co-latitude, while in their original work Jokipii and Kota (1989) fixed the value of $\delta$ between $10^{-3}$ and $3\times 10^{-3}$. In this work we use $\delta_{m}=2\times 10^{-5}$, tuned by comparison with Earth orbit observations during solar cycles 23-24 (Boschini et al., 2017). One has to note that the modified IMF introduced an additional magnetic component in the latitudinal direction that was not present in the pure Parker field. This should be considered when the diffusion tensor is generalized to heliocentric coordinates system using relationship reported e.g. in Burger et al. (2008, and reference there-in).

As already reported in the review by Owens and Forsyth (2013), while some basic descriptions – such as the Parker spiral (Parker, 1958) – are fully developed and already included into standard textbooks on space physics, other topics are still under development at the time of writing. A correct description of IMF is fundamental in a proper description of GCRs propagation (see also Raath et al., 2016, for a discussion about modified Parker’s magnetic field). As matter of fact, observations performed out-of-ecliptic by instruments on-board of Ulysses spacecraft (e.g., see Sanderson et al., 1995; Marsden, 2001; Balogh et al., 2001) showed the limit of the Parker field approach towards the polar regions and opened the way to a more complex description (e.g., see Jokipii and Kota, 1989; Smith and Bieber, 1991; Fisk, 1971; Hitge and Burger, 2010; Burger et al., 2008). Of great importance in the IMF description is the heliospheric current sheet (HCS) that divide the heliosphere in regions of inward or outward directed magnetic field lines (e.g., see section 9.1 of Solanki et al., 2006). The IMF is not a stable structure, but evolves with time following a 11-year cycle usually defined using sunspot numbers (Hathaway, 2015). At each cycle, the IMF reverses the magnetic polarity, usually defined by the sign of polar magnetic field strength (e.g., see Svalgaard et al., 1978), thus the Hale cycle or magnetic cycle of the IMF lasts $\approx 22$ years.

3 Effective Heliosphere Parameters and LIS’s

One of the success of HelMod model with respect to other available Monte Carlo codes for heliospheric propagation  – like, e.g., SolarProp–Gaggero et al., 2014 and HelioProp–Kappl, 2016)– is its reduced number of free parameters necessary for the description of modulation mechanisms. These parameters – as described later along the present section – are related to quantities determined from observations and need, in turn, to be tuned in order to obtain a comprehensive set of modulated spectra, like those discussed in Sect. 5. Free parameters tuned in this work are: a) the ratio between parallel and perpendicular diffusion coefficients and b) the $g_{low}$ parameter (both presented in Sect. 2.2). $\rho_{i}$ modifies the absolute scale of modulation intensity up to high rigidities; since this parameter also influences latitudinal gradients, its value is constrained in such a way that the obtained latitudinal gradients are in agreement with those found by Ulysses and presented in Sect. 5.4. The maximum value of $g_{low}$ tunes modulated differential intensity in the low rigidity range ($<$ 3 GV), mainly during low activity periods. In present approach other parameters, i.e., $\delta_{m}$ (Sect. 2.4) and drift suppression factor(Sect. 2.1) where not additionally tuned with respect to those previously estimated (e.g., see Bobik et al., 2013b), since their variation has a minor impact in the considered rigidity range for HelMod ($>$1 GV). Furthermore, the values of the diffusion parameter (described in Sect. 2.3) were evaluated using published data up to the end of solar cycle 23, and applied along the full solar cycle 24. However, one has to note that during solar cycle 22 for protons at low solar activity and positive polarity of the solar magnetic field, Ulysses probe has shown the presence latitudinal gradients with respect to the ecliptic plane. These gradients might affect the actual value of $K_{0}$ which is applied, mainly, as an overall diffusive scale factor for the whole heliosphere, although being evaluated on the ecliptic plane only.

Since any perturbation caused by Sun propagates into the heliosphere carried by SW, it takes typically $\sim 15$ Carrington rotations to reach the heliosphere boundary. In the meanwhile, the properties transported by SW change, e.g., new perturbations arise from Sun and propagate outwards, thus creating local areas, within the solar cavity, which have to be differently described. Therefore, in HelMod model the heliosphere is described as an effective heliosphere (Bobik et al., 2012) with a radius of 100 AU and is subdivided in 15 radially equally-spaced regions. Thus, in present bi-dimensional approximation, the model assumes that the propagation of GCR occurs within a static and spherical heliosphere with the TS located at 100 AU (for a further discussion see (Bobik et al., 2012)). Each i-th region experienced by CR particles is characterized by the heliospheric parameters evaluated at i-Carrington rotations back-in-time, corresponding to the time necessary to SW to reach it. Furthermore, as discussed in Della Torre et al. (2017a); Boschini et al. (2017), HelMod describes modifications of heliosphere dimensions as changes affecting only the effective distance between Earth and Sun. For example, in the re-scaled HelMod effective heliosphere if – because of the variation of SW ram pressure – TS is moved by 10 AU (e.g., see Washimi et al., 2011), Earth’s effective location is only shifted by 0.1 AU; Monte Carlo simulations (see Sect. 4) indicate that such variations have negligible effects on the modulated spectra when, instead, are determined at 1 AU (i.e., at Earth), thus for the purposes of the current work, we can properly assume that real and effective Earth locations are both at 1 AU from Sun.

It is worthwhile to remark that $\approx 100\,$AU is the average value of TS locations, which can be obtained from Table 2 of Whang et al. (2003). Furthermore (e.g., see Stone et al., 2005, 2008), Voyager 1 and 2 reached the TS in 2004 and 2007 at about 94.0 AU and 83.7 AU, respectively, from Sun, in agreement with the predictions from Whang and collaborators. Langner and Potgieter (2005) treated symmetric and asymmetric TS models and concluded that for $A>0$ cycles at solar minimum no significant difference occurs; for $A<0$ cycles variations remain negligible in nose direction while, approaching the tail direction, some differences can be appreciated at proton energies below 1–1.5 GeV. However, Langner and Potgieter (2005) and Potgieter (2008) suggested that, in general, a symmetric TS with a radial distance of $\approx 100\,$AU is still a reasonable assumption. In addition, it should be noted that, in literature, the heliospheric structure is considered latitudinally asymmetric (particularly) during solar minimum conditions mostly because the SW speed depends on the latitude and solar activity (e.g., see McComas et al., 1998). A more complete asymmetrical structure of the heliosphere can be evaluated using magneto-hydrodynamic (MHD) models (e.g., see Florinski and Pogorelov, 2009; Guo and Florinski, 2014) that include transport description in outer heliosphere, heliosheath and co-rotating interaction regions. Moreover, Voyager 1, 2 observations point to a dynamic TS that is moving inward/outward in the heliosphere (Stone et al., 2005; Richardson and Wang, 2011), while numerical models indicate that this TS movement could be as large as $\sim$20 AU over a complete solar cycle (see the discussion in Manuel et al., 2015, and references therein). Finally, even though variations of the real size of the heliosphere may be important for the analysis of CR propagation near the TS, we do not consider them in this work.

The diffusion parameter $K_{0}$ is determined (as discussed in Sect. 2.3) using the values of modulation strength, SSN444http://www.sidc.oma.be/sunspot-data/ values, NMCR555http://www.nmdb.eu/nest values at McMurdo station and the effective heliosphere radius (see discussion above). Other parameters (which depend on the solar activity) are the tilt angle $\alpha_{\rm t}$ of HCS, the magnetic field polarity [related to the sign of the coefficient $A$ in Eq. (9)], the magnetic field amplitude ($B_{\oplus}$) and, finally, the solar wind velocity ($V_{\rm sw}$). The latter two parameters are measured at Earth’s orbit and provided by OMNIWeb666http://omniweb.gsfc.nasa.gov/. The polarity of the magnetic field (obtained from Wilcox Solar Observatory Polar Field Observations777http://wso.stanford.edu/) and $B_{\oplus}$ determine the IMF described by means of Eq. (8). In the current model we use the so-called “line-of-sight” model (L-model) for the tilt angle $\alpha_{\rm t}$ of HCS (Hoeksema, 1995), that provides an overall general agreement with experimental data (e.g., see discussion in Bobik et al., 2012, 2013b; Della Torre et al., 2017a; Boschini et al., 2017). $\alpha_{\rm t}$ and the field polarity, obtained from Wilcox Solar Observatory, are used to deal with the drift velocity (as discussed in Sect. 2.4), which contributes to the overall convection velocity in Eq. (1). Drift contribution is relevant during low solar activity – e.g., for $\alpha_{\rm t}<30^{\circ}$ – and decreases with increasing solar activity. Since during the high activity period the heliospheric magnetic field is far from being considered regular, in analogy with other works (see e.g. Potgieter, 2008), we introduced a correction factor that suppress any drift velocity at solar maximum. Finally, the latitudinal dependence of the SW speed is the one described in (Bobik et al., 2012).

The local interstellar spectra (LIS) are input cosmic ray intensities for any modulation models. Fluxes are assumed isotropically distributed at heliosphere boundary, in a steady-state configuration. Recently, Della Torre et al. (2017a); Boschini et al. (2017) deduced LIS’s for protons, helium and antiprotons using the most recent experimental results combined with the state-of-the-art models (i.e., GALPROP888http://galprop.stanford.edu and HelMod) for propagation in galaxy and heliosphere. In fact, HelMod and GALPROP (Moskalenko and Strong, 1998; Strong and Moskalenko, 1998) were combined to provide a single framework and run to reproduce a comprehensive set of observations of CR species collected in different time periods, from 1997 up to 2015. The authors proposed an analytical expression for proton and helium nuclei (Della Torre et al., 2017a; Boschini et al., 2017) LIS’s as a function of the rigidity expressed in unit of [m2 s sr GV]-1 (see also Fig. 3):

[TABLE]

where $a_{i},b,c,d_{i},e_{i},f_{i},g$ are the numerical coefficients summarized in Table 3 (from Table 5 of Boschini et al. (2017)).

4 The Monte Carlo Code

For most applications Parker’s transport equation [Eqs. 1 and 3] has been solved using numerical methods, because its intrinsic complexity.

The traditional approach to solve multi-dimensional partial differential equations makes use of numerical integration methods such as the finite difference technique (e.g., see Jokipii and Kopriva, 1979; Kota and Jokipii, 1983; Potgieter and Moraal, 1985; Burger and Hattingh, 1995) or as the standard implicit difference technique (e.g., see Fisk, 1971; Kota and Jokipii, 1991). These methods have several disadvantages, mainly numerical instability problems when solving differential equation in higher dimensions (Pei et al., 2010; Kopp et al., 2012). A modern approach that in recent years has been used more frequently to solve numerically a variety of problems in space-physics is based on Monte Carlo methods (e.g., see Kruells and Achterberg, 1994; Fichtner et al., 1996; Yamada et al., 1998; Gervasi et al., 1999; Zhang, 1999; Alanko-Huotari et al., 2007; Pei et al., 2010; Strauss et al., 2011; Bobik et al., 2012; Gaggero et al., 2014; Kappl, 2016). A diffusion process described by a FPE can be written as well with a set of stochastic differential equations (SDE) (e.g., see Chapter 1.6-1.7 of Klöden and Platen, 1999). As reported in Bobik et al. (2016, and reference therein) this approach allows to get more flexibility in model implementation, more stability of numerical results and more possibility to explore physical results that are hard to handle with “classical” numerical methods. With the stochastic approach, the solution can be evaluated computing the SDEs both “forward-in-time” or “backward-in-time”. A comparison between the two approaches with a numerical estimation of systematic uncertainties can be found in Bobik et al. (2016). In “forward-in-time” approach quasi-particle objects were traced from the heliosphere boundary down to the inner part of heliosphere. In “backward-in-time” approach the numerical process starts from the target and trace-back quasi-particle objects till the heliosphere boundary. The “backward-in-time” method is widely used, due to faster evaluation of spectra at single points inside the heliosphere and is presented in this section.

A stochastic motion does not allow for a single particle study, but it is only possible to explore how the system evolves in average, considering all the particles as an ensemble. When the particles stochastic behavior is studied from probabilistic point of view, the theory of Markov stochastic process may provide very powerful mathematical tools (Zhang, 1999). As pointed out by Pei et al. (2010), it is interesting to note that, since each random process have to be independent of all others, one major advantage of the stochastic method is that it is very easy to parallelize the computation. Therefore, it can run in the same time on many CPUs of a local cluster reducing significantly the computation time and the hardware costs.

The equivalent set of SDEs, for the 2D approximation, is derived from the 3-D transport equations (reported in A) by integrating over the azimuthal component:

[TABLE]

where $K_{ij}$ is the symmetric part of diffusion tensor in heliocentric spherical coordinates, $v_{drift,i}$ follows the definition of Eq. (2), $V_{\textrm{sw}}$ is the magnitude of solar wind speed and $\omega_{i}$ is a random number following a Gaussian distribution with zero average and unit standard deviation. For stochastic integration with $L\neq 0$ the diffusion/convection process alone is not able to proper describe the stochastic evolution. Furthermore, $L$ has the meaning of an additional sources and losses term in FPE (see A) distributed inside the phase-space, such that it contributes to the solution with exponential corrections integrated along the stochastic path (Bobik et al., 2016).

The vector $\vec{q}=[r,\theta,T]$ represents a phase space density element, that in literature is usually labeled as a pseudo-particle, which time evolution is simulated by means of Eq. (11) from the inner part heliosphere up to the outer boundary. As presented in Bobik et al. (2016, and reference therein) the differential intensity $J$ can be obtained from the density of pseudo-particles by averaging over many realizations of the SDEs.

The procedure used to integrate the SDEs backward-in-time and to evaluate the solution at Earth (at $r=1$ AU and $\theta=\frac{\pi}{2}$) is described in Section 4.1.2 of Bobik et al. (2016, and reference therein), i.e.,

(a) pseudo-particles are generated at Earth with initial Energy $T_{i}$;

(b) each event is integrated over the time evolution of a pseudo-particle following the path described in Eq. (11), the integrated “loss” term ($L_{Int}$) is incremented by the quantity $L\Delta t$;

(c) when a pseudo-particle reaches the outer border of the effective heliosphere located at 100 AU ($r_{b}$) the value of differential intensity at boundary for the reconstructed energy $T_{r}$ is saved and, then, weighted for the “loss” factor $\exp{(-L_{Int})}$;

(d) the modulated differential intensity at Earth for initial energy $T_{i}$ is evaluated over $N$ realizations by mean of Equation 22 in Bobik et al. (2016):

[TABLE]

Equations (11, 12) allow one to evaluate modulated spectrum as solution of Eq. (1). The differential intensity evaluated using Eq. (3) follows the same procedure here described taking into account the relationship reported in Eq. (4). In this case the pseudo particle evolves in a phase-space that includes particle momentum instead of kinetic energy per nucleon (see A for details on how to derive 2-D SDE equations from generic 3-D solution).

An alternative approach is to use Monte Carlo integration to evaluate the normalized probability function $G(P_{0}|P)$ that gives a probability for a particle to be observed at Earth with a rigidity $P_{0}$ having a rigidity $P$ at the heliospheric boundary. Once $G(P_{0}|P)$ is evaluated it is possible to obtain the modulated spectrum directly from an arbitrary $J_{\rm LIS}$ provided, e.g., by GALPROP.  The modulated spectrum at specific rigidity $P_{0}$ is, then, obtained by (e.g., see Pei et al., 2010; Della Torre et al., 2017b; Boschini et al., 2017):

[TABLE]

This approach allows one to reduce the amount of simulations, when testing several LIS’s. Actually, in previous approach, a new LIS meant Monte Carlo realizations to be re-run for the same heliosphere parameters, while the latter approach allows one to use the same simulations for all LIS’s under test. We implemented an on-line calculator999http://www.helmod.org which, using a python script, reads the GALPROP outputs and provides the modulated spectra for periods of selected experiments for comparison with published data. The calculation of heliospheric propagation is replaced by the integration of Eq. (13) using the normalized probability functions, which are pre-evaluated by the HelMod code as described in the previous section. This method dramatically accelerates the modulation calculations while provides the same accuracy of the full-scale simulation.

5 Comparison with Observations During Solar Cycles 23-24

The current HelMod model provided modulated differential intensity for protons, helium nuclei and antiproton for low and high solar activities (as discussed in Sect. 1). In this article, we focus on HelMod results regarding protons and helium, whose LIS’s were recently investigated in Della Torre et al. (2017a); Boschini et al. (2017). The current parametrization is also suited to reproduce the high energy behavior of the measured spectra (e.g., see Figures 3 and 5 in Della Torre et al. 2017a regarding data up to 1 TV).

5.1 Low Solar Activity

HelMod was initially developed for periods of low solar activity, when the contribution of charge dependent transport processes – resulting in a magnetic drift convection, i.e., that expressed by $\vec{v}_{drift}$ in Eq. (1) – is so important enough to introduce peculiar features, for instance, in the time dependence observed for the positron fraction (e.g., see Della Torre et al., 2012) at rigidities lower than (10-20) GV and heliospheric latitudinal gradients of galactic cosmic rays distribution (e.g., see Bobik et al., 2013b).

The large-scale structure of the IMF is strongly affected by solar activity. More the IMF assumes a regular structure, more GCR particles experience the effects of magnetic drift transport. As already mentioned, the relevance of magnetic drift during such periods was widely recognized in literature (see e.g. Jokipii et al., 1977; Jokipii and Kopriva, 1979; Potgieter and Moraal, 1985; Boella et al., 2001; Strauss et al., 2011; Della Torre et al., 2012; Bobik et al., 2013b, a).

The latest low solar activity period was investigated, in particular, using data taken by PAMELA (e.g., see Adriani et al., 2013). Previously, AMS-01 mission (June 1998, Aguilar et al., 2002) on the space shuttle and BESS (Shikaze et al., 2007; Abe et al., 2016) on board of stratospheric balloons, sample few short time periods, during solar cycle 23. In Fig. 4, we show the comparison between experimental data from PAMELA and modulated spectra from HelMod for both protons and helium nuclei. It has to be remarked that, at rigidities lower than 1.5 GV, the observed helium spectrum is slightly lower than that expected by the modulated one.

Furthermore, in Fig. 5, we show the time evolution of a few single rigidity bins observed by PAMELA from 2006 to 2010 during the solar minimum (see also Potgieter et al., 2014; Vos and Potgieter, 2016, and reference therein). The HelMod modulated intensity evaluated for the same period and rigidity is superimposed with solid line. As described in Sect. 2.3, in HelMod the diffusion coefficient scales along the time using a practical relationship with smoothed sunspot numbers. The good agreement observed in Fig. 5 confirms that, for quiet periods, even if sporadic solar events can perturb the interplanetary particle transport, the latter is mostly regulated by its average properties.

5.2 High Solar Activity

The description of the transport properties during high solar activities is more challenging in comparison with those at low. The higher rate of solar events makes the interplanetary medium properties more complex and, in turn, the IMF far from being modeled with a regular shape.

The lack of systematic observations of GCR fluxes in high solar activity makes AMS-02 the first detector, which is able to obtain unprecedentedly precise and continuous measurements of GCRs, under such solar conditions. Previously, only BESS balloons (Shikaze et al., 2007) provided proton spectra for short-time periods during the peak of solar cycle 23. On other hands, AMS-02 already provided an unique data-set integrated over 3 year (Aguilar et al., 2015b, a) of data taking, during the solar activity peak of solar cycle 24.

Recently (see, Della Torre et al., 2017b, a; Boschini et al., 2017) HelMod model extended his results to reproduce AMS-02 flux of protons, helium nuclei and antiprotons integrated from 2011 to 2015, allowing one to begin unveiling the solar maximum period with an unprecedented detailed treatment. As discussed in Sect. 2.2, HelMod propagation model for high activity periods needed some additional refinement. First of all, the higher rate of solar particle emission makes hard to describe the interplanetary medium by means of its average properties only. In order to overcome this difficulty, NMCR was used instead of SSN as a proxy for scaling diffusion parameter’s time variation (see discussion in Sect. 2.2). Moreover, as discussed in Sect. 3, we introduced an additional correction factor that suppresses drift velocity at the solar maximum. These improvements allowed HelMod to reproduce the average proton and helium nuclei fluxes measured by AMS-02 (see Fig. 6) and by BESS (Della Torre et al., 2017a; Boschini et al., 2017), within the experimental error bars and the simulation uncertainties.

So far, it is important to remark that using AMS-02 data we could explore the solar maximum as a whole. When AMS-02 time dependent data will be available a deeper understanding of the particle transport in solar maximum conditions will be possible.

5.3 Helium over Proton Ratio

The helium over proton ratio can provide useful insights for understanding the propagation of cosmic rays.

In Fig. 7, we made a prediction for helium over proton ratios at different rigidities, as they should be observed by AMS-02 detector from June 2011 to May 2016. From an inspection of Fig. 7, one can observe that the ratio is almost constant with time for all considered rigidities. Helium has double charge and approximately four time the mass of protons, thus processes described in Eq. (1) produce a different level of solar modulation for similar kinetic energy per nucleon. In fact, diffusion process [see Eq. (5)], magnetic drift [see Eq. (2)] and adiabatic energy loss (latter term in Eq. 3) are naturally expressed in term of particle rigidity ($P=\frac{pc}{Ze}$). As a consequence, solar modulation acts in the same way for particle with same rigidity and charge sign. Furthermore, this should occur in spite of the large intensity variations as function of time for the proton fluxes – in particular at rigidities lower than about 5–6 GV –, observed by AMS-02 (e.g., see Della Torre and AMS-02 Collaboration, 2016). Thus, particle rigidity is the natural quantity for studying GCR propagation.

Finally, it is important to remark that, in case of particles with opposite charge sign, processes like the magnetic drift transport result in differently propagating positive charged and negative charged particles. Therefore, a different time behavior is expected to be observed and should be accounted for by accurate propagation models (e.g., see Della Torre et al., 2012; Adriani et al., 2016).

5.4 Dependence on Heliospheric Latitude

A unique set of measurements of heliosphere properties out of the ecliptic plane were performed from 1990 to 2009 by the Ulysses mission (see Heber, 2011 or Heber and Potgieter, 2006 for recent reviews). The COSPIN suite of instruments on board of such a spacecraft – described by Simpson et al. (1992) – includes the Kiel electron telescope (KET), designed to measure intensities and energy spectra of energetic particles separating electrons, protons and helium nuclei. The instrument covered the energy range from $\sim$5 MeV/nucleon to above 2 GeV/nucleon for protons and helium nuclei, and from $\sim 3$ MeV to above 300 MeV for electrons. Ulysses orbited around Sun reaching a maximum solar latitude of about $80^{\circ}$ in both the northern and southern hemispheres, at solar distances ranging from $\sim 1$ AU to $\sim 5$ AU and with an orbital period of approximately 5.5 years. Ulysses performed three “fast latitude scans” (FLS) of Sun. The first FLS in 1994/1995 took place near solar minimum with positive polarity ($A>0$). A non-symmetric galactic proton intensity with respect to the heliographic equator was observed. The minimal intensity for protons with energies $>0.1$ GeV was observed to be displaced by about 7–10 degrees towards the southern hemisphere (McKibben et al., 1996; Simpson, 1996; Heber et al., 1996). Moreover latitudinal gradients of  0.3%/degree for protons at low energies ($<0.1\,$GeV) and 0.22%/degree at higher energies ($>2\,$GeV) were observed.

The second FLS in 2000/2001 took place close to solar maximum with negative polarity ($\rm A<0$) and no latitudinal gradient was observed for any GCR species, indicating that drift effects are negligible at solar maximum (McKibben et al., 2003).

The third FLS in 2007/2008 was performed close to solar minimum, as the first FLS, but with reversed polarity ($A<0$). Observations made by PAMELA starting from 2006 (Casolino et al. 2008) allowed a comparative determination of the radial and latitudinal gradient during an $A<0$ solar magnetic epoch. For protons in the rigidity interval 1.6–1.8 GV the measured latitudinal gradient found was ($-0.024\pm 0.005$)%/degree and the radial gradient ($2.7\pm 0.2$)%/AU (De Simone et al., 2011). These measurements, performed during low solar activity, allowed a deeper insight into drift mechanisms and the structure of the heliospheric magnetic field in the polar regions. In this way, it was possible to test IMF models like the one proposed by Fisk (1996).

Detailed studies with HelMod model using Ulysses results explored the solar modulation outside the ecliptic plane. In Bobik et al. (2013b), it was proven that, using a propagation model as the one described in Sect. 2, HelMod code is able to reproduce qualitatively and quantitatively the latitudinal profile of the GCR intensity, and the latitudinal dip shift with respect to the ecliptic plane as observed in the inner part of heliosphere by the Ulysses spacecraft during the first FLS. As an illustrative example, in Fig. 8 we show the comparison between the measured helium normalized counting rate in the energy range (0.250-2.1) GeV and the modulated spectrum calculated using HelMod for 1 GeV. Both experimental data and simulations are normalized to the mean value to allow a relative comparison along the solar cycle. It is important to remark that, the aim of Fig. 8 is to show the qualitative agreement found between the HelMod spectra and observation data; in fact, HelMod calculations were performed for a mono-energetic bin, while KET observations are integrated over a large energy interval (e.g., see the discussion in De Simone et al., 2011). A more quantitative comparison with Ulysses data might need to combine together simulations for several energy bins weighted with Ulysses response function. The first FLS in 1994/1995 and the third in 2007/2008 are almost fully reproduced. Besides, the second FLS in 2000/2001 during the solar maximum still shows discrepancies, which deserve to be better understood.

6 Conclusions

The CR propagation inside the heliosphere was initially treated – almost 60 years ago – by Parker (1965) who provided a general theoretical framework for the heliospheric modulation through the interplanetary medium. Since then, continuous advances allowed a deeper and deeper understanding of the general properties of the IMF affecting particle motion. However, the description of the transport mechanisms occurring still needs further refinements for allowing a comprehensive treatment along the complete solar cycle. Moreover, since the heliosphere cannot be considered a steady-state environment, a model should be able to reproduce in a unique description both low and high activity periods. The idea behind HelMod is, in fact, that since all GCR species propagating in the heliosphere experience the same mechanisms, then the modulated spectra can be derived using the same time-dependent heliospheric parameters, i.e., the entire set of modulated observed spectra has to be reproduced with those parameters from the corresponding LIS’s.

In the present work, we have treated the heliospheric parameters of HelMod model and how they were implemented in the code for the version 3, as well the Monte Carlo integration technique with a full description of the SDEs, which finally allow the numerical solution of the FPE proposed by Parker (1965).

A relevant step forward in this field is currently due to the availability of data sets, regarding the time evolution of modulated CR spectra – provided by PAMELA and AMS-02 spectrometers – from the final part of the solar cycle 23 up to present. The numerical approach of HelMod model version 3 – discussed in this work and in the recent article by Della Torre et al. (2017b, a); Boschini et al. (2017) – is able to satisfactorily reproduce the full set of observations obtained, during the two latest solar cycles, for instance by BESS, AMS-01, PAMELA and AMS-02. It has to be pointed out that the unprecedented accuracy provided by AMS-02 proton data with experimental error down to 1–2% represents a challenge for a modulation model. The agreement found between proton and helium nuclei modulated spectra from AMS-02 and those calculated from HelMod is within the AMS-02 errors.

Although diffusion is the dominant modulation process for GCRs, it was shown how charge sign dependent processes, i.e., those resulting in a drift velocity convection, are fundamental elements for dealing with a complete solution for the particle transport through the heliosphere. Moreover, it was remarked that when the treatment accounts only for spectra observed at Earth, it can lead not to appropriately describe the complexity of the modulation phenomenon. Probes, like Ulysses, allow one to better investigate how to deal with a bi-dimensionally modeled heliosphere, i.e., the modulation effects occurring also out of the ecliptic plane as a function of solar latitude. For instance, the measured helium normalized counting rate from Ulysses and the modulated spectrum calculated using HelMod were found to agree, thus indicating the capabilities of the present model to investigate the CR transport up to large solar latitudes and distances from Sun up to 5 AU.

Finally, a proper description of heliospheric modulation, like that HelMod can provide, relates not only computed LIS’s to observed modulated differential intensities, but also may result in constraining the parameters of galactic production models, like GALPROP (as discussed in Della Torre et al., 2017a; Boschini et al., 2017).

Acknowledgements

This work is supported by ASI (Agenzia Spaziale Italiana) under contract ASI-INFN I/002/13/0 and ESA (European Space Agency) contract 4000116146/16/NL/HK.

We acknowledge the NMDB database (www.nmdb.eu), supported under the European Union’s FP7 programme (contract no. 213007) for providing data. The data from McMurdo were provided by the University of Delaware with support from the U.S. National Science Foundation under grant ANT-0739620. Finally, we acknowledge the use of NASA/GSFC’s Space Physics Data Facility’s OMNIWeb service, and OMNI data.

We wish to specially thank Pavol Bobik, Giuliano Boella, Davide Grandi, Karel Kudela, Simonetta Pensotti, Marian Putis, Davide Rozza, Mauro Tacconi and Mario Zannoni for their support to HelMod code and suggestions.

Appendix A Diffusion Tensor and Stochastic Differential Equations

One of the key point to provide a proper stochastic integration is to evaluate the correct set of SDEs. Although the procedure is relatively simple and straightforward, there are few points that are important to be remarked. In this section, we evaluate SDEs from Parker’s equation in the form of particle density for unit space and energy and omni-directional distribution function. More details on the procedure can be found in Zhang (1999); Pei et al. (2010); Kopp et al. (2012); Bobik et al. (2016) and reference therein.

A.1 SDE in Kinetic Energy

Equation (1), that controls cosmic rays modulation, can be always written by a multi-dimensional diffusion equation in the form of backward-in-time Kolmogorov equation (Zhang, 1999):

[TABLE]

The corresponding SDE set is, then, obtained as:

[TABLE]

It is important to note that Eq. (14) includes an additional parameter $L$, that is not directly taken into account for the evaluation of the stochastic path. Indeed, it represents an additional process by allowing the stochastic realization to be created at an exponential rate of $L$ as function of time (see Eq. 12). Equation (1), in a spherical heliocentric coordinate system, can be rewritten as:

[TABLE]

where the divergence of drift velocity is zero by definition:

[TABLE]

From a comparison between Eq. (14) and Eq. (A.1), it is possible to obtain:

[TABLE]

The matrix $B$ can be derived using the Cholesky–Banachiewicz or Cholesky–Crout algorithms. As pointed out in Appendix B of Pei et al. (2010), the so obtained matrix is not unique, but its solutions are stochastically equivalent (see e.g. Kopp et al., 2012).

A.2 FPE in momentum

To obtain SDE from Eq. (3), the latter should be rewritten in the form of Eq. (14) as described in A.1:

[TABLE]

One should note that in this case a) the linear term $L$ is equal to zero, thus the exponential weight in Eq. (12) can be neglected (as done, e.g., in Strauss et al., 2011), b) both spatial convection and diffusion matrix are the same as for SDE in kinetic energy, and c) the difference between SDE in kinetic energy and momentum is only in the energy loss terms: $A^{T}_{B}$, $A^{p}_{B}$ and $L$.