Numerical modeling of cosmic-ray transport in the heliosphere and interpretation of the proton-to-helium ratio in Solar Cycle 24

TL;DR

This paper develops a numerical model to understand cosmic-ray transport in the heliosphere, explaining the observed proton-to-helium ratio variations during Solar Cycle 24 through diffusion properties.

Contribution

It introduces a detailed numerical model linking cosmic-ray diffusion characteristics to the temporal evolution of the proton-to-helium ratio, incorporating analytical solutions.

Findings

Long-term proton-to-helium ratio behavior explained by diffusion dependence.

Model successfully reproduces observed temporal variations.

Highlights importance of mass/charge dependent diffusion in CR transport.

Abstract

Thanks to space-borne experiments of cosmic-ray (CR) detection, such as the AMS and PAMELA missions in low-Earth orbit, or the Voyager-1 spacecraft in the interstellar space, a large collection of multi-channel and time-resolved CR data has become available. Recently, the AMS experiment has released new precision data, on the proton and helium fluxes in CRs, measured on monthly basis during its first six years of mission. The AMS data reveal a remarkable long-term behavior in the temporal evolution of the proton-to-helium ratio at rigidity $R = p/Z <$ 3 GV. As we have argued in a recent work, such a behavior may reflect the transport properties of low-rigidity CRs in the inteplanetary space. In particular, it can be caused by mass/charge dependence of the CR diffusion coefficient. In this paper, we present our developments in the numerical modeling of CR transport in the Milky Way and in the heliosphere. Within our model, and with the help of approximated analytical solutions, we describe in details the relations between the properties of CR diffusion and the time-dependent evolution of the proton-to-helium ratio.

Equations

$${\partial }_{t}N=S+\vec{\nabla }\cdot (D\vec{\nabla }N)-N\Gamma +{\partial }_E(\dot{E}N),$$

$$\frac{\partial \psi }{\partial t}=\nabla \cdot \left(K\cdot \nabla \psi \right)-V\cdot \nabla \psi +\frac{1}{3}\nabla V\frac{\partial \psi }{\partial \ln \mathcal{R}}+Q$$

$$\frac{\partial \psi }{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^{2}K\frac{\partial \psi }{\partial r}\right)-V\frac{\partial \psi }{\partial r}+\frac{1}{3r^2}\frac{\partial }{\partial \mathcal{R}}{r}^{2}V\frac{\partial \psi }{\partial \ln \mathcal{R}}+Q$$

$$K\frac{\partial {}^2\psi }{\partial r^2}+\left(\frac{\partial K}{\partial r}+\frac{2K}{r}-V\right)\frac{\partial \psi }{\partial r}+\left(\frac{2V}{3r}+\frac{1}{3}\frac{\partial V}{\partial r}\right)\frac{\partial \psi }{\partial \ln \mathcal{R}}=0$$

$$\frac{\partial \psi }{\partial r}=\frac{{\psi }_{i+1}^{j+1}-{\psi }_{i-1}^{j+1}}{4{\Delta }_r}+\frac{{\psi }_{i+1}^j-{\psi }_{i-1}^j}{4{\Delta }_r}$$

$$\frac{{\partial }^2\psi }{\partial r^2}=\frac{{\psi }_{i+1}^{j+1}-2{\psi }_i^{j+1}+{\psi }_{i-1}^{j+1}}{2{\Delta }_r^2}+\frac{{\psi }_{i+1}^j-2{\psi }_i^j+{\psi }_{i-1}^j}{2{\Delta }_r^2}$$

$$\frac{\partial \psi }{\partial \ln \mathcal{R}}=\frac{{\psi }_i^{j+1}-{\psi }_i^j}{{\Delta }_{\ln \mathcal{R}}}$$

$${\psi }_{i-1}^j\left(\frac{{\phi }_1}{2{\Delta }_r^2}-\frac{{\phi }_2}{4{\Delta }_r}\right)+{\psi }_i^j\left(-\frac{{\phi }_1}{{\Delta }_r^2}-\frac{{\phi }_3}{{\Delta }_{\ln \mathcal{R}}}\right)+$$

$${\psi }_{i+1}^j\left(\frac{{\phi }_1}{2{\Delta }_r^2}+\frac{{\phi }_2}{4{\Delta }_r}\right)={\psi }_{i-1}^{j+1}\left(-\frac{{\phi }_1}{2{\Delta }_r^2}+\frac{{\phi }_2}{4{\Delta }_r}\right)+$$

$${\psi }_i^{j+1}\left(\frac{{\phi }_1}{{\Delta }_r^2}-\frac{{\phi }_3}{{\Delta }_{\ln \mathcal{R}}}\right)+{\psi }_{i+1}^{j+1}\left(-\frac{{\phi }_1}{2{\Delta }_r^2}-\frac{{\phi }_2}{4{\Delta }_r}\right)$$

$$\left[\begin{array}{ccccc}{\eta }_2& {\eta }_3& 0& \dots & 0\\ {\eta }_1& {\eta }_2& {\eta }_3& 0& \dots \\ 0& {\eta }_1& {\eta }_2& {\eta }_3& \dots \\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& \dots & \dots & {\eta }_2& {\eta }_3\end{array}\right]\left[\begin{array}{c}{\psi }_2^j\\ {\psi }_3^j\\ :\\ ⋮\\ {\psi }_{N-1}^j\end{array}\right]=\left[\begin{array}{c}{\eta }_4\\ {\eta }_4\\ :\\ ⋮\\ {\eta }_4-{\eta }_3{\psi }_N^j\end{array}\right]$$

$${\chi }^2(k_0)=\sum _j{\sigma }_j^{-2}{\left[{\hat{J}}_{\mathrm{p}}^j-\frac{1}{\Delta {\mathcal{R}}_j}{\int }_{{\mathcal{R}}_j}^{{\mathcal{R}}_{j+1}}{J}_{\mathrm{p}}(\mathcal{R})d\mathcal{R}\right]}^2$$

$$\beta (\mathcal{R})=\mathcal{R}/\sqrt{{\mathcal{R}}^2+(m_{p}A/Z{)}^2}$$

$$\partial \psi /\partial t=K\partial \psi /\partial r-V\psi$$

$$\mathcal{M}={\int }_r^d\frac{V}{K}dr\approx \frac{Vd}{K}\approx \frac{Vd}{K_0\beta \mathcal{R}/\mathrm{GV}}$$

$$\text{p/He}\approx \left(\frac{J_{\mathrm{p}}^0}{J_{\mathrm{He}}^0}\right)\frac{e^{-Vd/K_{\mathrm{p}}(\mathcal{R})}}{e^{-Vd/K_{\mathrm{He}}(\mathcal{R})}}$$

$$\text{p/He}\approx \frac{J_{\mathrm{p}}^0}{J_{\mathrm{He}}^0}\left\{1-\frac{\mu (t)}{\mathcal{R}/\mathrm{GV}}\left[\frac{1}{{\beta }_{\mathrm{p}}(\mathcal{R})}-\frac{1}{{\beta }_{\mathrm{He}}(\mathcal{R})}\right]\right\},$$

$${\Gamma }_{\text{p/He}}(\mathcal{R})\approx \frac{1-\frac{\mu (t_2)}{\mathcal{R}/\mathrm{GV}}\left[\frac{1}{{\beta }_{\mathrm{p}}(\mathcal{R})}-\frac{1}{{\beta }_{\mathrm{He}}(\mathcal{R})}\right]}{1-\frac{\mu (t_1)}{\mathcal{R}/\mathrm{GV}}\left[\frac{1}{{\beta }_{\mathrm{p}}(\mathcal{R})}-\frac{1}{{\beta }_{\mathrm{He}}(\mathcal{R})}\right]}$$

$$r(t,\mathcal{R})=\{\begin{array}{ll}A(\mathcal{R})& (t<t_0(\mathcal{R}))\\ A(\mathcal{R})+B(\mathcal{R})\times (t-t_0(\mathcal{R}))& (t\ge t_0(\mathcal{R}))\end{array}$$

Full text

††journal: Advances in Space Research

1 Introduction

The vast majority of cosmic rays (CRs) observed near-Earth originate well outside the solar system, from spectacular astrophysical phenomena such as supernova remnants or stellar winds. When entering the heliosphere, these particles are decelerated by the solar wind and deflected by the interplanetary magnetic fields. This makes the energy spectra of CRs in the inner heliosphere significantly different from those in the interstellar space. Moreover, the effect is known to change periodically, with the changing number of sunspots appearing on the solar corona.

The observed change of the Galactic CR flux over the variable Solar Cycle is known as solar modulation of CRs in the heliosphere. Along with its implications in solar or plasma astrophysics, solar modulation is an important phenomenon for Galactic CR physics studies, as it limits our ability to identify the CR sources or to search for dark matter annihilation signals (Grenier et al., 2015). Moreover, predicting the variations of the CR flux in the heliosphere represents a significant challenge for interplanetary space missions (Kudela et al., 2000; Tomassetti et al., 2017).

The solar modulation effect can be studied using various types of CR data such as counting rates from neutron monitors or direct measurements of CR fluxes (Usoskin et al., 2005, 2011; Bindi, 2017; Wiedenbeck et al., 2013). In the past decades, precious information have been gained thanks to space missions CRIS/ACE (Wiedenbeck et al., 2009), IMP-7/8 (Garcia-Munoz et al., 1997), Ulysses (Heber et al., 2009), and Voyager-1 (Cummings et al., 2016). The recent data from EPHIN/SOHO (Kühl et al., 2016), PAMELA (Adriani et al., 2013; Martucci et al., 2018), and AMS (Aguilar et al., 2018a, b) are giving a deep characterization of effect in the current Solar Cycle 24 (Bindi, 2017). The solar modulation effect is known to be caused by a combination of basic processes such as magnetic diffusion, drift motion, convection, and adiabatic energy losses (Parker, 1965; Krymsky, 1964; Fisk, 1971; Jokipii, 1967). All these processes are linked to the nature of the heliospheric medium and on the dynamics of charged particles in magnetic fields. The diffusive propagation arises from the erratic motion of CRs through magnetic turbulence, i.e., charged particle scattering off the small-scale irregularities of the magnetic fields. In the heliosphere, these irregularities are frozen in the expanding solar wind, and continuously transported radially towards the outer heliosphere. This gives rise to a global convective motion, and adiabatic energy losses associated with it. On the other hand, the large-scale global structure of the heliospheric magnetic field, and in particular the so-called heliospheric current sheet, give rise to charge-sign dependent effects in CR transport such as drift. All these processes are well correlated with the 11-year quasiperiodical cycle of solar activity (tout court Solar Cycle) and, for charge-sign dependencies, with with the 22-year cycle of magnetic polarity (Potgieter, 2013, 2014). At microscopic level, all these processes are governed by the rigidity (momentum/charge ratio) $\mathcal{R}=p/Z$, i.e., the quantity that best describes dynamics of particle (gyro)motion in magnetic fields. In principle, one can connect the spectrum of heliospheric turbulence with the rigidity dependence of the diffusion mean free path $\lambda$ of CR particles (Teufel & Schlickeiser, 2002).

In the recent measurements released by the AMS experiment (Aguilar et al., 2018a), the temporal dependence of the hadronic fluxes shows new details of CR modulation that can be used to test CR transport. More precisely, AMS has measured the temporal variation of the proton and helium fluxes in CRs, between 2011 and 2017, with a 27-day time resolution and a precision level of 1 %. When measured at the same value of rigidity $\mathcal{R}$, the CR fluxes of protons and helium are found to have very similar fine-structures in terms of time and relative amplitude. However, on yearly time scales, little differences have been observed between the two species. Such a difference between the proton and helium fluxes can be observed in the temporal evolution of their ratio, p/He, evaluated at rigidity $R\lesssim$ 3 GV. From the AMS data, the p/He ratio shows a remarkable decrease between March 2015 and May 2017, corresponding to the descending phase of solar activity, after the 2014 solar maximum and toward the 2019 solar minimum. This period is also called recovery phase for the Galactic CRs, because during this phase the absolute fluxes of proton and helium are seen to increase steadily. Possible explanations for the long-term p/He feature involve mass/charge dependent-effects in CR transport at low-rigidity, either differences in the local interstellar spectra (LIS) of CR proton and helium. In both cases, the p/He ratio can be also influenced by the 3He–4He isotopic mixing (Aguilar et al., 2018a; Jokipii, 1967; Gloeckler & Jokipii, 1967; Biswas et al., 1967; Herbst et al., 2017; Gieseler et al., 2017). In a recent paper (Tomassetti et al., 2018), we have show that the long-term p/He behavior observed by AMS represents a remarkable signature of the rigidity dependence of CR diffusion in the heliosphere that is described by a diffusion coefficient of the type $K(\mathcal{R})\propto\beta(\mathcal{R})\lambda(\mathcal{R})$. As we have argued, the velocity factor $\beta$ gives rise to little mass/charge dependence of CR diffusion in the low-rigidity limit, i.e., $\mathcal{R}\lesssim\,3$ GV. Similar results were also reported by other authors (Corti et al., 2019; Boschini et al., 2019; Luo et al., 2019). In all these works, albeit based on different modeling, the conclusions are driven by the availability of new important measurements: the Voyager-1 data collected outside the heliosphere, and the AMS multi-channel and monthly-resolved data, collected locally and over a significant fraction of the Solar Cycle. The improved modeling is presented in details in the rest of this paper. It consists in an accurate calculation of the proton and helium LIS’s (beyond the heliosphere), along with their corresponding modulated fluxes near-Earth (in the inner heliosphere). This paper is organized as follows. In Sect. 2, we recapitulate all the data we used in our analysis. In Sect. 3, we present the calculations of the LIS, i.e., our improved modeling in Galactic CR propagation in the ISM. In Sect. 4, we describe in details our numerical model of solar modulation and our data-driven approach that we have followed in order to test CR diffusion. In Sect. 5 we present the results and discuss the interpretation of the data using the numerical model and approximated calculations. We also discuss the main limitations of our modeling approaches. Conclusions are drawn in Sect. 6.

2 Observations

The model developments presented in this paper are motivated by the recent release of a large sample of CR data from the AMS collaboration. The data published by the AMS collaboration include lepton, proton, antiparticle spectra, primary nuclei, secondary to primary ratios and, more recently, flux time dependencies. All these data are available for download at the Cosmic-Ray Database of the Space Science Data Center (SSDC) at Italian Space Agency (Di Felice et al., 2017). In particular, the collaboration reported time-resolved measurements of CR proton and helium fluxes based on a total sample of 846 and 112 billion collected events, respectively (Aguilar et al., 2018a, b). Both proton and helium fluxes have been measured for the 79 Bartel’s rotations (BR), i.e., 27-day time intervals, ranging from May 2011 to May 2017. All proton fluxes range from 1 GV to 60 GV of rigidity. The helium fluxes range from 1.9 GV to 60 GV. This sample represents a simultaneous measurement of an uninterrupted time-series of CR proton and helium fluxes over an extended period of time. Overall, the variations of the proton and helium fluxes show similar fine structures, on monthly time scales. It was also noted that these structures are well correlated with the monthly number of sunspot (Corti et al., 2019), (i.e., with the Sun’s magnetic activity) and that their amplitudes decrease with increasing rigidity. With the precision of the AMS measurements, time-variations in the individual p-He fluxes are appreciated up to nearly 40 GV of rigidity. Most of these variations are canceled out in the flux ratio p/He of proton and helium at the same value of rigidity, The p/He ratio is observed to be constant at rigidity above 3 GV. Below 3 GV, however, the ratio show an interesting long-term time dependence. In particular, it decreases steadily in the period between Mar 2015 and May 2017, in the descending phase of the Sun’s activity and during which the absolute flux intensity increases rapidly.

To interpret these measurements within a data-driven approach, other sets of data are needed such as combined data from AMS and Voyager-1. We used these data to constrain the interstellar CR spectra for proton and helium, and to estimate their uncertainties. From the Voyager-1 experiment, we used the available interstellar data (i.e., collected after August 2012) down to 100 MeV/n energies. In particular, hydrogen flux data are available in the energy range 140-320 MeV. Data on CR helium and heavy ions data are available at 110-600 MeV/n (Cummings et al., 2016). These data have been collected by Voyager-1 in the period 2012-2016, i.e., while the spacecraft were moving in the interstellar space. Thus they provide direct and valuable constraints to CR propagation models. From AMS, we used data on the primary CR fluxes p-He-C-O at $\mathcal{R}>$ 60 GV (Aguilar et al., 2015a, b), and secondary-to-primary B/C ratio at $\mathcal{R}>$ 4 GV (Aguilar et al., 2018c). The minimal rigidity values are chosen to ensure that the data are representative of the interstellar CR flux (i.e., solar modulation free). We note that between the Voyager-1 flux data (below $\sim$ 1 GV) and the AMS data (above $\sim$60 GV), there is a large energy gap where no information is directly available on the LIS flux. Recovering the LIS fluxes inside the gap is possible, however, by means of a robust modeling of the Galactic CR propagation processes. The AMS data on primary fluxes and secondary-to-primary ratios are used for the determination of the model parameters, and eventually for the determination of the LIS shape (Sect. 3).

3 Galactic cosmic-ray propagation

The precise knowledge of the Galactic CR fluxes beyond the heliosphere is of crucial importance for the study of solar modulation, and in particular the interpretation of the tiny features observed in the p/He ratio. Traditionally, in solar modulation studies, the LIS fluxes constitute a mere physical input that is assumed to be known, so they are often provided as parametric input functions. The parametric approach was followed in all past studies where, in fact, large discrepancies were noted among the different works (Fisk, 1971; Gleeson & Urch, 1971; Moraal, 2013; Manuel et al., 2014; Herbst et al., 2017). A strong improvement was brought by the release of Voyager-1 data (Corti et al., 2016, 2019; Herbst et al., 2017) To build a LIS parametric model from the data, however, one has to deal with a large observational gap between the Voyager-1 LIS data (available only at $\mathcal{R}\lesssim$ 1 GV) and the AMS data (modulation-free only at $R\gtrsim$50 GV) (Corti et al., 2016). Hence, the use of a physically plausible shape for the LIS model is essential (Boschini et al., 2017). For this reason, we opted for a full calculation approach based on a physically motivated model of CR propagation in the Milky Way. The advantage of our approach is that it can provide a robust evaluation of the LIS uncertainties at all relevant energies.

To compute the Galactic CR fluxes of all particles, including proton and helium isotopes, we made use of a spatial dependent model of diffusive CR propagation in two halos (Tomassetti, 2012a, 2015b). In this model, the CR diffusion near the Galactic disk is characterized by a smaller rigidity dependence in comparison to CR diffusion in the halo away from the disk. The master equation describing CR transport in the Milky Way can be written as:

[TABLE]

where $N=N(E,r,z)$ is the particle number density as a function of energy and 2D spatial coordinates, $\Gamma=\beta cn\sigma$ is the destruction rate for collisions off gas nuclei, with density $n$, at velocity $\beta c$ and cross section $\sigma$. The source term $S$ is split into primary and secondary terms. The former describes the CR acceleration spectrum from Galactic sources (e.g., supernova remnants) $S_{\rm pri}\propto(\mathcal{R}/{\rm GV})^{-\nu}$, where $\nu=$ 2.28$\pm$0.12 for protons and $\nu=$ 2.35$\pm$0.13 for all heavier nuclei. The secondary production term for a given CR species is computed as $S_{\rm sec}=\sum_{\rm h}\Gamma_{h}^{\rm sp}N_{\rm j}$, describing the fragmentation of all heavier $h$-nuclei occurring at rate $\Gamma_{h}^{\rm sp}$. The term $\dot{E}=\--\frac{dE}{dt}$ describes standard processes of ionization and Coulomb losses. While all sources and gas are located in the disk ($z=0$), the CR transport occurs in the halo up to a vertical extent $z=\pm\,L$. The $L$ parameter sets the boundary conditions for the steady-state solution of Eq. 1, $\partial_{t}{N}\equiv{0}$. The CR transport is governed by a two-zone diffusion coefficient $D\equiv\beta D_{0}(\mathcal{R}/GV)^{\delta_{i/o}}$, where $D_{0}/L=0.01\pm$0.002 kpc/Myr. The two diffusion regimes are described by the index ${\delta_{i/o}}$, that takes two values: $\delta_{i}=0.18\pm$0.05 in the inner halo (near-disk region $|z|<\xi\,L$), and $\delta_{o}=\delta_{i}+\Delta$ in the outer halo ($|z|>\xi\,L$), where $\xi=0.12\pm$0.03 and $\Delta=0.55\pm$0.11. Once Eq. 1 is resolved for all CR species, giving $N=N(E,r,z)$, the local interstellar fluxes are computed as $J=\frac{\beta c}{4\pi}N_{\odot}$. Here the CR density is evaluated at the position of the Solar System $z_{\odot}\cong$0 and $r_{\odot}\cong$8.3 kpc.

In all the above mentioned parameters, the best-fit values and their corresponding uncertainties have been estimated by means of a Monte-Carlo Markov Chain sampling algorithm Feng et al. (2016). In addition to our benchmark setting, the effect of CR reacceleration in Alfvénic waves was also tested. An advantage of reacceleration models is their explanation for the GeV peak in the B/C ratio. We found, however, that our model favors only modest values for the Alfvénic speed, in agreement with the basic expectations for the interstellar medium, from which $v_{A}\approx\,B/\sqrt{4\pi\rho}\lesssim$ 6 km/s. In many past studies based on homogeneous diffusion, however, large Alfvénic speeds were reported (up to $v_{A}\sim\,30-40$ km s*-1*) (Trotta et al., 2011; Vladimirov et al., 2011). These models, however, require the introduction of low-energy breaks in the injection spectra (at $E\sim$ 10 GeV) in order to match the primary CR fluxes, while additional high-energy breaks are needed to describe the observed spectral hardening at $E\sim$ 200 GeV. On the other hand, in our model, the full shape of primary CR spectra is interpreted naturally in terms of CR propagation in two diffusive regions. Within this scenario, the GeV peak of the B/C ratio should be interpreted as a low-energy change of CR diffusion, due e.g., to magnetic damping (Ptuskin et al., 2006). A critical discussion to CR reacceleration can be found in Drury & Strong (2017).

Having modulation-free CR data is essential, in this phase of the work, to obtain unbiased LIS’s. Thus, we have examined the presence and influence of potential biases caused by residual modulation. In the LIS fitting procedure, the Voyager-1 data were used in combination with the high-rigidity AMS data on primary CR fluxes (Aguilar et al., 2015a, b), along with the newly released data on the B/C ratio (Aguilar et al., 2018c). With the use of minimal rigidity thresholds (60 GV for the absolute fluxes and 4 GV for the B/C ratio (Tomassetti, 2017a)) the effects of solar modulation are negligible at the $\lesssim$ 1 % level in the selected data. Nonetheless, as in (Feng et al., 2016), we accounted for residual modulation using a force-field modulation potential $\phi$ as nuisance parameter.

The comparison between pre-AMS (Feng et al., 2016) and post-AMS uncertainties is illustrated in Fig. 1. As shown in the figure, the new AMS data on the B/C ratio pose tight constraints to the CR propagation parameters, that are now comparable to the uncertainties arising from the fragmentation cross-sections (Tomassetti, 2015c, 2017a). From the tight B/C-driven constraints, we obtained improved calculations of the secondary 2H and 3He isotopes. Their production is computed using an improved set of fragmentation cross-sections (Tomassetti, 2012b; Tomassetti & Feng, 2017). In practice, the deuteron (2H) contribution turned out to be fully negligible for all purpose, being 2H/1H$\lesssim$ 1 %. On the other hand, the 3He nuclei contribute by $\sim$ 15 % to the total He flux in the GV region. Accounting for the 3He is essential for a precise description of the CR helium data (Biswas et al., 1967; Boschini et al., 2017). In agreement with Corti et al. (2019), we found that it causes a moderate impact on the temporal dependence of the p/He ratio.

In Fig. 2 the final proton and helium LIS’s resulting from calculations are shown. In the figure, the pink shaded bands represent the 1-$\sigma$ uncertainties estimated from the fit. It can be seen that the highest uncertainties lie in the $\sim$ 1-10 GeV energy region, where direct LIS data are not available. This region is also sensitive to CR transport parameters such as $\delta_{i/o}$, $D_{0}/L$, or $v_{A}$. However, proton and helium fluxes show very similar dependencies upon these parameters, so that the associated uncertainties are correlated. As a result, the p/He ratio of Fig. 2c show small uncertainties at all rigidities. In particular, it is interesting to note that the p/He ratio measured inside the heliosphere is very similar to its interstellar value, when it is computed as a function of rigidity. This explains why, before AMS, no time variations were observed for the p/He ratio (Bindi, 2017).

In Fig. 3, we plot some recent LIS derivations in comparison with the data (Tomassetti et al., 2018, 2017; Corti et al., 2016; Tomassetti, 2015a, 2017a, 2017b). These models are in agreement with our calculations within uncertainties, and they give consistent results. The high degree of agreement among the different CR propagation models is mainly ascribed to the tight constraints of the recent data. In particular, the Voyager-1 data have driven the CR LIS models to a new standard of precision in the low-energy region. From the pre-Voyager-1 studies, it can be noted that the proposed LIS were highly discrepant in the GeV and sub-GeV energy region (Herbst et al., 2017).

4 Solar modulation calculations

The transport of cosmic rays (CRs) in the heliosphere is described by the Krymsky-Parker equation (Krymsky, 1964; Parker, 1965). For a given CR particle type, the evolution of its phase-space density $\psi=\psi({\bf r},\mathcal{R},t)$ is governed by:

[TABLE]

where $\mathcal{R}$ is the particle rigidity. The various terms represent convection with the solar wind of speed ${V}$, spatial diffusion with tensor ${K}$, adiabatic momentum losses, and CR sources. Analytic solutions of this equation rely in simplifying approximations pertaining diffusion or losses. In spherical symmetry, for a wind radially flowing at speed $V$ and isotropic diffusion tensor, manageable numerical solutions can be found (Fisk, 1971; Gleeson & Urch, 1971; Moraal, 2013). The equation can be rewritten as:

[TABLE]

For steady-state conditions ($\partial\psi/\partial{t}=0$) and with zero source terms ($Q=0$), it becomes a parabolic equation of second order in space and first order in rigidity:

[TABLE]

The boundary conditions can be set at $r=r_{\odot}$, the Sun’s radius, and at $r=d$, the heliosphere’s boundaries. At $r=r_{\odot}$, the absence of any sources or sinks is assumed, so that ${\partial\psi}/{\partial{r}}|_{r_{\odot}}=0$ and the wind vanishes. At the heliopause $r=d$, the condition is stated by requiring that the CR flux is equal to its local interstellar (LIS) value, so that $J_{T}|_{r=d}=\mathcal{R}^{2}\psi=J_{\rm LIS}$. The initial condition results from the fact that the LIS is unchanged inside the heliosphere at sufficiently high rigidity $\mathcal{R}_{\rm max}$: $J(r,\mathcal{R}_{\rm max})=J_{\rm LIS}(\mathcal{R}_{\rm max})$ for every $r$ between $r_{\odot}$ and $d$. It should be noted that this condition accounts for the LIS at the boundaries ($\psi_{0}\equiv{Q}$) without the explicit specification of a source term $Q$ in the transport equation. To obtain the steady-state solutions for $\psi(r,\mathcal{R})$, we make use of the Crank-Nicolson implicit scheme. In this method, the domain is divided in a rectangular grid with $r$ ranging from $r_{\odot}$ to $d$ and $\rm{ln}\mathcal{R}$ ranging from $\rm{ln}\mathcal{R}_{\rm min}$ to $\rm{ln}\mathcal{R}_{\rm max}$. The CR distribution function $\psi(r,\mathcal{R})$ is represented as $\psi^{j}_{i}$ in in the $N\times{M}$ space of coordinates $r_{i}=r_{\odot}+i\Delta_{r}$ and $\rm{ln}\mathcal{R}_{j}=\rm{ln}\mathcal{R}_{0}+j\Delta_{\rm{ln}\mathcal{R}}$, where $i=\,1\dots N$ and $j=\,1\dots M$. The method consists of taking the average between the derivative evaluated at node $j$ and at node $j+1$, leading to discrete first and second order derivatives:

[TABLE]

After defining the quantities $\varphi_{1}=K$, $\varphi_{2}=-V+\frac{2K}{r}$, and $\varphi_{3}=\frac{2}{3}\frac{V}{r}+\frac{1}{3}\frac{\partial{K}}{\partial{r}}$, the equation becomes:

[TABLE]

By renaming the LHS coefficients as $\eta_{1}$, $\eta_{2}$, $\eta_{3}$, and the RHS one as $\eta_{4}$, the system can be rewritten in a compact form $A{\bf x^{j+1}}=B{\bf x^{j}}$, where $A$ and $B$ are tridiagonal matrices:

[TABLE]

The partial differential equation is then transformed in a set of linear equation relating $\psi^{j}$ and $\psi^{j+1}$. Every time the matrix is inverted, $\psi^{j}$ progresses one step in rigidity toward $\rm{ln}(\mathcal{R}_{\rm min})$. The full solution is found in the whole domain by stepping down in rigidity and solving the tridiagonal system at each step. In this work, we have set $N=$ 610 radius nodes, $M=$ 500 rigidity nodes, $\mathcal{R}_{\rm min}=$ 50 MV and $\mathcal{R}_{\rm max}=$ 1 TV for all CR particles. The tridiagonal matrix equations are solved using an iterative method. The equation has been resolved for several isotopic species relevant for this work. In contrast to the force-field approximation, widely used in CR astrophysics, the numerical approach presented here is suitable to test arbitrary forms of the CR diffusion coefficient $K$, or to study mass/charge-dependent effects arising from its functional dependence on rigidity or velocity. In this work, we adopt a benchmark diffusion coefficient of the form $K=K_{0}\beta(\mathcal{R})\left(\mathcal{R}/{\mathcal{R}_{0}}\right)$ where, as we will discuss, the mass/charge dependence of $K(\mathcal{R})$ are contained in the $\beta(\mathcal{R})$ factor. The overall normalization is set by the $K_{0}$ parameter (at the reference rigidity $\mathcal{R}_{0}\equiv$ 1 GV) and it is of the order of $\sim{10^{22}}$ cm2s*-1*. In the follow, the main focus is given to the fluxes of total hydrogen (1H +2H) and total helium (3He +4He).

5 Results and discussion

5.1 Fitting the cosmic ray proton fluxes

Our fits to the proton fluxes are shown in if Fig. 4. In the 16 panels figure, the AMS measurements of the proton fluxes are shown as open squares for 16 BR’s out of 79 (Aguilar et al., 2018a). The light orange lines describe our calculations for the best-fit proton fluxes, $\hat{J}_{\rm p}$. The fits have been carried out independently for each BR using the least squares method. For a given BR, the minimizing function is defined as:

[TABLE]

In the equation, a rigidity integration of the model output is computed within each bin $\Delta\mathcal{R}_{j}$, while for the data the mean value is used. The $\sigma_{j}$ factors account for experimental uncertainties in the measured fluxes. The theoretical flux is $J_{p}(\mathcal{R})$, resulting from the Galactic LIS calculations of Sect. 3 and from solar modulation model of Sect. 4.

In the solar modulation calculations, we have adopted a benchmark diffusion coefficient $K$ of the form $K=K_{0}\beta\left(\mathcal{R}/{\rm{GV}}\right)$, where we express its normalization as $K_{0}\equiv{10^{22}}k_{0}$ cm2s*-1*, so that the adimensional factor $k_{0}$ is determined from the data. Because the fits are repeated for each BR, the procedure returns a time-series of best-fit values $\hat{k}_{0}$ describing the temporal evolution of CR diffusion coefficient. Such an evolution is at the basis of the time-dependent nature of CR modulation (Manuel et al., 2014; Bobik et al., 2016), At this point, we note that $k_{0}$ is degenerated with the wind speed $V$. In particular, the data are only sensitive to the combination $k_{0}/V$, so that the use of different (or variable) $V$-values can be compensated by a proper rescaling of $k_{0}$. Hence, the wind speed $V$ was kept fixed to the reference value of 400 km s*-1*. In our model, the problem is modeled under the usual quasi-steady approach, i.e., by providing a time series of steady-state solutions for $J_{\rm p}$ corresponding to a time series of input parameters $k_{0}$. Such an approach is justified by the different timescales between the CR transport in the heliosphere and the changing solar activity. In principle, the CR modulation phenomena cover large spectrum of timescales, from hours to years, in relation to the physics process involved. In particular, the CR diffusion timescale decreases rapidly with increasing rigidity (Strauss et al., 2011). At the typical rigidity of the AMS data ($\mathcal{R}\gtrsim\,2$ GV), the characteristic timescale of CR modulation does not exceed the monthly resolution of the data. Moreover, the analyzed long-term p/He behavior is manifested at yearly timescales. Nonetheless, we stress that each element of the best-fit time series $k_{0}(t)$ should be regarded as an effective scale value for the diffusion coefficient (i.e., averaged over the CR propagation histories), not necessarily representing the instantaneous conditions of the heliospheric turbulence (Strauss et al., 2011). To test the robustness of the fitting results, we also made specific checks for the possible influence of drift effects. For this purpose, we made dedicated runs within a formulation where the global effect of drift is captured by an additional free parameter $v_{D}$ in the convective term, representing the drift speed, without giving adiabatic cooling (Jokipii et al., 1977; Isenberg & Jokipii, 1978). This causes no appreciable effects in the phenomenology of the p/He ratio while, for the absolute fluxes in the $\mathcal{R}\gtrsim$ 2 GV region, the $v_{D}$ parameter is highly degenerated with the diffusion parameters. Degeneracy between drift and diffusion is quite unavoidable in a 1D formulations, and can be removed only in more spatially complex models. Finally, to test the presence of undetected systematic bias in the physical inputs (or in the data), we repeated the fitting procedure after leaving a global normalization factor as nuisance parameter.

The best-fit time-series $\hat{k}_{0}$ is shown in Fig. 5(a), along with the time profile of the CR proton fluxes at $\mathcal{R}=2$ GV, in Fig. 5(b), in comparison with the AMS data. In the figure, the green shaded band represents the level of uncertainty on the proton LIS, if it is directly propagated to the modulated spectra. However, since $k_{0}$ is a free parameter, global variations in the proton LIS are partially re-absorbed in the fit, i.e., by a proper scaling of the $k_{0}$-values and the nuisance parameters. The relevant uncertainties in the proton flux are shown by the orange shaded band. This band accounts for the standard uncertainties from the fitting method, i.e., corresponding to one-sigma variations of $\chi^{2}$ curve as function of the free/nuisance parameters, where the procedure is repeated for a large number of LIS functions. In practice, the LIS-functions are Monte-Carlo generated according to the probability density functions of the underlying key parameters (see Sect. 3), i.e., all input LIS’s that produce the band of Fig. 2 (Feng et al., 2016). With this procedure, uncertainties associated with the LIS modeling are included in the final error band, with the proper account for their bin-to-bin correlation. The orange band can also be regarded as the level of precision to which our LIS knowledge could be reached, potentially, if the solar modulation is well understood and precisely modeled.

5.2 Predicting the helium flux and the p/He ratio

Once the $k_{0}$ time-series is determined from the AMS proton data, we have used our proton-calibrated model to predict the flux of CR helium. This includes the calculation of 3He and 4He isotopes. Helium flux calculations are shown in Fig. 4 for 16 selected epochs. The AMS data on CR helium are shown as open circles. They are well matched by the blue line, showing the model prediction for the total helium flux (i.e., 3He plus 4He). The green band shows the sub-dominant 3He component and its uncertainty. Because 3He is purely secondary in CRs (i.e., it is produced by fragmentation of primary CR particles such as 4He or 12C), its calculations relies strongly on the CR propagation parameters which, in turn, are constrained by the B/C ratio. Thus, an improvement on the LIS helium calculation is due to the new B/C data from AMS (Aguilar et al., 2018c). From the proton-driven calculations of the helium flux, we made a corresponding prediction for the p/He ratio and in particular for its temporal dependence. The temporal dependence of CR helium at $\mathcal{R}\approx$ 2 GV is shown in Fig. 5(c), where our predictions for the total flux are compared with the AMS data. The 3He isotopic flux and its temporal dependence are also shown in the panel, as green dashed line. Similarly, Fig. 5(d) reports predictions and data for the p/He ratio. Along with the benchmark model, in Fig. 5(d) we also plot the predictions from a diffusion coefficient of the form $K{\propto}k_{0}\mathcal{R}$ (pink dashed line), i.e., without any $\beta$-dependence. In this case there are no $A/Z$ terms in $K(\mathcal{R})$, which is a sole function of rigidity, and it can be seen that the p/He ratio is remarkably constant. Within this model, the diffusion of proton and helium is identical, at a given rigidity, so that any difference in their modulation is only caused by their different LIS shapes. This shows that the different spectral shapes play a minor effect in the p/He time dependence. On the other hand, before the recent measurements from AMS and PAMELA, the same spectral shape was assumed for the protons and helium LIS. The impact of p-He LIS differences in the solar modular is discussed in Gieseler et al. (2017); Corti et al. (2019); Webber (2015).

The p/He predictions from our benchmark diffusion coefficient ($K\propto\beta\mathcal{R}$) are shown in the blue solid line. It can be seen that, from this scenario, the predicted ratio agrees very well the observations. The relevant uncertainties of the model predictions are shown as shaded bands. They arise primarily from uncertainties in the p-He LIS models (as discussed above) and from uncertainties on the 3He components. These contributions are linked to the Galactic CR propagation parameters. A minor contribution to the uncertainty band arises from the statistical errors of the fitting procedure. Within the estimated uncertainties, our predictions of the temporal evolution of the p/He evolution agree well with the AMS data.

From the data, it can be noted that the p/He behavior is correlated with the intensity of the individual p-He fluxes Based on such a correlation, one may expect that the p/He ratio is subjected to a cyclic modulation behavior, This is indeed the case, as shown in Boschini et al. (2019). As results of simulations, the authors found that the p/He ratio is subjected to a quasi-periodical 11-year behavior. The absolute variation of the ratio is related to the level of solar activity. With the present data, the variation in the p/He ratio is observed during the recovery phase which begins approximately one year after the maximum of sunspots in Solar Cycle 24. During this phase the CR fluxes increase to nearly a factor of $2$ in about two years. The one year delay with respect to the 2014 solar maximum can be related to the CR modulation time lag. Such a lag can be interpreted as the time spent by the solar wind to transport the Sun’s induced magnetic perturbations, from the corona to the outer heliosphere (Tomassetti et al., 2017).

5.3 On the interpretation of the p/He behavior

An important results of our work in (Tomassetti et al., 2018) is that the time evolution of the low-rigidity p/He ratio (as measured by AMS at $\mathcal{R}\approx$ 3 GV) is described very well with a $\beta\times\lambda(\mathcal{R})$ dependence of the CR diffusion coefficient in the heliosphere. This implies that the measured long-term dependence of the p/He ratio is a signature of the universal rigidity dependence of the mean free path $\lambda(\mathcal{R})$. Given that the diffusion coefficient is $K(\mathcal{R})=\beta\lambda(\mathcal{R})/3$, the measured differences between the proton and helium flux evolution are caused by their difference in the $A/Z$ ratio between charge and mass number. In the non-relativistic limit, proton and helium nuclei at the same value of rigidity, must travel at different speed with $\beta_{p}>\beta_{He}$. This is apparent from the expression of the CR speed as function of rigidity:

[TABLE]

Thus, at fixed rigidity, CR protons ($Z/A=$1) experience a slightly faster diffusion than helium ($Z/A\approx$1/2). On the other hand, in the relativistic limit one has $\beta\rightarrow$1, so that the diffusion of the two particles becomes identical.

The rigidity dependence is illustrated more clearly in Fig. 6. In the left panel, the p/He time profile is plotted for several rigidity values from $\sim$ 2 GV to $\sim$ 20 GV (from top to bottom). From this figure, the overall trend of the AMS data is well described by the model. It can also be seen that, for increasing rigidity, the p/He ratio decreases rapidly from about $8$ to about $5$ (on average), following its interstellar rigidity dependence. The p/He dependence upon rigidity is also shown in Fig. 2(c). Similarly, the temporal dependence of the p/He ratio weakens as well at increasing rigidity, so that at $R\gtrsim$ 3 GV the p/He is essentially constant. An illustration of the rigidity dependence is given in the right panel of Fig. 6. The plot shows the quantity $\Gamma_{\rm p/He}\equiv{[\textsf{p/He}]_{t_{2}}}/{[\textsf{p/He}]_{t_{1}}}$ computed at two epochs $t_{1}$ and $t_{2}$. As reference epochs we have chosen February 2014 ($t_{1}$, where the p/He ratio is maximum) and May 2017 ($t_{2}$, where the ratio is minimum). As shown by the solid line, the $\Gamma_{\rm p/He}$ function is predicted to increases monotonically with rigidity, but it reaches a plateau as soon as the rigidity is relativistic (i.e., a few GV). This prediction agrees very well with the AMS measurements.

5.4 The convection-diffusion approximation

A simple and elegant interpretation of the p/He behavior can be obtained within the convection-diffusion approximation (CDA). The CDA approach represents a low-order approximated solution of Eq. 2, where all spatial dependencies are neglected (0D). In the CDA, the key hypothesis is that the diffusive current of CR particles ($\sim{K}\partial\psi/\partial{r}$) is exactly countered by the convective current ($\sim{V}\psi$). Any other process (such as magnetic drift or loss terms $\partial/\partial\rm{ln}\mathcal{R}$) are neglected. Since the convective current moves always outward, following the solar wind, the key CDA hypothesis implies that the diffusive current acts inward, which is not strictly true. From the simple balance of the two terms the transport equation becomes:

[TABLE]

In steady-state conditions, this equation has a simple analytical solution $\psi=\psi_{0}e^{-\mathcal{M}}$, where $\psi_{0}$ is the LIS phase-space density for the considered CR particle type. The changing conditions of the background plasma are captured by the adimensional parameter $\mathcal{M}$:

[TABLE]

For typical values $V\sim$ 400 km/s, $d\sim$ 100 AU, and $K_{0}\sim$ 1022 cm2/s, the parameter $\mathcal{M}$ at the GV scale is of the order of unity. In practice, $\mathcal{M}$ is determined experimentally from the CR data for a given LIS hypothesis. For a diffusion coefficient of the type $K\propto\,\beta\mathcal{R}$, it follows that to first approximation the modulation strength is inversely proportional to $\beta\mathcal{R}$. Within the CDA, the p/He ratio is of the type:

[TABLE]

where $K_{\rm p}$ and $K_{\rm He}$ are the diffusion coefficients of CR protons and helium nuclei, the difference of which arises from the $A/Z$ factors contained in the CR velocity of Eq. 12. The CDA solution shows explicitly that the modulation depends upon the parameter combination $\mu{\equiv}Vd/K_{0}$, which captures the effects of changing conditions in solar activity and the underlying parameter degeneracy. We have also verified that such a degeneracy is respect even by the numerical solution, as long as drift and latitudinal dependence are not considered. In the limit $\mathcal{M}\lesssim$ 1, the p/He ratio of Eq. 15 can be rewritten in the following form

[TABLE]

that shows clearly the dependence of the ratio upon $\beta(\mathcal{R})$. Similarly, the $\Gamma_{\textsf{p/He}}(\mathcal{R})$ function can be readily calculated as:

[TABLE]

From the above relations, one can easily understand the dependence of the p/He ratio upon the variable strength of CR modulation, expressed by the effective parameter $\mu(t)$,

Our calculations within the CDA framework are shown in Fig. 7. Again, the CR diffusion coefficient is of the type $K\propto{k_{0}}\beta\mathcal{R}$, where we used $A/Z\cong$ 1 for CR protons and $A/Z\cong$ 3.8 for helium, to account for their average isotopic composition. As discussed, the time dependence of the combination $Vd/k_{0}$ is represented by a time-series of $k_{0}$ parameters. In the CDA calculations of Fig. 7, the input time-profile $k_{0}=k_{0}(t)$ is taken from the fit results obtained in Sect. 5.1. From the figure, the temporal evolution of the p/He ratio is described fairly well within the CDA framework. Similarly, a reasonable rigidity dependence for $\Gamma_{\textsf{p/He}}(\mathcal{R})$ is obtained.

From Eq. 17, one may note that the $\Gamma_{\textsf{p/He}}(\mathcal{R})$ function appears as LIS-independent within the CDA approach, due to the fact that the LIS part is factorized out in in Eq. 15. We emphasize, however, that the LIS dependence is implicit from the data-driven (e.g., proton driven) determination of the parameter $\mathcal{M}$, and thus the CDA calculations for the p/He evolution do depend on the assumed LIS. On the other hand, the different LIS shapes for protons and helium has an impact in the p/He time evolution that cannot be accounted within the CDA framework. Since energy losses are neglected, the CDA is also known to systematically underestimate the low-energy CR fluxes in the inner heliosphere. It follows that this approximation is far inadequate for describing the data at the accuracy level of AMS. The limitations of the CDA approximation are reviewed in details in Refs. (Moraal, 2013; Gleeson & Urch, 1971).

5.5 The parametric description

A parametric description for the evolution of the p/He ratio can be obtained by means of the simple broken-line model (Aguilar et al., 2018a):

[TABLE]

In this function, $A=A(\mathcal{R})$ is the average p/He ratio value from May 2011 to $t_{0}$, being $t_{0}=t_{0}(\mathcal{R})$ the time when the p/He ratio deviates from the average $A(\mathcal{R})$-value, and $B(\mathcal{R})$ is the slope of its temporal variation. Using Eq. 18, we made fits to the AMS data for every bin of rigidity between 2 and 60 GV. At rigidity above 3.29 GV, the ratio becomes consistent with a constant value $A$ at the 95 % confidence level, i.e., it becomes time independent. Below 3.29 GV, the observed ratio becomes steadily decreasing ($B<$ 0) after some epoch $t=t_{0}$. The combined best-fit value for $t_{0}$ is February 28 $\pm$ 42 days, 2015, at all rigidity values. This epoch coincides with the beginning of the recovery phase of the CR fluxes.

The best-fit of the p/He data using Eq. 18 is shown in Fig. 8. In the left panel, the best-fit broken-line description of the p/He behavior is shown in comparison with the data. Our numerical calculations are also shown for comparison. Overall, in spite of its lack of physical motivations, Eq. 18 provides a satisfactory parametric description of the data at a fixed rigidity. In fact, the broken-line description was aimed at capturing the gross feature of the p/He evolution, i.e., its decreasing behavior during the recovery phase. On the other hand, since the fits are made independently for each rigidity bin, the $\Gamma_{\textsf{p/He}}(\mathcal{R})$ behavior shown in the right panel appears dominated by fluctuations of the same order of the data.

In this work, we have provided a simple explanation for such a feature. According to our calculations, the expected feature for the p/He ratio is a progressive variation, the amplitude of which depends on the changes in the CR diffusion coefficient. In the period between March 2015 and May 2017, the variation is appreciable in the AMS data. This period is the so-called “recovery phase” of the CR flux modulation, i.e., when solar activity decreases rapidly toward the end of Solar Cycle 24.

6 Conclusions

In this paper we discussed in details the solar modulation calculations that we have developed for interpreting the recent CR data released from AMS. To describe the data at precision level demanded by AMS, we carried out improved calculations for the interstellar fluxes of Galactic CR nuclei, along with their uncertainties, and numerical calculations of CR propagation in the heliosphere tightly calibrated against the new AMS proton data. Our data-driven approach enabled us to predict the temporal dependence of the p/He ratio and to test specific forms for the diffusion coefficient of CRs in the heliosphere. In the calculations discussed here there are numerous simplifications. For instance, we have neglected the tensor nature of CR diffusion which is essential for describing the effects of anisotropic diffusion or curvature drift. We have also introduced simplified description of the heliosphere, by neglecting the latitudinal dependence of the solar wind, the termination shock, the heliospheric current sheet. Accounting for all these features would inevitably increase the model complexity, thus requiring more parameters. These simplification should be kept on mind in the interpretation of the basic parameters of the model, as they should be regarded as effective values. For example, the neglection of magnetic drift may affect the best-fit value for the diffusion parameter $k_{0}$. We expect that any drift-induced bias, however, would be the same for CR proton and helium particles. As we have discussed, the monthly-resolved measurements of the p/He ratio are well suited to test CR diffusion in the heliosphere. Being a ratio of positively charged particles, it cancels out potential biases arising from, e.g., magnetic drift or anisotropic diffusion.

As we have show throughout this paper, the temporal behavior of the p/He ratio is caused by the $A/Z$-dependence in the CR diffusion coefficient. When expressed as function of rigidity, it can be written as $K(\mathcal{R})=\beta(\mathcal{R})\times\lambda(\mathcal{R})/3$, where the mean free path $\lambda$ is a universal function of rigidity, so that the $A/Z$-terms are contained in the equation for the CR speed as function of rigidity $\beta(\mathcal{R})$, Eq. 12. The remarkable agreement between our prediction and the AMS data supports the concept of universality in the CR propagation histories in the heliosphere. In the future, further tests will be made using data over a larger fraction of the Solar Cycle, or measurements of different nuclear ratios.

Acknowledgements

B.B., E.F., and N.T. acknowledge support of Italian Space Agency (ASI) under agreement ASI-UniPG 2019-2-HH.0. F.B. and M.O. acknowledge support of Portuguese Science Foundation (FCT) under agreement CERN/FIS-PAR/0020/2017. The data used in this work were obtained from the SSDC/ASI Cosmic-Ray Database (https://tools.ssdc.asi.it/CosmicRays).