Revisit cosmic ray propagation by using $^{1}$H, $^{2}$H, $^{3}$He and $^{4}$He

The secondary-to-primary ratios are unique tools to investigate cosmic ray propagation mechanisms. In this work, we use the latest data of deuteron-to-helium~4 ratio and helium~3-to-helium~4 ratio measured by PAMELA combined with other Z$\leq$2 primary fluxes measured by PAMELA and Voyager-1, to constrain the cosmic ray acceleration and propagation models. The analysis is performed by interfacing statistical tools with the GALPROP propagation package. To better fit both the modulated and unmodulated low energy cosmic ray data, we find that a time-, charge- and rigidity-dependent solar modulation model is better than the force-field approximation. Among all the studied cosmic ray propagation models, the diffusion-reacceleration-convection model is strongly supported by the derived Bayesian evidence. The robustness of the estimated diffusion slope $\delta$ is cross-checked by another low-mass secondary-to-primary ratio, i.e. the antiproton-to-proton ratio. It is shown that the diffusion-reacceleration-convection model can reconcile well with the high energy antiproton-to-proton ratio. This indicates that the estimated value of $\delta$ is reliable. The well constraint $\delta$ from the `best' model is found to be close to 1/3, inferring a Kolmogorov-type interstellar magnetic turbulence.


I. INTRODUCTION
Understanding cosmic ray (CR) acceleration and propagation mechanisms is an attractive topic over the last century. The standard CR theory reveals that Galactic CRs are accelerated via Fermi acceleration in supernova shocks, consequently diffuse on turbulence in Galactic magnetic field and meanwhile possibly experience reacceleration and convection [1]. Though this paradigm is generally accepted, the details of the involved physical processes remain unclear. For instance, the spectral index of the diffusion coefficient presented in recent literatures differs from 0.2 to 0.6 [2][3][4], which could indicate a big difference on the properties of Galactic magnetic turbulence. Recently, benefiting from the launches of various space missions such as PAMELA [5], AMS-02 [6] and DAMPE [7], CR physics has entered a data-driven epoch. These new instruments circling around Earth for a long duration and suffering no residual atmosphere interference, eliminate the statistical and systematic uncertainties of the measurements. Besides, a spacecraft Voyager-1 launched at 1977, has reached the boundary of Heliosphere since 2012 and can directly observe interstellar CRs. By using the unprecedented precise data measured by these experiments, we are allowed to investigate the acceleration and propagation phenomena elaborately and reliably.
In order to probe CR propagation models, the most commonly used secondary-to-primary ratio in pervious theoretical studies was the boron-to-carbon (B/C) ratio. The diffusion-reacceleration (DR) model was usually proposed which could naturally explain the peak around 1 GeV in the B/C ratio (e.g. see [2,[8][9][10]). However, a discrepancy between the antiproton (p) data and the DR model prediction (see [10] and references therein) can't be ignored. It was suggested in [2] that the Z ≤ 2 CR species might suffer different propagation properties with respect to the heavier elements. Consequently, an evident problem is that using only the B/C ratio may bias our evaluation on CR transport models. Some studies argued that the "p discrepancy" could be caused by the inaccuratep production cross-sections and new cross-sections were proposed [11][12][13][14]. Incorporating different hadronic models can lead to an uncertainty onp production achieving 50% [15]. Since thep problem is far from certain, alternatively, it is very important and necessary to employ other low-mass secondaries to investigate CR propagation mechanisms. The deuterium-to-helium 4 ( 2 H/ 4 He) ratio and the helium 3-to-helium 4 ( 3 He/ 4 He) ratio measured by PAMELA give us this opportunity. In this paper, our main goal is using the 2 H/ 4 He and 3 He/ 4 He ratios as well as other Z ≤ 2 primary spectra to constrain CR acceleration and propagation models.

II. DATA
Combining the secondary-to-primary ratios and the primary spectra allows us to explore both propagation and acceleration processes. PAMELA has recently published the 2 H/ 4 He ratio from 0.1 GeV/n to 1.4 GeV/n and the 3 He/ 4 He ratio from 0.1 GeV/n to 1.1 GeV/n with excellent precision [16]. These data are used in our fittings. The measurements of the 2 H/ 4 He and 3 He/ 4 He ratios at higher energies are only measured by IMAX92 [17,18], but with very limited accuracies. They are not used here. For the primary spectra, PAMELA has measured the proton (p) flux from 400 MeV to 1.2 TeV and the helium (He) flux from 100 MeV/n to 600 GeV/n [19]. These data will also be included in this work. We note that AMS-02 has released the p and He fluxes with better statistics [20,21]. But these data were taken from May of 2011 to November of 2013. A reversal of Heliospheric magnetic field (HMF) polarity from A < 0 to A > 0 occurred between October of 2012 and June of 2013. During this period, the HMF varied dramatically and the solar effects could be very complicated. In order to facilitate the solar modulation calculation, the AMS-02 p and He data are not utilized in this analysis.
Apart from the PAMELA data, we also use the CR spectra measured by Voyager-1 [22]. It is for the first time that scientists have observed interstellar CR spectra. Elemental intensities for the Z ≤ 2 particles provided by Voyager-1 include the p flux from 3 MeV/n to 346 MeV/n and the He flux from 3 MeV/n to 661 MeV/n [23,24]. The prominent roles of the Voyager-1 data can be stressed as following. First, CRs measured by Voyager-1 suffer no burden from Solar wind, which can solely constrain CR acceleration and propagation models without considering the influence of solar modulation. Second, by comparing the modulated data to the unmodulated data, the solar modulation effect can be probed straightforwardly.

III. SOLAR MODULATION
Solar modulation is one of the main obstacles for us to extract message from CR data observed on Earth. In literatures, a most frequently used model for heliospheric modulation was the force-field approximation [25]. This model depends on one free parameter, the so-called Fisk potential Φ. Though this model has been widely used due to its easiness in usage, it does not give a reliable description on the heliospheric transport. Recently, some studies treated this phenomenon in a more reasonable way [26][27][28][29][30]. For example, in [30], a solar propagation code HelMod [31,32] was incorporated with the CR transport code GALPROP [33,34] to calculate the CR spectra. Their procedure was divided into two steps. First, with fixed HelMod parameters, the data with rigidities larger than 2 GV were used to evaluate the best-fit GALPROP parameters. Second, a delicate inter-calibration between GALPROP and HelMod was made. This approach allowed a physical calculation of solar modulation. However, the calibrating step was complicated and time-consumed.
Another heliospheric model was proposed in [27], based on several well measured time-dependent solar observables. Considering a CR particle with charge q, rigidity R and velocity β = υ/c, this model constructs the potential Φ as: where A and |B tot | are the polarity and magnitude of the HMF measured at Earth, α (t) is the tilt angle of the heliospheric current sheet (HCS), R 0 is a reference rigidity, φ 0 and φ 1 are two normalization factors. It introduces a Heaviside step function H (−qA (t)) which equals to zero or unity, depending on the sign of qA (t). The second term of Equation 1 represents a drift movement along the HCS, which is the main process suffered by particles with qA (t) < 0. For qA (t) > 0, CRs travel rather directly from the polar regions of the heliopause to Earth, and H (−qA (t)) = 0. This model, referred as CM model in the following paragraphs, is physically motivated and can be interfaced with GALPROP as easy as the force-field approximation.
In this work, we attempt to interface the CM model with GALPROP to perform a simultaneous scan on all the involved parameters. By comparing the results evaluated separately from the CM model and from the force-field model, we can understand better the uncertainties caused by solar modulation. The PAMELA measurements employed in this work were derived from data collected between July of 2006 and December of 2008. As observed, |B tot | and α vary with time [35,36]. To determine each CR spectrum over the whole period, we calculate the modulated spectrum for each individual six-month duration and then average them. We set R 0 to 0.5 GV [27]. It was found in [37] that the estimated values of φ 1 varies significantly by considering different CR propagation models. The estimation of φ 0 may also be relevant with the CR source and propagation parameters chosen in the analysis. Therefore, the factors φ 0 and φ 1 are both allowed to be free in our fittings.

A. Model parameter description
To solve the CR propagation equation, we use the version r2766 [55] of GALPROP. We focus our study on the transport and source parameters characterising different processes described as following. Diffusion is a consequence of CRs interacting with magnetohydrodynamic turbulence. The diffusion coefficient D xx is estimated to be a power law in rigidity ρ [1,38], and is usually assumed as: where D 0 is the normalization of diffusion coefficient at a reference rigidity ρ 0 , β is the particle velocity, η is a low energy dependence factor which could possibly caused by the turbulence dissipation effect [38] and δ is the diffusion spectral index. The free parameters linked to diffusion include D 0 , δ and η. Diffusion happens not only in position space, but may also in momentum space, which results in reacceleration of CR particles. The associated diffusion coefficient in momentum space D pp is correlated to D xx , and is expressed as [39]: where v A is the velocity of fluctuations in the hydrodynamical plasma, called the Alfvén velocity. v A is the main free parameter related to reacceleration. Besides, CR particles may transport in bulk away from Galactic plane by the convection process. The convection velocity V c (z) is usually assumed to vary linearly with the distance from Galactic plane: here V (0) and dV /dz are the main free parameters. According to the diffusive shock acceleration theory, the injected density q for CR species i is assumed to follow a rigidity power law. The general form is given as: where f (R, z) is the spatial distribution of sources and N i is the normalization abundance for CR species i. Some studies suggested a broken power law of injection spectrum [10,24], which employed different injection indices ν 1 and ν 2 above and below a reference rigidity ρ br . Meanwhile, heavier elements may have different indices with protons. We assume that ν 1 and ν 2 are the injection indices for protons, and ν 1He and ν 2He are those for helium nuclei. Their differences are set as: The parameters ν 1 , ν 2 and ∆ He are allowed to vary freely. In GALPROP, the source abundance of protons N p is normalized based on the propagated proton spectrum at 100 GeV, and N i for other species are scaled by their abundances relative to that of protons. In our previous study [40], we found that for different models, N p varied slightly. Thus we set N p = 4.69 × 10 −9 cm −2 sr −1 s −1 MeV −1 to fit the PAMELA proton data. The helium abundance relative to proton is used as a free parameter, i.e. X He . As suggested in [41], the production cross-sections of 2 H and 3 He provided in [42] are adopted in this work to give better description of 2 H and 3 He production.
To summarize, our free parameters include the propagation parameters: D 0 , δ, η, v A , V (0), dV /dz; the source parameters: ν 1 , ν 2 , ∆ He , X He ; and the solar modulation parameters: φ 0 and φ 1 (for the CM model), or Φ (for the force-field model). A degeneracy exists between D 0 and the halo size of Galaxy z h . It can only be broken by unstable species. In this work, z h is fixed at 4 kpc to be consistent with earlier GALPROP studies [10,30,43] and to ease the comparison of results.

B. Statistical methods
To perform a quick understanding on how well a model reflects the observed data, the χ 2 "goodness of fit" test is used. The minimization library MINUIT [44] is interfaced with GALPROP for the analysis. This allows a direct recognization of improper models. For those models which can generally reproduce the data, a further Bayesian approach is used do a comparison on selected models. Given the data set D and the vector of free parameters Θ assuming a model H, Bayes' theorem states that: where P (Θ|D, H) represents the posterior probability density function (p.d.f.) of the parameters, P (D|Θ, H) = L (Θ) is the likelihood function, P (Θ|H) is the prior distribution. The Bayesian evidence Z = P (D|H) is a key ingredient to comparing alternative models. A larger value of Z will be achieved by a model which depends on fewer free parameters but fits better the data. A comparison between two hypotheses H 1 and H 2 can be performed by comparing their respective evidences. Assuming Z 2 /Z 1 > 1, model H 2 is favored versus model H 1 , and vice versa. If the logarithm of this ratio is larger than an empirical thresholds 5.0, it represents a strong evidence [45]. Furthermore, the p.d.f.s of the parameters can be derived by the Bayesian inference, which will help us to understand the uncertainties and correlations of parameters. To implement sufficiently fast Bayesian analysis, a package MultiNest [46][47][48] is utilized to study CR models.

V. RESULTS AND DISCUSSIONS
Different types of propagation models are explored in this work: (1) the plain diffusion model (PD), (2) the plain diffusion model with an ad hoc break in the diffusion coefficient (PDbr), (3) the diffusion-reacceleration model (DR), (4) the diffusion-convection model (DC); (5) the diffusion-reacceleration-convection model (DRC). All models studied here consider a break in the injection spectrum. The data employed in the analysis only contain Z≤2 data. Since 2 H and 3 He are mainly produced via interactions of 1 H and 4 He with interstellar medium, we let the nuclear chain start from 4 He to economize the calculation consumption. A test on DRC model was carried out by setting nuclear chain begin from 30 Si. It was found that the best-fit parameters only vary a few percent. This implies that ignoring the influence of heavy elements is an acceptable choice.

A. A comparison between different solar modulation models
The CM model and the force-field model are interfaced separately with GALPROP. All the related models add a suffix "-CM" and a suffix "-FF" respectively. The results show that by incorporating the force-field approximation, none of the propagation models give a value of reduced χ 2 close to 1. Tracing back to the contribution of each data set to the overall χ 2 , it seems that all the models are not able to simultaneously fit the modulated and unmodulated fluxes. On the other hand, for the same propagation configuration, assuming a CM scenario leads to a much smaller χ 2 value. This suggests that compared with the force-field approximation, the rigidity-dependent CM model can describe the data better.
Since the secondary-to-primary ratios used in our fitting procedure only cover low energy data, using an incorrect solar modulation approximation may result in a bias on the evaluation of CR acceleration and propagation models, as well as the involved parameters. For example, for the force-field scenario, among all the studied models, the DR-FF model gives the smallest χ 2 = 2.14. It is found that an inclusion of a convection effect does not improve the goodness of fit. But this is not the case for the CM scenario, in which the DRC-CM model gives the smallest χ 2 = 1.21, as shown in Table I. This model provides the best goodness of fit and prefers a strong convection with dV /dz = 30 ± 5 km s −1 kpc −1 .
Considering a typical value of φ 1 around 8∼12 GV given in Table I, it is clear that the rigidity dependency of the CM model substantially influences the CR particles with rigidities below a few GV. At larger rigidities, since the second term of Equation 1 only contributes a few percent to the total value of Φ (R), the CM scenario is almost equivalent with the force-field model. As we known, the index of CR spectrum at high energies only depends on the value of ν 2 + δ 2 . But it is unclear which mechanisms shape the low energy spectra. Using the CM solar modulation model may help us to clarify this uncertain question.     [1]. Unlike the PDbr-CM and DC-CM models, the DR-CM and DRC-CM models can generally reproduce all the data. It implies that the reacceleration process, which is usually suggested to explain the behavior of the low energy B/C ratio, is also preferred by the Z ≤ 2 data. Compared with the DR-CM model, the DRC-CM model with one more parameter achieves a smaller χ 2 value equal to 1.21, allowing a more satisfactory description of data. The DR-CM model and the DRC-CM model are further studied in the Bayesian analysis. Priors are specified on the free parameters to restrict them in physically reasonable regions, as shown in Table II. The results are summarized in Table III, including the logarithm of Bayesian evidence for each model and the posterior mean with standard deviation on parameters. The differences in log-evidence between these two models is ln (Z 2 /Z 1 ) → 24. An inclusion of convection significantly increases the evidence, suggesting that the DRC-CM model is the "best" one among all   Parameter evaluations from the Bayesian method are consistent with the results obtained from the χ 2 analysis. The "best" DRC-CM model yields a diffusion spectral index δ = 0.29 ± 0.04, encouraging a Kolmogorov-type magnetic turbulence. We test a different fit using default GALPROP 2 H and 3 He production cross-sections. The value of δ can be upshifted by 0.07, but with much worse goodness of fit. Even in this case, it still favors a Kolmogorov-type turbulence. The Alfvén velocity v A allowed for the DRC-CM model is around 30.6 ± 1.4 km/s. A strong reacceleration level is needed. In our analysis, V 0 is fixed at 0 km/s for convection models. This is not an arbitrary choice. If we free V 0 in the fit, it converges at zero. It hints that adopting a constant convection velocity may rest on an incorrect assumption. Nevertheless, in our work we fix V 0 to reduce the number of free parameters and diminish the calculation time. The estimated dV /dz value is equal to 31 ± 4 km s −1 kpc −1 , suggesting a strong convection velocity. Constraints are also placed on the injection parameters. The results show that a softening on injection spectrum at 7.3 ± 0.4 GV is necessary for the DRC-CM model to recover the data. For all the models, without a break on the spectral index, the derived χ 2 values are highly enlarged. This implies that a single power law is insufficient to reproduce the data. The posterior mean values of ν 1 and ν 2 are 1.766 ± 0.026 and 2.544 ± 0.027, respectively. A sum of δ + ν 2 is around 2.8 to hold the measured high energy proton flux.
The 1D and 2D marginalized posterior p.d.f.s of free propagation and source parameters for the DRC-CM model are given in FIG. 2. As shown, a negative correlation exists between η and δ. This is expected from Equation 2. In order to reproduce the low energy secondary-to-primary ratios, a larger η requires a smaller δ. The convection velocity dV /dz is found to be negatively correlated with D 0 . For CR particles, a fast diffusion needs a slow convection to maintain a sufficient transport duration in Galactic halo. The relationships between other transport parameters do not present arresting correlations. Considering the source parameters, a noticeable negative relation appears between ν 2 and X He . This is because a flatter injection spectrum will deviate more from the measured helium flux if it adopts a lower value of helium abundance normalization. The injection indices ν 2 and ν 1 exhibit a positive correlation. This can be understood by taking into account their relations with δ. To achieve the observed propagated CR slope, a larger    Table III). The dark/light purple area represents the 68%/95% credible interval.
δ leads to a flatter injection spectrum both at low rigidities and at high rigidities. This causes a positive correlation between ν 2 and ν 1 . The degeneracy between δ and ν 2 can be broken by using high energy secondary-to-primary ratios. In this work, to clarify the reliability of estimated δ, we use thep/p ratio to do a cross-check in Sect. V C.
The theoretical calculations in 68% and 95% credible intervals for the fitted CR spectra are shown in FIG. 3. It can be found that the DRC-CM model can satisfactorily recover the PAMELA data including the 2 H/ 4 He ratio, the 3 H/ 4 He ratio, the proton flux and the helium flux. Meanwhile, this model also accommodates the interstellar CR spectra released by Voyager-1. Our results are compared with a diffusive reacceleration model proposed in an earlier paper [2]. That model was claimed as the best-fit model for the low-mass CR elements. However, it shows dramatic disagreements with the observed 2 H/ 4 He, 3 H/ 4 He ratios and the unmodulated fluxes. The bias on their parameters is due to the employment of an inappropriate description on solar modulation, as discussed in Sect. V A.

C. Cross-check with thep/p ratio
The prediction for thep/p ratio is calculated. The defaultp production cross-sections embedded in GALPROP [8,49,50] is used here. As shown in FIG. 4, compared with the conventional DR model suggested in [10], which adopted a force-field approximation, the DRC-CM model (see Table I) yields 30%∼60% higherp production around a few GeV. It seems that an inclusion of a convection mechanism and incorporating a charge-dependent heliospheric modulation model both influence substantially to thep/p ratio. The DRC-CM model agrees fairly well with the PAMELA data below 1 GeV and above 10 GeV. At moderate energies, the DRC-CM model yields a slightly lower p/p ratio. This discrepancy might be caused by the inaccuracy ofp production cross-sections [2,15]. Nevertheless, at energies larger than 20 GeV, the calculatedp/p ratio for the DRC-CM is in an excellent agreement with the observations by both PAMELA and AMS-02 [51]. This is an indication that the derived value of δ remains robust for all the Z≤2 particles.

D. Comparison with other analyses
Our results agree fairly well with the estimates given in [3]. By using thep/p ratio and the Z≤2 primaries, they found δ = 0.30 ± 0.03 0.02 (stat) ± 0.10 0.04 (sys) and v A = 25 ± 2 km/s. Compared with the values of δ obtained by using the B/C ratio, for example δ = 0.395 ± 0.025 and ∼ 0.35 as presented in [30] and [52], our estimate is slightly lower. Taking into account the possible systematic uncertainties in production cross-sections and/or in solar modulation, light and heavier nuclei may still have compatible diffusion process. Concerning the convection process, a constant convection velocity was assumed in [3]. And in [30], they found the convection velocity with V 0 = 12.4 ± 0.8 km/s and dV /dz = 10.2 ± 0.7 km s −1 kpc −1 . In contrast, we find that V 0 converges at 0 km/s and a larger dV /dz is required. In [30], we note that their constraints on transport parameters depended on the approximate solar modulation parameters. As a consequence, a different choice on modulation parameters might lead to a different estimation on CR propagation parameters. In our analysis, we free simultaneously the solar modulation parameters. The difference  Table III. The dark/light purple area represents the 68%/95% credible interval. The blue line is computed using the best-fit parameters derived in [2] by using only low-mass data (p, He,p). In bottom panels, the upper purple band and the dashed line represent the calculated unmodulated spectra.
in the fitting procedure may be responsible for the inconsistencies on convection. Anyway, our result on V 0 is more appealing for the sake of the continuity on convection velocity at Galactic plane.

VI. CONCLUSION
For the first time, we use the PAMELA measurements of hydrogen and helium isotopes to perform an elaborate study and a comparison on cosmic ray propagation models. Accurate and reliable constraints are generated by using the 2 H/ 4 He and 3 He/ 4 He ratios combined with the Z≤2 primaries. Serval hints can be inferred from our results. First, a global fit on both the modulate and unmodulated data disfavors the widely used force-field approximation. Instead, the tension between CR propagation models and the observations can be released by incorporating the time-, charge-, rigidity-dependent CM solar modulation model. Second, various models are comprehensively studied and compared by using statistical methods. From a quick χ 2 minimization analysis, it is found that only models including the reacceleration effect can accommodate individual ratios and fluxes. A Bayesian analysis provides evidences which strongly support the DRC-CM model over the DR-CM model to explain the Z≤2 data. Thep/p ratio not involving in the fitting procedure further demonstrates the plausibility of our DRC-CM model and the robustness of the estimated diffusion spectral index δ. Third, the value of δ is estimated to be close to 1/3, which indicates a Kolmogorov-type turbulence. Even taking into account the possible systematic errors might introduced in 2 H and 3 He production cross-sections, the Kraichnan type of turbulence is still disfavored.
This work improves our knowledge on CR acceleration and propagation mechanisms. A further study utilizing heavier elements will be implemented in future. AMS-02 has published the B/C ratio [53] and the carbon flux [54] with much better precision than PAMELA. But we need to realize that the AMS-02 measurements covered a complicated  Table III. The dark/light purple area represents the 68%/95% credible interval for the DRC-CM model. The red line is derived by using the best-fit parameters for the DRC-CM model. The blue line presents the predictions for the DR-CM. The green line shows the theoretical calculation for the best-fit DR model derived in [10].
period during which HMF polarity traversed from A<0 to A>0. The solar modulation effect for this period will be explored by combining the forthcoming AMS-02 periodic fluxes.