Propagation of cosmic rays in the AMS-02 era

TL;DR

This study uses recent cosmic ray data from AMS-02 and PAMELA to constrain propagation models, finding that models with reacceleration best fit the data but reveal differences in propagation for nuclei and leptons.

Contribution

It introduces a comprehensive comparison of propagation models using new cosmic ray data, highlighting the importance of reacceleration in fitting proton and B/C observations.

Findings

Reacceleration models fit proton and B/C data well.

Diffusion-reacceleration models predict too many low-energy positrons.

Non-reacceleration models better fit low-energy antiproton fluxes.

Abstract

In this work we use the newly reported Boron-to-Carbon ratio (B/C) from AMS-02 and the time-dependent proton fluxes from PAMELA and AMS-02 to constrain the source and propagation parameters of cosmic rays in the Milky Way. A linear correlation of the solar modulation parameter with solar activities is assumed to account for the time-varying cosmic ray fluxes. A comprehensive set of propagation models, with/without reacceleration or convection, have been discussed and compared. We find that only the models with reacceleration can self-consistently fit both the proton and B/C data. The rigidity dependence slope of the diffusion coefficient, $\delta$, is found to be about $0.38-0.50$ for the diffusion-reacceleration models. The plain diffusion and diffusion-convection models fit the data poorly. We compare different model predictions of the positron and antiproton fluxes with the data. We find that the diffusion-reacceleration models over-produce low energy positrons, while non-reacceleration models give better fit to the data. As for antiprotons, reacceleration models tend to under-predict low energy antiproton fluxes, unless a phenomenological modification of the velocity-dependence of the diffusion coefficient is applied. Our results suggest that there could be important differences of the propagation for nuclei and leptons, in either the Milky Way or the solar heliosphere.

Tables

Table 1. Table 1: Data taking time of various measurements and the average modelled sunspot numbers one year before the data taking time.

	Time	$\overline{N}$
ACE(10Be/9Be)	08/1997-04/1999	23.5
ACE(B/C)	05/2011-05/2016	54.3
PAMELA-2006($p$)	11/2006	17.4
PAMELA-2007($p$)	12/2007	7.3
PAMELA-2008($p$)	12/2008	3.0
PAMELA-2009($p$)	12/2009	1.0
AMS-02($p$)	05/2011-11/2013	40.8
AMS-02(B/C)	05/2011-05/2016	54.3
PAMELA($\overline{p}$)	07/2006-12/2008	10.0
AMS-02($e^{+}$)	05/2011-11/2013	40.8
AMS-02($\overline{p}$)	05/2011-05/2015	51.5

Table 2. Table 2: Posterior mean and 68 % percent 68 68\% credible uncertainties of the model parameters

	Unit	PD	DC	DC2	DR	DR2	DRC
$D_0$	$\left({10}^{28}{\mathrm{cm}}^2{s}^{-1}\right)$	$5.29\pm 0.51$	$4.20\pm 0.30$	$4.95\pm 0.35$	$7.24\pm 0.97$	$4.16\pm 0.57$	$6.14\pm 0.45$
$\delta$		$0.471\pm 0.006$	$0.588\pm 0.013$	$0.591\pm 0.011$	$0.380\pm 0.007$	$0.500\pm 0.012$	$0.478\pm 0.013$
$z_h$	$\left(\mathrm{kpc}\right)$	$6.61\pm 0.98$	$10.90\pm 1.60$	$10.80\pm 1.30$	$5.93\pm 1.13$	$5.02\pm 0.86$	$12.70\pm 1.40$
$v_A$	$\left(\mathrm{km}{s}^{-1}\right)$	—	—	—	$38.5\pm 1.3$	$18.4\pm 2.0$	$43.2\pm 1.2$
$d{V}_c/dz$	$\left(\mathrm{km}{s}^{-1}{\mathrm{kpc}}^{-1}\right)$	—	$5.36\pm 0.64$	$5.02\pm 0.55$	—	—	$11.99\pm 1.26$
$R_0$	$\left(\mathrm{GV}\right)$	—	—	$5.29\pm 0.23$	—	—	—
$\eta$		—	—	—	—	$-1.28\pm 0.22$	—
$\log \left(A_p\right)$666Propagated flux normalization at 100 GeV in unit of ${\mathrm{cm}}^{-2}{s}^{-1}{\mathrm{sr}}^{-1}{\mathrm{MeV}}^{-1}$		$-8.334\pm 0.003$	$-8.334\pm 0.003$	$-8.336\pm 0.003$	$-8.347\pm 0.002$	$-8.334\pm 0.002$	$-8.345\pm 0.002$
${\nu }_1$		$2.44\pm 0.01$	$2.45\pm 0.01$	$2.43\pm 0.01$	$1.69\pm 0.02$	$2.04\pm 0.03$	$1.82\pm 0.02$
${\nu }_2$		$2.34\pm 0.03$	$2.30\pm 0.01$	$2.30\pm 0.01$	$2.37\pm 0.01$	$2.33\pm 0.01$	$2.37\pm 0.01$
$\log \left(R_{br}\right)$777Break rigidity of proton injection spectrum in unit of $\mathrm{MV}$		$5.06\pm 0.13$	$4.82\pm 0.05$	$4.78\pm 0.06$	$4.11\pm 0.02$	$4.03\pm 0.03$	$4.22\pm 0.03$
${\Phi }_0$	$\left(\mathrm{GV}\right)$	$0.595\pm 0.005$	$0.537\pm 0.006$	$0.419\pm 0.005$	$0.180\pm 0.008$	$0.290\pm 0.014$	$0.220\pm 0.008$
${\Phi }_1$	$\left(\mathrm{GV}\right)$	$0.495\pm 0.011$	$0.485\pm 0.011$	$0.472\pm 0.012$	$0.487\pm 0.011$	$0.485\pm 0.011$	$0.482\pm 0.013$
${\chi }^2/\mathrm{dof}$		748.6/463	591.0/462	494.6/461	438.8/462	341.0/461	380.5/461

Table 3. Table 3: Summary of different propagation models versus the data

	B/C & protons	positrons	antiprotons
PD	Poor	Too few	Fair
DC	Poor	Too few	Fair
DC2	Poor	Too few	Good
DR	Good	Too many	Slightly few
DR2	Good	Too many	Fair
DRC	Good	Too many	Slightly few

Equations

$$\begin{array}{rl}\frac{\partial \psi }{\partial t}=& Q(\mathbf{x},p)+\nabla \cdot \left(D_{xx}\nabla \psi -{\mathbf{V}}_c\psi \right)+\frac{\partial }{\partial p}{p}^2{D}_{pp}\frac{\partial }{\partial p}\frac{1}{p^2}\psi \\ & -\frac{\partial }{\partial p}\left[\dot{p}\psi -\frac{p}{3}\left(\nabla \cdot {\mathbf{V}}_c\psi \right)\right]-\frac{\psi }{{\tau }_f}-\frac{\psi }{{\tau }_r},\end{array}$$

$$D_{xx}=D_0{\beta }^{\eta }{\left(\frac{R}{R_0}\right)}^{\delta },$$

$$D_{pp}{D}_{xx}=\frac{4p^2{v}_A^2}{3\delta (4-{\delta }^2)(4-\delta )\omega },$$

$$f(r,z)={\left(\frac{r}{r_{\odot }}\right)}^{1.25}\exp \left(-3.56\cdot \frac{r-r_{\odot }}{r_{\odot }}\right)\exp \left(-\frac{|z|}{z_s}\right),$$

$$q(R)\propto\left\{\begin{array}[]{ll}\left(R/R_{\mathrm{br}}\right)^{-\nu_{1}},&R<R_{\mathrm{br}}\\ \left(R/R_{\mathrm{br}}\right)^{-\nu_{2}},&R\geq R_{\mathrm{br}}\end{array}.\right.$$

$$\Phi ={\Phi }_0+{\Phi }_1\times \frac{N(t)}{N_{\max }},$$

Full text

I Introduction

The propagation of cosmic rays (CRs) in the Milky Way is a fundamental question to understand the origin and interactions of Galactic CRs. It also provides us a useful tool to probe the properties of the interstellar medium (ISM). It is well known that the charged CRs will propagate diffusively in the Galactic magnetic field, experiencing possibly the reacceleration, convection, spallation and energy loss processes Gaisser (1990); Strong et al. (2007). The propagation process can be described with the diffusive transport equation Berezinskii et al. (1990); Gaisser (1990). Depending on different simplifications, the transport equation can be solved analytically Webber et al. (1992); Bloemen et al. (1993); Maurin et al. (2001, 2002); Shibata et al. (2004). Also there were efforts to include most of the relevant processes and the observation-based astrophysical inputs, and to solve the propagation equation numerically, e.g., GALPROP Strong and Moskalenko (1998); Moskalenko and Strong (1998) and DRAGON Evoli et al. (2008).

To understand the propagation of CRs is not only important for the CR physics itself, but also the basis of searching for the exotic signal from particle dark matter. The propagation of CRs couples closely with the production, leading to the entanglement between source parameters and propagation parameters. Fortunately, the spallation of the CR nuclei when colliding with the ISM will produce secondary nuclei (with kinetic energy per nucleon unchanged). The ratio between those secondary nuclei and the parent nuclei will cancel out the source information, leaving basically the propagation effect. Widely used are the Boron-to-Carbon (B/C) and sub-Iron-to-Iron ((Sc+Ti+V)/Fe) ratios. The unstable-to-stable ratio of the secondary isotopes plays another important role to constrain the CR propagation. The unstable nuclei with lifetimes comparable to the diffusion time of the CRs, such as 10Be ($\tau=1.39\times 10^{6}$ yr) and 26Al ($\tau=7.17\times 10^{5}$ yr), can be used as the clocks to measure the residual time of CRs in the Milky Way halo.

Many works have been dedicated to using the secondary-to-primary ratios and the unstable-to-stable isotope ratios to constrain the CR propagation parameters (see e.g., Swordy et al. (1990); Mueller et al. (1991); Strong and Moskalenko (1998); Yanasak et al. (2001); Maurin et al. (2001); Lionetto et al. (2005); Ave et al. (2009); Pato et al. (2010); di Bernardo et al. (2010); Obermeier et al. (2012)). However, due to the large number of the model parameters and the degeneracy between different parameters, the investigation of the parameter space is incomplete and the conclusion might be biased. In addition, more and more data have been accumulated nowadays. It is necessary to combine different data sets in a statistical way. Recently several works employed the Markov Chain Monte Carlo (MCMC) method to try to take a full scan of the parameter space with large samples of the data Putze et al. (2009, 2010); Trotta et al. (2011); Jin et al. (2015); Jóhannesson et al. (2016); Korsmeier and Cuoco (2016); Feng et al. (2016). The MCMC method is known to be efficient for the minimization of high-dimensional problem and is widely used in different areas.

We have developed a tool, CosRayMC, through embeding the CR propagation code in the MCMC sampler Lewis and Bridle (2002), which have already been applied to the study of the CR lepton excesses Liu et al. (2010, 2012a, 2012b); Yuan et al. (2015); Yuan and Bi (2013, 2015). In light of the newly reported CR nuclei and B/C data by PAMELA and AMS-02, we apply this tool to re-visit the CR propagation and constrain the propagation parameters in this work. Compared with previous studies Trotta et al. (2011); Jin et al. (2015); Jóhannesson et al. (2016); Korsmeier and Cuoco (2016); Feng et al. (2016), we will present an extensive study of different propagation models, including the plain diffusion scenario, the diffusion reacceleration scenario and the diffusion convection scenario. Furthermore, we will employ a phenomenological treatment of the time-dependent solar modulation based on the solar activities. Finally, the predicted positron and antiproton fluxes of different propagation models will be compared with the data as a consistency check.

This paper is organized as follows. In Sec. II we define the propagation model configurations. In Sec. III we describe the fitting procedure. The fitting results and expectations of secondary positron and antiproton fluxes are presented in Sec. IV. We discuss our results in Sec. V, and finally conclude in Sec. VI.

II Propagation models

Galactic CRs are accelerated in cosmic accelerators such as supernova remnants and pulsars before they are injected into the ISM. During their propagation in the Galaxy, secondary particles can be produced by the collisions between primary CRs and the ISM. The propagation of CRs in the Galaxy is usually described by the diffusive transport equation

[TABLE]

where $\psi$ is the differential density of CR particles per momentum interval, $Q$ is the source term, $D_{xx}$ is the spatial diffusion coefficient, $\mathbf{V}_{c}$ is the convective velocity, $D_{pp}$ is the diffusion coefficient in the momentum space describing the reacceleration effect, $\dot{p}\equiv\mathrm{d}p/\mathrm{d}t$ is the momentum loss rate, $\tau_{f}$ and $\tau_{r}$ are correspondingly the time scales for nulear fragmentation and radioactive decay.

The diffusion coefficient is usually assumed to vary with rigidity by a power-law form

[TABLE]

where $D_{0}$ is the normalization factor, $R_{0}$ is a reference rigidity, $\delta$ is the power-law index which depends on the property of turbulence in the interstellar medium (ISM), $\beta$ is the velocity in unit of light speed, and $\eta$ is a phenomenological parameter describing the velocity dependence of the diffusion coefficient at low energies, which is generally to be 1. For single power-law form of $D_{xx}$, we fix $R_{0}$ to be 4 GV. For the broken power-law case (see below), $R_{0}$ is left to be free in the fitting.

We assume the convection velocity linearly and continuously vary from the Galactic disk to halo, $\mathbf{V}_{c}=\mathbf{z}\cdot\mathrm{d}V_{c}/\mathrm{d}z$, where $\mathbf{z}$ is the position vector in the vertical direction to the Galactic disk. Such a form can avoid the discontinuity at the Galactic plane.

The reacceleration effect would lead to a diffusion in the momentum space. Its diffusion coefficient in momentum space, $D_{pp}$, is related with the spatial diffusion coefficient as Seo and Ptuskin (1994)

[TABLE]

where the $v_{A}$ is the Alfven velocity and $\omega$ is the ratio of magnetohydrodynamic wave energy density to magnetic field energy density. Since $\omega$ can be effectively absorbed in $v_{A}$, we assume it to be 1.

The source function $Q(\mathbf{x},p)$ is expressed as $f(\mathbf{x})q(p)$, where $f(\mathbf{x})$ is the spatial distribution and $q(p)$ is the injection energy spectrum of CR sources. The spatial distribution is assumed to follow that of supernova remnants

[TABLE]

with parameters slightly adjusted to match the Galactic diffuse $\gamma$-ray emission and the ratio of H2 to CO Strong et al. (2007); Trotta et al. (2011); $r_{\odot}=8.5$ kpc is the distance from the Sun to the Galactic center; $z_{s}\approx 0.2$ kpc is the characteristic height of the Galactic disk. The nuclei injection spectrum is assumed to be a broken power-law function of rigidity

[TABLE]

Power-law form of particle spectrum is expected from the simple shock acceleration mechanism. However, it has been found that single power-law spectrum is somehow not enough to describe the observational data, especially when there is strong reacceleration of CRs Trotta et al. (2011). The observations of $\gamma$-ray emission from a few supernova remnants which are interacting with molecular clouds also suggest a broken-power law of CRs in/around the source Ackermann et al. (2013). Note that we neglect the potential second break at hundreds of GV of the CR nuclei Panov et al. (2007); Ahn et al. (2010); Adriani et al. (2011); Aguilar et al. (2015a, b), which is beyond the energy range we are interested in. Since we focus on the B/C ratio, the small difference between the spectra of protons and heavier nuclei Adriani et al. (2011); Aguilar et al. (2015b) is also neglected.

The diffusive nature of charged particles in the Milky Way has been well established Gaisser (1990). However, whether the reacceleration and/or convection plays significant roles in regulating the propagation of CRs is unclear. The widely existed galactic winds suggest that convective transport of CRs may be relevant Jokipii (1976). On the other hand, the observed peak of the B/C around $\sim 1$ GeV/n by HEAO-3 Engelmann et al. (1990) may require an effective reacceleration Moskalenko et al. (2002). While the reacceleration model can fit the B/C data, it would under-predict antiprotons Moskalenko et al. (2002). An adjustment of the $\eta$ parameter in the diffusion coefficient was introduced to solve such a discrepancy di Bernardo et al. (2010). The modification of the low energy diffusion coefficient is also physically motivated from the potential resonant interaction of CR particles and the magnetohydrodynamic (MHD) waves which results in dissipation of such waves Ptuskin et al. (2006).

In this work we will test all these kinds of models with the new observational data. Specifically, the propagation models include:

1. the plain diffusion (PD) model without reacceleration and convection,
2. the diffusion convection (DC) model, 3) the diffusion convection model with a break of the rigidity-dependence of the diffusion coefficient (with $\delta=0$ below the break rigidity $R_{0}$ Moskalenko et al. (2002); DC2), 4) the diffusion reacceleration (DR) model, 5) the diffusion reacceleration model with $\eta$ left free to fit (DR2), and 6) the diffusion reacceleration convection (DRC) model. The relevant propagation parameters are $(D_{0},\,\delta,\,z_{h},\,v_{A},\,dV_{c}/dz,\,R_{0},\,\eta)$.

We keep in mind that the above described propagation framework is actually simplfied. The diffusion coefficient may vary in the Milky Way due to different magnetic field distributions in the disk and halo (e.g., Tomassetti (2012); Guo et al. (2016)). In particular, CRs may be confined much longer around the sources than expected due to non-linear self-generation of MHD waves via streaming instability D’Angelo et al. (2016). These complications are less clear and beyond the scope of the current work. The caveat is that considering these effects may result in different results from the adopted framework (see e.g., Feng et al. (2016)).

III Fitting procedure

III.1 CosRayMC

The CosRayMC code is a combination of the numerical propagation code GALPROP333http://galprop.stanford.edu/ Strong and Moskalenko (1998); Moskalenko and Strong (1998) and the MCMC sampler (adapted from CosmoMC Lewis and Bridle (2002)). The MCMC technique is widely applied in astrophysics and cosmology to investigate the high-dimensional parameter space from observational data. It works in the Bayesian framework. The posterior probability of model parameters $\boldsymbol{\theta}$ in light of the observational data $D$ is ${\mathcal{P}}(\boldsymbol{\theta}|D)\propto{\mathcal{P}}(D|\boldsymbol{\theta}){\mathcal{P}}(\boldsymbol{\theta})$, where ${\mathcal{P}}(D|\boldsymbol{\theta})$ is the likelihood and ${\mathcal{P}}(\boldsymbol{\theta})$ is the prior probability of $\boldsymbol{\theta}$.

The Markov chain is generated following the Metropolis-Hastings algorithm. The general process is as follows. One first proposes a random step in the parameter space. Then the acceptance probability is calculated by the ratio of the target probabilities of this proposed point to the former one. If the proposed point is accepted, then repeat this procedure. Otherwise, go back to the former point and have another trial. The stationary distribution of the chain samples will approaches the target probability distribution ${\mathcal{P}}(\boldsymbol{\theta}|D)$. For more details, one can refer to Neal (1993); Gamerman (1997).

III.2 Data sets

We adopt the most recently available accurate data sets of CRs by PAMELA and AMS-02 in our fittings. For the B/C ratio, we employ the just-released data by AMS-02 which cover an energy range of hundreds of MeV/n to TeV/n Aguilar et al. (2016a). In order to have better constraints on the low energy behavior of the B/C ratio, we also employ the data from ACE-CRIS444http://www.srl.caltech.edu/ACE/ASC/level2/lvl2DATA_CRIS.html with the same period as that of AMS-02. To constrain the lifetime of CRs in the Galaxy, we also use the 10Be/9Be data from some old measurements: Ulysses Connell (1998), ACE Yanasak et al. (2001), Voyager Lukasiak (1999), IMP Simpson and Garcia-Munoz (1988), ISEE-3 Simpson and Garcia-Munoz (1988), and ISOMAX Hams et al. (2004). The proton fluxes are employed to constrain the injection parameters of CRs. As will be discussed in the next subsection, we will try to give a more reasonable treatment of the solar modulation effect, the time-dependent proton fluxes from 2006 to 2009 measured by PAMELA Adriani et al. (2013) and the average flux from 2011 to 2013 by AMS-02 Aguilar et al. (2015a) are used. Table 1 summarizes the observational time of each data sets.

III.3 Solar modulation

In this work we use the force-field approximation to account for the solar modulation of low energy CRs when propagating in the heliosphere Gleeson and Axford (1968). However, since the various data sets in our work cover a wide time window in which solar activities varies much, they should not share a common modulation potential. Fig. 1 shows the sunspot numbers of different time from 1995 to present555https://solarscience.msfc.nasa.gov/SunspotCycle.shtml. The data we use are basically from the end of solar cycle 23 to the begining of solar cycle 24, except for the 10Be/9Be data. More importantly they are roughly in the period that the polarity of the solar magnetic field is in the same $A^{-}$ cycle. This enables us to have a relatively simple approach of the solar modulation with a correlation with solar activities.

Here we employ a linear evolution behavior of the modulation potential with respect to the evolution of the sunspot number

[TABLE]

where $N_{\rm max}\approx 72.2$ is the model predicted maximum sunspot number in solar cycle 24 (shown by the solid line in Fig. 1; Hathaway et al. (1999)), $N(t)$ is the sunspot number during which the data were collected, $\Phi_{0}$ and $\Phi_{1}$ are free parameters which will be derived through fitting to the CR data. The average sunspot numbers for various CR data taking time are given in Table 1. Note that we always count the sunspot number for the time one year before the actual data taking time, due to the possible delay of the modulation effect compared with solar activity. This treatment is consistent with the fact that the PAMELA proton flux in 12/2009 is higher than that in 12/2008, while the solar minimum of cycle 23 ended at the beginning of 2009. Given a typical speed of $\sim 500$ km/s, solar winds need about one year to fill the heliosphere with a scale of $\sim$100 astronomical units, which further supports our treatment.

IV Results

IV.1 Fitting results of various models

We use the MCMC algorithm to determine the model parameters of the six models as described in Sec. II through fitting to the data. The posterior mean and $68\%$ credible uncertainties of the model parameters are given in Table 2. Since the data are precise enough, we obtain statistically good constraints on the model parameters. Some of the model parameters, such as the injection spectral indices, are constrained to a level of $\lesssim 1\%$. The propagation parameters are constrained to be about $10\%-20\%$, which are relatively large due to the degeneracy among some of them. For the rigidity-dependence slope of the diffusion coefficient, $\delta$, the statistical error is only a few percent. Compared with previous studies Trotta et al. (2011); Jin et al. (2015); Jóhannesson et al. (2016), our results are widely improved. The one-dimensional (1-d) probability distributions and two-dimensional (2-d) confidence regions of the major propagation parameters are summarized in Figs. 2-7. We also show explicitly the comparison of the data with the fitting results (with $95\%$ credible bands) in Figs. 8-10.

The fittings show that the models with reacceleration (DR, DR2, and DRC) can fit the B/C and proton data well, while the other three non-reacceleration models fit the data relatively poorly888Note, however, the study of CR electrons and positrons results in a different conclusion, i.e., the convection models are more favored than the reacceleration models Lin et al. (2015).. The reduced chi-squared values are all smaller than 1 for the three reacceleration models. For the non-reacceleration models, the $\chi^{2}$ values indicate $p-$values of $\sim 7.8\times 10^{-16}$, $4.3\times 10^{-5}$, and $0.14$ for the PD, DC, and DC2 models, respectively. From Fig. 8 we can see that the predicted B/C ratios for non-reacceleration models do not match the low energy ($E_{k}\lesssim 1$ GeV/n) data well. This is perhaps due to larger solar modulation potentials for non-reacceleration models, which are required by the proton data. These results illustrate the importance of including the low energy ACE data of B/C and the primary CR flux data when studying the propagation of CRs.

There is a clear degeneracy between $D_{0}$ and $z_{h}$. This is because the B/C data can only constrain $D_{0}/z_{h}$ effectively Maurin et al. (2001); Jin et al. (2015). The unstable-to-stable secondary ratio is expected to break such a degeneracy. However, the current 10Be/9Be ratio data are of relatively poor quality. The $95\%$ credible region of $D_{0}$ is $[5.2,9.2]\times 10^{28}$ cm2 s*-1* for the DR model, and the cresponding value of $z_{h}$ is $[3.7,8.2]$ kpc. As a comparison, they are $[5.45,11.20]\times 10^{28}$ cm2 s*-1*, and $[3.2,8.6]$ kpc in Ref. Trotta et al. (2011). Our results improve moderately compared with that of Ref. Trotta et al. (2011). Through analyzing the synchrotron radiation and the electron/positron fluxes, Di Bernardo et al. also found a relatively large propagation halo height ($z_{h}>2$ kpc; Di Bernardo et al. (2013)), which is consistent with our results.

There are some other correlations among the propagation parameters. For example, for the DC and DC2 scenario, an anti-correlation between $D_{0}$ and $dV_{c}/dz$ can be found (Figs. 3 and 4). This can be understood that, a larger convection velocity tends to blow the particles away from the disk, resulting in a lower flux, which can be compensated by a longer propagation time (hence a smaller $D_{0}$). A positive correlation between $\delta$ and $dV_{c}/dz$ can be understood similarly. Since the convection is only important for low energy particles, a larger convection velocity will lead to harder spectra of the CR fluxes and B/C ratio, which can be compensated by a larger value of $\delta$. For the DR2 scenario, we find anti-correlations between $v_{A}$ and $\delta$, $\eta$ and $\delta$, and positive correlation between $v_{A}$ and $\eta$. A larger $v_{A}$ value gives softer spectra of the CR fluxes and B/C ratio, and hence suggesting a smaller $\delta$. The anti-correlation between $\eta$ and $\delta$ can be understood as: a smaller $\eta$ (note that $\eta<0$) gives a larger diffusion coefficient at low energies, and results in harder spectra after the propagation. A larger value of $\delta$ is then able to compensate such an effect.

The slope $\delta$ of the diffusion coefficient is well constrained (with statistical uncertainty of a few percents) given the model setting. However, there are relative large differences among different model configurations. For the reacceleration models, $\delta$ is about $0.38$ for the DR model, and about $0.5$ for the DR2/DRC models. For the DC/DC2 models, $\delta$ is even larger (about $0.6$). These results can be understood via the correlations between $\delta$ and other parameters as described above. The fitting to the B/C ratio above 65 GV gives a slope of $-0.333$ Aguilar et al. (2016a). Our results show that in specific models the value of $\delta$ may differ from that directly inferred from the data. This is because, on one hand, the low energy spectrum of the B/C ratio depends on propagation models, and on the other hand, the uncertainties of high energy data are relatively large. It is currently difficult to distinguish the Kolmogrov ($\delta=1/3$; Kolmogorov (1941)) and the Kraichnan ($\delta=1/2$; Kraichnan (1965)) type of interstellar turbulence. Nevertheless, we find that for some of the propagation model settings, such as the DR2 and DRC models, the Kraichnan type of turbulence is favored. For the DR model, the fitting value of $\delta$ is closer to, but still different from, that predicted by the Kolmogrov theory.

For reacceleration models, the Alfven velocity $v_{A}$ is about 38 km s*-1* for the DR model, which decreases (increases) to about 18 (43) km s*-1* for the DR2 (DRC) model. The major effect of reacceleration is to produce a “GeV bump” of the CR flux and B/C ratio. For the DR2 model, a larger $\delta$ gives higher B/C ratio at lower energies, and hence a smaller reacceleration effect is needed. This can also be seen from the anti-correlation between $v_{A}$ and $\delta$ (Fig. 6). The effect of convection is, however, opposite from that of reacceleration. Therefore for the DRC model, a larger value of $v_{A}$ is favored given a non-zero value of $dV_{c}/dz$.

A break of the injection spectrum around $10-20$ GV is favored in the reacceleration models. Such a break is required to fit the proton fluxes, in order to reduce the “GeV bump” produced by the reacceleration. Such a break is not necessary for the non-reacceleration models. Nevertheless, we find that a spectral hardening with a change of the slope of $\sim 0.10-0.15$ is favored by the fitting. Such a break enables a better fit to the high energy proton flux by the AMS-02 which shows a spectral hardening above $\sim 330$ GV. The break rigidity is not exactly the same as that obtained directly from the data, because the low energy spectral behavior also enters in the fitting.

As for the solar modulation, we find that the time-dependent term of the modulation potential, $\Phi_{1}$, is similar for all models. It reflects the differences of the proton fluxes at different time. The platform term $\Phi_{0}$ differ from each other. In general, non-reacceleration models need remarkably larger $\Phi_{0}$ to accommodate the low energy data of protons.

IV.2 Positrons

The fluxes of secondary positrons can be calculated self-consistently given the fitting propagation and source parameters. Fig. 11 shows the expected $2\sigma$ bands of positron fluxes, compared with the AMS-02 data Aguilar et al. (2014). We find that the reacceleration models which fit the B/C and proton data well would result in a remarkable bump at $\sim$GeV energies and exceed the data significantly. This is consistent with that found in earlier studies Moskalenko et al. (2002); Trotta et al. (2011). For the non-reacceleration models, on the other hand, the expected positron fluxes are lower than the data by a factor of $\sim 2-3$. These results indicate that the production and propagation of positrons may be significantly different from that of the CR nuclei.

For all these models, the predicted positron spectra at high energies ($\gtrsim 10$ GeV) are much softer than that of the data, which indicate the existence of primary positron sources, e.g., pulsars Shen (1970); Harding and Ramaty (1987); Zhang and Cheng (2001).

IV.3 Antiprotons

Fig. 12 shows the results of antiprotons from the models, compared with the PAMELA Adriani et al. (2010) and AMS-02 Aguilar et al. (2016b) measurements. We find that the model predictions are roughly consistent with the data. More detailed comparison shows that in general the non-reacceleration model predictions match the data better than the reacceleration models. For the DR and DRC models, there are slight deficits of low energy ($\lesssim 10$ GeV) antiprotons compared with the data. The DR2 model can marginally fit the data. The prediction of the DC2 model is consistent with the data. For the PD and DC models, however, they slightly under-predict antiprotons around 10 GeV and over-predict lower energy antiprotons. At the high energy end ($E\gtrsim 100$ GeV), there might be excesses of the data (see also Huang et al. (2016); Li (2016); Feng and Zhang (2017); Cholis et al. (2017)). For models with larger $\delta$ values such as the DC, DC2, and DR2 models, the excesses are remarkable. For the other three models with relatively smaller $\delta$ values such excesses are less significant.

V Discussion

V.1 The discrepancy between non-reacceleration models and

CR nuclei data

It seems that the non-reacceleration models have difficulty to fit the proton fluxes and the B/C ratio simultaneously. We find that for the non-reacceleration models the required solar modulation potential ($\Phi_{0}$) is significantly higher than that of reacceleration models, which results in poor fittings to the low energy B/C data of ACE. To test that whether such a discrepancy is due to the difference of solar modulation between protons and heavy nuclei, we do similar fittings using the preliminary Carbon flux by AMS-02 AMS-02 collaboration (2016) instead of the proton fluxes. We find similar conclusion as above, which means that the difference of solar modulation between protons and heavy nuclei is not the major reason of this discrepancy.

Another possible reason is the injection spectrum of CRs. For the three non-reacceleration models the injection spectrum at low rigidities is proportional to $R^{-(2.4-2.5)}$, which is quite soft compared with that of the reacceleration models $R^{-(1.7-2.0)}$. Even though we enable a break of the low energy spectrum of all the models, the fitting results turn out to favor a high energy hardening instead of a low energy break for the non-reacceleration models. We have added another break in the injection spectrum of Eq. (4), and redo the fittings. Still no effective improvement of the fittings is found.

V.2 The discrepancy between all models and the positron data

Our results show that the reacceleration models would over-predict low energy ($\sim$GeV) positrons compared with the measurements, while the non-reacceleration models tend to under-predict positrons. Similar results for reacceleration models have also been obtained in Ref. Di Bernardo et al. (2013). One kind of uncertainty is the hadronic $pp$-interaction. In this work we use the parameterization of positron production in $pp$-interaction of Ref. Kamae et al. (2006). As illustrated in Ref. Delahaye et al. (2009), some other parameterizations would give a positron yield spectrum differing by a factor of $\lesssim 2$ in a certain energy range. However, the uncertainty of the hadronic interaction may not be able to fully solve this discrepancy, especially for the reacceleration models. The other models adopted in Ref. Delahaye et al. (2009) give even more positrons between GeV and TeV, which makes the reacceleration models exceed the data even more. Therefore our results indicate that the propagation of CR nuclei and leptons, either in the Galaxy or in the heliosphere, might be different. Given the very efficient energy losses of leptons, they may experience large fluctuations in the Galaxy Pohl and Esposito (1998). The solar modulation effects may also be different between nuclei and leptons due to their distinct mass-to-charge ratios. The charge-sign dependent solar modulation may take effect either Clem et al. (1996); Della Torre et al. (2012); Maccione (2013); Potgieter et al. (2014); Kappl (2016).

V.3 The Voyager-1 measurements in outer heliosphere

The Voyager-1 spacecraft has traveled by more than 100 astronomical units from the Earth. It has been thought to approach the edge of the heliosphere since a sudden drop of the intensity of low energy ions and an abrupt increase of the CR intensity from outside the heliosphere were observed Stone et al. (2013). The measured CR flux by Voyager-1 can thus be believed to be a direct measurement of the local interstellar CRs. The Voyager-1 data would be helpful in better constraining the source injection parameters as well as the solar modulation parameters. However, as shown in Ref. Cummings et al. (2016), the very low energy ($\lesssim 50$ MeV/n) B/C spectrum measured by Voyager-1 is difficult to be modelled in various models. Further tuning of the modelling and/or better understanding about the measurements may be necessary. The Voyager-1 data will be included in future studies.

V.4 Reacceleration models and antiprotons

The reacceleration models would generally under-estimate the low energy antiproton fluxes. Several kinds of scenarios were proposed to explain this. In Ref. Moskalenko et al. (2003) it was proposed that a local and fresh source, probably associated with the Local Bubble, might produce additional low energy primaries and hence decrease the measured secondary-to-primary nuclei ratio. The annihilation of several tens of GeV dark matter particles may also be responsible for the low energy excess of antiprotons Hooper et al. (2015); Cui et al. (2016); Cuoco et al. (2016). Alternatively, an empirical adjustement of the velocity-dependence of the diffusion coefficient with a $\beta^{\eta}$ term, i.e., the DR2 model in this work, was suggested to be able to explain the B/C and antiproton data di Bernardo et al. (2010). In this treatment a larger $\delta$ value and a weaker reacceleration effect is required, which enables more production of low energy secondary particles (both Boron and antiprotons). As shown in Fig. 12, the DR2 model does improve the fitting. However, the physical motivation for such a term is not well justified. Finally, the uncertainties of the production cross section of antiprotons make this problem still inconclusive Donato et al. (2001); Giesen et al. (2015); Lin et al. (2016).

VI Conclusion

In this work we adopt the precise measurements of the B/C ratio and the time-dependent proton fluxes by AMS-02 and PAMELA to constrain the injection and propagation parameters of Galactic CRs. We employ a self-consistent treatment of the solar modulation by means of a linear correlation of the modulation potentials with solar activities. We have carried out a comprehensive study of a series of CR propagation models, including the PD, DR, DC, DRC, and two variants of the DR and DC models. The predictions of secondary positrons and antiprotons based on the fitting parameters are calculated and compared with the data.

We summarize the comparison of various models with different data sets in Table 3. It is shown that no model can match all these data simultaneously, which suggests that the actual case for the origin, propagation, and interaction of CRs is more complicated than our current understanding. For the CR nuclei only, we find that the DR2 model may give the best match to all the data. However, the phenomenological modification of the diffusion coefficient (the $\beta^{\eta}$ term) may need to be understood further di Bernardo et al. (2010).

We list our main conclusion as follows.

• •

The reacceleration models (DR, DR2, and DRC) can fit both the B/C and proton fluxes well, while non-reacceleration models (PD, DC, and DC2) can not. The failure of non-reacceleration models can not be simply ascribed to the differences of solar modulation or the source injection spectra between protons and heavier nuclei.

• •

The statistical uncertainties of the propagation parameters are constrained to a level of $10\%-20\%$, thanks to the precise measurements of CR data by AMS-02. However, there are relatively large differences (up to a factor of $\sim 2$) among different model settings.

• •

For reacceleration models, the value of $\delta$ is found to be about $0.38-0.50$, which slightly favor the Kraichnan type of interstellar turbulence.

• •

The reacceleration models will over-produce positrons but under-produce (except DR2) antiprotons in general. The non-reacceleration models, on the other hand, predict fewer positrons and (marginally) consistent antiprotons when compared with the measurements.

• •

Our results suggest that there are significant differences of the propagation in either the Milky Way or the heliosphere between nuclei and leptons.

With more and more precise data available, we are able to investigate the CR-related problems in great detail. It turns out that the problem seems to be more complicated than what we expected based on the rough measurements in the past. The final understanding of the propagation of CRs may need not only the CR data themselves but also the full improvements of the understanding of the astrophysical ingradients of the Milky Way, as well as the nuclear and hadronic interactions.

Acknowledgements.

We thank the ACE CRIS instrument team and the ACE Science Center for providing the ACE data. This work is supported by the National Key Research and Development Program of China (No. 2016YFA0400200), the National Natural Science Foundation of China (No. 11475191), and the 100 Talents program of Chinese Academy of Sciences.