Nuclear de-excitation lines as a probe of low-energy cosmic rays

The energy density of Galactic cosmic rays (CRs), $w_{\rm CR}\approx 1\leavevmode\nobreak\ \rm eV/cm^{3}$, is comparable to the energy density contributed by the interstellar magnetic fields and thermal gas. CRs play an essential role in the process of star formation through the heating and ionization, initiating several crucial chemical reactions in the dense cores of molecular clouds (Papadopoulos, 2010).

The flux and spectra of CRs inside the Solar System have been measured with unprecedented accuracy by space detectors such as PAMELA¹¹1A Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (Adriani et al., 2013, 2014), AMS-02²²2The Alpha Magnetic Spectrometer (Aguilar et al., 2015a, b), CREAM³³3The Cosmic Ray Energetics and Mass experiment (Yoon et al., 2017), DAMPE⁴⁴4The Dark Matter Particle Explorer (An et al., 2019), and CALET⁵⁵5The CALorimetric Electron Telescope (Maestro et al., 2020). Most of these direct measurements from Earth focus on the CR spectra above a few GeV/nuc, below which the fluxes are strongly affected by solar modulation. A few years ago, the Voyager 1 satellite passed through the termination shock and measured the low-energy CR (LECR) spectra (from several MeV/nuc up to hundreds of MeV/nuc) beyond the Solar System (Cummings et al., 2016). The measurements outside the heliosphere are thought to be free of solar modulation and thus may represent the LECR spectra in the local interstellar medium. However, this information cannot be extrapolated to other parts of the Galaxy, in which the CRs are not homogeneously distributed (Aharonian et al., 2018; Jóhannesson et al., 2018; Baghmanyan et al., 2020). Furthermore, the flux of LECRs can vary dramatically on a smaller scale because of higher energy losses and propagation effects, and the flux is usually much higher around possible CR acceleration sites. This is indirectly supported by studies of the interstellar ionization rates (in particular, see, e.g., Indriolo et al., 2009; Indriolo & McCall, 2012; Indriolo et al., 2015). These researches also argued for an LECR component in addition to the standard contribution by supernova remnants (SNRs), which is also supported by the observation of primary Be (Tatischeff & Kiener, 2011).

In this regard, $\gamma$-rays produced in interactions of CRs with the ambient gas can be used as a unique tool to study the spectral and spatial distributions of CRs throughout the Milky Way. For CRs whose kinetic energy exceeds the kinematic threshold of $\pi$-meson production, $E_{\rm th}\simeq 280$ MeV/nuc, the best $\gamma$-ray energy band for exploration is 0.1 - 100 GeV because of the copious production of $\pi^{0}$-decay $\gamma$-rays and the potential of the currently most sensitive detector, Fermi LAT (Atwood et al., 2009). At energies below this kinematic threshold, the nuclear de-excitation lines provide the most straightforward information about the LECR protons and nuclei (e.g., Ramaty et al., 1979; Murphy et al., 2009).

In this paper, we treat the production of nuclear de-excitation lines by combining the recent advances in the modeling of nuclear reactions with the propagation and energy losses of subrelativistic and transrelativistic CRs. In Sec. 2 we calculate the distributions of LECRs near the sources, taking the processes of diffusion and energy losses into account. In Sec. 3 we describe the method for calculating the emissivity of de-excitation $\gamma$-ray lines using the code TALYS (Koning et al., 2008). We then present the results and discuss their implications. Finally, we summarize the main results in Sec. 4.

The energy losses of CR protons with kinetic energies above $\sim$1 GeV are dominated by inelastic nuclear reactions (see Fig.1), and the corresponding cooling weakly depends on energy. As a result, the propagation has very little effect on the initially injected CR spectral shape. At lower energies, the losses are dominated by the ionization and heating of the ambient medium. In this regime, as shown in Fig.1, the cooling time is energy dependent, thus the propagation can significantly distort the initial spectrum. On the other hand, because of the energy-dependent particle diffusion, LECRs propagate more slowly than relativistic CRs. Consequently, the spatial distribution of LECRs should be much more inhomogeneous than that of high-energy CRs. LECRs are expected to be concentrated in the vicinity of their accelerators with energy distributions that strongly depend on the distance to their production sites. The propagation of LECRs in the proximity of their sources should therefore be treated with great care and depth.

As introduced in Padovani et al. (2009), $L=-dE/dN_{\rm H_{2}}$, represents the energy loss-rate per column density $dN_{\rm H_{2}}$, then the ionization cooling time for protons colliding with a medium atomic hydrogen density $n$ can be estimated by $\tau\sim\frac{E}{L\,n\,\beta\,c}$, where $\beta=v/c$. For protons with kinetic energy 10 MeV that collide with ${\rm H_{2}}$, $L\simeq 2\times 10^{-16}\,\rm eV\,cm^{2}$ (see Fig.7 in Padovani et al. 2009). The cooling time of 10 MeV protons in a dense environment with $n\geq 100\,\rm cm^{-3}$ therefore is about several thousand years, much shorter than the duration of typical CR accelerators such as SNRs. Therefore we adopt the steady-state solution for continuous injection to estimate the spectra of LECRs. To be more specific, we assume a continuous injection of protons from a stationary point source. Then the steady-state energy and radial distribution distribution of the LECR protons is obtained by application of the analytical solution derived by Atoyan et al. (1995),

Here $Q(E)$ represents the CR injection rate, $P(E)$ is the energy loss-rate of CR protons summarized in Padovani et al. (2009), $r$ represents the radial distance to the source, and $\Delta u(E,x)=\int^{x}_{E}\frac{D(E^{\prime})}{P(E^{\prime})}dE^{\prime}$, where $D(E)$ is the energy-dependent diffusion coefficient of CRs. In general, this is a quite uncertain parameter that depends on the level of turbulence in the environment. In the Galaxy, it is derived from observations of secondary CRs. We adopted the diffusion coefficient derived for Galactic CRs at energies above 1 GeV (Strong et al., 2007) and assumed that it can be extrapolated to low energies down to 1 MeV: $D(E)=4\times 10^{28}\chi(E/1\rm\ GeV)^{0.5}\ {\rm cm^{2}\,s^{-1}}$. Because the turbulence level near the CR sources is high, the diffusion could be much slower (Malkov et al., 2011; D’Angelo et al., 2018). Effects like this have previously been observed in the $\gamma$-ray band. The electron halos near pulsars reveal a small diffusion coefficient (Abeysekara et al., 2017; Di Mauro et al., 2019). To take this effect into account, we considered a broad range of the parameter $\chi=1,0.1,0.01$. Here the total proton injection rate $Q$, integrated from 1 MeV to 100 MeV, and the distance of the source $d$, are fixed by the value of parameter $Q/d^{2}=10^{38}\,\rm erg\,s^{-1}\,kpc^{-2}$.

The injection spectrum of CRs depends on the acceleration mechanism and the conditions inside the accelerator such as the turbulence level and the magnetic field. Typically (although not always), in the case of diffusive shock acceleration, for instance, the distribution of accelerated particles can be presented as a power law in momentum $p$, $Q(p)=Q_{0}p^{-s}$ (see, e.g., Amato, 2014). In this paper, it is more convenient to write the distribution in terms of kinetic energy $E$, $Q(E)=Q_{0}p^{-s}/\beta$, where $Q_{0}$ is the normalization parameter derived from the power of the CR source. By solving Eq.1, we obtain the steady-state distribution of protons at different radial distances $r$ to the source.

The energy spectra of CR protons diffusing from the site of their acceleration are shown in Fig.2. In the left panel of Fig.2, the curves are calculated for four injection spectra of protons with power-law index $s=$1.0, 2.0, 3.0, and 4.0, in which the radial distance $r=10\,\rm pc$, the medium density $n=1\,\rm cm^{-3}$ , and the parameter $\chi=1$. In the right panel of Fig.2, we present the CR spectra for the injection spectrum with $s=2.0$ at various distances, in which $r=$1 pc, 10 pc, and 100 pc. By changing the diffusion coefficient, that is, the parameter $\chi=$1, 0.1, and 0.01, we find that for a small diffusion coefficient, the energy spectra of protons become very hard, especially at large distances from the source.

In Fig. 3 we show the radial dependence of proton fluxes. In the left panel, we show the effect of the diffusion coefficient on the flux of 10 MeV protons. At small distances from the source, their proton flux is significantly higher for slow diffusion, but at large distances, the flux drops as the cooling becomes an important factor. The contrast of the spectra under different assumptions of $n$ (the red and black lines in the left panel of Fig.3) also shows the effect of the ambient gas density on the CR spectrum. As expected, because of enhanced energy losses, the radial distribution of the CR fluxes become sharper than $1/r$ at larger distances, $r\geq 10$ pc. In the right panel, we show the $r$-profiles for different CR energies calculated for the parameter $\chi=1$. For the radial profile of the CR, the flux at 1 GeV is close to $1/r$. This is consistent with predictions for the continuous injection in Atoyan et al. (1995), as long as the energy losses are negligible. At low energies, $E\leq 1\,\rm GeV$ because of the high energy losses and slow propagation, we see a stronger radial dependence of the CR flux. Namely, at small distances, the flux follows the $1/r$ profile, but at larger distances the cooling starts to play significant role and causes the radial distribution to become as steep as $1/r^{3}$. For 1 MeV protons, the transition between $1/r$ and $1/r^{3}$ takes place at just several pc. Last but not least, we note that the hardening of the spectrum of CR protons at low energies in the local interstellar medium, as shown by the gray line in Fig. 2, fitted by Phan et al. (2018) using data from Voyager 1 and AMS-02 detectors), can be naturally explained by the propagation and energy losses.

In the derivation of Eq.(1), only the diffusion and energy loss of protons are taken into account. Meanwhile, the advection may also play a non-negligible role in the vicinity of CR sources. In this case, an advection term $V\frac{\partial F_{\rm p}(r,E)}{\partial r}$ should be added to the propagation equation, where $V$ is the advection velocity. By dimensional analysis, the advection dominates diffusion in the region $r>D/V\sim 30\,{\rm pc}\frac{D}{10^{27}\leavevmode\nobreak\ \rm cm^{2}s^{-1}}/\frac{V}{100\leavevmode\nobreak\ \rm km/s}$. The Alfv$\acute{e}$nic velocity in the interstellar medium is estimated as several $\rm km/s$ assuming the density of $\sim 1\,\rm cm^{-3}$ and magnetic fields of about $3\,\rm\mu G$ (see, e.g., Han, 2017). In this case the advection can be neglected close to the source. However, for a strong outflow of gas or a stellar wind near the source, with a speed as as high as $1000\leavevmode\nobreak\ \rm km/s$ (see, e.g., Domingo-Santamaría & Torres, 2006), at distances $r\geq 0.1$ pc, the advection can dominate the diffusion, resulting in the steady-state solution

where $E_{0}$ is obtained from the equation $\frac{r}{V}=\int_{E}^{E_{0}}\frac{dE^{\prime}}{P(E^{\prime})}$. The corresponding CR fluxes are shown in Fig.4. The fast advection results in the strong suppression of CR flux even compared with the fast diffusion case ($\chi=1$) in the vicinity of the source.

Moreover, in the calculation above we assumed that the power-law diffusion coefficient extends to low energy. Recent analyses show indications of a low-energy break in the diffusion coefficient in the ISM (Vittino et al., 2019; Weinrich et al., 2020), which was also predicted in Ptuskin (2006) because the damping of CRs terminates the cascade of turbulence and induces faster diffusion for LECRs. However, it is not straightforward to adopt this scenario in the vicinity of the CR sources, where the external turbulence is also stronger by an order of magnitude. The detailed calculation requires a self-consistent treatment of the CR propagation and magnetic turbulence development and damping near the CR sources, which need further investigations.

Together with magnetic fields and turbulent gas motions, Galactic CRs play a dominant role in the energy balance of the interstellar medium. Although subrelativistic (suprathermal) particles supply a substantial fraction of the CR pressure, the real contribution of LECRs remains uncertain. Because of the slow propagation and severe energy losses, the local LECRs make negligible contributions to the fluxes beyond 100 pc. For the same reason, LECRs are expected to be inhomogeneously distributed in the Galactic disk; subrelativistic protons and nuclei are expected to be concentrated around their acceleration sites. Because all potential CR source populations (SNRs, stellar clusters, individual stars, etc.) are linked in one way or another to the star-forming region, the effective confinement of LECRs in these regions is expected to produce feedback that stimulates the star formation through the ionization of the nearby molecular clouds. LECRs play a significant role also in the chemistry of the interstellar medium. Thus an unbiased, observation based information about LECRs at the sites of their concentrations in the Milky Way is crucial for understanding the fundamental processes linked to the dynamics and chemistry of the interstellar medium, star formation, etc. The most direct channel of information is provided by nuclear de-excitation $\gamma$-ray lines resulting from the interactions of protons and nuclei of LECRs with the ambient gas.

We explored the effect of the initial spectra shape, energy loss, and diffusion coefficient on the spatial and energy distributions of LECRs around the source that continuously injects accelerated protons and nuclei into the surrounding medium. LECRs, lose much energy during propagation, especially at energies below 100 MeV/nuc. At large distances, depending on the diffusion coefficient, their radial profiles become steeper than $1/r$ (see Fig.3), which is the typical radial distribution of particles at higher energies at which the energy losses can be neglected. We calculated that characteristics of the nuclear $\gamma$-ray line emission initiated by interactions of LECRs with the ambient gas using the cross sections gathered from experimental data and theoretical modeling with the code TALYS-1.95. We found that for the standard diffusion coefficient characterizing the propagation of CRs in the interstellar medium, the $\gamma$-ray emission is mainly produced in regions around the source within $\sim 25$ pc. The slower diffusion with the parameter $\chi<1$ makes the $\gamma$-ray source more compact and brighter. The expected $\gamma$-ray fluxes are unfortunately well below the sensitivity of current $\gamma$-ray detectors, however. With the arrival of the proposed $\gamma$-ray detectors that are dedicated to low energies, such as e-ASTROGAM (de Angelis et al., 2018) and AEMGO (McEnery et al., 2019), the detection of nuclear $\gamma$-rays, in particular, the lines at 0.847, 1.63, 4.44, and 6.13 MeV of nearby (within a few kpc) accelerators of LECRs would become feasible provided that the LECR accelerators are surrounded by dense ($\geq 100\,\rm cm^{-3}$) gaseous regions in which LECRs propagate significantly more slowly than in the interstellar medium, $\chi\ll 1$.