The Curious Case of High-energy Deuterons in Galactic Cosmic Rays

Deuteron isotopes ${}^{2}$ H are rare particles in Galactic cosmic rays (CRs). They are destroyed rather than formed in thermonuclear reactions in stellar interior so that, from supernova remnants (SNRs) as sources of CRs, no significant amount of deuteron is expected to be released by diffusive shock acceleration (DSA). An important process by which high-energy deuterons can be created in the Galaxy is nuclear fragmentation of CR nuclei with the gas of the interstellar medium (ISM). The main source of deuteron production is fragmentation of ${}^{4}$ He and ${}^{3}$ He isotopes, along with C–${\mathsf{N}}$–O nuclei, in collisions with interstellar hydrogen and helium. An important contribution comes from the reaction ${\mathsf{p}}+{\mathsf{p}}\to {\mathsf{d}}+\pi$. From nuclear fragmentation, interactions of CRs with the ISM are also known to generate Li–Be–B nuclei and antiparticles, that are otherwise rare from stellar nucleosynthesis processes. The measured abundances of these elements in the cosmic radiation enables us to pose tight constraints on the astrophysical models of CR propagation in the Galaxy and, in the context of indirect searches of dark matter, to asses the astrophysical background of CRs antiparticles (Grenier et al. 2015).

Very recently, the Alpha Magnetic Spectrometer (AMS-02) experiment has reported a new precise measurement of the B/C ratio in CRs at energy between ∼ 0.5 and 1000 GeV/n (Aguilar et al. 2016). At kinetic energies above ∼10 GeV/n, the ratio is found to decrease steadily with increasing energy. Since boron nuclei are mainly produced by C–${\mathsf{N}}$–O fragmentation in the ISM, the observed trend of the B/C ratio reflects the conception that CRs diffusively propagate in the Galactic magnetic fields with an average diffusion coefficient that increases with energy. According to this picture, the ratio between ${}^{2}$ H and its main progenitor ${}^{4}$ He must follow the same behavior in the GeV–TeV energy range. Recent calculations have shown, indeed, that the observed boron and deuteron abundances at GeV/n energies can both be self-consistently described by diffusion models of CR propagation (Coste et al. 2012; Tomassetti 2012a). In the GeV–TeV energy region, the d/He ratio is expected to decrease rapidly, as fast as the B/C ratio does, but CR deuterons have never been detected at these energies.

Quite unexpectedly, a new analysis of the data collected by the satellite mission SOKOL has determined the d/He ratio at 0.5–2 TeV/n energy (Turundaevskiy & Podorozhnyi 2016). In this measurement, the discrimination between ^1,2H and ^3,4He isotopes was performed by means of neural network analyses of the topology of hadronic showers developed by these particles. By means of two different Monte Carlo simulations, consistent results have been obtained: 0.114 ± 0.023 and 0.099 ± 0.021 for the ${}^{2}$ H/(${}^{1}$ H+${}^{2}$ H) ratio, 1.64 ± 0.30 and 1.43 ± 0.27 for the ${}^{2}$ H/${}^{4}$ He ratio.

In this Letter, we show that the above results represent a striking anomaly in CR physics that cannot be explained by standard models of CR propagation and secondary production in the ISM. Then, we argue that a deuteron excess can be explained in terms of hadronic production occurring inside SNRs, which was proposed in Blasi (2009) to explain the positron excess in CRs (see also Kachelrießet al. 2011; Serpico 2012). Using new evaluations of fragmentation cross-sections, we calculate for the first time the high-energy production of CR deuterons in SNRs and in the ISM, demonstrating that this mechanism can account for the new d/He data. Along with the d/He ratio, we also compute the B/C ratio under the same framework, showing that there are conflicting results in the model predictions for the two observables. We therefore conclude that this tension can be resolved if the deuterons progenitors are not accelerated in the same sources of boron progenitors. We discuss our results and their implications for the interpretation of antiproton data.

We compute the spectrum of CR nuclei accelerated in SNRs within the linear DSA theory and including the production of secondary fragments. Similar calculations are done in earlier works (Blasi 2009; Blasi & Serpico 2009; Mertsch & Sarkar 2009, 2014; Tomassetti & Donato 2015; Herms et al. 2016). We follow closely the derivation of Tomassetti & Donato (2012). In the shock rest-frame (x = 0), the upstream plasma flows in from $x\lt 0$ with speed u₁ (density n₁) and the downstream plasma flows out to $x\gt 0$ with speed u₂ (density n₂). The compression ratio is $r={u}_{1}/{u}_{2}={n}_{2}/{n}_{1}$. For a nucleus with charge Z and mass number A, the equation describing diffusion and convection at the shock reads

$$\begin{eqnarray}&&u\displaystyle \frac{\partial f}{\partial x}=D\displaystyle \frac{{\partial }^{2}f}{\partial {x}^{2}}+\displaystyle \frac{1}{3}\displaystyle \frac{{du}}{{dx}}p\displaystyle \frac{\partial f}{\partial p}-{{\rm{\Gamma }}}^{\mathrm{tot}}\ f+Q,\end{eqnarray} \tag{ 1 }$$

where f is the phase space density, D(p) is the diffusion coefficient near the shock, u is the fluid speed, and ${{\rm{\Gamma }}}^{\mathrm{tot}}={cn}{\sigma }^{\mathrm{tot}}$ is the total fragmentation rate for cross-sections ${\sigma }^{\mathrm{tot}}$ and background density n, which is assumed to be composed of H and He like the average ISM. The source term includes particle injection at the shock, $Q=Y\delta (x)\delta (p-{p}^{\mathrm{inj}})$, with ${p}^{\mathrm{inj}}\equiv {{ZR}}^{\mathrm{inj}}$ and ${R}^{\mathrm{inj}}\cong 0.5\,$ GV for all nuclei. The Y-constants set the normalization of each species. We assume strong shocks ($r\approx 4$) and a diffusion coefficient $D={\kappa }_{B}\tfrac{p/Z}{3B}$, where ${\kappa }_{B}$ parameterizes the deviation of D(p) from the Bohm value due to magnetic damping. The resulting acceleration rate at momentum p is ${{\rm{\Gamma }}}^{\mathrm{acc}}\sim {u}_{1}^{2}/20\,D$. For an SNR of age ${\tau }_{\mathrm{snr}}$, the condition ${{\rm{\Gamma }}}^{\mathrm{acc}}={\tau }_{\mathrm{snr}}^{-1}$ defines the maximum momentum scale attainable by DSA. In the presence of hadronic interactions, the additional requirement ${{\rm{\Gamma }}}^{\mathrm{tot}}\ll {{\rm{\Gamma }}}^{\mathrm{acc}}$ must be fulfilled. The downstream solution reads

$$\begin{eqnarray}&&{f}_{2}(x,p)={f}_{0}(p)\left(1-\displaystyle \frac{{{\rm{\Gamma }}}_{2}^{\mathrm{tot}}}{{u}_{2}}x\right)+\displaystyle \frac{{q}_{2}}{{u}_{2}}x,\end{eqnarray} \tag{ 2 }$$

where ${f}_{0}(p)$ is the distribution function at the shock. As found in Mertsch & Sarkar (2009), ${f}_{0}(p)$ is given by

$$\begin{eqnarray}{f}_{0}(p) & = & \alpha {\displaystyle \int }_{0}^{p}{\left(\displaystyle \frac{p^{\prime} }{p}\right)}^{\alpha }Y\delta (p^{\prime} -{p}^{\mathrm{inj}}){e}^{-\chi (p,p^{\prime} )}\displaystyle \frac{{dp}^{\prime} }{p^{\prime} }\\ & & +\alpha {\displaystyle \int }_{0}^{p}{\left(\displaystyle \frac{p^{\prime} }{p}\right)}^{\alpha }\displaystyle \frac{{q}_{1}D}{{u}_{1}^{2}}(1+{r}^{2}){e}^{-\chi (p,p^{\prime} )}\displaystyle \frac{{dp}^{\prime} }{p^{\prime} },\end{eqnarray} \tag{ 3 }$$

with $\alpha =3r/(r-1)$, $\chi \approx \alpha ({{\rm{\Gamma }}}_{1}^{\mathrm{tot}}/{{\rm{\Gamma }}}^{\mathrm{acc}})[D(p)-D(p^{\prime} )]$, and the subscript i = 1 (i = 2) indicates the upstream (downstream) region. The first term of Equation (3) describes primary particles injected at the shock and it is of the form $\sim {p}^{-\alpha }{e}^{-\chi }$. The second term of Equation (3) describes the production and acceleration of CR fragments from heavier progenitors. For each $k\to j$ process, the q-term of Equation (3) is given by ${q}_{{kj}}^{\sec }(p)={\xi }_{{kj}}^{-3}{f}_{j}(x,p/{\xi }_{{kj}}){{\rm{\Gamma }}}_{{kj}}^{\mathrm{fr}}$, where ${{\rm{\Gamma }}}_{{kj}}^{\mathrm{fr}}={\sigma }_{{kj}}n\beta c$ is the secondary production rate, and ${\sigma }_{{kj}}$ is the corresponding cross-section. The momentum inelasticity factor ${\xi }_{{kj}}={A}_{j}/{A}_{k}$ expresses the conservation of kinetic energy per nucleon between progenitor and fragment.

Equation (1) is solved for all relevant species. We considered primary nuclei (with $Y\gt 0$) ${\mathsf{p}}$, ${}^{4}$ He, C, ${\mathsf{N}}$, O tuned to recent data as in Mertsch & Sarkar (2014), and secondary B production from C–${\mathsf{N}}$–O collisions with hydrogen and helium gas. The adopted fragmentation cross-sections are those re-evaluated in Tomassetti (2015b). We account for deuteron production from collisions of ${\mathsf{p}}$, ${}^{3}$ He, and ${}^{4}$ He off hydrogen and helium. Measurements and calculations for these reactions are available only below 10 GeV/n of energy. Thus, we have performed new calculations at $E=10\mbox{--}1000$ GeV/n using the hadronic Monte Carlo generator QGSJET-II-04 (Ostapchenko 2011). The cross-sections for the dominant channels are shown in Figure 1. The solid lines are our parameterizations adapted from Tomassetti (2012a) and extended to higher energy. Cross-sections for collisions off the He target are larger by nearly a factor of two and contribute to ∼20 %. The contribution from fragmentation of heavier nuclei is negligible. In all reactions, kinetic energy per nucleon is approximately conserved except for the ${\mathsf{p}}$+${\mathsf{p}}$ reaction, which peaks at a proton energy of ∼600 MeV/n and produces deuterons within a narrow energy range around 150 MeV/n (Meyer 1972). Hence, these deuterons experience DSA acceleration similarly to primary components with non-zero Y-factor, i.e., to a power-law spectrum $\sim {p}^{-\alpha }$. In contrast, deuterons ejected by helium fragmentation have an initial spectrum ${p}^{-\alpha }$ so that, for those produced within a distance $d\lesssim \,D/{u}_{1}$ from the shock and undergo DSA, their final spectrum reaches the form ${f}_{0}^{\sec }\sim \,{f}_{0}^{\mathrm{pri}}D(p)\propto {p}^{-\alpha +1}$ (see Equation (3)). This means that, for Bohm-type diffusion, the secondary deuteron spectrum is one power harder than that of its progenitors. The total CR flux produced by SNRs is computed as

$$\begin{eqnarray}&&{{ \mathcal S }}^{\mathrm{snr}}(p)=4\pi {p}^{2}{{ \mathcal R }}_{\mathrm{SN}}{\int }_{0}^{{\tau }_{\mathrm{snr}}{u}_{2}}4\pi {x}^{2}{f}_{2}(x,p){dx}\end{eqnarray} \tag{ 4 }$$

where ${{ \mathcal R }}_{\mathrm{SN}}\cong 25$ Myr⁻¹ kpc⁻² is the explosion rate per unit volume and ${\tau }_{\mathrm{snr}}\cong 40$ kyr is the age of the SNR.

Zoom In Zoom Out Reset image size

To model the subsequent propagation of CRs in the ISM, we adopt a two-halo model of CR diffusion and nuclear interactions (Tomassetti 2012b, 2015a). The Galaxy is modeled as a disk of half-thickness $h\cong 100$ pc containing SNRs and gas with number density $\tilde{n}\cong 1$ cm⁻³. The disk is surrounded by a diffusive halo of half-thickness L and zero matter density. We give a one-dimensional description in the thin disk limit ($h\ll L$). For each CR nucleus, the transport equation reads

$$\begin{eqnarray}&&\displaystyle \frac{\partial { \mathcal N }}{\partial t}=\displaystyle \frac{\partial }{\partial z}\left(K(z)\displaystyle \frac{\partial { \mathcal N }}{\partial z}\right)-2h\hat{\delta }(z){\tilde{{\rm{\Gamma }}}}^{\mathrm{tot}}{ \mathcal N }+2h\hat{\delta }(z){{ \mathcal S }}^{\mathrm{tot}},\end{eqnarray} \tag{ 5 }$$

where ${ \mathcal N }$ is its density, $\hat{\delta }$ is a Dirac function of the z-coordinate, K(z) is the diffusion coefficient of CRs in the Galaxy, and ${\tilde{{\rm{\Gamma }}}}^{\mathrm{tot}}=\beta c\tilde{n}{\sigma }^{\mathrm{tot}}$ is the destruction rate in the ISM at velocity $\beta c$ and cross-section ${\sigma }^{\mathrm{tot}}$. The source term ${{ \mathcal S }}^{\mathrm{tot}}$ is split into a primary term ${{ \mathcal S }}^{\mathrm{snr}}$, obtained from Equation (4) as solution of the DSA equation, and a secondary production term ${{ \mathcal S }}^{\sec }={\sum }_{k}{\tilde{{\rm{\Gamma }}}}_{k}^{\mathrm{fr}}{{ \mathcal N }}_{k}$, from fragmentation of k-type nuclei in the ISM with rate ${\tilde{{\rm{\Gamma }}}}_{k}^{\mathrm{fr}}$. To compute the interaction rates in the ISM ${\tilde{{\rm{\Gamma }}}}^{\mathrm{in}/\mathrm{fr}}$, we adopt the same fragmentation network (and same cross-sections) as occurring inside SNRs. Equation (5) is solved in steady-state conditions $\partial { \mathcal N }/\partial t=0$. The derivation of the full solution is in Tomassetti (2012b). The diffusion coefficient is taken of the type $K{(p,z)=\beta {K}_{0}(({pc}/{Ze})/{GV})}^{\delta (z)}$ where K₀ expresses its normalization. For the scaling index, $\delta (z)$, we adopt $\delta ={\delta }_{0}$ in the region of $| z| \lt \xi \,L$ (inner halo) and $\delta ={\delta }_{0}+{\rm{\Delta }}$ for $| z| \gt \xi \,L$ (outer halo). Our default parameters are set as $\xi \,L\cong \,0.1$ ${\delta }_{0}\cong 1/3$, ${\rm{\Delta }}\cong 0.55$, and ${K}_{0}/L\cong 0.01$ kpc Myr⁻¹. The differential energy fluxes of each species are given by $J(E)=\tfrac{\beta c}{4\pi }{ \mathcal N }$. Solar modulation is described in force-field approximation using the parameter $\phi =500$ MV for a medium-level modulation strength. The d/He and B/C ratios as function of kinetic energy per nucleon are eventually calculated as ${J}_{{}^{2}{\rm{H}}}/{J}_{{}^{4}\mathrm{He}}$ and $({J}_{{}^{10}{\rm{B}}}+{J}_{{}^{11}{\rm{B}}})/({J}_{{}^{12}{\rm{C}}}+{J}_{{}^{13}{\rm{C}}})$, respectively. In the following, we consider two model implementations representing two alternative scenarios:

• 1.  
  Scenario #1 (B/C-driven, conservative)—standard model without interactions in sources, which is the case of a CR flux released by young SNRs with amplified magnetic fields ($B\gtrsim 100$ $\mu {\rm{G}}$) and/or low background density (${n}_{1}\approx \,{10}^{-3}$ cm⁻³). In this model, secondary production of CRs deuterons or Li–Be–B occurs only in the ISM. This model is tuned to match the new B/C data from AMS-02 at GeV/n–TeV/n energies.
• 2.  
  Scenario #2 (d/He-driven, speculative)—model with copious production and acceleration of secondary particles in an SNR shockwave, which is the case for a GeV–TeV flux provided by old SNRs with damped magnetic fields, slow shock speed, or dense ambient medium, i.e., with the combination ${n}_{1}{\kappa }_{B}{B}^{-1}{u}_{8}^{-2}\sim 400$. This model is tuned against the d/He data including the new SOKOL measurement at TeV/n energies.

Model calculations are shown in Figure 2 for the B/C ratio the d/He ratio at energies between ∼0.5 GeV and 2 TeV per nucleon. Scenario #1 is plotted as green solid lines. In this model, secondary nuclei such as ${\mathsf{d}}$ or B are entirely generated in the ISM, i.e., without SNR components, and thus secondary/primary ratios decrease steadily as ${J}_{s}/{J}_{p}\propto (L/{K}_{0})[\xi +(1-\xi ){\rho }^{-{\rm{\Delta }}}]{\rho }^{-{\delta }_{0}}$ where $\rho =({pc}/{Ze})/{GV}$. In the high-energy limit one has ${J}_{s}/{J}_{p}\propto {E}^{-1/3}$. It can be seen that this model fits remarkably well the new AMS-02 data on the B/C ratio and dictates a similar trend for the d/He ratio, which is predicted to reach the level of $\sim {10}^{-2}$ in the TeV/n energy scale. We therefore conclude that the SOKOL measurement of the d/He ratio is at least two orders of magnitude higher than that expected from standard models where CR deuterons are produced by fragmentation in the ISM.

Zoom In Zoom Out Reset image size

In Scenario #2, shown as red dashed lines, hadronic interaction processes inside SNRs generate a source component of secondary nuclei, which is harder than that arising from CR collisions with the ISM and, as discussed, even harder than that of primary ${\mathsf{p}}$–He spectra. It is then possible, with fragmentation inside SNRs, having secondary/primary ratios that increase with energy. Figure 2 shows that Scenario #2 matches fairly well the d/He ratio measurements at GeV/n and at TeV/n energies, therefore providing an explanation for the new SOKOL data. Under this model, however, the B/C ratio is also predicted to increase, at energies above ∼50 GeV/n, in remarkable contrast with the new AMS-02 data. While interactions inside SNRs seem to be the only mechanism capable of explaining a rise in the d/He ratio, it is apparent that the observed decrease of the B/C ratio conflicts with this mechanism. We also note that, at the ∼1 GeV/n energy region where secondary CR production in the ISM dominates, the two ratios are consistently described by both models #1 and #2, at least within the precision of the current data.

As we see it, the only solution to this tension is a situation where the connection between d/He ratio and B/C ratio is broken. This situation is realized if the sources accelerating helium and protons (and producing deuterons) are not the same as those accelerating heavier C–${\mathsf{N}}$–O nuclei and, in particular, deuterons must be accelerated by a low-metallicity source, which may be the case for SNRs expanding over H-dominated molecular clouds. Such a possibility was also discussed in Cholis & Hooper (2014, see Section VI), and proposed in other works (Fujita et al. 2009; Kohri et al. 2016), all focused on antiparticle excesses in CRs. In such a scenario, the B/C ratio can no longer be used to place constraints on antimatter spectra. In this respect, it is important to note that the connection between d/He ratio and antimatter/matter ratios would still be preserved because, in contrast to Li–Be–B nuclei, secondary deuterons share their progenitors with positrons and antiprotons. It is then interesting to calculate the antiproton/proton ratio in light of new data released by AMS-02 (Aguilar et al. 2016). Calculations were performed in Blasi & Serpico (2009) and subsequent works (e.g., Mertsch & Sarkar 2014). In this work, the antiproton distribution at the shock ${f}_{0}(p)$ is calculated numerically, as done in Herms et al. (2016), in order to drop the "inelasticity approximation" that links the antiproton momentum p to the primary proton momentum p_p through an assumedly constant factor $\xi \equiv \langle p/{p}_{p}\rangle$. It was noted that such an approximation leads to an overestimate of the high-energy antiproton production in SNRs (Kachelrießet al. 2011). Along with antiproton production from ${\mathsf{p}}$–He collisions with hydrogen and helium gas, we also account for tertiary reactions (such as $\bar{p}+p\to {\bar{p}}^{\prime }+X$) and for destruction processes. All these processes are implemented in both acceleration and propagation. The corresponding cross-sections are taken from Feng et al. (2016). The resulting $\bar{p}$/$p$ ratio is shown in Figure 3. It can be seen that Scenario #2 is preferred by the AMS-02 data. Reducing nuclear uncertainties in antiproton production is clearly essential for a complete discrimination between the different models (Feng et al. 2016).

Zoom In Zoom Out Reset image size

Given the implications of these new data on the phenomenology of CR propagation, we believe that the situation deserves more clarification on the experimental side. The SOKOL analysis relies on unconventional techniques of deuteron/proton mass separation, which is always a very challenging task. For instance the d/He measurement might be overestimated due to undetected background arising, e.g., from the mass distribution tails of CR protons or from ${}^{4}$ He nuclei fragmenting in the top of the instrument. On the other hand, it is very unlikely for such a background to affect the results by two orders of magnitude. The deuteron spectrum and the d/He ratio are being precisely measured by the AMS-02 experiment at $E\sim 0.1\mbox{--}10$ GeV/n with standard spectrometric techniques. In addition, AMS-02 is also equipped with a Transition Radiation Detector (TRD), designed for lepton/hadron mass separation, which can provide direct measurements of the Lorentz factor γ at TeV/n energies (Obermeier & Korsmeyer 2015). In standard magnetic spectrometers, the CR mass is derived from velocity and momentum measurements, $M\propto \,p/(\gamma \beta c)=\tfrac{p}{\beta c}\sqrt{1-{\beta }^{2}}$, so that its corresponding resolution $\delta M/M$ is given by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta M}{M}\right)}^{2}={\left(\displaystyle \frac{\delta p}{p}\right)}^{2}+{\gamma }^{4}{\left(\displaystyle \frac{\delta \beta }{\beta }\right)}^{2},\end{eqnarray} \tag{ 6 }$$

showing that the mass resolution degradates rapidly, at relativistic energies, due to the ${\gamma }^{4}$-factor. In contrast, with the opportunity of performing direct TRD-based γ-measurements, AMS-02 may have the capability to detect CR deuterons at the ${ \mathcal O }(\mathrm{TeV})$ energy scale.

This work is aimed at interpreting new data, registered by the SOKOL experiment in space, that revealed a surprisingly high abundance of CR deuterons in the TeV/n region. In contrast to antiparticle excesses that can be explained, e.g., by pulsar models or dark matter annihilation (Serpico 2012; Feng & Zhang 2016), the SOKOL data demand an enhanced high-energy production of CR deuterons from hadronic interactions. We found that no explanation for this measurement can be provided in terms of standard collisions of CRs with the gas of the ISM. As we have shown, a viable solution for this puzzle is the occurrence of nuclear fragmentation inside SNRs, but this mechanism conflicts with the new AMS-02 data on the B/C ratio.

Thus, if the SOKOL measurement is taken as face value, we conclude that the sources accelerating helium and protons (thereby producing deuteron) may not be the same as those accelerating C–${\mathsf{N}}$–O nuclei (otherwise producing Li–Be–B nuclei), and that the former are more efficient in the production and acceleration of secondary particles. Under such a scenario, the connection between B/C ratio and antiparticle/particle ratios would also be broken, and thus the B/C ratio should not be used to place constraints on the astrophysical antimatter background. On the other hand, the d/He ratio would still represent a direct diagnostic tool for assessing this background.

We are grateful to Tanguy Pierog, Colin Baus, and Ralf Ulrich for the CRMC interface to MC generators. N.T. acknowledges support from the H2020-MSCA-IF-2015 action, grant No. 707543-MAtISSE. J.F. acknowledges support from China Scholarship Council and the Taiwanese Ministry of Science and Technology under grant Nos. 104-2112-M-001-024 and 105-2112-M-001-003.