New solar modulation modeling of the galactic proton ﬂux measured by the AMS02 and PAMELA experiments

. A thorough understanding of solar e ﬀ ects on the galactic cosmic rays is relevant both to infer the local interstellar spectrum characteristics and to investigate the dynamics of charged particles in the heliosphere. We present a newly developed numerical modulation model to study the transport of galactic protons in the heliosphere. The model was applied to the 27-day averaged galactic proton ﬂux recently released by the PAMELA and AMS02 experiments, covering an extended time period from mid-2006 to mid-2017.


Introduction
Nearly all direct observations of Cosmic Rays (CRs) are made inside the heliosphere. As a result, the CR spectrum observed at Earth can differ significantly from the Local Interstellar Spectrum (LIS) due to many physical processes that modify the spectrum. This is particularly true for CRs with kinetic energies below ∼20 GeV, which can be efficiently diffused, curved and de-accelerated adiabatically during their propagation from heliosphere to the Earth. The changes in the Galactic CR (GCR) intensities and energy spectra are related to the solar activity, and are referred to as CR solar modulation. A thorough understanding of solar effects on the GCR propagation is therefore relevant both to infer the LIS and to investigate the dynamics of charged particles in the heliosphere. It is also important for astronauts and the electronic components radiation hazard during long-duration missions.
The new precise data from AMS02 [1] and PAMELA [2] experiments offer a unique possibility to study the solar modulation over a long period of time.

The proton model
Propagation in the heliosphere is described by the Parker's equation [3]; Here t is the time, p is the CR particle momentum, V is the solar wind speed, V D is the averaged guiding center speed for a pitch angle-averaged distribution function f and K is the diffusion tensor. In a e-mail: emanuele.fiandrini@pg.infn.it this work, a 2D numerical integration with ∂/∂t = 0 is made using a Stochastic Differential Equation approach with a backward time integration, based on [4]. The heliosphere is a dynamic void in the local interstellar medium where Sun activity dominates. In this work the heliosphere has a spherical symmetry with the boundary located at r = 122 AU, where local interstellar spectrum is given as boundary condition, with the termination shock at r = 85 AU and the Earth position at r = 1 AU in the equatorial plane. The Heliospheric Magnetic Field (HMF) is wounded up in a modified Parker's spiral [5,6], as where B 0 is the HMF value at Earth position, r is the Sun radius and δ(θ) = 8.7 × 10 −5 / sin(θ) if θ < 1.7 o or θ > 178 o and = 0 otherwise. A major corotating structure of importance to CR modulation is the wavy Heliospheric Current Sheet (HCS), which divides the solar magnetic field into hemispheres of opposite polarity. The HCS is parametrized by using its tilt angle α and is well correlated to solar activity. The waviness of the HCS is one of the best proxy for solar activity relevant to CR modulation and is widely used in numerical approaches. The tilt angle has been observed at the Wilcox solar observatory since 1976. Solar wind velocity is radially directed from close to the Sun, up to the inner heliosheath. During periods of minimum solar activity the flow is latitude dependent, changing from an average of 450 kms −1 in the equatorial plane to 800 kms −1 in the polar regions. This significant effect disappears with increasing solar activity, when the average speed is ∼400 kms −1 at all the latitudes. The latitudinal dependence of the solar wind speed is parametrized as in [7] and was made dependent on the solar activity by increasing the polar angle at which the solar wind speed changes from slow to fast as a function of the tilt angle.
Diffusion is described by a diagonal diffusion tensor with three diffusion coefficients, namely parallel diffusion K || , the transverse radial, K ⊥r , and transverse polar diffusion coefficient, K ⊥θ . The parallel diffusion can be computed using quasi linear theory [8]. The diffusion coefficient parallel to the background field, K || , was parametrized as in [10] with β = v/c, the ratio of the particle's speed to the speed of light. K 0 is a constant in units of 10 23 cm 2 s −1 , with R 0 = 1 GV, B the HMF magnitude and B 0 the local field value, such that the units are in K 0 . Here a and b are a power indexes that determine the slope of the rigidity dependence, respectively, below and above a rigidity with the value R k , whereas h = 3.0 determines the smoothness of the transition. The value R k determines the rigidity where the break in the power law occurs. The parameters K 0 , a, b depend on the solar activity, therefore vary in time and are free parameters that can be used to fit measured CR flux at a given epoch. Perpendicular diffusion in the radial direction is parametrized as K ⊥r = 0.02 × K || , which a widely used assumption [9], while the polar perpendicular diffusion was parametrized as K ⊥θ = 0.02 × g(θ) × K || , where g(θ) is a function that enhances the diffusion by a factor d = 3 near the poles, as defined by Eqn. (8) in [10]. The motion in the large-scale average HMF is given by pitch-angle averaged velocity V D . In this 2D modeling, the HCS is described as a region with an extension that depends on the rigidity of CR particles and on the value of the tilt angle. In this approach, the drift velocity is given by The first term describes the gradient and curvature drifts, while the second one describes the motion in the region affected by the HCS. f (θ) is a function that describes the transition between the regions inside and outside the HCS, given in [6], e θ is the unit vector along the polar direction and K A is the anti-symmetric part of the diffusion tensor [11], given by Here Q is the CR charge, R Ak is the break in K A and K A0 =1. The local interstellar spectrum (LIS) of galactic CR is an input for any modulation models. Fluxes are assumed isotropically distributed at heliosphere boundary, in a steady-state configuration. In the current model the proton LIS flux that was used is taken from [13].

Application to AMS02 and PAMELA data
The numerical model was applied to the Bartel-rotation averaged AMS02 data from May 2011 to May 2017 [1] and to the Bartel-rotation averaged PAMELA data from June 2006 to January 2014 [2]. The model has six free parameters that can be varied to fit the measured data sets. They are: tilt angle α(t), HMF at Earth B 0 (t), the HMF polarity A(t), the parallel diffusion normalization coefficient K 0 (t), the two spectral indexes of the diffusion coefficient, a(t) and b(t). The first three define the average status of heliosphere at a given epoch, the last three fit the data sets and define the modulation of CR at the same epoch. The break R k is fixed at R k = 3 GV and R Ak =0.3 GV, since it was seen that can be treated as constant in time. As a first approximation, to take into account the time-varying status of the heliosphere, we take the averaged tilt angle and HMF over a time period of 1 year preceding the selected BR, using a Backward Moving Average (BMA). To fit the measured data sets of CR of AMS02 and PAMELA, a six dimensional discrete grid of the model parameters vector q = (α, B 0 , A, K 0 , a, b) was built with the parameters values shown in the Table 1. For each of the 1.215 × 10 6 points in the grid, a proton flux J m (E, q) was simulated. Then, for each data set J d (E, t), a χ 2 was defined as The errors are given by at a given epoch, defined as the time t of the beginning of the corresponding BR, the first three parameters, α(t), B 0 (t) and A(t) fix the average status of the heliosphere. For x = K 0 (t), a(t) and b(t)), the minimum chi-squared curve, χ 2 min (x), is found, regardless of the values of all the other parameters, as shown in Figure 1a for the spectral index below the break a for one data set, where the χ 2 min (x) is represented by the black dots. Then, an interpolation is done on the χ 2 min (x) to find its miminum and the value of the parameter x at minimum, 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 , as shown in 1b, where the solid line indicates x best and the dashed lines x ± . The model was applied to the Bartel-Rotation averaged proton fluxes measured by the PAMELA and AMS02 experiments [1,2]. The data time-series cover the solar minimum from 2006 to 2009, the ascending phase to solar maximum, occurred on ∼ February 2013, when the Sun polarity A reversed, and the following descending phase until May 2017. The comparison to the data, |J d −J best | σ , shows that the difference between any data set and the corresponding best flux model are all well below 1. In Figure 2 are shown the model parameters K 0 , a and b as a function of time for the data sets of AMS02 and PAMELA. It is clearly visible a correlation to the solar activity measured by the averaged Sun Spot Number (SSN) with a time lag between the SSN and the model parameters. To find the time lag between the model parameters and the SSN we compare K 0 (t) with SSN(t -∆t) varying ∆t and find the ∆t that maximizes the correlation. It turns out to be ∆t = 1.04 ± 0.1 years. K 0 is anti-correlated with SSN when time lag is taken into account, the slopes a and b are less clearly correlated with SSN but a variation can be seen between solar max and solar min.

Conclusions
A new modeling of the proton flux has been implemented. It successfully describes AMS02 and PAMELA data covering 11 years of CR measurements. A clear correlation between solar activity and CR modulation is found. The model array can be used to find the modulation parameters at any epoch.