Cosmic-ray transport in the heliosphere with HelioProp
Andrea Vittino Carmelo Evoli Daniele Gaggero
astro-ph.HE hep-ph physics.space-ph





§ INTRODUCTION


Before reaching the site of their detection at Earth, Cosmic Rays (CRs) propagate across the Solar System, where they are subject to the interaction with the turbulent solar wind and with the Heliospheric Magnetic Field (HMF). The ensemble of effects that arise as the result of these interactions can significantly alter CR intensity, in particular at low energies. This mechanism goes typically under the name of solar modulation (we address the reader to <cit.> for a review). 

An accurate description of solar modulation represents a necessary ingredient in the study of CR properties, also in connection with the search for exotic CR sources, such as Dark Matter (DM). This is the reason why we have developed HelioProp, a tool to model the transport of CRs in the heliosphere by means of the realistic three-dimensional model based on the Stochastic Differential Equation (SDE) technique described in <cit.>. 

Results based on an earlier version of the code have been reported in <cit.>. The analyses described here have been conducted with a novel version of HelioProp which is currently under development and which will be released soon as a completely open-source package. In addition, HelioProp will be part of the Dragon2 project <cit.>. Together, the two codes aim at describing CR transport from the source to the observer under conditions that are as general and realistic as possible. 

This proceeding is organised as follows: in Section <ref> we briefly describe the most important ingredients that characterise the physics of solar modulation, while in Section <ref> we illustrate the main features of the numerical method that is adopted in HelioProp to model CR transport in the heliosphere. In Section <ref> we present some results obtained with the code, in particular for what concerns the modelling of CR proton and electron fluxes and the impact of solar modulation on DM searches. Lastly, in Section <ref> we report our conclusions.
 



§ THE PHYSICS OF SOLAR MODULATION


Following the model proposed by Parker <cit.>, the HMF is typically assumed to possess a spiral structure. If we consider a spherical coordinate system {r,θ,ϕ} with origin at the Sun location, the HMF can be written as:

    B⃗(r,θ,ϕ) = A_c B_0 ( r_0/r)^2 [ r̂ - ( Ω r  sinθ/V_sw) ϕ̂],
 
where r_0 = 1 AU and Ω is the differential rotation rate of the Sun, while V_sw represents the velocity of the solar wind, which is a radial flow continuously directed outward from the Sun. All the results that are shown in this work have been derived by assuming for the radial and latitudinal dependence of V_sw the expression given in <cit.>.  Concerning the parameters that define the normalization of the HMF, B_0 can be related to the intensity of the HMF at Earth B_e,




while A_c is the HMF polarity, defined as A_c = ± H (θ - θ') where the sign is known to shift with a period of roughly 11 years, while H is the Heaviside step function and θ' is the angular extent of the wavy Heliospheric Current Sheet (HCS). This quantity is defined as follows:

    θ' = π/2 + sin^-1[ sinα sin( ϕ - ϕ_0 + Ω r/V_sw)]
 
where ϕ_0 is an (arbitrary) constant azimuthal phase and α∈ [0,π/2] is the HCS tilt angle <cit.>, which is a quantity that varies according to solar activity and reaches its lowest values (α≈ 10^∘) in solar minima periods.
 
The transport of Cosmic-Rays (CRs) across the heliosphere is described in terms of a transport equation <cit.>:

    ∂ f/∂ t = - V⃗_sw·∇ f  + ⟨v⃗_d⟩·∇ f + ∇· (K⃗_⃗s⃗·∇ f) + 1/3 (∇·V⃗_sw ) ∂ f/∂ln R

where R=pc/eZ is the rigidity of the CR particle under consideration while f is its distribution function, related to the CR differential intensity j by f = R^2j. The four terms on the r.h.s. are used to model, respectively, convection under the influence of the solar wind, drifts, spatial diffusion and adiabatic energy changes. 

Concerning the process of CR drifts, it is modelled through the drift velocity (averaged with respect to the pitch angle) ⟨v⃗_d⟩ that appears in Eq. (<ref>). Since drifts are caused both by the interaction of CRs with the gradient and curvature of the HMF and by the change in the field polarity at the crossing of the HCS, the drift velocity can be written as the sum of two components: ⟨v⃗_d⟩ = ⟨v⃗_d⟩_GC + ⟨v⃗_d⟩_HCS. These two components are implemented in HelioProp according to the prescriptions given in <cit.>, to which we address the reader for additional details. It is important to remark that, following <cit.>, we include in our model a progressive weakening of drifts towards small rigidities. This is done by multiplying the drift velocity by the reduction factor:


    f(R)_drifts = (R/R_d)^2/1+(R/R_d)^2


Concerning the process of CR spatial diffusion, it is related to the tensor K⃗_⃗s⃗, which represents the symmetric part of the diffusion tensor. More precisely, one can work in the reference frame of the average HMF (frame that is denoted here with the superscript B) and decompose the total diffusion tensor into a symmetric and an antisymmetric component as follows:   

    K^B_ij = K^B,S_ij + K^B,A_ij =
    [   k_∥     0     0;     0 k_⊥,r     0;     0     0 k_⊥,θ;       ]
    +
    [    0    0    0;    0    0  k_A;    0 -k_A    0;      ],

where k_∥ is the diffusion coefficient along a direction parallel to the one of the HMF, while k_⊥,r and k_⊥,θ are the diffusion coefficients along the radial and polar perpendicular directions. The antisymmetric part of the diffusion tensor K^B,A_ij accounts for drifts and is already included in Eq. (<ref>) in terms of the drift velocity v⃗_ d = ∇×(K_A B⃗/|B|). On the other hand, as already said, the symmetric tensor K^B,S_ij, once translated into its counterpart in spherical coordinates K^S_ij, enters the transport equation to model spatial diffusion. As one can from Eq. (<ref>), diffusion is totally described once that the coefficients k_∥, k_⊥,r and k_⊥,θ are determined. As typical for the most realistic descriptions of solar modulation, HelioProp features a fully anisotropic diffusion, which means that it works under the assumption that k_∥ k_⊥,r k_⊥,θ. We express the three diffusion coefficients in terms of the Mean Free Paths (MFPs) in the associated direction k_i = v/3 λ_i with v being the CR velocity. We follow a common approach (see <cit.> for example) and define the parallel mean-free-path as a broken power-law with break at rigidity R_k and indices a and b below and above the break: 

    λ_∥ = λ_0 (1 nT/B) ( R/1 GV)^a(( R_k/1 GV)^c + ( R_k/1 GV)^c /1 + ( R_k/1 GV)^c)^b-a/c.

As for the MFPs along the perpendicular directions, HelioProp allow the user to implement any generic function. In this work, we assume: 

    λ_⊥,r =f_⊥λ_∥ ,     λ_⊥,θ =f_⊥λ_∥ H(θ) ,


with the function H(θ) being defined as in <cit.>. 




§ NUMERICAL SOLUTION OF THE TRANSPORT EQUATION



HelioProp models solar modulation by numerically solving Eq. (<ref>). This is done by following the SDE technique described in <cit.>, which was originally proposed in the context of CR solar modulation in <cit.>. In brief, this method is realised by writing the transport equation as a backward Kolmogorov equation: 

    -∂ j/∂ s = ∑_i ( A_i∂ j/∂ x_i) + 1/2∑_i,k( C_ik∂^2 j/∂ x_i ∂ x_k )
 
where j is the CR differential intensity, while the parameter s is a backward time, related to the standard time t by t = t_fin - s, where t_fin is the final time. The Itō's lemma <cit.> states that the above equation is equivalent to a set of stochastic processes: 

    dx_i = A_i(x_i)ds + ∑_j B_ij (x_i) dW_i

where we have used C = B^T B, while x_i = {r,θ,ϕ,E} and W_i is a Wiener process, i.e. a stochastic variable for which the increments Δ W_i = w_i - w_i-1 are independent one from the other and each of them has a gaussian probability distribution with zero mean and standard deviation equal to Δ s_i = t_i - s_s-1. We address the reader to <cit.> for the explicit form of the stochastic processes described by Eq. (<ref>).

Under a practical point of view, when following this approach one models solar modulation by combining the backward trajectories in the phase space of N pseudo-particles. More precisely, these pseudo-particles are injected at Earth at s=0 with phase space coordinates x_Earth and then they are back-propagated through the heliosphere by following the stochastic processes described by Eq. (<ref>) until they reach the boundary of the heliopause (HP). 

The CR flux at Earth is then obtained by averaging the Local Interstellar (LIS) spectra evaluated at the coordinates x_HP that the pseudo-particles have at the HP: 

    j_Earth (x_Earth) = 1/N∑_k=1^N j_LIS(x_HP)
  






§ RESULTS




 §.§ Proton and electron fluxes





The first test that we perform consists in comparing the results of HelioProp with the proton and electron fluxes measured by PAMELA over short time intervals (one month for the protons, six months in the case of electrons) along the period of minimal solar activity and negative polarity of the HMF that goes from 2006 to 2009 <cit.>[Because of space limitations, here we restrict ourselves to 4 electron datasets out of the 7 released by the PAMELA Collaboration. ]. 

For each data-taking period that we consider, we fix the tilt angle α and the HMF at Earth B_e to the values reported in <cit.>. We are thus left with a set of free parameters that include the quantities λ_0, a, b, c and R_k that appear in Eq. (<ref>), the factor f_⊥ of Eq. (<ref>) and the rigidity R_d of Eq. (<ref>). For simplicity, we impose that electrons and protons share the same values for both R_d and f_⊥. In particular, we find that good fits to data can be found by assuming R_d = 0.32 GV, which is the same value used in <cit.>, and f_⊥ = 0.021, that is very close to the commonly used value 0.02 (see <cit.>).  We assume the parallel MFP λ_∥ of protons to have a linear dependence on the CR rigidity as in the model of <cit.>. This means that, for protons, a = b = 1 (and therefore c and R_k do not play any role, as clear from Eq. (<ref>)). In the case of electrons, we find that we are able to reproduce the data behaviour by assuming a = 0, b = 1.55, c=3.5. A rigidity-dependence of the electron MFP in terms of a  broken power-law of this kind (i.e. with a flattening at small rigidities) appears to be in a qualitative agreement with previous analyses of electron data, as the one reported in <cit.>. With this considered, the only parameters that are left free to vary from one dataset to the other are λ_0 (for both protons and electrons) and R_k (only for electrons). Concerning the proton and electron LIS fluxes, we adopt the parameterizations based on Voyager 1 and PAMELA data presented in <cit.>.  



A summary of the relevant parameters used in this analysis is reported in Table <ref>. Results are reported in Fig. <ref>, In the top row we show PAMELA data compared with the best-fit configurations obtained within the framework of our model, while in the bottom row the corresponding parallel MFPs at Earth are shown as a function of the CR rigidity. 

As it can be seen, the model that we are using can reproduce the observed behaviour of both protons and electrons. It is important to remark that the approach followed here is data-driven, since we fix the values of the the parameters characterising CR diffusion and drifts by fitting PAMELA data; as discussed in Section <ref>, HelioProp allows for a large freedom in the definition of the solar modulation setup and therefore one could instead follow a more theoretically motivated strategy and adopt for these parameters the results obtained in the context of turbulence transport models (see, e.g.,  <cit.>). 









 §.§ Antiprotons and consequences for DM searches



In the context of indirect DM searches in charged CRs, it is customary to treat solar modulation in terms of the force field approximation <cit.>. Here we provide a simple example of what can be the impact of adopting a more realistic solar modulation setup, as the one implemented in HelioProp. We do this by comparing predictions of DM-generated antiproton fluxes obtained with HelioProp and with the force-field approach. Analyses along these lines have been already performed with the previous version of HelioProp and can be found in <cit.> (for the case of antiprotons) and in <cit.> (for the case of antideuterons). 

Since we are dealing with an exotic source of antiprotons contributing mostly at low energies, we cannot tune the solar modulation parameters by directly using antiproton observations and we must rely on proton data alone. Obviously, we have to consider a proton and an antiproton datasets derived from observations performed over the same time interval. This is the case of the PAMELA datasets presented in  <cit.>. We find that we can fit the PAMELA proton dataset by using the same setup described in Section <ref> (with λ_∥ = 0.016 AU at 100 MV at Earth). In the framework of the force field approximation, a comparably good fit is provided by assuming a force-field potential of 450 MV. 





As for antiprotons, we adopt a LIS derived from a DRAGON run performed within the KRA propagation model (see <cit.> for details about this propagation setup)[It is important to remark that the DRAGON run has been performed by tuning the proton injection spectrum in order to match the LIS used for protons in this work.]. The antiproton fluxes at Earth obtained by modulating this LIS with HelioProp and with the FF approach are shown in Fig. <ref>: as one can see, the difference between the two modulation setups is minimal (below 10%) in the case of the antiproton flux produced by spallation processes. On the contrary, if one considers the antiproton flux produced by the annihilation of light DM, as we do here by assuming a DM particle with mass 10 or 20 GeV that annihilates into bb̅, the difference is larger and can reach the 35%. This is shown in the right panel of Fig. <ref>, where it is also shown that such a distance can be larger than the experimental uncertainty. This proves that when deriving constraints to DM properties from low-energy CR data, a realistic treatment of solar modulation is necessary.    







§ CONCLUSIONS


We have presented here the new version of HelioProp the numerical code designed to model solar modulation in a three-dimensional and charge-dependent way by using the SDE approach described in <cit.>.  HelioProp is designed to be part of the Dragon2 project <cit.>, but will also be released soon as an independent and fully open source package. 

We have briefly illustrated the main features of the code, under a physical and numerical point of view. In addition, we have presented a couple of basic tests aimed at investigating its performances and potentialities. In particular, we have shown that HelioProp is able to reproduce the observed behaviour of low-energy CR electrons and protons and we have illustrated how it can be an important tool for a more accurate treatment of DM indirect detection.   








