New results in solar modulation modeling in light of recent cosmic-ray data from space

TL;DR

This paper introduces an improved, validated model of solar modulation of cosmic rays, leveraging recent space-borne data to better understand how solar activity influences cosmic ray flux within the heliosphere.

Contribution

The paper presents a novel, validated model of solar modulation that incorporates recent cosmic-ray data and detailed heliospheric physics, enhancing understanding of GCR transport and solar activity effects.

Findings

Enhanced model of solar modulation validated with recent data

Identified key parameters linking GCR flux to solar activity

Provided new insights into GCR transport mechanisms in the heliosphere

Abstract

Thanks to space-borne experiments such as the AMS-02 and PAMELA missions in low-Earth orbit, along with the Voyager spacecrafts in the interstellar space, a large collection of multi-channel and time-resolved Galactic cosmic ray (GCR) data has recently become available. Here we present an improved measured-validated model of the "solar modulation" effect, i.e., the temporal evolution of the GCR flux inside the heliosphere caused by the 11-year variability cycle of the Sun's magnetic activity. We present our improved modeling of the structure of the heliosphere, the physical mechanisms of diffusion, drift, and energy losses of GCR particles in the heliosphere. We present our results for the temporal dependence of the key model parameters and their relationship with solar activity proxies. We discuss implications for the GCR transport in magnetic turbulence, and new insights on our understanding of the solar modulation phenomenon.

Equations

$$\frac{\partial f}{\partial t}=\nabla \cdot [{\mathbf{K}}^S\cdot \nabla f]-({\vec{V}}_{sw}+{\vec{V}}_D)\cdot \nabla f+\frac{1}{3}(\nabla \cdot \vec{V})\frac{\partial f}{\partial (lnR)}$$

$${\lambda }_{\parallel }\sim r_L^2\frac{⟨B^2⟩}{w(k_{\mathrm{res}})}\sim R^2/w(k_{\mathrm{res}})$$

$$K_{\parallel }=\frac{K_0}{3}\beta \left(\frac{B_0}{B}\right){\left(\frac{R_0}{R}\right)}^a\times {\left[\frac{(R/R_0{)}^h+(R_k/R_0{)}^h}{1+(R_k/R_0{)}^h}\right]}^{\frac{b-a}{h}}$$

$${\chi }^2(\vec{q})=\sum _i\frac{{\left[J_d(E_i,t)-J_m(E_i,\vec{q})\right]}^2}{{\sigma }^2(E_i,t)}$$

Full text

1 Introduction

The spectrum of Galactic cosmic rays (GCRs) observed near-Earth is significantly different from the so-called Local Interstellar Spectrum (LIS), i.e., their spectrum beyond the boundaries of the heliosphere, because GCRs are subjected to several physical mechanisms during their propagation through the interplanetary space. More specifically, charged particles in the heliospheric plasma can be spatially diffused, magnetically drifted, advected and decelerated by the outflowing solar wind and its embedded magnetic field. The modifications of the GCR intensities and energy spectra are known to vary with time, in connection with the Sun’s magnetic activity. This phenomenon is referred to as solar modulation of GCR, and it is observed to change periodically, following the periodical change in the sunspot number (SSN) observed in the solar corona. The SSN is a widely used in solar physics, e.g., as proxy for solar activity or characterizing the Solar Cycle, as well as for studying its correlation with the GCR modulation properties [1, 2].

A thorough understanding of solar modulation is then important, in the physics of GCRs, either to infer their LIS’s or to investigate the dynamics of charged particles in the heliospheric turbulence. Understanding the evolution of the GCR radiation is also essential for astronauts and the electronic components radiation hazard during long-duration missions. Along with the Voyager-1 data beyond the heliosphere [3], the new precise data from AMS-02 [4, 5] and PAMELA [6, 7] experiments offer a unique possibility to study the solar modulation over a long period of time.

2 The Numerical Model

The propagation of GCRs in the heliosphere is described by the Parker-Krymsky equation [8]:

[TABLE]

The equation describes the time evolution of the GCR phase space density $f(t,R)$, where $R=p/Z$ is the GCR particle rigidity (momentum/charge ratio), $\vec{V}_{sw}$ is the solar wind speed, $\vec{V}_{D}$ is the guiding center drift speed, and $\mathbf{K}^{S}$ is the symmetric part of the GCR diffusion tensor [9, 10]. The GCR flux J=J(t,R) is related to $f$ by $J=R^{2}f$. In this work, the equation is solved in steady-state conditions ($\partial/\partial{t}=0$), using the stochastic approach [11], based an upgraded and customized version of the code Solarprop [12, 13]. We use a 2D model, with radius $r$ and heliolatitude $\theta$ as spatial coordinates. The heliosphere is regarded as a spherical cavity, placed around the Sun, where the wind flows radially. We place the termination shock (TS) at $r_{\rm TS}=85$ AU from the Sun. Here the wind drops to subsonic speed. The heliospheric boundary is taken at $r_{\rm HP}=122$ AU, the heliopause (HP), where the wind vanishes and the GCR flux joins its interstellar value. Within this structure, the Earth lies in the equatorial plane at $r_{\rm E}=$1 AU. Along with its radial direction, the speed of the solar wind is known to depend upon latitudinal and solar activity. The latitudinal dependence appears well pronounced during period of minimum activity (low SSN) while, in contrast, the wind profile becomes highly symmetric during solar maximum (high SSN). Furthermore, in the equatorial plane the speed is highly stable with time, so that it is often set to the constant value of $400$ km s*-1*. In our model, we adopt a global parameterization for $V_{sw}(r,\theta)$, as in [14]. Regarding its influence on the GCR modulation, the wind transports a frozen-in turbulent Heliospheric Magnetic Field (HMF), wounded up in a rotating spiral structure. Of great importance for GCR modulation is the Heliospheric Current Sheet (HCS), sketched in Fig. 1 The HCS is a co-rotating magnetic structure which divides the HMF into two hemispheres of opposite polarity.

The angular extension of the HCS is described by the so-called tilt angle $\alpha$. The tilt angle is the angle between the Sun’s rotational axis and its main magnetic axis. It sets the level of the HCS waviness, which is time dependent. The $\alpha$-angle is constantly monitored by the Wilcox Solar Observatory (WSO) on a 10-day basis [15, 16]. It ranges from a few degrees during solar minimum (flat HCS) to nearly 900 during maximum and reversal (wavy HCS). In particular, magnetic drift of GCRs is important across the wavy layer of the HCS. More in general, the motion of GCR particles in the HMF is decomposed in a regular gradient-curvature large-scale motion and a HCS drift motion, both occurring on the background average HMF, plus an erratic random-walk scattering off the small-scale HMF turbulence, from which GCR diffusion originates. All these effects can be included in a diffusion tensor $\mathbf{K}$, made of a symmetric part that describes the diffusion and a asymmetric part that describes drift: $\mathbf{K}=\mathbf{K}^{S}+\mathbf{K}^{A}$, with $K_{ij}^{S}=K_{ji}^{S}$ and $K_{ij}^{A}=-K_{ji}^{A}$. The elements of the symmetric tensor that describes the diffusion $K^{S}_{ij}$ can be derived from kinetic theory, following a microscopic approach from GCR pitch-angle scattering off the random HMF irregularities. A key input parameter for the cosmic ray transport is the parallel spatial diffusion coefficient $K_{\parallel}$, which is usually expressed in terms of the mean free path $\lambda_{\parallel}$ along the background HMF: as $K_{\parallel}=\beta c\lambda_{\parallel}/3$, where $\beta=v/c$ is the GCR speed. GCRs are also subjected to perpendicular diffusion, with an associated length $\lambda_{\perp}$, that is usually assumed to be proportional to the parallel diffusion length: $\lambda_{\perp}=\xi\lambda_{\parallel}$, with $\xi\cong\,0.02$ [17]. The scattering of the GCR on the HMF irregularities is maximum when the GCR Larmor radius $r_{L}$ is of the same order as the spatial scale size of the irregularities $\hat{\lambda}$. Thus, the scattering for GCR particles with rigidity $R$ will depends upon the HMF power spectrum around the wave number $k$, corresponding to the spatial scale $\sim r_{L}(R)$. The essential dependence of $\lambda_{\parallel}$ on the HMF power spectrum can be therefore expressed as:

[TABLE]

where $\langle B^{2}\rangle$ is the mean square value of the background field and $w(k_{\rm{res}})$ is the power spectral density of the HMF fluctuations. The spatial scale of the resonant component of the spectrum is $\hat{\lambda}=2\pi/k_{\rm{res}}$. From the cyclotron resonance condition, $r_{L}\sim\hat{\lambda}$, particles with rigidity $R$ will resonate at $k_{\rm{res}}\sim 1/R$. The $w(k)$ follows a power-law $w(k)\propto k^{-\eta}$, where the index $\eta$ depends on the type and on the spatial scales of the turbulence energy cascade [18]. Thus, $\lambda_{\parallel}$ depends on the turbulence spectral index as $\lambda_{\parallel}\sim R^{2-\eta}$. On a wide range of scale, at least three regimes (with different index $\eta$) can be distinguished for the HMF power spectrum [18] that, in turn, will correspond to three diffusion regimes for GCRs: the $1/f$ range at large scales (i.e., small $k$, high rigidity), the inertial range at intermediate scales, and sub-ion range at the smallest scales (large $k$, low rigidity). At rigidity above few hundreds MVs, relevant to this work, GCRs resonate in the $1/f$ and $inertial$ scale range. Thus, a good parameterization for the rigidity dependence of $\lambda_{\parallel}$ is the double power-law function, defined by two indices and a critical rigidity value. In this work, we have adopted the following parametric description for the parallel diffusion coefficient $K_{\parallel}$ [14]:

[TABLE]

where $K_{0}$ is a constant in units of $10^{23}$ $cm^{2}s^{-1}$, with $R_{0}\equiv$ 1 GV, $B$ the HMF magnitude and $B_{0}$ the local field value (at 1 AU). The parameters $a$ and $b$ are two power indices that set the slope of the rigidity dependence, respectively, below and above the rigidity $R_{k}$. The parameter $h$ determines the smoothness of the transition between the two regimes. The perpendicular diffusion follows from $\lambda_{\perp}=\xi\lambda_{\parallel}$, with the addition of polar corrections (based on the Ulysses data at high latitudes) made to $\lambda_{\perp}$ and to the HMF model [19, 20]. The parameter set ${K_{0},a,b}$ then specifies the diffusion of GCRs on the HMF turbulence. Because these parameters depend on the properties of the background plasma, they may be subjected to temporal evolution following the Solar Cycle [21]. To test the time evolution of these parameters, we use the recent data on GCR protons from AMS-02 and PAMELA [4, 6, 7].

At this point it is important to note that Eq. 1 requires the flux at the HP to be specified as boundary condition, for a given GCR particle species. Determining the GCR LIS requires the development of dedicated models of Galactic physics processes, e.g., from the distribution of GCR sources, particle acceleration mechnisms and diffusive propagation in the Galactic turbulence [22]. In the work presented here, the GCR proton LIS relies on improved calculations from recent works [23, 24, 25]. In these calculations, the proton LIS was tightly constrained with low-energy data (at $\sim$ 100 -500 MeV of kinetic energy) collected from Voyager-1 beyond the HP [3] and AMS-02 precision measurement performed in the high-energy regime ( $E\gtrsim$ 100 GeV), i.e., where solar modulation effects become asymptotically negligible [4, 26, 27]. The GCR proton LIS of this work agrees fairly well with other spectra proposed recently [28, 30, 29, 31, 32, 33]. Our LIS is shown in Fig. 1, where comparison is made with the data from Voyager-1, AMS-02, and PAMELA.

3 The parameter extraction

Our model has six free parameters, in principle time-dependent, that have to be determined using GCR or Solar data at a given epoch $t$: the tilt angle of the HCS, $\alpha(t)$, the HMF intensity at Earth’s location $B_{0}(t)$, the HMF polarity $A(t)$, the global normalization factor of the diffusion tensor $K_{0}(t)$, and two spectral indices of its rigidity dependence, $a(t)$ and $b(t)$, defined at rigidity below and above the break parameter $R_{k}$. In this work we follow the quasi steady-state approach, where the time-dependence of the problem is described as a continuous series of steady-state solution of Eq. 1. In practice, time-series of data are used to extract a corresponding time-series of best-fit parameters. For a given epoch defined by the time $t$, the parameters $\alpha(t)$, $B_{0}(t)$ and $A(t)$ can be directly determined with Solar observations from WSO. These parameters describe the average status of the HMF in a given epoch. Since $\alpha$ and $B_{0}$ are subjected to continuous variations, we have computed an average value at a given epochs $t$, performed over a time window $[t-\Delta{T},t]$, with $\Delta{T}=$ 12 months. The window is chosen so that the average values of $\hat{\alpha}$ and $\hat{B}_{0}$ reflect the average HMF conditions sampled by GCRs when propagating from the heliopause to Earth [31]. In fact, the model assumes a static heliosphere, where its status is frozen during the travel of GCRs toward Earths.

The remaining parameters $K_{0}$, $a$, and $b$ are those describing the properties of GCR diffusion on the heliospheric plasma. They have been determined using GCR proton measurements from AMS-02 and PAMELA. In practice, to fit the GCR data, we have built a six-dimensional discrete grid of the model parameters vector $\vec{q}=$ ($\alpha$, $B_{0}$, $A$, $K_{0}$, $a$, $b$), for about $10^{6}$ grid nodes, i.e., parameter configurations. For each configuration, a differential energy spectrum of GCR protons $J_{m}(E,\vec{q})$ was calculated. Each calculation consist in the Monte-Carlo simulation of about 6$\times\,10^{4}$ pseudo-particles, backwardly-propagated from Earth to the HP, and then re-weighted for their LIS, as to obtain the locally modulated GCR flux. In practice, $2\times 10^{3}$ pseudo-particles were simulated for each of the 60 energy bins, ranging from 20 MeV to 200 GeV and following a logaritmically increasing step. Once the proton grid has established, the time-resolved GCR proton data has been utilized. From the data at a given epoch (time $t$), the corresponding quantities $A$, $\hat{\alpha}$, and $\hat{B_{0}}$ are determined in the first place. Then, for each GCR data set $J_{d}(E,t)$, a $\chi^{2}$ statistics was computed for all the modulated fluxes in the model array:

[TABLE]

the model $J_{m}$ energy flux at a given energy was computed by interpolated along on the grid, so that the resulting $\chi^{2}(\vec{q})$ function of Eq. 4 is a continuous function of the parameters (with the exception of the polarity parameter $A$, that takes two values). The errors are given by $\sigma^{2}(E_{i},t)=\sigma_{d}^{2}(E_{i},t)+\sigma_{mod}^{2}(E_{i},t)$, i.e., by the sum in quadrature of two contributions: errors in the experimental measurements, for the $i$-th energy bin, and theoretical uncertainties in the calculated flux $J_{m}$, modulated, related to LIS and to the Monte-Carlo generation statistics. Physically, this is due to adiabatic losses, that prevent a large fraction of low energy to reach the Earth. To reduce fluctuations, a high number of pseudo-particle trajectories has to be simulated. In particular, we generated about $2\times 10^{3}$ pseudo-particles per energy bin per run, which demanded a considerable CPU time. The evaluation of the flux fluctuations was done as follows. The ratio between the modulated flux and the corresponding LIS $J_{mod}/J_{\rm{LIS}}$, is approximately equal to the ratio $\sim N_{m}/N_{G}$, where $N_{m}$ is the number of pseudo-particles that reach the HP with energy $E$, and $N_{G}$ is the number of pseudo-particles generated at the same energy. Since the propagation process is stochastic, we expect that the relative error of the modulated flux decreases as ${\delta}J_{mod}/J_{mod}=1/\sqrt{N_{m}}$ with $N_{m}=\left(J_{mod}/J_{LIS}\right)N_{G}$. The resulting errors amount to $\sim$ 10$\div$20 % for GCR energy of 20 MeV, and become constant at the level of $\sim\,{2\%}$ at GCR energy above a few GeVs. The determination of the remaining GCR diffusion parameters ($K_{0}(t)$, $a(t)$ and $b(t)$) proceeds as follows. Given a data set $J_{d}(E,t)$, for every parameter $q$ = $K_{0}(t)$, $a(t)$ and $b(t)$, the corresponding $\chi^{2}(q)$ distribution is evaluated, its absolute minimum $\chi^{2}$ is found, and the corresponding parameter $q_{\rm best}$ is then taken as best-fit. To find the minimum two methods were implemented, both based on a multidimensional grid interpolation, and the associated systematic errors were accounted for.

4 Results and discussion

For each free parameter, our modified least-square minimization procedure returns a time-series of its best-fit values and their uncertainties. The results of the fitting are summarized in Fig. 2, showing the diffusion scaling $K_{0}$ (top panel) and its spectral indices $a$ (middle) and $b$ (bottom). The colors represent the data sets used to perform the fits, i.e., time-resolved proton measurements from AMS-02 (red full dots) and PAMELA (black full squares). It is important to note that the period covered by our study spans across a significant fraction of Solar Cycle, including the interesting phase of magnetic reversal, i.e., when the HMF polarity inversion takes place. We indicate by $T_{inv}$ the epochs of reversal, occurred between 2012 and 2013, which is shown in the figure by vertical blue line, along with its uncertainty band [34]. We can consider three relevant epochs: (i) $t\ll T_{\rm rev}$, before HMF reversal and with negative polarity ($A<$0), (ii) $t\gg T_{\rm rev}$, after reversal and under positive polarity ($A>$0), (iii) $t\sim T_{\rm rev}$ during the reversal, when the HMF is disorganized so that its polarity $A$ is not well defined. In the top panel of Fig. 2, the smoothed SSN is also shown, as dotted blue line, renormalized as to fit the scale. The SSN is a good proxy for solar activity and for studying its correlation with GCR modulation properties [2]. From the figure, it can be seen that the diffusion normalization parameter $K_{0}$ evolves with time, following a well defined pattern. As clear from the figure, the change in $K_{0}$ appears well anti-correlated with the smoothed SSN, i.e., with solar activity. In particular, GCR diffusion is faster during solar minimum, and slower during maximum. This behavior is at the basis of the time-dependent modulation phenomenon. In fact, a change in diffusion produces a higher (lower) intensity during epochs of solar minimum (maximum) and vice-versa. In the ideal limit of infinitely fast diffusion, the modulation effect would disappear (e.g., without Sun, or in the limit of high-energy GCRs). During the HMF reversal, it can also be seen that the best-fit $K_{0}$ curve has a little discontinuity, following the discontinuity of the HMF polarity that switches from $A=-1$ to $A=+1$, across $t=T_{\rm rev}$. In practice, magnetic reversal occurs gradually, across the region of the shaded band. In fact it involves the magnetic flip of the north/south poles and the ri-organization of the HMF. This process has a duration of about one year, as suggested by HMF observations [34] and found by AMS-02 data on the $e^{+}/e^{-}$ ratio evolution [5]. For the diffusion parameter values, one may expect a similar transition occurring during this phase. In Fig. 2, the best-fit $K_{0}$-values are shown under both polarities in the reversal phase. A smooth transition between the two fits is also shown.

Interestingly, we also found that the diffusion spectral indices $a$ and $b$ show clear temporal dependence. In our work, the index $a$ is related to the power spectrum in the $1/f$ range, while $b$ is related to the inertial range of the turbulent energy cascade of HMF. In both ranges, our parameters are in agreement with the measured slopes of the HMF power spectrum on Jan-Feb 2007 [18]. In most of numerical models of solar modulation, these parameters are usually assumed to be time-independent. Variations in these parameters imply unexpected changes in the HMF turbulence spectrum. This is also suggested in other recent works [35, 36, 37].

5 Acknowledgements

We acknowledge the support of ASI under agreement ASI-UniPG 2019-2-HH.0. In this work, we made use of the ASI-SSDC Cosmic-Ray Database at ASI, the SILSO/SIDC SSN Database at the Royal Observatory of Belgium, and real-time data from the Wilkox Solar Observatory. We thank Todd Hoeksema for support with the WSO solar data data.