Evidence for a Time Lag in Solar Modulation of Galactic Cosmic Rays

In recent years, new-generation experiments of cosmic-ray (CR) detection have reached an unmatched level of precision that is bringing transformative advances in astroparticle physics (Grenier et al. 2015; Amato & Blasi 2017). Along with calculations of CR propagation in the galaxy, the interpretation of the data requires a detailed modeling of the so-called solar modulation effect. Solar modulation is experienced by all CR particles that enter the heliosphere to reach our detectors near Earth. Inside the heliosphere, CRs travel through a turbulent magnetized plasma, the solar wind, which significantly reshapes their energy spectra. This effect is known to change with time, in connection with the quasi-periodical 11 year evolution of the solar activity and to provoke different effects on CR particles and antiparticles (Potgieter 2013, 2014).

Observationally, an inverse relationship between solar activity (often monitored by the number of solar sunspots) and CR flux intensity has been known about for a long time. The effect of solar modulation in the low-energy CR spectra ($E\,\lesssim$ GeV) is measured by several experiments (Valdés-Galicia & González 2016; Bindi et al. 2017). Theoretically, the paradigm of CR propagation in the heliosphere has been developed soundly over the past decades. Solar modulation is caused by a combination of basic transport processes such as diffusion, convection, adiabatic cooling, and drift motion, yet the underlying physical mechanisms and their associated parameters remain under active investigation. In CR astrophysics, solar modulation is often treated using simplified models where the key parameter(s) are degenerated with the parameters of CR propagation in the galaxy. The development of solar modulation models can be ensured by two crucial factors: (i) the precise knowledge of the interstellar spectra (LIS) of CRs outside the heliosphere, and (ii) the availability of time series of CR data on different species. Recent accomplishments from strategic space missions have enabled us to make significant progress in this field. The entrance of Voyager 1 into interstellar space provided us with the very first LIS data on CR protons and electrons (Stone et al. 2013; Cummings et al. 2016). Long-duration space experiments PAMELA (in orbit since 2006) and AMS (since 2011) have started releasing a continuous stream of time-resolved data on CR particles and antiparticles (Adriani et al. 2013, 2016; Bindi et al. 2017). These measurements add to a large wealth of low-energy data collected in the past decades by space missions CRIS/ACE (Wiedenbeck et al. 2009), IMP-7/8 (Garcia-Munoz et al. 1987), Ulysses (Heber et al. 2009), EPHIN/SOHO (Kühl et al. 2016), and from ground data provided continuously by the neutron monitor (NM) worldwide network (Mavromichalaki et al. 2011; Steigies 2015).

In this Letter, we present new calculations of CR fluxes near Earth that account for the dynamics of CR modulation in the expanding heliosphere. Using a large collection of modulated and interstellar CR data collected in space, we have constructed a predictive and measurement-validated model of solar modulation that depends only on direct solar-activity observables: the sunspot number (SSN) and the tilt angle of the heliospheric current sheet (HCS). In our calculations, all relevant processes are formally expressed using "retarded" relations in order to account for a time lag, ${\rm{\Delta }}T$, between solar-activity observations and effective conditions of the modulation region. Regarding this region as a bubble with a radius $d\sim$ 100–120 au and radially flowing wind of speed $V\sim$ 300–700 km s⁻¹, we expect a time lag of the order of ${\rm{\Delta }}T$ ~ 0.5–1 year. Such a lag was reported in a number of empirical correlation studies between NM rates and solar-activity indices (Mavromichalaki & Petropoulos 1984; Nymmik 2000; Belov et al. 2005; Lantos 2005; Badruddin et al. 2007; Aslam & Badruddin 2012, 2015), although it is usually ignored in CR modulation models. As we will show, recent direct measurements of CRs reveal the existence of an eight-month lag. This result sets the timescale of the changing conditions of the heliosphere, enabling us to forecast the near-Earth CR flux well in advance. Along with the intrinsic interest in plasma or solar astrophysics, our results address a prerequisite for modeling space weather effects, which is an increasing concern for space missions and air travelers.

The transport of CRs in heliosphere is described by the Krymsky–Parker equation for the omnidirectional phase space density $\psi (t,p,{\boldsymbol{r}})$ expressed as function of time t, momentum p, and position ${\boldsymbol{r}}$ (Potgieter 2013):

$$\begin{eqnarray}&&\displaystyle \frac{\partial \psi }{\partial t}=-({\boldsymbol{V}}+{{\boldsymbol{v}}}_{d})\cdot {\rm{\nabla }}\psi +{\rm{\nabla }}\cdot ({\boldsymbol{K}}\cdot {\rm{\nabla }}\psi )+\displaystyle \frac{1}{3}({\rm{\nabla }}\cdot {\boldsymbol{V}})\displaystyle \frac{\partial \psi }{\partial \mathrm{ln}p}.\end{eqnarray} \tag{ 1 }$$

The various terms represent convection with the solar wind of speed $| {\boldsymbol{V}}| \,\cong$ 400 km s⁻¹, drift motion with an average speed ${{\boldsymbol{v}}}_{d}$, spatial diffusion with tensor ${\boldsymbol{K}}$, and adiabatic momentum losses. We set up a minimal 2D description, with ${\boldsymbol{r}}=(r,\theta )$, radius, and heliolatitude (Bobik et al. 2012). The parallel component of the diffusion tensor is ${K}_{\parallel }={\kappa }^{0}\tfrac{{10}^{22}\beta p/\mathrm{GeV}}{3B/{B}_{0}}$, in units of cm² s⁻¹, where we have factorized an adimensional scaling factor, ${\kappa }^{0}$, of the order of unity. The perpendicular component is ${K}_{\perp }\cong 0.02\,{K}_{\parallel }$ (Giacalone & Jokipii 1999). The regular magnetic field (HMF) is modeled with the usual Parker structure, $B=\tfrac{{AB}_{0}{r}_{0}^{2}}{{r}^{2}}\sqrt{1+{{\rm{\Gamma }}}^{2}}$, with ${\rm{\Gamma }}=({\rm{\Omega }}r/V)\sin \theta$, where ${\rm{\Omega }}=2.866\cdot \,{10}^{-6}$ rad s⁻¹ is the angular rotation of the Sun, ${B}_{0}{r}_{0}^{2}\,\cong$ 3.4 nT au² sets the field intensity at ${r}_{0}=1$ au, and $A=\pm \,1$ sets the magnetic polarity cycle of the Sun. The polarity is positive (negative) when the HMF points outward (inward) in the northern hemisphere. This model accounts for gradient and curvature drift effects. In particular, drift is important across the wavy layer of the HCS, i.e., the surface where polarity changes from north to south. The angular extension of the HCS is described by the tilt angle α. The drift velocity components v_r and ${v}_{\theta }$ are both proportional to $\tfrac{2r\beta p}{{qA}}$, so that the sign of ${{\boldsymbol{v}}}_{d}$ depends on the product qA (Potgieter 2014). We account for the presence of the termination shock (TS), placed at a default distance ${r}_{\mathrm{ts}}\cong 85$ au, and for its time dependence over the solar cycle (Manuel et al. 2015; Vos & Potgieter 2016). Beyond the TS, plasma density and B increase by a factor ${s}_{\mathrm{ts}}\,=$ 3, the TS compression ratio, while V and K decrease by $1/{s}_{\mathrm{ts}}$. These quantities are changed gradually, between their pre- and post-shock values, across a TS thickness of $L=1.2\,\mathrm{au}$. For instance, the wind speed is modeled as $V(r)=\tfrac{{s}_{\mathrm{ts}}+1}{2{s}_{\mathrm{ts}}}{V}_{0}-\tfrac{{s}_{\mathrm{ts}}-1}{2{s}_{\mathrm{ts}}}{V}_{0}{\tan }^{-1}(\tfrac{r-{r}_{\mathrm{ts}}}{L})$ (from Equation (3) in Langner et al. 2003), with ${V}_{0}\,\cong$ 400 km s⁻¹. A similar procedure, based on an $\arctan$-like function of θ, is also applied to regularize the drift velocity near the HCS (from Equation (22) and below in Bobik et al. 2012). The modulation boundaries are set at the heliopause (HP) at ${r}_{\mathrm{hp}}=122\,\mathrm{au}$. The CR distribution in the heliosphere is computed using the stochastic differential equation approach of Kappl (2016). This method consists of the backward-in-time propagation of a large number of pseudo-particles from Earth to the HP boundaries (Alanko-Huotari et al. 2007; Strauss et al. 2012; Raath et al. 2016). For a given particle type, the steady-state solution of Equation (1) ($\partial \psi /\partial t=0$) is obtained by sampling. The LIS fluxes ${J}^{\mathrm{IS}}$ are calculated within an improved model of CR acceleration and propagation (Tomassetti 2015; Tomassetti & Donato 2015; Feng et al. 2016). Proton and electron LISs are well constrained by the Voyager 1 and AMS data (Cummings et al. 2016; Aguilar et al. 2015a, 2015b). Antiproton and positron LISs rely on secondary production calculations and are affected by larger uncertainties (Tomassetti 2012; Feng et al. 2016). The CR number density is given by $N(p){dp}\,=4\pi {p}^{2}\psi {dp}$. The resulting flux is $J=\tfrac{\beta c}{4\pi }N$, which we express in units of GeV⁻¹ m⁻² s⁻¹ sr⁻¹. Due to drift, CR particles and antiparticles sample different parts of the heliosphere. During negative polarity, protons and positrons ($q\gt 0$) travel near the HCS and drift across that plane, while electrons and antiprotons ($q\lt 0$) propagate preferentially through the polar regions, with faster diffusion and smaller losses. The role of positive and negative particles interchanges with polarity. Along with qA, the interplay of the various physics processes depends on the levels of HCS waviness and HMF irregularities that are contained in α and ${\kappa }^{0}$. Tilt-angle measurements $\hat{\alpha }$ are provided on a 10 day basis from the Wilcox Solar Observatory using two reconstruction procedures: the "classic" model-L and the improved model-R (Hoeksema 1995; Ferreira & Potgieter 2004).³ A smoothed interpolation of model-R is adopted as the default. The diffusion coefficient is also time-dependent, due to changes in the HMF turbulence (Manuel et al. 2014). A basic diagnostic for the HMF turbulence is the manifestation of solar sunspots (Plunian et al. 2009, Bobik et al. 2012; Boschini et al. 2017). Here, we adopt a simple two-coefficient relation ${\kappa }^{0}\equiv a+b\mathrm{log}(\hat{S})$, where the $\hat{S}$-function smoothly interpolates the monthly series of SSNs provided by the SIDC—Royal Observatory of Belgium (Clette & Lefèvre 2016).⁴

The temporal behavior of these quantities, shown in Figure 1 from 2000 to 2017, is at the basis of the time-dependent nature of solar modulation. In practice, the problem is modeled under a quasi-steady fashion, i.e., by providing a time series of steady-state solutions corresponding to a time series of input parameters. While this approach is justified by the different timescales of CR transport and solar activity, the correspondence between solar-activity observations and effective conditions of the heliosphere may require that dynamics considerations come into play. A finite amount of time is needed, in fact, for the properties observed in the solar corona to be transported in the outer heliosphere by the plasma. To tackle this issue, we introduce a parameter ${\rm{\Delta }}T$ in our calculation, describing a time lag between the solar-activity indices of Figure 1 and the medium properties of the modulation region, i.e., the spatial region effectively sampled by CRs. As discussed, evidence for a lag of ~6–12 months has been reported in NM-based empirical studies (Mavromichalaki & Petropoulos 1984; Nymmik 2000; Belov et al. 2005; Lantos 2005; Badruddin et al. 2007; Chowdhury et al. 2016; Mishra & Mishra 2016).

Zoom In Zoom Out Reset image size

Our model is specified by three free parameters only, $a,b$, and ${\rm{\Delta }}T$, that we constrain using a large amount of data. We use monthly resolved proton data from the PAMELA experiment (Adriani et al. 2013) collected in solar minimum conditions between 2006 July and 2010 January, and new data from the EPHIN/SOHO space detector (Kühl et al. 2016), resolved yearly between 2000 and 2016. We also include data from the BESS Polar-I (Polar-II) mission from 2004 December 13 to 21 (2007 December 23 to 2008 January 16; Abe et al. 2008, 2012). These measurements are given in terms of time series of energy spectra, ${\hat{J}}_{j,k}=\hat{J}({t}_{j},{E}_{k})$, each one representing a snapshot of the CR flux near Earth at epoch t_j. Calculations $J({t}_{j},{E}_{k})$ are performed using retarded functions of the physics inputs ${\alpha }_{j}=\hat{\alpha }({t}_{j}-{\rm{\Delta }}T)$ and ${\kappa }_{j}^{0}=a+b\mathrm{log}(\hat{S}({t}_{j}-{\rm{\Delta }}T))$. We then build a global ${\chi }^{2}$-estimator:

$$\begin{eqnarray}&&{\chi }^{2}(a,b,{\rm{\Delta }}T)=\displaystyle \sum _{j,k}{\left[\displaystyle \frac{J({t}_{j},{E}_{k};a,b,{\rm{\Delta }}T)-{\hat{J}}_{j,k}}{{\sigma }_{j,k}}\right]}^{2}.\end{eqnarray} \tag{ 2 }$$

The quantity ${\sigma }_{j,k}$ includes experimental errors in the data and model uncertainties due to finite statistics of the simulation. The following sources of errors are also accounted for: (i) LIS uncertainties from the constraints provided by Voyager 1 and AMS data; (ii) uncertainties on $\hat{\alpha }(t)$ from the L–R model discrepancy, and from the different smoothing procedures on the HCS drift (see Burger & Potgieter 1989); (iii) uncertainties on $\hat{S}(t)$ from the smoothed SSN variance (see Figure 1); and (iv) uncertainties on the TS strength and position, including its time variations. The free parameters are estimated by means of standard minimization techniques. In practice, for both polarities, calculations are performed by interpolation over 4D grids of α-values ranging from 0° to 80° with 1° steps, ${\kappa }^{0}$-values ranging from 0.1 to 6 with 0.1 of average (but nonequidistant) resolution steps, E-values consisting of 40 log-spaced points between 0.08 GeV and 50 GeV, and 7 TS position values, ${r}_{\mathrm{ts}}$, ranging from 70 to 100 au. This task required the simulation of about 14 billion proton trajectories, corresponding to several months of CPU time. Along with protons, simulations of $\bar{p}$ and ${e}^{\pm }$ particles have been carried out.

The global fit has been performed on 3993 proton data points collected between 2000 and 2012 (in $A\lt 0$ conditions) at kinetic energy between 0.08 and 50 GeV. The best-fit parameters are $\hat{a}=4.07\pm \,0.95,\hat{b}=-1.39\pm 0.34$, and $\widehat{{\rm{\Delta }}T}=8.1\,\pm 1.2$ months, giving ${\chi }^{2}/{df}=2643/3990$. The fit was repeated after fixing ${\rm{\Delta }}T\equiv \,0$, i.e., under the conventional "unretarded" scenario, returning $\hat{a}=3.52\pm 0.91$, $\hat{b}=-1.16\pm 0.28$, and ${\chi }^{2}/{df}=4979/3991$. Results are shown in Figures 2 and 3, which illustrate the time and energy dependence of calculations in comparison with the data. The fits give satisfactory results at all energies and epochs. From Figure 2, it can be seen that the model reproduces very well the time evolution of the proton flux, at $E=1\mbox{--}1.5\,\mathrm{GeV}$. In the $A\gt 0$ period, i.e., from the 2013 reversal and up to 2018, the proton flux is predicted to increase with time. The post-reversal time evolution of the proton flux is currently being measured by the AMS experiment (Bindi et al. 2017). It is also clear from the figure that the retarded scenario (with ${\rm{\Delta }}T$ ≡ 8.1 months; thick red line) allows for a much better description of the time evolution of the proton flux. The results illustrated here are obtained using a static TS position at the default value of ${r}_{\mathrm{ts}}=85$ au. Extreme values of ${r}_{\mathrm{ts}}^{\pm }=95/75\,\mathrm{au}$ have also been tested. They lead to a ±10% shift in the near-Earth flux intensity, for fixed input parameters, in agreement with Manuel et al. (2015). However, this global shift is re-absorbed by different best-fit values for the parameters $\hat{a}$ and $\hat{b}$, while consistent time-lag values are inferred in all cases. Finally, we have tested the case of a dynamic TS by means of a time-dependent ${r}_{\mathrm{ts}}(t)$-function. The TS time dependence is poorly known and must rely on calculations. We adopt a reference function based on Richardson & Wang (2011, 2012) with a rescaling in order to match the average TS positions determined from calculations based on Voyager 1 and 2 (Washimi et al. 2011; Webber & Intriligator 2011; Provornikova et al. 2014). The robustness of the results is studied using a parametric form of the type ${c}_{1}{r}_{\mathrm{ts}}(t+{c}_{2})$ with ${c}_{1}=1\pm \,0.12$ and ${c}_{2}=\pm \,3\,\mathrm{years}$, along with ${s}_{\mathrm{ts}}=3\pm 1$. We estimated a ±0.8 month uncertainty for the inferred ${\rm{\Delta }}T$ value, which is accounted for as a systematic error in the quoted results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To inspect our results in greater depth, we have also performed a time series of fits to the single energy spectra ${\hat{J}}_{j}$ by directly using the diffusion scaling as a free parameter. This provided a time series of 62 ${\hat{\kappa }}_{j}^{0}$-values corresponding to various epochs t_j. Fits have been done after fixing ${\rm{\Delta }}T\equiv \,\widehat{{\rm{\Delta }}T}=8.1$ months and ${\rm{\Delta }}T\,\equiv$ 0, respectively. The reconstructed time evolution of ${\hat{\kappa }}^{0}(t)$ is shown in Figure 4(a). At this point one can inspect the correlation between ${\hat{\kappa }}^{0}$ and SSN, which is shown in Figures 4(b) and 4(c) for the two scenarios, i.e., with and without time lag. In the first case, the relation between ${\hat{\kappa }}_{j}^{0}$ and the "delayed" SSN is well described by a simple universal relation between ${\kappa }^{0}$ and SSN, shown as a yellow line together with its uncertainty band. The resulting linear correlation coefficient is ${\rho }_{{\rm{\Delta }}T}=-0.89$. In contrast, under the standard ${\rm{\Delta }}T=0$ scenario, a simple one-to-one correspondence between diffusion scaling and SSN observations cannot be established along the entire variation range of SSNs. In particular, any ${\kappa }^{0}(\hat{S})$ relation would fail in describing the solar minimum of 2009, i.e., in proximity to the transition between descending and ascending phases of solar activity. In this case, the correlation coefficient is ${\rho }_{0}=-0.64$. Similar considerations can be drawn from the comparison of the PAMELA data in Figure 2 with the SSN evolution of Figure 1(b): the CR proton flux keeps increasing for a few months after the solar activity minimum. These findings explain why other authors, when proposing simple relations between ${\kappa }^{0}$ and SSN had to adopt different coefficients for descending and ascending phases (Bobik et al. 2012; Boschini et al. 2017). Remarkably, this problem is naturally resolved in our model: once the time lag is properly accounted for, the ${\kappa }^{0}$–SSN relations can be captured by a simple universal function. Similar results are also found from the correlations between ${\kappa }^{0}$ and the tilt angle. The two quantities are anticorrelated, with coefficients ${\rho }_{{\rm{\Delta }}T}=-0.78$ and ${\rho }_{0}=-0.57$ corresponding to the two scenarios. The inferred lag is also found to be highly stable to changes in the fitting energy range (moving the minimal energy from 0.08 to 8 GeV) and in the input spectrum (changing the LIS intensity within a factor of two at $E\,\approx$ 1 GeV).

Zoom In Zoom Out Reset image size

Finally, our model is used to predict the time evolution of antimatter-to-matter ratios such as e⁺/e⁻ and $\bar{p}/p$. Calculations are shown in Figure 5 under both scenarios ${\rm{\Delta }}T$ = 8.1 months and ${\rm{\Delta }}T$ = 0. Measurements of the relative variation of the ratios are from PAMELA (Adriani et al. 2016), AESOP (Clem & Evenson 2009), BESS Polar (Abe et al. 2008, 2012), BESS TeV (Haino et al. 2005), and AMS (Aguilar et al. 2016). For the e⁺/e⁻ ratio, we note that calculations within ${\rm{\Delta }}T=8.1$ months are favored, although the data do not permit a resolute discrimination. Across the magnetic reversal, shown in the figure as shaded bars, a remarkable increase (decrease) of the e⁺/e⁻ ($\bar{p}/p$) ratio is predicted. It should be noted, however, that the dynamics of the transition cannot be modeled during reversal because the HMF polarity is not well defined. One may also expect that drift motion is ineffective in this phase. Nonetheless, a rise in the e⁺/e⁻ ratio profile has been detected in the PAMELA data (Adriani et al. 2016), and this rise is found to occur a few months after completion of the Sun's polarity reversal. A corresponding decrease of $\bar{p}/p$ ratio is predicted to occur at the same timescale. Crucial tests can be performed by AMS measurements of the two ratios, or even better, by measurements of individual particle fluxes for $p$, $\bar{p}$, ${e}^{+}$, and ${e}^{-}$, under both polarity conditions and across the reversal.

Zoom In Zoom Out Reset image size

New experiments of CR detection in space have determined the evolution of CR fluxes near Earth with an unmatched time resolution. These data allow for detailed studies of the solar modulation effect and its dynamical connection with the evolving solar activity. In this Letter, we have reported new calculations of CR modulation based on a physically consistent numerical model that accounts for particle diffusion, drift, convection, and adiabatic cooling. We have adopted a simple formulation where the time-dependent physics inputs of the model consist in SSN and HCS tilt angle. We have shown that this model reproduces well the time evolution of the CR proton spectra measured by PAMELA, EPHIN/SOHO, and BESS experiments. Our model is highly predictive once the correspondence between modulation parameters and solar-activity indices is established. Our study revealed an interesting aspect of the dynamics of CR modulation in the expanding wind, that is, the presence of time lag ${\rm{\Delta }}T$ between solar data and the condition of the heliosphere. Using a large compilation of CR proton data we found ${\rm{\Delta }}T=8.1\pm 1.2$ months, which is in agreement with basic expectations and with recent NM-based analysis (Aslam & Badruddin 2015; Chowdhury et al. 2016; Mishra & Mishra 2016). An interesting consequence of this result is that the galactic CR flux at Earth can be predicted, at any epoch t, using solar-activity indices observed at the time t−${\rm{\Delta }}T$. This result is of great interest for real-time space weather forecast, which is an important concern for human spaceflight.

In this work, the parameter ${\rm{\Delta }}T$ has been determined using CR protons during negative polarity. This parameter has to be viewed as an effective quantity representing the average of several CR trajectories in the heliosphere during $A\lt 0$ conditions. Further elaborations may include the use of NM data, for larger observation periods, or the accounting for a latitudinal dependence in the wind profile or in the diffusion coefficient. Since CR particles and antiparticles sample different regions of the heliosphere, we expect slightly different time lags, ${\rm{\Delta }}{T}_{\pm }$, depending on the sign of qA (and, in particular, ${\rm{\Delta }}{T}_{-}\lesssim \,{\rm{\Delta }}{T}_{+}$). With the precision of the existing data, we were unable to test this hypothesis. A detailed re-analysis of our model, in this direction, will be possible after the release of the monthly resolved data from AMS on CR particle and antiparticle fluxes.

All calculations presented here will be made available and kept updated at the SSDC Cosmic-Ray Database hosted by the Space Science Data Center of the Italian Space Agency (Di Felice et al. 2017).⁵

We thank our colleagues of the AMS Collaboration for valuable discussions. N.T. and B.B. acknowledge the European Commission for support under the H2020-MSCA-IF-2015 action, grant No. 707543 MAtISSE—"Multichannel Investigation of Solar Modulation Effects in Galactic Cosmic Rays." M.O. acknowledges support from FCT—Portugal, Fundação para a Ciência e a Tecnologia, under grant SFRH/BD/104462/2014.