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

The solar modulation effect of cosmic rays in the heliosphere is an energy-, time-, and particle-dependent phenomenon which arises from a combination of basic particle transport processes such as diffusion, convection, adiabatic cooling, and drift motion. Making use of a large collection of time-resolved cosmic-ray data from recent space missions, we construct a simple predictive model of solar modulation which depends on direct solar-physics inputs: the number of solar sunspots and the tilt angle of the heliospheric current sheet. Under this framework, we present calculations of cosmic-ray proton spectra, positron/electron and antiproton/proton ratios and their time dependence in connection with the evolving solar activity. We report evidence for a time-lag $\Delta{T}=8.1\pm\,1.2$ months, between solar activity data and cosmic-ray flux measurements in space, which reflects the dynamics of the formation of the modulation region. This result enables us to forecast the cosmic-ray flux near Earth well in advance by monitoring solar activity


Introduction
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(Potgieter , 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 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(Adriani et al. , 2016Bindi 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, DT , between solaractivity observations and effective conditions of the modulation region. Regarding this region as a bubble with a radiusd 100-120 au and radially flowing wind of speedṼ 300-700 km s −1 , we expect a time lag of the order of DT ∼0.5-1 year. Such a lag was reported in a number of empirical correlation studies between NM rates and solaractivity indices (Mavromichalaki & Petropoulos 1984;Nymmik 2000;Belov et al. 2005;Lantos 2005;Badruddin et al. 2007;Aslam & Badruddin 2012, 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.

Methodology
The transport of CRs in heliosphere is described by the Krymsky-Parker equation for the omnidirectional phase space density y ( ) r t p , , expressed as function of time t, momentum p, and position r (Potgieter 2013): The various terms represent convection with the solar wind of speed @ | | V 400 km s −1 , drift motion with an average speed v d , spatial diffusion with tensor K, and adiabatic momentum losses. We set up a minimal 2D description, with q = ( ) r r, , radius, and heliolatitude (Bobik et al. 2012). The parallel component of the diffusion tensor is GeV 3 22 0 , in units of cm 2 s −1 , where we have factorized an adimensional scaling factor, k 0 , of the order of unity. The perpendicular component is @  K K 0.02 (Giacalone & Jokipii 1999). The regular magnetic field (HMF) is modeled with the usual Parker 2.866 10 6 rad s −1 is the angular rotation of the Sun, @ B r 0 0 2 3.4 nT au 2 sets the field intensity at = r 1 0 au, and =  A 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 q v are both proportional to b r p qA 2 , so that the sign of 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 85 ts 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 ts 3, the TS compression ratio, while V and K decrease by s 1 ts . These quantities are changed gradually, between their pre-and post-shock values, across a TS thickness of = L 1.2 au. For instance, the wind speed is modeled as = - . 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 fluxes J 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. 2015aAguilar et al. , 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 p y p dp , which we express in units of GeV −1 m −2 s −1 sr −1 . Due to drift, CR particles and antiparticles sample different parts of the heliosphere. During negative polarity, protons and positrons ( > q 0) travel near the HCS and drift across that plane, while electrons and antiprotons ( < q 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 k 0 . Tilt-angle measurements â 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). 3 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. 2012Boschini et al. 2017). Here, we adopt a simple two-coefficient relation where theŜ-function smoothly interpolates the monthly series of SSNs provided by the SIDC-Royal Observatory of Belgium (Clette & Lefèvre 2016). 4 The temporal behavior of these quantities, shown in Figure   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 DT 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). Our model is specified by three free parameters only, a b , , and DT , that we constrain using a large amount of data. We use monthly resolved proton data from the PAMELA experiment (Adriani et al. 2013)   The quantity s 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 â ( ) t from the L-R model discrepancy, and from the different smoothing procedures on the HCS drift (see Burger & Potgieter 1989); (iii) uncertainties onˆ( ) 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, k 0 -values ranging from 0.1 to 6 with 0.1 of average (but nonequidistant) resolution steps, E-values consisting of 40 logspaced points between 0.08 GeV and 50 GeV, and 7 TS position values, r 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 ofp and  e particles have been carried out.

Results
The global fit has been performed on 3993 proton data points collected between 2000 and 2012 (in 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 parameterŝ a andb, 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 t ts -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 . We estimated a±0.8 month uncertainty for the inferred DT value, which is accounted for as a systematic error in the quoted results.
To inspect our results in greater depth, we have also performed a time series of fits to the single energy spectraĴ j by directly using the diffusion scaling as a free parameter. This provided a time series of 62 k j 0 -values corresponding to various epochs t j . Fits have been done after fixing D º D =  T T 8.1 months and D º T 0, respectively. The reconstructed time evolution of k ( ) t 0 is shown in Figure 4(a). At this point one can inspect the correlation between k 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 k j 0 and the "delayed" SSN is well described by a simple universal relation between k 0 and SSN, shown as a yellow line together with its uncertainty band. The resulting linear correlation coefficient is r = -D 0.89 T . In contrast, under the standard D = 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 k (ˆ) S 0 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 r = -0.64 0 . 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 k 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 k 0 -SSN relations can be captured by a simple universal function. Similar results are also found from the correlations between k 0 and the tilt angle. The two quantities are anticorrelated, with coefficients r = -D 0.78 T and r = -0.57 0 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 1 GeV). Finally, our model is used to predict the time evolution of antimatter-to-matter ratios such as e + /e − andp p. Calculations are shown in Figure 5 under both scenarios DT = 8.1 months and DT = 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(Abe et al. , 2012, BESS TeV (Haino et al. 2005), and AMS (Aguilar et al. 2016). For the e + /e − ratio, we note that calculations within D = 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  (Abe et al. 2008(Abe et al. , 2012Adriani et al. 2013;Aguilar et al. 2015a;Kühl et al. 2016). Calculations for D = T 0 are shown as thin dashed lines. The shaded bars indicate the magnetic reversals of the Sun's polarity (Sun et al. 2015).  et al. 2008, 2012Adriani et al. 2013;Kühl et al. 2016). . Scatter plots representing the best-fit k 0 -values from various data sets (Abe et al. 2008(Abe et al. , 2012Adriani et al. 2013;Kühl et al. 2016) vs. time (a) and vs. SSNs after accounting (b) and not accounting (c) for the time lag. The yellow lines represent the adopted parameterization and its uncertainty band. e + /e − (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 ofp 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,p, + e , ande , under both polarity conditions and across the reversal.

Conclusions
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 solaractivity 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 DT between solar data and the condition of the heliosphere. Using a large compilation of CR proton data we found D =  T 8.1 1.2 months, which is in agreement with basic expectations and with recent NMbased 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−DT . 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 DT 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 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, D  T , depending on the sign of qA (and, in particular,  D D -+ T 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   (Haino et al. 2005;Abe et al. 2008Abe et al. , 2012Clem & Evenson 2009;Adriani et al. 2016;Aguilar et al. 2016). The shaded bars indicate the magnetic reversals of the Sun's polarity (Sun et al. 2015).