Numerical modeling of cosmic rays in the heliosphere: Analysis of proton data from AMS-02 and PAMELA

Galactic cosmic rays (CRs) inside the heliosphere are affected by solar modulation. To investigate this phenomenon and its underlying physical mechanisms, we have performed a data-driven analysis of the temporal dependence of the CR proton flux over the solar cycle. The modulation effect was modeled by means of stochastic simulations of cosmic particles in the heliosphere. The model were constrained using measurements of CR protons made by AMS-02 and PAMELA experiments on monthly basis from 2006 to 2017. With a global statistical analysis of these data, we have determined the key model parameters governing CR diffusion, its dependence on the particle rigidity, and its evolution over the solar cycle. Our results span over epochs of solar minimum, solar maximum, as well as epochs with magnetic reversal and opposite polarities. Along with the evolution of the CR transport parameters, we study their relationship with solar activity proxies and interplanetary parameters. We find that the rigidity dependence of the parallel mean free path of CR diffusion shows a remarkable time dependence, indicating a long-term variability in the interplanetary turbulence that interchanges across different regimes over the solar cycle. The evolution of the diffusion parameters show a delayed correlation with solar activity proxies, reflecting the dynamics of the heliospheric plasma, and distinct dependencies for opposite states of magnetic polarity, reflecting the influence of charge-sign dependent drift in the CR modulation.


I. INTRODUCTION
Galactic cosmic rays (CR) are high-energy charged particles produced by astrophysical sources, distributed in our galaxy, which travel through the interstellar medium and finally arrive at the boundary of the nearby region to Earth where the Sun's activity dominates: the so called heliosphere. When entering the heliosphere, CRs travel against the expanding solar wind (SW) and interact with the turbulent heliospheric magnetic field (HMF) [1]. They are subjected to basic transport processes such as convection, diffusion and adiabatic energy losses. They are also subjected to the gradient-curvature drifts in the large-scale HMF and to the effects of the heliospheric current sheet (HCS). Magnetic drift depends on the charge-sign of the particles and on the polarity of the HMF; CRs drift along different trajectories according to the polarity of the HMF. The cumulative effects of these processes are behind the so-called solar modulation phenomenon of CRs, that is, the modification of the energy spectra of CRs in the heliosphere, which is driven by the Sun's magnetic activity. Due to solar modulation, the CR flux observed at Earth is significantly different from that in interstellar space, known as Local Interstellar Spectrum (LIS). Solar modulation depends on the CR particle species, its energy, and its charge sign. It is also a time-dependent and space-dependent phenomenon, i.e., it depends on where and when the CR flux is measured inside the heliosphere. The solar modulation effect decreases with increasing energy of the CR particles. With the precision of the new CR data from AMS-02, the modulation effect is appreciable at kinetic energies up to dozens GeV. Solar activity shows a 11-year cycle, from its minimum when the Sun is quiet and the CR intensity is at its largest, to its maximum of solar activity when the CR flux is minimum. The intensity and the energy spectra of the CR flux are therefore anticorrelated with solar activity, in relation with its varying proxies such as the number of sunspot (SSN) or the tilt angle of the solar magnetic axis with respect to the rotation axis α [2][3][4]. Along with the 11-year solar cycle, the HMF polarity shows a remarkable 22-year periodicity, with the magnetic reversal occurring during each maximum of solar activity. This periodicity is important for CR modulation, and in particular to study the effects of particle drifts in the large-scale HMF.
Since CR modulation is a manifestation of the CR propagation through the heliosphere, CR data can be used to investigate the fundamental physics processes governing the transport of charged particles through the heliospheric plasma. In particular, precise measurements of the energy and time dependence of the CR fluxes may help to disentangle the interplay of the different physics mechanisms at work. In this respect, the physical understanding of CR modulation in the heliosphere is one of the main objectives of many theoretical and observational studies [5][6][7][8]. Besides, modeling the CR modulation is essential for the search of new physics signatures in the fluxes of CR antimatter such as positrons or antiprotons. An antimatter excess in CRs may suggest the occurrence of dark matter annihilation processes or the existence of new astrophysical sources of antimat-ter. Since the low-energy spectra of CRs are influenced by solar modulation, any interpretation about the origin of antiparticles requires an accurate modeling of the charge-sign and energy dependent effects of CR modulation [9]. Understanding the evolution of the CR fluxes in the heliosphere is also important for assessing the radiation hazard of astronauts, electronics, and communication systems for low-Earth-orbit satellites or deep space missions [10,11]. In fact, the Galactic CR flux constitutes a significant dose of ionizing radiation for human bodies and electronics, and thus an accurate knowledge of the temporal and spatial variation of the CR in the heliosphere will reduce the uncertainties in the radiation dose evaluation [12]. An important challenge, in this context, is to establish a predictive model for solar modulation that is able to forecast the CR flux evolution using solar activity proxies.
From the observational point of view, a substantial progress has been made with the new measurements of the proton flux from the Alpha Magnetic Spectrometer (AMS-02) experiment in the International Space Station [13,14] and the PAMELA mission onboard the Resurs-DK1 satellite [15,16], along with the data provided by the Voyager-1 spacecraft beyond the heliosphere [17]. In particular, AMS-02 and PAMELA have recently released accurate measurements of CR proton spectra over Bartels' rotation basis (BR, 27 days), over extended energy range and for extended time periods, covering the long solar minimum of 2006-2009 (cycle 23/24), the ascending phase of cycle 24, the solar maximum and HMF reversal of 2013-2014, and the subsequent descending phase towards the new minimum until May 2017. Therefore, the data allows for the study of the CR propagation in the heliosphere under very different conditions of solar activity and epochs of opposite HMF polarities, which may bring a substantial advance in the understanding of the solar modulation phenomenon.
In this paper, we present a data-driven analysis of the temporal dependence of the flux of CR protons, which constitute the most abundant species of the Galactic cosmic radiation. The analysis has been conducted using a stochastic model of CR propagation, i.e., a Monte Carlo based approach in which the solar modulation effect is computed by statistical sampling. Using the recent timeand energy-resolved measurements of CR proton fluxes on BR basis, by means of a procedure of statistical inference, we determine the temporal and rigidity dependencies of the mean free path of CRs propagating through the heliosphere, along with the corresponding uncertainties. The rest of this paper is organized as follows. In Sect. II, we describe in details the numerical implementation of the CR modulation model, which is based on known and conventional mechanisms of particle transport in the heliosphere. In Sect. III we present the procedure for the data-driven determination of the key model parameters and their uncertainty, which is based on a grid sampling over a multidimensional parameter space. In Sect. IV we present the fit results and discuss their interpretation, in terms of physical mechanisms of CR transport, in relation with the properties of heliospheric environment or with known proxies of solar activity. We then conclude, in Sect. V, with a summary of our study and a discussion on its future developments.

II. THE NUMERICAL MODEL
To get a realistic description of CR modulation phenomenon, one needs to capture the essential features of CR transport in the heliosphere. The diffusive propagation of the charged particles in the turbulent heliospheric plasma is described by the Parker's equation [18]: The equation, along with its boundary conditions, describes the evolution of the distribution function f (t, r, R) for a given particle species, where t is the time, and R is the particle rigidity, i.e., the momentum per charge units R=p/Z. In this paper, we will focus on cosmic protons, so that R ≡ p. The quantity K is the drift-diffusion tensor of the CR particles in the turbulent HMF of the heliosphere. Because of the complexity of the transport equation, analytical solutions can be found only for very simplified situations such as in the Force-Field or the Diffusion-Convection approximations [19,20]. The full solution of Eq.(1) can be obtained numerically. Here we employ the stochastic method, that has become widely implemented in recent years thanks to the enormous progress in computing speed and resources [8,21,22]. The method consists of transforming the Parker's equation into a set of Stochastic Differential Equations (SDE) and then using Monte Carlo simulations to sample the solution, i.e., the differential CR intensity for a given species, at a given position in heliosphere [23,24].
In general, the flux of CRs inside the heliosphere is time-dependent, reflecting the varying conditions of the medium over which they propagate [25]. A common practice is to follow a quasi steady-state approximation where the time-dependent CR modulation is described as a succession of steady-state solutions (∂/∂t = 0) and the effective status of the heliospheric plasma during the CR propagation is defined in a suitable way. The approximate way of taking into account the varying status of the heliosphere during the CR propagation is described in Sect. II. Furthermore, in the SDE method, pseudoparticles are propagated backward in time from the Earth position to the heliospheric boundaries. The numerical engine for handling the Monte Carlo generation and the trajectory tracing is extracted from the publicly available code SolarProp [21]. Based on the SolarProp simulation framework, we have implemented a customized model that is described in the following.

A. The modulation region
The heliosphere is a dynamic void in the ISM generated by the SW and regulated by Sun's activity. The relevant boundary for the CR modulation phenomenon is the heliopause (HP), which separates the heliospheric plasma from the local ISM. The HP is usually modeled as a spherical structure of radius r HP ≈ 122 AU, where the Sun lies at its center. Within the heliosphere, the termination shock (TS) is located at r TS ∼ = 85 AU, while the Earth position is at r 0 ≡ 1 AU placed in the equatorial plane. The large-scale HMF -The outward flowing SW embeds a frozen-in HMF which is wounded up in a modified Parker spiral [26]. The ideal Parker's field is given by: where r and θ are helioradius and colatitude, B 0 is the HMF value at Earth position, A = ± 1 is the field polarity, and H is the Heavyside step function. The winding angle ψ of the field line is defined as tan ψ = Ω(r − r ) sin θ/V sw ; the angle Θ determines the position of the wavy HCS, given by Θ = π/2 + sin −1 [sin α sin (Ωr/V w )] [27]. Here the quantity Ω is the average equatorial rotation speed ≈ 2.73×10 −6 rad s −1 , α is the HCS tilt angle and r = 696.000 km is the radius of the Sun. The Parker's model overwounds by several degrees beyond the value of the winding angle ψ, determined by the model at the polar regions. To avoid this, one has to consider that solar wind disturbances and plasma waves propagating along the open field lines modify the magnetic field at the polar regions, so that it does not degenerate to a straight line along the polar axis. Here we adopt the modification of Jokipii & Kota [28]: where δ(θ) = 8.7 × 10 −5 / sin(θ) if 1.7 • < θ < 178.3 • and 3×10 −3 otherwise [29]. The winding angle ψ is then modified as: The term involving the dimensionless constant δ reflects the fact that the random field is equivalent to a small latitudinal component B θ ∼ δ(θ)r/r . In this way, modifications on HMF and winding angle are effective only near the polar regions, as shown in Fig. 1 where the two quantities are shown as function of colatitude. It is worth noticing that the definitions of B θ and δ(θ) imply ∇ · B = 0. Polarity and Tilt Angle -An important characteristic for the CR solar modulation is that the HMF follows a ∼ 22-year cycle, known as magnetic polarity cycle, characterized by a N/S reversal about every ∼ 11 years, during the maximum of solar activity. The period when B is directed outwards in the northern hemisphere of the Sun is known as positive polarity epoch(A > 0), while when it has the opposite direction are known as (A < 0) cycle. In practice the quantity A is a dichotomous variable that expresses the sign of B-field projection in the outward direction from the northern hemisphere, A ≡ B N /|B N | (or the inward projection of B S in the southern hemisphere). In practice it can be determined using observations of the polar HMF in proximity of the Sun (Sect. III B). The relevance of magnetic polarity in the context of solar modulation arises from CR drift motion: it can be seen (Sect. II B) that the equations ruling CR drift in the HMF depend upon the sign of the product between A andq = Q/|Q|, where Q is the CR electric charge. Thus, opposite drift directions are expected for oppositeqA conditions. A major co-rotating structure relevant to CR modulation is the HCS, which divides the HMF into hemispheres of opposite (N/S) polarity and where B = 0. Due to the tilt of the solar magnetic axis, the HCS is wavy. The level of the HCS wavyness changes with time and it is set by the tilt angle α(t). Typically, it varies from α ∼ 5 • during solar minimum to α ∼ 70 • during solar maximum. The tilt angle is reconstructed by the Wilcox Solar Observatory using two different models for the polar magnetic field: the socalled L-model and R-model. In this work the classical L-model reconstruction is used as default.
The Wind -The SW speed V sw is taken as radially directed outward. However, the wind field exhibits a radial, latitudinal, and temporal dependence, where the latter is related to the solar cycle. During periods of solar minimum, the flow becomes distinctively latitude dependent, changing from ∼400 km s −1 in the equatorial plane (slow speed region) to ∼ 800 km s −1 in the polar regions (high speed region), as observed by Ulysses [30]. This effect is mitigated during epochs of solar maximum, when the angular extension of the slow-speed region increases to higher latitudes. Beyond the TS, the SW slows down by a factor 1/S, where S = 2.5 is the shock compression ratio, as measured by the Voyager probes [31]. In this region, the wind is slowed down to subsonic speed. To incorporate such features in our model, we adopt the parametric expression given in [32]: where V 0 = 400 km s −1 , and L = 1.2 AU is the scale thickness of the TS. The top and bottom signs correspond to the northern (0 ≤ θ ≤ π/2) and southern hemisphere (π/2 ≤ θ ≤ π) of the heliosphere, respectively. The angle θ T determines the polar angle at which the SW speed changes from a slow to a fast region. It is defined as θ T = α+δα, where α is the tilt angle of the HCS and δα = 10 • is the width of the transition. With this approach, the angular extension θ T of the SW profile changes in time and it is linked to the level of solar activity, using the angle α as proxy. The expression is valid for r r , i.e., away from the Sun. Beyond the TS, the real SW speed is expected to decrease as r −2 , so that ∇ · V sw = 0 and CR particles do not experience adiabatic cooling. The radial and latitudinal SW profile is shown in Fig. 2 for two values of α corresponding to solar minimum (α ∼ = 10 • ) and solar maximum (α ∼ = 60 • ) conditions.

B. The particle transport
The Parker's equation for the particle transport contains all physical processes experienced by a given species of CR particles traveling in the interplanetary space. In Eq.(1), the drift-diffusion tensor can be written as: in a reference system with the third coordinate along the average magnetic field. The symbol K denotes the diffusion coefficient along the field direction, while K θ⊥ and K r⊥ the diffusion coefficients along the perpendicular and radial direction, respectively. K A expresses the value of the antisymmetric part of the diffusion tensor, where its explicit form results from the effects on the motion of CR particles due to drift. V sw is the SW speed and V D is the guiding center speed for a pitch angle-averaged nearly isotropic distribution function. The equation can be then re-written as: The motion of the CR particles in the HMF is usually decomposed in a regular gradient-curvature and HCS drift motion on the background average HMF and a diffusion due to the random motion on the small-scale fluctuations of the turbulent HMF. All these effects are included in the diffusion tensor K of Eq.(6), which can be decomposed in a symmetric part that describes the diffusion and an antisymmetric one that describes the drifts, i.e., Particle moving in a magnetic turbulence are pitch-angle scattered by the random HMF irregularities. This process is captured by the symmetric part of the diffusion tensor K S , which is diagonal if the z-coordinate is aligned with the background HMF. Three diffusion coefficients are therefore needed, namely, parallel diffusion K , transverse radial, K ⊥r , and transverse polar diffusion coefficient K ⊥θ . The coefficients can also be expressed in terms of mean free path λ along the background HMF, e.g., K = βcλ /3 (with β = v/c). The determination of the diffusion coefficients is a key ingredient to study the propagation of charged particles in turbulent magnetic fields like the HMF and is the subject of many theoretical and computational studies. The Quasi Linear Theory (QLT) has been successful at describing parallel diffusion, especially in its time-dependent and non-linear extensions [33]. Regarding perpendicular diffusion, the QLT provides upper limits within the field line random walk description [33,34], while the best approaches follow the nonlinear guiding center theory [35][36][37].
From a microscopic point of view, CR diffusion is linked to the resonant scattering of particles with rigidity R with the HMF irregularities around the wave number k res ∼ 2π/r L , where r L = R/B. The essential dependence of λ on the HMF power spectrum can be expressed as λ ∼ r 2 L B 2 /w(k res ) ∼ R 2 /w(k res ), where B 2 is mean square value of the background field and w(k res ) is the power spectrum of the random fluctuations of the HMF around the resonant wave number. The power spectral density follows a power-law as w(k) ∼ k −ν , where the index ν depends on the type and on the spatial scales of the turbulence energy cascade [38,39]. Therefore, λ depends on the turbulence spectral index as λ ∼ R 2−ν In this work, for the rigidity and spatial dependence of the parallel diffusion coefficient, we adopt a double power-law rigidity dependence and an inverse proportionality with the local HMF magnitude, following Ref. [32]: In this expression, K 0 is a constant of the order of 10 23 cm 2 s −1 , R 0 = 1 GV to set the rigidity units, B the HMF magnitude and B 0 the field value at Earth and written in a way such that the units are in K 0 . Here a and b are power indices that determine the slope of the rigidity dependence, respectively, below and above a rigidity R k , whereas h determines the smoothness of the transition. The perpendicular diffusion in the radial direction is calculated as K ⊥r = ξ ⊥r ×K || , while the polar perpendicular diffusion was parameterized as is a function that enhances K ⊥θ by a factor d near the poles, defined as [32]: Here The enhancement in the latitude direction of K ⊥θ , together with the anisotropy between the perpendicular diffusion coefficients and HMF modification at the polar regions, is needed to account for the very small latitudinal dependence of the CR intensity, as it was observed in the Ulysses data [30,40]. The adoption of constant ξ ⊥factors implies that K ⊥ and K follow the same rigidity dependence, which may be a simplification in the high-R domain [36,41]. Nonetheless, QLT-based simulations agree for nearly rigidity-independent ξ, with the typical value of 0.02-0.04 [34,43]. In this work, the parameters ξ ⊥r and ξ ⊥θ are fixed to the value 0.02. We now turn on drift effects, that account for the charge-sign and polarity dependence of CR transport in the HMF [27,44]. The regular motion of CRs on the large-scale HMF is given by the pitch-angle averaged guiding center drift speed V D . It can be related to the antisymmetric part of the diffusion tensor [45]: where the antisymmetric part of the tensor has the form: Here ijk is the Levi-Civita symbol, u(θ) is a function that describes the transition between the region influenced by the HCS and the regions outside of it and ζ(R) is a function of rigidity that suppresses drifts at low rigidity. To determine the value of K A , we note that the small value of the ratio K ⊥ /K suggests that CR particles move over many gyro-orbits in a mean free path, therefore the drift motion is weakly affected by scattering. In the weak scattering approximation, one has: where Q is the CR particle charge and K 0 A is a normalization factor ≤ 1. Drift motion is relevant close the HCS, where CRs cross many times regions of opposite HMF polarity. A 2D description of HCS drift is given in Burger & Hattingh [45]. In this approach, the drift velocity is given by: where the two vectors are defined as follows: The G-term in Eq.(14) describes the gradient-curvature drifts, the H-term describes the particle motion across the region affected by the HCS, e θ is the unit vector along the polar direction, and u(θ) is given by: with H the Heaviside step function, and The angle 2r L /r depends on the maximum distance that a particle can be away from the HCS while drifting. Finally, the function u(θ) is such that u(π/2) = 0, u(c h ) = 0.5 and ∂u(π/2)/∂θ = 1. CR drift coefficients are expected to be reduced in presence of turbulence as results theoretically and from numerical test-particle simulations [46,47]. In this work, we use a simple approach to incorporate drift reduction. Following Ref. [47], we adopt a reduction factor of the type: where the reduction occurs at rigidity below the cutoff value R A = λ ⊥ δB T , which depends on the perpendicular diffusion length and total variance of the HMF. The reduction is effective at R R A , when ζ ≈ (R/R A ) 2 1, while in the high-R limit one has ζ ≈ 1. The cut-off value R A depends on the HMF turbulence through λ ⊥ and δB T . With typical values of λ ⊥ ≈ 1.5 × 10 −3 AU and δB T ≈ 3.5 nT for the considered epochs, one can estimate R A ≈ 0.3 0.6 GV. In this work we have fixed it at 0.5 GV, corresponding to a proton kinetic energy of 125 MeV. The normalization K 0 A factor is fixed to 1, so that the whole drift reduction is regulated by ζ.
The most relevant feature of magnetic drift is that its direction depends on the sign of the charge,q = Q/|Q|, and on the HMF polarity A, via the productqA, so that particles with oppositeqA will drift in opposite directions and will follow different trajectories in the heliosphere. This characteristic is expected to give observable charge-sign dependence in the CR modulation. Finally, in a reference frame with the z coordinate along the average magnetic field, the diffusion tensor is given by Eq. (6). The effective diffusion tensor in heliocentric polar coordinates is obtained by a coordinate transformation in the modified Parker's field. In our 2D approach, the relevant components are K rr = K cos 2 ψ + K ⊥r sin 2 ψ, K θθ = K ⊥θ and K θr = K A sin ψ = −K rθ .

C. The proton LIS
To resolve the modulation equation for cosmic protons, their LIS must be specified as boundary condition. The determination of the CR proton LIS requires a dedicated modeling effort, starting from the distribution of Galactic CR sources and accounting for all the relevant physical processes that occur in the interstellar medium. In this work, we adopt an input LIS for CR protons that relies on a two-halo model of CR propagation in the Galaxy [48,49]. In this model, the injection of primary CRs in the ISM is described by rigidity-dependent source terms S ∝ (R/GV) −γ with γ = 2.28±0.12 for protons. The diffusive transport in the L-sized Galactic halo is described by an effective diffusion coefficient D = βD 0 (R/GV ) δ i/o with D 0 /L = 0.01±0.002 kpc/Myr [9,49]. The two spectral indices δ i/o describe two different diffusion regimes in the inner/outer halo, with δ i = 0.18±0.05 for |z| < ξ L (inner halo), and δ o = δ i + ∆ for |z| > ξ L (outer halo), with ∆ = 0.55±0.11. The z variable here is the vertical spatial coordinate. The half-thickness of the halo is L ∼ = 5 kpc and the near-disk region (inner halo) is set by ξ = 0.12±0.03. Finally, we considered the impact of diffusive reacceleration. Within the two-halo model, the interstellar Alfvénic speed is constrained from the data to lie between 0 and 6 km s −1 . Calculations of the proton LIS were constrained by various sets of measurements: low-energy proton data (at 140 -320 MeV) collected by Voyager-1 beyond the HP, high-energy proton measurements (E 60 GeV) made by AMS-02 in low Earth orbit, along with measurements of the B/C ratio from both experiments. The latter were essential to constrain the diffusion parameters of the LIS model [9]. Details on this model are provided elsewhere [49,50]. The resulting proton LIS is shown in Fig. 3 in comparison with the data from Voyager-1, along with PAMELA and AMS-02 measurements made in March 2009 and April 2014, respectively. The uncertainty band associated with the calculations is also shown in the figure. This model is in good agreement with other recently proposed LIS models [5,22,[52][53][54].

III. DATA ANALYSIS
In this section, we present the analysis method by which we extract knowledge and insights from the data using the mathematical framework described Sect. II. In practice, we defined a set of physics observables, to be computed as model predictions, and a set of model parameters to be determined by statistical inference.
A. The cosmic ray data The data used in this work consist in time-resolved and energy-resolved measurements of CR proton fluxes, in the kinetic energy range from ∼ 80 MeV to ∼ 60 GeV. Specifically, we use the 79 BR-averaged fluxes measured by the AMS-02 experiment in the International Space Station from May 2011 to May 2017 [13], and the 47+36 BRaveraged fluxes observed by the PAMELA instrument in the satellite Resurs-DK1 from June 2006 to January 2014 [15,16]. The data sample corresponds to a total of 10,101 data points collected over a time range of about 11 years, from the solar minimum from 2006 to 2009, the ascending phase to solar maximum, when the HMF polarity A reversed from A <0 to A >0, and the following descending phase until May 2017. These data have been retrieved by the ASI-SSDC Cosmic Ray Data Base [55].
The intensity of the CR proton fluxes in the energy range between 0.49 -0.62 GeV are shown in Fig. 4 as a function of time for both the PAMELA and AMS-02 data sets. From the figure, the complementarity of the two experiments is apparent. It can be seen that the highest intensity of the CR is reached during ∼ December 2009, i.e., under the solar minimum, while the lowest intensity occurs in ∼ February 2014, around solar maximum. The vertical dashed line of the figure shows the HMF reversal epoch T rev , along with the transition region shown as a shaded area where the HMF is disorganized and the polarity is not defined. The determination of T rev and the transition region are presented later on.

B. The parameters
The numerical model presented in Sect. II makes use of several physics input to be determined with the help of observations. Inputs include solar parameters, characterizing the conditions of the Sun or the interplanetary plasma, and transport parameters that describe the physical mechanisms of CR propagation through the plasma. Solar and transport parameters are inter-connected each other and they may show temporal variations related to the solar cycle. For instance, solar parameters such the magnetic field magnitude, its variance and its polarity are transported from the Sun into the outer heliosphere, therefore provoking time-dependence CR diffusion and drift.
We identified, in our model, a set of six time-dependent key parameters that are of relevance for the phenomenology of CR modulation. They are the tilt angle of the HCS α(t), the strength of the HMF at the Earth's location B 0 (t), the HMF polarity A(t), and the three diffusion parameters appearing in Eq. (8): the normalization factor of the parallel diffusion tensor, K 0 (t), and the two spectral indices of the rigidity-dependence of CR diffusion, a(t) and b(t), below and above the break R k , as seen in Eq. (8). Note that all key parameters are expressed as continuous functions of time t, but in practice, they have been determined for the epochs corresponding to the CR flux measurements.
The three solar parameters α, B 0 , A can be determined from solar observatories: data of HMF polarity and tilt are provided by the Wilcox Solar Observatory on 10-day or BR basis. Measurements of the HMF B 0 at 1 AU are done in-situ on daily basis, since 1997, by the Advanced Composition Explorer (ACE) on a Lissajous orbit around L1 [56]. It is important to notice that, in this study, our aim is to capture the effective status of the large-scale heliosphere sampled by CRs detected at a given epoch t, and this is connected to solar-activity parameters that are precedent to that epoch. In fact, several studies have reported a time lag of a few months between the solar activity and the varying CR fluxes [53,57], reflecting the fact that the perturbations induced by the Sun's magnetic activity take a finite amount of time to establish their effect in the heliosphere. To tackle this issue, for each epoch t associated to a given CR flux measurement, we perform a Backward Moving Average (BMA) for α and B 0 , and A, i.e., a time-average of these quantities calculated over a time window [t − τ, t]. The window extent τ is the time needed by the SW plasma to transport the magnetic perturbations from the Sun to the HP boundary, which ranges between ∼ 8 months (fast SW speed) and ∼ 16 months (slow SW speed). In the case of α, the window is large because the HCS is always mostly confined in the slow (equatorial) SW region. In the case of B 0 , the BMA has to be computed by an integration over the latitudinal profile of the SW speed at a given epoch. Our estimations are consistent with the lag reported in other studies [53,57] and supported by correlative analysis that we made a posteriori. Figure 5 shows the reference parameters B 0 , α calculated for for each reference epoch t corresponding to a BR-averaged CR measurement. A similar estimate is done for the polar magnetic field and for the resulting polarity A, in Fig. 5d. The latter can be regarded as a "smoothed" definition for the magnetic polarity A, otherwise dichotomous (A=±1). When the HMF is in a defined polarity state, one has A = ±1. During the HMF reversal transition epoch (shaded area in the figures), as the polarity is not well defined, the estimate of A takes a floating value between −1 and +1.
At this point, we also recall that several parameters entering the model that have been kept constant in the simulation, i.e., assumed to be known or time-independent.  [15,16] and AMS-02 (filled circles) [13,14]. The vertical dashed line shows the epoch of the HMF polarity inversion, along with the shaded area indicating the reversal epoch.
The HP and TS positions were fixed at r HP =122 AU and r TS =85 AU, deduced from the Voyager-1 observations. The data suggest that the TS may vary over the solar cycle of the order of a few AU, but its impact in the CR fluxes is not negligible [53].
The h parameter of Eq. (8), describing the smoothness of the transition between the two diffusion regimes below and above R k , was kept constant at h = 3. Within the precision of the data, the h parameter has no appreciable impact on the CR fluxes. Similarly, the rigidity break R k for K was kept fixed at the value 3 GV. This parameter represents the scale rigidity value where the CR Larmor radius matches the correlation length of the HMF power spectrum, which is at the GV scale. Regarding the value of R k , we found that time variations on this quantity do not give appreciable variations in the CR fluxes [see,e.g. 32]. The ξ ⊥i coefficients for the diffusion tensor, for which the values used here represent a widely used assumption [e.g., 40]. The polar enhancement factor of Eq.(9) is kept constant at d = 3 for ξ ⊥θ so that the condition K ⊥ /K 1 is still fulfilled at the polar regions. Regarding magnetic drift, the critical rigidity R A of Eq.(18) is kept constant at 0.5 GV following previous studies and independent observations on the CR latitudinal gradient [32,59]. This choice could be tested only with low-rigidity CR data (R R A ), as our results are insensitive to the exact value of R A . The normalization factor for drifts speeds K 0 A was chosen to be unity such to set "full drift" speeds in the propagation model for all the periods, and this the drift reduction is entirely given by Eq. (18). Reductions in the K 0 A -value may occur during periods of strong magnetic turbulence, e.g., during solar maximum [25,59].

C. The statistical inference
The parameter grid -The transport parameters K 0 (t), a(t) and b(t) have been determined from the AMS-02 and PAMELA data by means of a global fitting procedure. For this purpose a six dimensional discrete grid of the model parameters vector q = (α, B 0 , A, K 0 , a, b) was built, i.e., the model was run for every node of the grid such to produce a theoretical calculation for the CR proton flux. In the grid, the parameter α ranges from 5 • to 75 • with steps of 10 • , B 0 from 3 to 8 nT with steps of 1 nT, and the polarity A takes the two values A = +1 and A = −1. The parameter K 0 ranges from 0.16 to 1.5 × 10 23 cm 2 s −1 , with steps of 0.08 × 10 23 cm 2 s −1 , the indices a and b range from 0.45 to 1.65 with steps of 0.05. The total number of grid nodes amounts to 938,400. For each node of the parameter grid, a theoretical prediction for the modulated proton flux J m (E, q) was evaluated, as function of kinetic energy, over 120 energy bins ranging from 20 MeV to 200 GeV with log-uniform step. Using the SDE technique, 2 × 10 3 pseudo-particles were Monte Carlo generated and retro-propagated for each energy bin. This task required the simulation of about 14 billion trajectories of pseudo-protons, corresponding to several months of CPU time. Once the full grid was completed, the output flux was tabulated and properly interfaced with the data. For each data set J d (E, t), representing a set flux measurements as function of energy for a given epoch t, a χ 2 estimator was evaluated as: Similarly to J m , the χ 2 estimator is built such to be a continuous function of the parameters q, except for the variable A that is treated as discrete. From the χ 2 estimator, the transport parameters {K 0 , a, b} can be determined by minimization at any epoch, while the solar parameters {B 0 , α, A} can be considered as "fixed inputs", as they are determined by the epoch t using the BMA reconstruction presented above. For a given set of BMA inputs such as B 0 and α, the flux J m (E, q) can be ex- pressed as a continuous function of the parameters by means of a multilinear interpolation over the grid nodes.
In the α − B 0 plane, one has α j < α(t) < α j+1 and B 0k < B 0 (t) < B 0,k+1 , where α j and B 0,k are the closest values of the grid corresponding to their BMA averages. Regarding polarity A, both ±1 evaluations were done under the assumption that the polarity is known. The flux model dependence upon energy should also be handled. In Eq. (19), E i are the mean measured energies reported from the experiments (coming from binned histograms). In general, the E i array does not correspond to the energy grid of the model. The model evaluation of J m (E, q) at the energy E i was done by log-linear interpolation.
The uncertainties -The σ factors appearing in Eq. (19) represent the total uncertainties associated with the flux. They can be written as are the experimental errors associated to the flux measurement of the i-th energy bin around E i , while σ 2 m (E i , t) are the theoretical uncertainties of the flux calculations evaluated at the same value of energy. Uncertainties in experimental data are of the order of 10 % in the PAMELA data and ∼ 2% in the AMS-02 data, although they depend on kinetic energy. Theoretical uncertainties include statistical fluctuations of the finite SDE generation of pseudo-particle trajectories. Uncertainties are relevant at low energy where, due to the heavy adiabatic energy losses, the Monte Carlo sampling suffers from a smaller statistics. Thus, after repeating many times the simulation with the same modulation parameters, the modulated flux will fluctuate around an average value because of the random process of pseudoparticles propagation with the SDE approach. These fluctuations can be arbitrarily reduced with the increase of the pseudo-particle generation, but at the expense of a large CPU time. The evaluation of these uncertainties can be done as follows. Given N m as the number of pseudo-particles that reach the boundary with energy E, and N G as the number of pseudo-particles generated at the same energy, the ratio of the modulated flux to the LIS flux is J m /J LIS ≈ N m /N G . Since the propagation process is stochastic in nature, the relative error of the modulated flux scales as δJ m /J m = 1/ √ N m , where N m = N G (J m /J LIS ). We found that the generation of N ∼ = 2 × 10 3 pseudo-particles for each energy bin is sufficient for being not dominated by SDE-related uncertainties. The relative uncertainties as function of kinetic energy are shown in Fig. 6. The errors are about ∼ 10 − 20% at 20 MeV of energy and decrease with increasing energy. They become constant at ∼ 2% above few GeVs. A minor source of systematic error comes from the multilinear interpolation of the parameter and energy grid, i.e., from the method we used to evaluate the flux at any arbitrary set of parameters and energy. From dedicated runs, we have estimated that the uncertainty introduced by the interpolation, rather than the direct simulation with of J( q, E), is always of the order of 1 %. An important source of systematic error is the uncertainty coming from the input LIS of CR protons, see Sect. II C. The LIS uncertainties are highly energydependent. They are significant in the energy region of ∼ 1-10 GeV (up to 30 % and more), where direct interstellar data are not available but the modulation effect is still considerable. However, in this energy region, the Galactic transport parameters regulating the LIS intensity are in degeneracy with the free parameters of CR diffusion (Sect. III B) and in particular with K 0 [50]. Such a degeneracy translates into a correlation between the best-fit K 0 values and the LIS intensity at the GeV scale which, in turn, determines the absolute scale of the the modulated CR flux J 0 at the GeV scale. The K 0 − J 0 correlation is also discussed in Sect. IV A. To estimate the impact of the LIS uncertainty on the temporal dependence of the best-fit parameters of CR diffusion in heliosphere, we proceeded as in Ref. [50,51]. We performed dedicated runs of fitting procedure for a large number of randomly generated LIS functions where, for each input LIS, the time-series of the diffusion parameters were determined. In practice, the LIS functions were generated using the Monte Carlo framework in Ref. [49], i.e., according to the probability density function of the Galactic CR transport parameters. With this procedure, the systematic uncertainties associated with the LIS modeling are included in the final errors with a proper account for their correlations.
The reversal phase -The parameter T rev marks the epoch of the 2013 magnetic reversal, where the HMF flipped from negative to positive polarity states The polarity of the HMF, however, is well defined only for t T rev and t T rev , where the large-scale HMF structure follows a dipole-like Parker's field to a good approximation. During reversal, the polarity of the field is less sharply defined and the HMF field follows a more complex dynamic [e.g., 60]. A way to account for this situation is to use a generalized definition of polarity, such as the BMA reconstruction A of Fig. 5 which ranges from -1 to +1. For any given parameter configuration q, the flux model J m (E, q) can be built as a linear combination of fluxes with defined polarities, weighted by a transition function P≡(1 − A)/2: where q (±) = {α, B 0 , A ± , K 0 , a, b} is a vector of parameters with fixed polarity A = ±1, and J (±) m are the corresponding modulated fluxes. The weight P ranges from 1 to 0, for floating polarity A ranging from -1 to 1. The time-dependence of the P(t)-function associated to the polarity A(t) of Fig. 5 can be expressed as follows: where δT ∼ = 3 months. The transition function P(t) is such that P ∼ = 0 (P ∼ = 1) for t 3 T rev (t 3 T rev ) within 1 % level of precision, i.e., when t = T rev ± 3δT , the flux is 99% made of a fixed polarity, while the maximum mixing is for t = T rev when P(t) = 1/2. It is worth noticing that Eq.(20) relies on the implicit assumption that, during HMF reversal, the modulated flux of CRs can be regarded as a superposition of fluxes with positive and negative polarity states. We also note that this  [60,61]. The parameter extraction -Our determination of the diffusion parameters K 0 (t), a(t) and b(t) is based on the least squares method. In practice, we proceeded as follows. Given a set of CR proton flux measurements J d (E, t), for each parameter x= K 0 (t), a(t), and b(t), the corresponding χ 2 (x) distribution, defined as in Eq. (19), is evaluated. The evaluation is done for all values of the other parameters y = x, marginalized over the hidden dimensions. This returns a curve χ 2 min (x) as function of the parameter x and minimized over all hidden dimensions. From the minimization of χ 2 min (x), the best-fit parameterx and its corresponding uncertainty are estimated. For the minimization, we tested two approaches. One method consisted in the interpolation with a cubic spline of the whole χ 2 min (x) curve. A second method, similar to Corti et al. [5], consisted in the determination of the minimum x i,min point from a parameter scan over the grid, and then by making a parabolic re-fitting of the χ 2 min (x) curve around the x i,min and its adjacent points. The position of the minimum and its uncertainty can be calculated as estimation of x best . The errors on the parameters are estimated as , where x ± is the parameter value such that χ 2 min (x ± ) = χ 2 min (x best ) + 1 above and below x best , which is the standard error estimation of the least squares method. The little discrepancy of the two methods was used as a systematic errors which, however, turned out to be negligible in comparison with the standard errors of the fit. The shapes of the χ 2 min projections as function of the diffusion parameters is illustrated in Fig. 7 for two distinct epoch March 2009 (BR 2379 during solar minimum) and April 2014 (BR 2466, during solar maximum). For each curve, the best-fit parameterx is shown (vertical line) along with its associated uncertainty σ x (shaded band). In the two considered epochs, the data come from PAMELA and AMS-02 experiment, respectively. As seen from the figure, AMS-02 gives in general large χ 2 -values in comparison with PAMELA. In both time series the convergence of the fit is good and the parameters are well constrained. It can be seen that the AMS-02 data provide tight constraints on the K 0 and b parameters, while the parameter a is more sensitive to low-rigidity data and thus it is better constrained by PAMELA. After the best-fit pa- rameters have been determined for a give set of data, the best model flux J best (E) is recalculated using a multilinear interpolation over the 5-dimensional grid such that x i ≤ x best < x i+1 where x = α, B 0 , K 0 , a, and b. In this procedure the polarity A is not involved, because it is regarded as fixed parameter. The flux determination done under both A + /A − hypotheses gives the two J ± fluxes of Eq. (20). The best model is shown in Fig. 8

IV. RESULTS AND DISCUSSION
Here we present the results of the fitting procedure described in Sect. III C and implemented using the considered data set on CR protons of Sect. III A. We found that the agreement between best-fit model and the measurements on the fluxes of CR protons was in general very good for all the data sets and over the whole rigidity range. In Fig. 9 the best-fit models for the proton fluxes are shown as colored lines for some selected epochs, along with the CR proton LIS. The calculations are compared with the data from experiments PAMELA and AMS-02 at the corresponding epochs. The long-dashed line represents the proton LIS model used in this work and presented in Sect. II C.

A. Temporal dependencies
The main results on the parameter determination procedure are illustrated in Fig. 10. The figure shows the best-fit model parameters K 0 , a, and b as function of the epoch corresponding to the measurements of AMS-02 (filled circles) and PAMELA (open squares). The vertical dashed line and the shaded area around it represent the reversal phase, as in the previous figures. As a proxy for solar activity, Fig. 10d shows the monthly SSN data. The solid line shows the smoothed SSN values, obtained with a moving average within a time window of 13 months, along with its uncertainty band. It can be seen that the diffusion parameters show a remarkable temporal dependence, and such a dependence is well correlated with solar activity. From the figure, it can be seen that the normalization of the parallel diffusion coefficient K 0 shows a clear temporal dependence. The diffusion normalization appears to be maximum in the A < 0 epoch before reversal (t T rev ), and in particular during the unusually long solar minimum of 2009-2010. The minimum of K 0 is reached during solar maximum in 2014, about one year after polarity reversal. From the comparison between panel (a) and panel (d), the K 0 parameter appears anti-correlated with the monthly SSN. Physically, larger values of K 0 imply faster CR diffusion inside the heliosphere, thereby causing a milder attenuation of the LIS, i.e., giving a higher flux of cosmic protons in the GeV energy region. In contrast, lower K 0 values imply slower CR diffusion which is typical in epochs of high solar activity where the modulation effect is significant. Qualitatively, this behavior can be interpreted within the Force-Field approximation where, in fact, positive correlation is expected between SSN and the modulation potential φ ∝ 1/K 0 [9]. Within the framework of the Force-Field model, the parameter φ is interpreted as the average kinetic energy loss of CR protons inside the heliosphere. For similar reasons, a positive correlation between the best-fit K 0 -value and the CR flux intensity J 0 at a given energy as can be noticed, in particular, from the comparison of Fig. 10a with Fig. 4. Our finding are in agreement with earlier works [5,53,62]. During the reversal phase, the temporal evolution of the model parameters in Fig. 10 is obtained using the weighted linear combination of model fluxes with opposite polarities given by Eq. (20). During this epoch, the diffusion of CRs is slow and the tilt angle α reaches large values, typically higher than 65 • .
The inferred K 0 -values and their temporal evolution are related to the level of magnetic turbulence in the heliospheric plasma. As clear from the figure, the diffusion is faster when the Sun is quiet with low turbulence levels and vice-versa. From Eq.(8), the CR diffusion coefficients are linked to the HMF intensity and its temporal evolution which, however, from Fig. 5, appears to be quite shallow in the epoch considered. As recently suggested in Ref. [63], the relation between the diffusion coefficient and the magnitude of the local HMF can be described by a power-law, but the two quantities obey to different relationships for ascending and descending phases of the Solar Cycle. Physical explanation for these behaviors may involve temporal variations in the spectrum of heliospheric turbulence during the solar cycle [64,65], that we discuss in the following. Investigations on the correlations between solar and diffusion parameters are made in Sect. IV C.

B. The evolving turbulence
The a and b parameters shown in Fig. 10 describe the rigidity dependence of CR diffusion tensor K below and above the break value R k . These parameter can test how the Sun variability affects the spectrum of magnetic irregularities of the heliospheric plasma, that is, its turbulence spectrum. From figure, it can be noted that both parameters show a characteristic temporal dependence in the epoch considered. In the negative polarity epoch of t T rev , and in particular during solar activity minimum, the spectral indices of CR diffusion are seen to vary smoothly and slowly with time.
The two spectral indices show a different temporal dependence. The index a is found to be essentially time independent, with an average value of a = 1.21±0.06, while the index b shows a distinct long-term evolution in the considered period. During the long unusual minimum from 2006 to 2009, b remains constant at a value of b = 0.74±0.03, as long as the solar activity is quiet and the corresponding number of monthly sunspots is below ∼ 50. Subsequently, in ∼ 2010-2011, when the ascending phase of the solar cycle sets in, b starts to increase steadily. During this period, the CR flux decreases steadily as well. The increase keeps going during the whole reversal phase, i.e., at full maximum solar activity. Here the b parameter reaches an average maximum value of 1.3 ± 0.07. After this phase and during the flux recovery phase in the positive polarity epoch, the index b decreases steadily during the descending phase of the solar cycle, until it recovers the values of the previous solar minimum. Instead, the index a shows no prominent features over the whole descending phase.
It should be noted, however, that the a parameter is poorly constrained in the A > 0 phase, because the AMS-02 data are available only above 1 GV of rigidity, and thus they are not highly sensitive to this parameter. From the figure, it can be seen that the index b is negatively correlated with the diffusion normalization parameter K 0 : during minimum, where K 0 is large and the CR diffusion is therefore fast, its rigidity dependence is shallow (b ≈ 0.8) in comparison to solar maximum, where diffusion is slow and its rigidity dependence is more pronounced (b ≈ 1.3). Since the two indices are related to the power spectrum of the heliospheric turbulence, they could be used to infer the spectral index ν of the power spectrum density of HMF irregularities (see Sect. II B). Keeping in mind that λ ∝ R 2−ν , the index a is related to the power spectrum density in the energy-containing range, while the index b is related to the power spectrum in the inertial range of the turbulent energy cascade of HMF. The results indicate that the diffusion spectrum in the energy-containing regime does not depend on the solar activity, while, in the inertial range, the spectrum appears to evolve as a function of the solar activity, with a clear delayed peak at the solar maximum. The spectral index of the turbulence in the energy-containing range is ν ec = 0.79±0.13 over all the period examined in this work, while in the inertial range the spectral index evolves from ν in = 0.74±0.08 at solar minimum to ≈1.3±0.15 during the solar maximum.
The temporal and rigidity dependence of the CR mean free path λ (t, R) can be determined from Eq.(8) using our best-fit parameters. At the R≈1 GV rigidity scale, our λ is found to range between 0.05 AU and 0.3 AU, depending on solar activity. This result is in excellent agreement with the large collection made in Ref. [66] of observational measurements on the scattering mean free path [42]. In addition, our result show that the CR variability involves the rigidity dependence of the diffusion tensor, in particular via the spectral indices a = a(t) and b = b(t). An important implication of this finding is that the parallel diffusion coefficient cannot be write as a product K (t, R) = f (t)×g(R), where a universal rigidity dependence g(R) is modulated in amplitude by means of a factorized function f (t) [50,62]. Mathematically, this makes the K (t, R) function of Eq.(8) a non separable function of rigidity and time variables. Physically, it indicates that the HMF turbulence spectrum varies significantly over the solar cycle, depending on the cycle phase. In particular, the power spectrum is observed to be steeper around solar maximum and flatter during solar minimum, with a quasi-periodical pattern. The temporal variability of HMF turbulence is also studied from the analysis of neutron monitor data [64]. These findings suggest that during epochs of quiet activity, kinetic self-organized turbulence dominates the CR spectrum, such as, e.g., a Kolmogorov-type cascade, while random processes and transient events in the heliosphere play a key role during high-activity epochs of the solar cycle. The use of wider sets of data may allow to provide better clarification on such a behavior.

C. Cross-correlations
We now inspect the running cross-correlation between solar and transport parameters. Figure 11 displays the scatter diagrams of the best-fit diffusion parameters against the BMA reconstruction of the local HMF value, B 0 (left column) and the HCS tilt angle α (right column). In panel (a), the diffusion normalization parameter K 0 is shown. The different markers are used to indicate the reconstructions obtained during epochs of positive (blue circles) and negative polarity (pink squares), as well as during reversal phase (green triangles). This behavior can be compared with the one found by Wang et al. [63] where, from an analysis of the ascending and descending phases of the solar cycle (both during negative polarity) two distinct power-law relations were observed between diffusion coefficient and local HMF magnitude. Our results confirm the relationship between K 0 and B 0 becomes complex when the examination is done over a large fraction of the solar cycle that include polarity changes.
In particular, two distinct relationships can be observed for A < 0 and A > 0 polarity conditions. Regarding the correlation between the spectral index parameters a and b with the HMF magnitude B 0 , smoother relationships were found. The index a is nearly constant with time, while the index b increases slowly during solar maximum, i.e., during the reversal phase. Both parameters are seen to depend only weakly on the polarity phase, and no particular cross-correlation is observed between two spectral indices. The scatter plot of K 0 versus tilt angle is also shown, in Fig. 12 where, again, the different style of the markers refer to the different phases of solar activity. The dependence is similar to that observed with the HMF intensity, showing a pronounced negative correlation and a characteristic modulation loop. The correlation between the flux intensity J 0 and the diffusion normalization K 0 is shown in Fig. 13. In this figure, the flux intensity J 0 is extracted from the data at the reference kinetic energy E 0 = 0.49 − 0.62 GeV, as in Fig. 4, while K 0 is the best-fit value at the corresponding epoch. From the figure, the CR flux intensity appears in general well correlated to the normalization factor of the diffusion coefficient, which appears to be the driving parameter of the modulation model. It can also be seen that relationship between J 0 and K 0 is remarkably linear during epochs of well-defined polarity. We describe it with the following empirical relation: By making separate fits for the two polarity epochs, we obtained η + = (2212 ± 250) × 10 −23 for A > 0, and η − = (1929 ± 260) × 10 −23 cm −4 GeV −1 sr −1 for A < 0. The best-fit offset are J + off = −46 ± 21 for positive polarity, and J − off = −286 ± 68 m −2 s −1 GeV −1 sr −1 for negative polarity. The two fits are shown in Fig. 13 as dashed line. It is interesting to note that, within the fitting errors, the two slopes η + and η − turned out to be consistent each other, i.e., the slope of J 0 (K 0 ) is polarity and charge-sign independent. Polarity-effect results into different offsets J ± of f for the two phases. This result may help to quantify the effects of drift motions to the CR modulation. The diffusion coefficient appears to be independent upon thê qA sign product, as indicated by the consistency between η + and η − values from the fit. For a given K 0 value, the resulting difference in the fluxes is only due to the opposite directions of the net drift and convective flux for epochs of opposite polarities. The quantity ∆J ≡ J + off − J − off can be used as a measurement of the net effect of drift on the total CR flux, for a given level of CR diffusion.
We also note that in the figure, the fit results obtained under periods of undefined polarity (green triangles) connect smoothly the two regimes. In this epoch the role of drift is not well understood, but the flux J 0 remains correlated with K 0 . To close the loop, it may take an entire cycle of magnetic polarity.

D. Lags and loops
From Fig. 10, it can be noticed that a time shift of a few month is present between the smoothed SSN (the S(t) function) and the best-fit modulation parameters K 0 (t), a(t) and b(t). For instance, the highest CR flux intensity was reached around October 2009, with J max = 2289±220 m −2 s −1 GeV −1 sr −1 , i.e., about eight months after the SSN minimum of February 2009. Similarly, the minimum flux intensity was observed around February 2014, J min = 498± 23 m −2 s −1 GeV −1 sr −1 , while solar maximum occurred in April 2013. To estimate the average time lag between K 0 (t) and the smoothed SSN S(t), we compare the correlation between K 0 (t) and S(t − ∆T lag ). The best-value for the lag ∆T lag can be obtained by a scan of ∆T lag , in order to determine the Pearson linear correlation coefficient ρ as function of ∆T lag . The ∆T lag parameter which maximizes ρ is then taken as best estimate of the average time lag between the SSN and CR modulation parameters. For the analyzed period, we obtain ∆T lag = 11.4 ± 1.4 months. Thus, on average, the modulation of CRs observed at the epoch t is related to manifestations of solar activity at the epoch t − ∆T lag . The correlation between diffusion parameters and smoothed SSN is shown in Fig. 14, where the model parameters at the epoch t are shown as a function of the SSN at the same epoch (left column) and at the epoch t−∆T lag (right column). In general, when the time lag is not taken into account, the diffusion normalization K 0 (t) appears as a multivalued function of SSN, showing a characteristic hysteresis structure over the different phases of the solar cycle. When the lag is taken into account, the curve of K 0 vs SSN shrinks, approaching a single-valued function. This would allow, in principle, to forecast the modulation parameters at the epoch t from observations of SSN made in advance by ∆T lag ). However, the a and b parameters versus the delayed SSN do not show clear one-to-one relationships, which suggests that the use of a single lag value may be a too simplistic approach. The calculated lag depends weakly on the BMA averages used to define the heliosphere status. On the other hand, the BMA procedure of Sect. III B is well motivated by the observation of such a lag. In this respect, an estimate of the uncertainty on ∆T lag can be done by varying the time window T BMA used to get the average conditions (B 0 and α) of the heliosphere. Our estimation of ∆T lag is fairly consistent with other recent works [3,53,67]. Nonetheless, there are some discrepancies with the reported values if one account for even/odd cycle dependence of the lag. Our estimation of the time lag lies in solar cycle 24, but it appears longer than that reported in previous even-numbered solar cycles, though it is comparable to the lag observed in odd-numbered solar cycles [68][69][70]. In this respect, as well as in other characteristics, cycle 24 is unusual when compared to previous even cycles. Other differences may be related to the rigidity of CR particles, as past studies are based on neutron monitors rates. The global dependence of the time lag upon the solar cycle and on the rigidity of the CR particles will be addressed in a forthcoming paper.

V. CONCLUSIONS AND DISCUSSION
Thanks to the recent availability of time-resolved data from space, the study of CRs in the heliosphere has become an active topic of investigation. In particular, the recent data released by AMS-02 and PAMELA on the monthly evolution of proton and helium permits new investigation of the solar modulation phenomenon over a large fraction the of solar cycle. These data have triggered new efforts at establishing advanced models of CR propagation in heliosphere [22,71,[73][74][75]. In particular, many recent studies were focused on specific aspects of the CR modulation such as, e.g., the particle dependence of CR diffusion [5,50], the relationship between modulation and solar activity proxies [63,72], the derivation of improved LIS evaluation [20,75], or the extraction of CR modulation parameters using statistical inference [5], which is also the main goal of the present paper. More specifically, in this paper, we have investigated the propagation of Galactic CRs in the heliosphere using a numerical model based on stochastic simulations and calibrated by means of a large set of experimental data. The data consist of time-series of CR proton fluxes reported by AMS-02 and PAMELA experiments in low Earth orbit. The measurements are made on 27-day basis, corresponding to a solar rotational period, and cover a time range of 11 years, corresponding to a solar cycle period. The sample include epochs of very different solar conditions such as solar minimum, solar maximum, ascending and descending phases, as well as positive and negative HMF polarity states. The time range and resolution of these data is therefore optimal for the study long-term modulation of Galactic CRs, and in particular, for investigating influence of solar variability in the diffusive propagation of CRs in the heliospheric turbulence.
In our calculations we have used, as time-dependent physical inputs, BMA values of the tilt angles α of the HCS, the local HMF strength at 1 AU B 0 , and the magnetic polarity A. These quantities constitute a very good proxies for solar activity. In this analysis, we have been focused on the parameters describing the temporal and rigidity dependence of CR diffusion. We have determined the time-series of the diffusion normalization, K 0 , and that of the spectral indices a and b that control the dependence of CR diffusion upon rigidity.
In practice, to perform a statistical inference using the data, and to account for the evolving conditions of the heliospheric plasma, we have built a large array of differential energy fluxes J(E), evaluated at Earth's location, corresponding to 938,400 parameter configurations. To sample such a 6-dimensional parameter space, we have simulated about 14 billions trajectories of cosmic protons in the interplanetary space. Each simulated particles was backwardly propagated from Earth's vicinity to the heliospheric boundaries. The array of models generated in this work can be used to estimate the modulation parameters of CR protons at any epoch and for any set of experimental data, ranging from 20 MeV to hundreds GeV of kinetic energy. We also note that in our model, the time dependence of the problem is treated by providing a time series of steady-state solutions for J p associated with a time series of input parameters k 0 , which is a simplification. Such an approach stands as long as the timescales between CR transport in the heliosphere does not exceed the analyzed changes in solar activity. To extend the analysis to smaller time-scale (e.g., daily) or to lower energies (e.g., MeV-scale), a time-dependent solution of the Parker's equation should be considered. Nonetheless, we also stress that the time-series of best-fit parameters derived in this work should be regarded as effective values, averaged over the CR propagation histories, not necessarily representing the instantaneous conditions of the heliospheric plasma.
Our approach is also simplified in several aspects, for example regarding the rigidity and spatial dependence of the diffusion tensor, or its perpendicular components. Nonetheless, in comparison to our earlier works, we have introduced several new recipes that capture most of the relevant features of CR propagation in the heliosphere. The agreement of our calculations with the CR flux data is very satisfactory. As we have shown, using CR proton data, it is possible to determine the detailed evolution of the rigidity dependence of the diffusion coefficient with the solar activity, ad thus, the physical nature of the turbulence embedded in the frozen-in HMF carried out by the SW. Our findings indicate that solar variability has an important effect on the turbulence spectrum of HMF irregularities, and an imprint of this mechanism can be observed in the rigidity dependence of the diffusion tensor. In particular, we have reported a remarkable long-term dependence for the two spectral indices a and b. These results show that the turbulence regime evolves with time, following the solar cycle, and thus the temporal and rigidity dependencies of CR diffusion coefficients cannot be described by a separable function of the type K (t, R)≡K 0 (t)×f (R). In this respect, we remark that the time-rigidity separability for CR diffusion is assumed by several models of solar modulation, although such an assumption is not supported by theoretical considerations [19,62,63]. Moreover, the study of the correlation between solar and diffusion parameters reveals charge-sign dependent features in the CR modulation effect, such as different patterns for the different phase of the HFM polarity cycle.
We remark that solar cycle 24 has been unusual when compared to the previous cycles, therefore also the CR modulation conditions were unusual. The solar minimum between cycles 23 and 24 was quite longer and deeper than expected [32,68]. while the maximum of cycle 24 was the smallest recorded in a century of standardized SSN observations, and with a double-peak structure [76]. In our analyzed data sample, the correlation between CR flux modulation and solar activity as measured by the SSN is apparent. The CR proton intensity modulation, in anti-phase with solar activity, in the considered period shows an average time lag of about 11 months. A next phase of this work is to study the dependence of the lag on solar activity parameters (such as SW speed or HMF polarity) and CR transport properties (such as diffusion or drift coefficients), in order to understand the dynamics of the physical mechanisms behind the solar modulation phenomenon. Further steps also include the implementation of a better description of the HMF, of the diffusion tensor and the drift reduction factor during solar maximum. In particular, we assumed "full drift" at any phase of the cycle, including the HMF reversal epoch where the modulated flux of CRs was modeled as superposition of fluxes with positive and negative polarity states. While our approach provided a good description of the flux evolution in the reversal region, one may ar-gue that large-scale drift may be suppressed during solar maximum due to the more chaotic structure of the HMF. This idea can in principle be tested using the data. In particular, the availability of time-dependent measurements on CR antiprotons will be precious to study the modulation effect across solar maximum. Data of the temporal dependence of CR antiprotons are still lacking, but the AMS-02 experiment has the capability of making such a measurement.