© L. Brahimi et al. 2020

1. Introduction

Cosmic rays (CRs) are a major component of the interstellar medium (ISM), along with the gas and magnetic field (Ferrière 2001), and have a non-negligible impact on the global dynamics of our galaxy (Grenier et al. 2015). Through their streaming they can drive galactic winds (Recchia et al. 2017; Girichidis et al. 2016). They contribute to the enrichment of the ISM in light elements by spallation and are also a strong source of ionization especially in the dense and cold ISM (Padovani et al. 2009). Furthermore, they may contribute to the turbulent dynamics of our galaxy through the generation of magnetic perturbations (Commerçon et al. 2019).

The origin of CRs is still under debate (Drury et al. 2001). Strong shocks in supernova remnants (SNRs) seem to be able to accelerate CRs via diffusive shock acceleration (DSA) up to a few hundred TeV (Lagage & Cesarsky 1983) or even to PeV at early times (when shock speeds v_sh are as high as 0.1 times the speed of light, e.g. Schure & Bell 2013). Once released into the ISM, CRs have a random walk induced by scattering off turbulent magnetic perturbations generated by large-scale plasma motions and/or by the waves that they self-generate while drifting at speeds greater than the local Alfvén speed (Wentzel 1968; Kulsrud & Pearce 1969; Skilling 1975; Yan & Lazarian 2004; Farmer & Goldreich 2004). However, currently there is no fully consistent theory capable of simultaneously explaining both CR acceleration at SNR shocks and their escape into the ISM (see Telezhinsky et al. 2012 for a discussion). Besides ionization studies, probes of ISM and CR interactions are essentially based on gamma-ray observations, such as for example the recent H.E.S.S. Tcherenkov telescope observations from the forward shock of the SNR RX J 1713.7-3946 (H.E.S.S. Collaboration et al. 2018) or W 28 (Aharonian et al. 2008) which provide new constraints on the CR escape process.

The ISM is a very complex medium and the way it is structured is still an open question. It is accepted that the ISM can be modeled in terms of distinct phases in rough pressure equilibrium. Depending on the degree of ionization, we find first the ionized phases: the hot ionized medium (HIM) and the warm ionized medium (WIM), and then the partially ionized phases: the warm neutral medium (WNM), the cold neutral medium (CNM), the diffuse molecular medium (DiM), and molecular clouds (MCs); see e.g. (Field et al. 1969; McKee & Ostriker 1977; Reynolds 1989; Gaensler et al. 2008). Partially ionized phases cover a minor fraction of the galactic disk volume but contain the majority of the mass. Most of the ISM matter is concentrated in MCs, the birthplace of stars. In partially ionized medium, low-energy CRs constitute an important source of ionization (Phan et al. 2018)¹, especially close to an accelerator where one can expect enhanced CR density for some amount of time (Vaupré et al. 2014; Nava et al. 2016). Cosmic rays can also exert a force on the gas through the production of magnetohydrodynamic (MHD) waves. The turbulent dynamics of WNM and CNM as well as MCs may be affected in the presence of an intense CR flux (Everett & Zweibel 2011; Commerçon et al. 2019).

In this paper, we particularly focus on the way CRs escape from SNR shocks into the ISM, the way they interact with the partially ionized ISM phases, and how their transport properties are modified in these media. Our model relies on the CR cloud model (CRC; see Sect. 2.2 for further details) proposed by Malkov et al. (2013), Ptuskin et al. (2008) and in Nava et al. (2016), referred to hereafter as N16. This study completes a parallel investigation of CR escape in the ionized ISM by Nava et al. (2019), referred to hereafter as N19. Cosmic rays at different energies escape from the remnant at different times depending on the properties of the ambient ISM. In this study we discuss three different escape models when radiative losses become important during SNR expansion. Depending on the local CR over pressure, CRs trigger a resonant streaming instability and produce magnetic fluctuations (Skilling 1971, 1975)². We assume in this study that the amplitude of the self-generated waves is weak and treat the problem according to the quasi-linear theory of CR transport (Schlickeiser 2002). In partially ionized phases, the waves triggered by CRs are damped by ion-neutral collisions and by nonlinear interactions with the large-scale turbulence (see N16 and Farmer & Goldreich 2004). This work provides some answers to different open questions: For how long/far does a SNR produce CR over-density in the neutral/molecular ISM? What is the space-time-dependent rate of self-generated waves around a SNR accelerator? How can we evaluate the impact of CR self-confinement over the grammage (see also D’Angelo et al. 2016)? Over what distance should one account for the effect of stronger local CR gradients?

In Sect. 2 we present our model: we describe the properties of neutral ISM phases adopted in this study, we discuss the different wave-damping processes, we discuss the way CRs escape, and we generalize the work proposed by Malkov et al. (2013) and N16. In particular we propose variants of the escape model based on the work of Sutherland & Dopita (2017). In Sect. 3 we present our results for each escape model for one particular medium at different energies. In Sect. 4 we calculate the grammage of CRs produced during their propagation close to their sources also considering a setup including multiple phases. Finally, we conclude in Sect. 5.

2. The model

2.1. Partially ionized interstellar medium phases

2.1.1. Phase properties

We first describe the main physical properties of weakly ionized ISM phases. We start from the assumption that the weakly ionized ISM is composed of hydrogen and helium with a density ratio n_H/n_He ≈ 13.28. The total plasma density is given by n_T = n_H + n_He. For each phase, we consider a thermal plasma using a two-fluid MHD model: each fluid (neutral and ionized) is composed of a dominant species described by its mean mass (μ_i = m_i/m_H and μ_n = m_n/m_H) and its density (n_n and n_i). We define the ionization fraction X = n_i/(n_i + n_n) of the plasma. We consider the plasma to be filled with a magnetic field composed of a regular and a turbulent component. The turbulent component results from an injection of free energy at large scales L ≃ 50 − 100 pc (see N16). The essential chemical and physical properties of cold ISM have been taken from the work of Snow & McCall (2006) and are summarized in Table 1. These values are subject to some uncertainties that are also included in the table.

Neutral atomic phases. The WNM and CNM are thermally stable solutions of the thermal instability together with the WIM (Field et al. 1969). Both phases are composed of atomic gas. The WNM and CNM are detected by HI emission and absorption line surveys. For both WNM and CNM, we assume that the dominant neutral species is atomic hydrogen (HI) while the dominant ion is ionized hydrogen (HII) in the WNM and ionized carbon (CII) in CNM, respectively (Snow & McCall 2006). The WNM is a relatively low density medium (n_T = 0.2 − 0.5 cm³) with a high equilibrium kinetic temperature (6 000 − 10 000 K) whereas CNM is denser (n_T = 20 − 50 cm⁻³) and colder (T = 30 − 100 K). The typical size of the phases are difficult to derive: 21 cm line surveys give empirical relations between the gas velocity dispersion σ and the typical size R of the region, suggestive of a turbulent cascade with a scaling σ ∝ R^k, with k ∼ 0.35 (Larson 1979; Roy et al. 2008). These analyses reconstruct σ from column density and line spin temperature measurements and have R ranging between 10⁻² and 100 pc. Another aspect is that these media are not isotropic in shape but rather show elongated structures (Heiles & Troland 2003). The aspect ratio between parallel and perpendicular length scales can be up to two orders of magnitude. We then decide to keep a maximum length scale for these two phases of ∼100 pc which is close to the typical size of large-scale turbulent motion injection. We fix the phase size interval to R ∼ 1 − 100 pc, smaller scales being of little interest in the framework of this study which involves CRC with sizes larger than 1 pc.

Diffuse molecular medium. A fundamental difference between cold atomic and cold molecular phases is the level of penetrating UV radiation (Snow & McCall 2006). As UV radiation decomposes molecules, molecular clouds need to be surrounded by a shielding medium in order to conserve their chemical equilibrium. This shield ensures a transition from an external part composed of atomic gas (HI dominated) to a molecular interior (H₂ dominated) and has a total minimal visual extinction magnitude of 0.2. For this medium, we assume a neutral fluid composed of a mixture of 50% neutral hydrogen (HI) and 50% molecular hydrogen (H₂). As enough interstellar radiation is present to photo-dissociate CO molecules, the dominant ionized species is CII as in the CNM phase. The DiM is also characterized by higher density (n_T = 100 − 500 cm⁻³); see Snow & McCall (2006). Its characteristic size is not well constrained; we consider typical sizes derived by Solomon et al. (1987) from CO surveys between 1 pc to 50 pc, the upper limit corresponding to giant molecular clouds.

We note that the properties of the phases and the dimension of our model constrain the number of cases we can study in our work. We use a model of CR propagation restricted to the propagation of CRs along a background magnetic field, hence a 1D model. This is the flux tube approximation (see Sect. 2.2), which imposes that the respective sizes of the phase and the CRC be R_phase >  R_CRC(t_esc), where t_esc is the time at which CRs at a particular energy escape from the CRC. We verify that this condition applies a fortiori to each phase in our study. Another aspect is that the partially ionized ISM is likely a relatively inhomogeneous medium, and we should not expect the SNR to propagate into a single phase over several tens of parsecs. To account for this inhomogeneity we need to investigate CR propagation in multi-phase ISM, which can require multi-dimensional simulations. This aspect is discussed in Sect. 5 and merits future dedicated study. However, in Sect. 4 we consider an example of multi-phase ISM in the context of the calculation of the CR grammage.

2.1.2. Magnetohydrodynamic wave damping in partially ionized media

Ion-neutral damping. Alfvén waves are supported by the motion of ions. In a two-fluid model energy transfer rate operates between neutrals and ions which produces a damping of the Alfvén waves. In the WNM phase, the ion-neutral collision rate ν_in expresses the energy transfer rate from ions to neutrals and is given in Eq. (19) of N16 by

$$\begin{array}{c}\hfill {\nu }_{\mathrm{in}}=n_{\mathrm{n}}\sigma v_{\mathrm{th}}=2n_{\mathrm{n}}(8.4\times {10}^{-9}{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}){\left(\frac{T}{{10}^4\mathrm{K}}\right)}^{0.4},\end{array}$$(1)

where σv_th is the fractional rate of change of the proton velocity v_th averaged over the velocity distribution, n_n (n_i) is the neutral (ion) volume density in cm⁻³, and T is the plasma temperature in K.

In the case of a colder phase (T ≤ 50 K), such as in CNM or DiM, the ion-neutral collision rate is expressed by (Jean et al. 2009)

$$\begin{array}{c}\hfill {\nu }_{\mathrm{in}}=n_n(2.1\times {10}^{-9}{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}).\end{array}$$(2)

The neutral-ion collision rate ν_ni expresses the energy transfer from neutrals to ions; it is linked with ν_in by the relation (Xu et al. 2015)

$$\begin{array}{c}\hfill {\nu }_{\mathrm{in}}=\chi {\nu }_{\mathrm{ni}},\end{array}$$(3)

where χ = (m_n/m_i)(X⁻¹ − 1) = ρ_n/ρ_i is the ratio between the mass density of neutrals and ions. Ion-ion ν_ii and neutral-neutral ν_nn collision rates express the internal energy exchanges within each fluid. In a weakly ionized plasma, ν_ii is negligible because of the low ion relative density. In the case of very low ionization rates, ν_nn may contribute to the damping of Alfvén waves (Xu et al. 2015). It is however found to be negligible in the cases under study here.

As shown by Eq. (3), the quantity of energy transferred between ions and neutrals depends on the ion/neutral plasma mass density. If X ≪ 1, then ν_in ≫ ν_ni controls the energy transfer between ions and neutrals. The Alfvén wave dispersion relation is calculated in a two-fluid approach by Soler et al. (2013) and has been solved using different approximations by Xu et al. (2016) and Lazarian et al. (2004). The imaginary part of the solution of the dispersion relation corresponds to the Alfvén wave damping rate; it is presented in Fig. 1 for each phase by a bold black line as a function of the energy of CRs. As in this work we only consider CR in resonance with waves, there is a one to one relationship between the particle Larmor radius r_L = γmc/qB₀ (or energy, in the relativistic domain) and the wave number k, such that kr_L ∼ 1. Here we introduce the CR charge q, mass m, and Lorentz factor γ = (1 − (v/c)²)^−1/2 where v/c is the ratio of the speed of the particle to the speed of light. The background magnetic field has a strength B₀. We observe two asymptotic behaviors. At low CR energy (E ≪ eB₀V_Ai/ν_in), ions and neutrals are weakly coupled and the energy transfer from ions to neutrals is maximal because their motions are uncorrelated. The asymptotic Alfvén wave damping rate is then given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{in}}=-\frac{{\nu }_{\mathrm{in}}}{2},\end{array}$$(4)

Fig. 1. Ion-neutral damping rate (${\mathrm{\Gamma }}_{\mathrm{d}}^{\mathrm{in}}$, black curves), Farmer and Goldreich damping rate (${\mathrm{\Gamma }}_{\mathrm{d}}^{\mathrm{FG}}$, red curves), and Lazarian damping rate (${\mathrm{\Gamma }}_{\mathrm{d}}^{\mathrm{L}}$, blue curves) as a function of resonant CR energy. The shaded region refers to the range of energy where the real pulsation of Alfvén waves is not defined, which is why the turbulent damping, which depends on the Alfvén, is not defined in this scale range.

and the Alfvén speed is the one produced by the ionic component defined by $V_{\mathrm{Ai}}=B_0/\sqrt{4\pi {\rho }_i}$. At high CR energy (E ≫ eB₀V_Aχ/ν_in), the motions of ions and neutrals become correlated and the energy transfer from ions to neutrals drops (McIvor 1977). The Alfvén waves damping rate is therefore given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{in}}=-\frac{{\xi }_{\mathrm{n}}{V}_{\mathrm{A}}^2{e}^2{B}_0^2}{2{\chi }^{-1}{\nu }_{\mathrm{in}}}{E}^{-2},\end{array}$$(5)

where ξ_n = ρ_n/(ρ_n + ρ_i) is the neutral fraction in the gas, and the Alfvén speed is given by $V_{\mathrm{A}}=B_0/\sqrt{{\rho }_{\mathrm{n}}+{\rho }_{\mathrm{i}}}$. Finally, in the energy range [eB₀V_Ai/ν_in, eB₀V_Aχ/ν_in] Alfvén waves do not propagate.

In Fig. 1 we see that the damping of Alfvén waves is stronger in the diffuse molecular phases (DiM) than in atomic phases (WNM, CNM). The cutoff band in the WNM is located around 10 TeV while in the other mediums it is located around 100 GeV. We can also observe that the forbidden propagation band gets wider as the ionization rate decreases.

Turbulent damping. Alfvén waves generated by CRs can also interact with large-scale injected Alfvénic turbulence. Although the sources of this turbulence are poorly known, among the main contributors we can cite supernovae explosions or galactic differential rotation (Mac Low & Klessen 2004). Interactions between wave packets propagating in opposite directions lead to a distortion of the wave packets. A CR-generated wave of wavelength λ ∼ r_L is “damped” when the size of the distortion produced by the interacting wave packet is δ ∼ λ (see Sect. 4.2 in L16). In this process, waves generated by CRs interact with background perturbations through a three-wave process (Ng & Bhattacharjee 1997) which leads to the production of larger k modes. This process is called turbulent damping (see Farmer & Goldreich 2004; FG04, Lazarian 2016; L16) and we use this term hereafter. It should be stressed that this process does not correspond to a real damping but produces a transfer of perturbations from the resonant scale to smaller scales.

In N16 we use the turbulent damping term proposed by FG04. The authors of L16 propose an extension of the FG04 model to other turbulent regimes, including magnetosonic turbulence. Their results are based on a particular approach to the phenomenon of MHD turbulence (see Lazarian & Vishniac 1999). By adopting this model we acknowledge its limitation. More involved discussions of wave interaction in either super- of sub-Alfvénic turbulence can be found in Matthaeus et al. (1994), Zank et al. (1996, 2012), Adhikari et al. (2015), albeit mostly in the context of solar wind turbulence. Adapted to ISM, these models may produce a different turbulent damping rate which could modify the generality of our conclusions especially at high CR energies where ions and neutrals are in the coupled regime of collisions. However, the turbulent damping rate has to increase to Fat least one order of magnitude to be competitive at 10 TeV (see Fig. 1 below).

We now reproduce the main results obtained by L16. In the sub-Alfvenic regime (with M_A = V_L/V_A <  1) the turbulence at large scales is weak mostly 2D; it cascades perpendicularly to the mean magnetic field direction. The velocity perturbations scale as v_k = V_L(kL)^−1/2. Below the scale $\mathcal{l}=L{M}_{\text{a}}^2$ the turbulence becomes strong and is described $v_{\mathrm{k}}=V_{\mathrm{L}}{(kL)}^{-1/3}{M}_{\mathrm{A}}^{4/3}$. In the above expression, M_A = V_L/V_A is the Alfvénic Mach number, V_L is the rms turbulence speed at the injection scale L, and k ≡ k_⊥ is the perpendicular wave number.

In the sub-Alfvenic, strong turbulence regime, (M_A <  1 and $k^{-1}<L{M}_A^4$) the damping term is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{L},\mathrm{sub},\mathrm{A},\mathrm{s}}\approx \frac{V_{\mathrm{A}}{M}_{\mathrm{A}}^2}{k^{-1/2}{L}^{1/2}}\mathrm{for}\frac{l_{\min }^{4/3}{M}_{\mathrm{A}}^{4/3}}{L^{1/3}}<k^{-1}<L{M}_{\mathrm{A}}^4,\end{array}$$(6)

where M_A = V_L/V_A is the Alfvénic Mach number, V_L is the rms turbulence speed at the injection scale L, and k corresponds to the turbulent scale. Here, l_min corresponds to the turbulence cutting scale and is given by the Eq. (13) below.

In a sub-Alfvenic, weak turbulence regime (M_A <  1 and $k^{-1}>L{M}_{\text{A}}^4$) the damping term is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{L},\mathrm{sub},\mathrm{A},\mathrm{w}}\approx \frac{V_{\mathrm{A}}{M}_{\mathrm{A}}^{8/3}}{k^{-2/3}{L}^{1/3}},\end{array}$$(7)

for $L{M}_{\text{A}}^4<k^{-1}<L{M}_{\text{A}}$. Alternatively, if k⁻¹ >  LM_A, then it is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{L},\mathrm{sub},\mathrm{A},\mathrm{w}}\approx M_{\mathrm{A}}^2\frac{V_{\mathrm{A}}}{L}\mathrm{for}{k}^{-1}<L,\end{array}$$(8)

and

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{L},\mathrm{sub},\mathrm{A},\mathrm{w}}\approx M_{\mathrm{A}}^2\frac{V_{\mathrm{A}}L}{k^{-2}}\mathrm{for}{k}^{-1}\gg L.\end{array}$$(9)

In the super-Alfvenic turbulence regime ($\mathcal{l}=L/M_{\text{a}}^3$), the turbulence at scales larger than $\mathcal{l}=L/M_{\text{a}}^3$ is isotropic and hydrodynamic and is described by a Kolmogorov model. Below this scale, strong turbulence applies and the cascade leads to anisotropic velocity perturbations, again following $v_{\mathrm{k}}=V_{\mathrm{L}}{(kL)}^{-1/3}{M}_{\mathrm{A}}^{4/3}$. The damping rate is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{L},\mathrm{super}}\approx \frac{V_{\mathrm{A}}{M}_{\mathrm{A}}^{3/2}}{L^{1/2}{k}^{-1/2}}\end{array}$$(10)

if $\frac{l_{\min }^{4/3}}{L^{1/3}}{M}_{\mathrm{A}}<k^{-1}<L{M}_{\mathrm{A}}^{-3}$, and

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{L},\mathrm{super}}\approx \frac{V_{\mathrm{A}}{M}_{\mathrm{A}}}{L^{1/3}{k}^{-2/3}}\end{array}$$(11)

if $k^{-1}>L{M}_{\mathrm{A}}^{-3}$.

Although not trivial, the choice of the injection scale L and the associated turbulent velocity V_L is essential to characterize the nature of the turbulence and the behavior of the turbulent damping. Following the same assumptions as N16, we chose L = 50 pc and V_L = max(V_A, c_s) where V_A corresponds to the Alfvén speed in the coupled regime and c_s is the sound speed. For each medium, this assumption leads to M_A ≈ 1 because V_A exceeds c_s. The MHD turbulence is in the trans-Alfvénic turbulent regime. However, in the DiM, values extracted from Table 1 lead to M_A ∼ 5. Increasing by a factor five however does not change the impact of turbulent damping as we limit our investigation to CRs with energies below 30 TeV. Even at these energies, ion-neutral collisions dominate the damping of slab waves. Acknowledging the uncertainties over the turbulence injection scale and the magnetic field strength and densities in these media, we decided to keep a trans-Alfvénic turbulence injection model in all of the phases under study. The choice for the injection scale L is rather arbitrary as very little is known about its value in the ISM. It is also very likely that L varies depending on the position in the disk and between the disk and the galactic halo. Typical scales of 1-100 pc are reported in Haverkorn et al. (2008). Furthermore, several sources of turbulent motion injection may compete at a given location, as may be the case in the local ISM (Burlaga et al. 2015). Turbulence can be injected by supernovae, massive star winds, galactic shear motions, HII region expansion, or by CR themselves; see e.g. Giacinti & Kirk (2017).

Figure 1 shows the two turbulent damping models. The red solid line represents the FG04 model and the bold blue line the L16 model. In a trans-Alfvénic regime, the turbulent damping rate is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{trans},\mathrm{L}04}\approx \frac{V_{\mathrm{A}}}{k^{-1/2}{L}^{1/2}}\mathrm{for}\frac{l_{\min }^{4/3}}{L^{1/3}}<k^{-1}<L.\end{array}$$(12)

The difference between the FG04 model and the L16 model appears at

$$\begin{array}{c}\hfill l_{\min }^{-1}={(2{\nu }_{\mathit{ni}})}^{3/2}{L}^{1/2}{V}_{\mathrm{A}}^{-3/2}\sqrt{1+\frac{V_{\mathrm{A}}}{2{\nu }_{\mathrm{ni}}L}},\end{array}$$(13)

corresponding to the low-energy cutting scale of the turbulence. Unlike in the FG04 model, the turbulent cascade in the L16 model is cut at energies below ∼500 GeV, 1.6 GeV, and 0.4 GeV in the WNM, CNM, and DiM phases, respectively. However, in practice, this difference is not relevant in our study because at low energy, the Alfvén wave damping rate is always dominated by ion-neutral collisions, and the turbulent damping only becomes dominant at high energies (∼100 TeV).

2.2. Escape of cosmic rays from supernova remnants

The process of CR leakage from SNRs is connected to the acceleration process in the source interior. As CRs leak out, they extract some energy flux which increases the compressibility of the shocked fluid (Berezhko & Ellison 1999). High-energy particles streaming ahead of the shock front produce magnetic fluctuations which ensure the confinement of lower energies at the shock (Bell et al. 2013). Until now however, no model has been built that describes in a self-consistent way the acceleration mechanism at the SNR strong shocks and the process of CR escape from the accelerator (see the discussion in Telezhinsky et al. 2012), so the transition between the following two locations.

The source interior where the magnetic perturbations are generated by streaming instabilities and the level of turbulence is high (δB ∼ B₀). Cosmic rays have a Bohm diffusion and their mean free path is small (λ ∼ r_L), meaning that they likely diffuse isotropically (D_∥ ∼ D_⊥). The CR pressure is in excess compared to the magnetic pressure.

The surrounding ISM where magnetic perturbations are injected at large scale (δB ≪ B₀). The mean free path of CRs is high (λ ≫ r_L) and they diffuse along the mean magnetic field (D_∥ ≫ D_⊥). Their pressure is in equipartition with the magnetic pressure.

Current CR leakage models do not explain this transition and induce an unavoidable jump in the mean free path of the CRs (Zirakashvili & Ptuskin 2012; Reville et al. 2009; Berezhko 1996). Ptuskin et al. (2008) and Malkov et al. (2013) propose a model to account for the escape of CRs from SNRs which is disconnected from the acceleration process. This model assumes that escaping CRs progressively leak out from the source which can be seen as a CR cloud, hence the name cosmic ray cloud (CRC) model. The model further considers that while escaping into the ISM, CRs propagate along the background magnetic field and that the cloud expansion is sufficiently slow to approximate the CR transport by a 1D process (see Ptuskin et al. 2008; Malkov et al. 2013). Cosmic rays, whilst propagating along the background magnetic field, trigger slab-type Alfvén waves through the resonant streaming instability. The generation of self-generated turbulence produces a reduction of the CR diffusion coefficient. Finally, the model is developed in the quasi-linear approximation, so it is valid only in the case of well-developed, weak turbulence, an assumption that requires verification a posteriori. The 1D approximation for the CR transport is valid for distances to the CRC that are shorter than the background magnetic turbulence coherence length.

In this paper, we use the kinetic code developed and presented in N16 and N19. We adapt the simulations to the case of partially ionized phases. We again use the setup proposed by Malkov et al. (2013) and N16. We describe the transport of CRs and waves along magnetic field lines by solving two coupled equations given by

$$\begin{array}{c}\hfill \frac{\partial P_{\mathrm{CR}}}{\partial t}+V_{\mathrm{A}}\frac{\partial P_{\mathrm{CR}}}{\partial z}=\frac{\partial }{\partial z}\left(D\frac{\partial P_{\mathrm{CR}}}{\partial z}\right),\end{array}$$(14)

for the CR pressure P_CR, and

$$\begin{array}{c}\hfill \frac{\partial I}{\partial t}+V_{\mathrm{A}}\frac{\partial I}{\partial z}=2({\mathrm{\Gamma }}_{\mathrm{growth}}-{\mathrm{\Gamma }}_{\mathrm{d}})I+Q,\end{array}$$(15)

for the waves energy density I. Equation (14) governs the evolution of the CR pressure P_CR in time and space. The space propagation is controlled by two processes: an advection with the scattering centers at a speed V_A and a random walk along the background field with a diffusion coefficient D. We note that in our case both V_A and D are energy dependent. The Alfvén speed depends on the energy regime (see Sect. 2.1.2): if ions and neutrals are coupled, the Alfvén speed is calculated using the total density, while in the decoupled regime only the ion density should be retained. The Alfvén speed can also be space-dependent in case the ISM is inhomogeneous. In that case, the above equations have to be modified (Skilling 1975; Zank et al. 2012). It should be stressed that the above equations are restricted to the quasi-linear theory. Any modelling of the back-reaction of the self-generated turbulence over the background turbulence is therefore beyond the scope of this simple model and requires the use of an MHD code properly coupled to CR kinetics. We return to this point in Sects. 4 and 5. The diffusion coefficient is D(E) = D_B(E)/I(E) where D_B(E) is the Bohm diffusion coefficient (see Bell 1978a,b) and I(E) is the energy density of the resonant Alfvén waves. The background diffusion coefficient is defined by D₀ = D_ISM(E/10 GeV)^0.5 and has been fixed to D_ISM = 10²⁸ cm² s⁻¹ which is of the order of the value inferred from direct measurements of the ratios of CR primary to secondary. Cosmic ray transport is linked to the turbulence by Eq. (15). In the latter, the growth of Alfvén waves is governed by the relaxation of the streaming instability which corresponds to the first RHS term. Skilling (1975) gives the expression for the growth rate:

$$\begin{array}{c}\hfill 2{\mathrm{\Gamma }}_{\mathrm{growth}}I=-\frac{V_{\mathrm{A}}}{W_0}\frac{\partial P_{\mathrm{CR}}}{\partial z},\end{array}$$(16)

where $W_0=B_0^2/8\pi$ corresponds to the magnetic energy of the medium. The second RHS term in Eq. (15) describes the damping of the waves due to ion-neutral collisions and turbulent damping Γ_d = Γ_IN + Γ_turb, L16. We also add the nonlinear Landau effect to handle cases where I(E) becomes sufficiently close to one (see N16 and N19). Finally, the last term allows to consider all external sources of turbulence, and we set Q = 2Γ_dI₀ where I₀ is the background turbulence level. Hereafter, the latter is assumed to be identical whatever the ISM phase. A detailed investigation of the large-scale-injected turbulence in partially ionized media is a difficult task and is postponed to a future study (see the discussion in Xu et al. 2016).

The physical CR pressure is defined by

$$\begin{array}{c}\hfill P_{\mathrm{CR}}=\frac{4\pi }{3c^3}{E}^{4}f(E),\end{array}$$(17)

where f(E) is the CR energy space phase function defined by

$$\begin{array}{c}\hfill f(E)=\frac{3c^3}{16{\pi }^2{R}_{\mathrm{s}}^3}\frac{1}{E^2}\frac{\mathrm{d}N}{\mathrm{d}E},\end{array}$$(18)

where R_s is the shock radius and the CR energy spectrum dN/dE is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=\frac{(2-\mathrm{\Gamma })W_{\mathrm{CR}}}{E_{\max }^{2-\mathrm{\Gamma }}-E_{\min }^{2-\mathrm{\Gamma }}}{E}^{-\mathrm{\Gamma }},\end{array}$$(19)

if Γ ≠ 2 and

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=\frac{W_{\mathrm{CR}}}{ln(E_{\max })-ln(E_{\min })}{E}^{-\mathrm{\Gamma }},\end{array}$$(20)

if Γ = 2. Here, W_CR = 10⁵⁰ erg is the amount of energy of the SN explosion imparted into the CR pressure, and E_min = 2 × 10⁻⁴ GeV and E_max = 2 × 10⁴ GeV are the energy limits of the CR spectrum.

Cosmic rays are initially confined in a region of size a = 2R_s/3 (see footnote 1 of N16 for an explanation) and an initial pressure $P_{\text{CR}}^0$. Outside, the initial CR pressure is negligible. At the beginning of the simulation, the energy of turbulent waves is equal to its background value I₀. Following Malkov et al. (2013), we introduce the partial pressure Π given by

$$\begin{array}{c}\hfill \mathrm{\Pi }=\frac{V_{\mathrm{A}}}{D_{\mathrm{B}}}{\mathrm{\Phi }}_{\mathrm{CR}},\end{array}$$(21)

where ${\mathrm{\Phi }}_{\text{CR}}={\int }_0^{\infty }\text{d}z{P}_{\text{CR}}=a{P}_{\text{CR}}^0$ corresponds to the CR flux across a tube section. The physical meaning of Π is discussed in detail in N16. We can identify two regimes of CR propagation depending on the value of Π:

If Π >  max(1, τ_diff/τ_damp) where τ_diff ≈ a²/D corresponds to the CR diffusion time-scale and τ_damp ≈ [2Γ_d]⁻¹ is the damping time scale of waves. Cosmic ray transport is nonlinear and the streaming instability is efficient enough to produce slab waves.

If Π <  max(1, τ_diff/τ_damp). Cosmic ray transport is linear, the CR flux is too low to produce waves and the diffusion proceeds with D = D_ISM.

If in that latter case in addition we neglect advection, the system of Eqs. (14) and (15) becomes

$$\begin{array}{c}\hfill \frac{\partial P_{\mathrm{CR}}}{\partial t}=D_{\mathrm{ISM}}\frac{{\partial }^2{P}_{\mathrm{CR}}}{\partial z^2},\end{array}$$(22)

where D_ISM = D_B/I_ISM and

$$\begin{array}{c}\hfill \frac{\partial I}{\partial t}=0,\end{array}$$(23)

which means that I = I_ISM. This situation is systematically achieved after a sufficiently long time as the flux of CRs across the tube section decreases and the condition Π <  max(1, τ_diff/τ_damp) is realized. In this case, the analytic solution is given by

$$\begin{array}{c}\hfill P_{\mathrm{CR}}(z,t)=\frac{1}{2}{P}_{\mathrm{CR}}^0[\mathrm{erf}\left(\frac{a-z}{2\sqrt{D_{\mathrm{ISM}}t}}\right)+\mathrm{erf}\left(\frac{a+z}{2\sqrt{D_{\mathrm{ISM}}t}}\right)].\end{array}$$(24)

In the above framework, Malkov et al. (2013) defined the CRC half-life time $t_{\frac{1}{2}}$ as the time it takes for the CR integrated pressure in the initial sphere to be divided by two. The CRC half time is defined by the relation

$$\begin{array}{c}\hfill {\displaystyle {\int }_{-a}^{+a}{P}_{\mathrm{CR}}(z,t_{\frac{1}{2}})\mathrm{d}z=2a\frac{P_{\mathrm{CR}}^0}{2},}\end{array}$$(25)

where $P_{\text{CR}}^0$ is the initial CR pressure in the cloud. To each CR energy corresponds a particular half-life time. In order to determine at which moment and at which SNR radius a given population of CR escapes, N16 proposes to intersect the half-life time evolution function resulting from the Eq. (25) $t_{\frac{1}{2}}(a)$ with the solution t(R_SNR) describing the time evolution of the radius of an SNR. The approach adopted here is more relevant for type Ia SNe which evolve in an unperturbed ISM contrary to the case of core-collapse SNe where the wind of the stellar progenitor shapes the ambient medium. However, if the stellar wind magnetic field has a relatively homogeneous structure over scales that are larger than R_SNR one may still treat the escape problem as 1D but now using the path length along magnetic field lines.

The evolution of SNRs has been described by Cioffi et al. (1988), Truelove & McKee (1999) and can be decomposed into four expansion stages :

Free expansion,

$$\begin{array}{c}\hfill t<t_{\mathrm{free}}\Rightarrow t\propto R_{\mathrm{SNR}};\end{array}$$

Sedov-Taylor (ST),

$$\begin{array}{c}\hfill t_{\mathrm{free}}<t<t_{\mathrm{PDS}}\Rightarrow t\propto R_{\mathrm{SNR}}^{5/2};\end{array}$$

pressure-driven snowplow (PDS),

$$\begin{array}{c}\hfill t_{\mathrm{PDS}}<t<t_{\mathrm{MCS}}\Rightarrow t\propto R_{\mathrm{SNR}}^{10/3};\end{array}$$

momentum-conserving snowplow (MCS),

$$\begin{array}{c}\hfill t_{\mathrm{MCS}}<t<t_{\mathrm{merge}}\Rightarrow t\propto R_{\mathrm{SNR}}^4.\end{array}$$

The constant of proportionality is defined by the radius of the SNR at the transition between the free and ST stages (see N16), at this stage the ambient medium starts to dominate the SNR dynamics.

$$\begin{array}{c}\hfill R_{\mathrm{SNR}}=5.0{\left(\frac{E_{\mathrm{SNR},51}}{n_T}\right)}^{\frac{1}{5}}{[1-\frac{0.05M_{\mathrm{ej},\odot }^{\frac{5}{6}}}{E_{\mathrm{SNR},51}^{\frac{1}{2}}{n}^{\frac{1}{3}}{t}_{\mathrm{kyr}}}]}^{\frac{2}{5}}{t}_{\mathrm{kyr}}^{\frac{2}{5}}[\mathrm{pc}].\end{array}$$(26)

The SNR radius depends on the SN mechanical energy E_SN in units of 10⁵¹ ergs, the ejecta mass M_ej,⊙ in units of solar mass, the ambient density n_T in units of cm⁻³, and the time in kiloyears.

In order to understand the properties of the CR escaping process, we determined for each phase (WNM, CNM and DiM) and for each population of CRs (10 GeV, 100 GeV, 1 TeV, 10 TeV) the time at which half of the initial CR pressure has escaped from the CRC and to what radius the shock extends. Our results are presented in Fig. 2.

Fig. 2. Properties of the expansion of the CRC in the ISM. Left: evolution of the half-life time of the CRC as a function of its initial radius a in parsec units. Green, blue, red, and black solid lines refer to particle energies at 10, 102, 103, and 104 GeV. The black dashed line refers to the time evolution of the radius of the SNR front shock in the considered phase: WNM for the top figure, CNM for the middle one, and DiM for the bottom one. These plots have been obtained for Γ = 2.2. Right: evolution of the half-life time of the CRC derived from the left plots as a function of CR energy. The blue solid line refers to a CR energy spectral index of 2.2 while the red dashed line refers to a spectral index of 2.0. The gray zone corresponds to the no-wave propagation zone; it does not lead to any (a, t1/2) solution.

Warm neutral medium. In the WNM phase, we recover the results obtained in N16 (see Fig. 2 upper row). We observe that low-energy CRs escape late, at 10 GeV R_esc ≈ 25 pc and $t_{\frac{1}{2}}\approx 42\mathrm{kyr}$, at the start of the radiative stage. As energy increases, the radius and escape time decrease linearly. We note that the escaping process is nonlinear (Π ≈ 7 × 10³,  2 × 10³,  0.5 × 10³ ≫ 1) at 10, 100, and 1000 GeV respectively meaning that the self-generated turbulence can –a priori– confine particles near the source. However, even if the growth rate is high it has to be compared with the damping rate, which can itself also be high. This comparison is the object of Sect. 3. We do not find any escaping solution for 10 TeV CRs. Figure 1 shows that Alfvén waves do not propagate, hence particles cannot trigger any waves and their propagation is mostly controlled by perturbations at scales larger or equal to their Larmor radius (Yan & Lazarian 2004).

Cold neutral medium. In the CNM phase, the description of the CRC expansion is more complex (see middle row of Fig. 2). Cosmic rays escape earlier, but the radiative stage of the SNR occurs also earlier. Cosmic rays at 10 GeV escape the latest (∼9 kyr) at a shock radius of ∼5 pc. The solutions have Π ≈ 31 × 10³,  4 × 10³,  25,  4 ≫ 1 or > 1 at 10, 100, 10³ and 10⁴ GeV, respectively. At energies of 100 GeV and 1 TeV, CRs escape at a similar time (∼1.2 kyr) and radius (∼2.7 pc). This is due to the fact that the damping rate of 1 TeV waves occurs in the coupled ion-neutral collision regime and hence drops by a factor 100 with respect to the damping of low CR energy. At 10 TeV, CRs are released at a small radius (<2 pc), hence their pressure gradient is still high.

Diffuse molecular medium. In the DiM phase, the self-generated turbulence has larger amplitudes than in the CNM and the WNM (see Fig. 2, lower row). This is essentially due to the small escaping CR radius which is a consequence of SNR dynamics in a medium with a high density. Generally, CRs in the DiM escape earlier than in the CNM and they all escape during the radiative phase, which also occurs earlier. At 1 TeV, CRs seem to escape slightly before 10 GeV CRs. Similarly to the case of CNM, this is due to the difference between the damping rates in the coupled and decoupled regimes. Figure 1 also shows that the turbulence in resonance with 100 GeV CRs is immediately damped explaining their linear diffusion. The half-life time of the CRC at 10 TeV is very small (∼0.07 kyr) and its radius is also quite small (<1 pc). This produces substantial turbulence; values of Π are about 97 × 10³, 16, and 9 at 10, 10³, and 10⁴ GeV, respectively.

It is also possible to study the CRC expansion properties by observing the cloud half-life time evolution as a function of the CR energy (see left parts of Fig. 2). In the WNM, the half-life time decreases linearly with the energy of CRs. Low-energy CRs are efficiently confined while high-energy CRs propagate quasi linearly. In the colder phases (CNM, DiM), we observe a gap around the energies corresponding to the nonpropagating turbulence range followed by a plateau due to the substantial level of self-generated turbulence connected with the drop of the damping rate. At very high energies, the half-life time decreases again because CRs are injected with an energy spectrum with power-law of negative index, this means that the highest energies account for less pressure. We also performed simulations with a harder CR spectrum with Γ = 2.0 but keeping the same power imported into CRs, the half-life time decreases at low energy because of the decreasing value of the pressure but this effect is not strong (see the right-hand panels of Fig. 2). At high energy the effect is inverted because of the generation of larger numbers of Alfvén waves due to the higher level of CR pressure.

In Table 2 we summarize four different turbulence regimes operating in the phases under study. The production of self-generated turbulence can be seen to depend on the wave growth rate and the wave-damping rate imposed by the density and the ionization rate of the external medium. Let us briefly describe them.

Regime 1 is that of strong nonlinear propagation. It represents the case where the wave growth rate is strong which is the case when the escaping radius is small, and the damping rate is reduced as in the coupled regime. This is the regime investigated by Malkov et al. (2013). The self-generated turbulence rate is growing almost exponentially with time leading to a strong CRs confinement. This regime is represented by a continuous decrease of the CR diffusion coefficient between t_1/2 and 4t_1/2 which is the typical limit time investigated in this study (see also N16 and N19).

Regimes 2 and 3 are transitory. Regime 2 corresponds to the case where the wave growth rate is high and the damping rate is strong, but still the difference between the two rates is small. Regime 3 corresponds to the case where the streaming instability is weak but produces waves weakly damped by interaction with the plasma. These regimes are represented by a decline followed by a growth of the CR diffusion coefficient between t_1/2 and 4t_1/2.

Regime 4 is the regime dominated by the damping process. It represents the case where the waves growth rate is small which is the case when the CRC radius is large, and the damping rate is strong. This regime is represented by a continuous growth of the CR diffusion coefficient between t_1/2 and 4t_1/2.

2.3. Alternative escape models

In the model presented above and in N16 it is necessary to check that the shock speed is always larger than 110 km s⁻¹ which corresponds to the limit speed below which the shocked medium does not generate UV radiation and cannot ionize the precursor medium (Shull & McKee 1979). In that case, ion-neutral collisions quickly damp any magnetic fluctuations (Sutherland & Dopita 2017). We have checked in our simulations that the speed of the shock at the escape time of the lowest CR energies is always greater than this limit.

Another aspect associated with the ionization radiation produced by the shock wave is the possibility to have extended ionization fronts in the upstream medium (Zirakashvili & Ptuskin 2018). This effect is possibly expected for fast adiabatic shocks (Allen et al. 2008). Extended ionization precursors if they exist should rather be expected in the WNM, whereas in denser media their extension should be reduced by a factor 1/n_H where n_H is the hydrogen density. Balmer dominated shocks however show limited extension of a few percent of the SNR size of the heated gas precursors (Lee et al. 2007).

However, the way CRs escape strongly depends on the evolution of the magnetic perturbations in the shock precursor, especially once the SNR has entered the radiative phase. We now consider two extreme cases of CR release in the ISM which reflect our ignorance of the process. Hereafter, we refer to model 𝒫 (for primary), the escape model presented in N16 and N19 and above. Model 𝒫 is the model default hereafter unless otherwise specified.

In order to differentiate the escaping solutions obtained with the alternative models from the ones obtained with the main model, we refer to the SNR escaping radius variable as a_esc and the CR escaping time as t_esc. These are analogous to the variables a and $t_{\frac{1}{2}}$ but are obtained from different physical assumptions, as described below (they are identical in model 𝒫).

Model ℱ (fast). Magnetic field fluctuations are rapidly damped in the precursor once the shock enters the radiative phase implying that all accelerated particles are released at a(E) = R_PDS if E <  E^*. Here, E^* is the maximal energy of particles escaping at the start of the radiative phase in model 𝒫, and depends on the properties of the external medium. As a consequence, we observe in Table 3 that the radii and escaping times of CRs drop with respect to the values obtained in model 𝒫. As escaping radii are smaller, the wave growth rate increases leading to stronger nonlinear propagation effects with respect to model 𝒫.

Model 𝒮 (slow). We suppose that the shock precursor is sufficiently ionized. Once in the radiative phase, the fluid is compressed and all particles that have not yet escaped at t_PDS stay confined until the shock velocity drops below 110 km s⁻¹. In this case we observe in Table 3 that the radii and the escaping times of CRs are large and the propagation is less nonlinear with respect to the model 𝒫.

3. Cosmic ray-propagation results

We performed simulations of CR propagation in homogeneous realizations of WNM, CNM, and DiM phases at different CR energies: 10 GeV, 100 GeV, 1 TeV, and 10 TeV. We investigated the time-dependent transition between the CRC and the CR background. We note that in order to avoid any nonphysical pressure gradients at small escape times ($<t_{\frac{1}{2}}$) we smooth the initial pressure step by a hyperbolic tangent that is 1% of the escape radius. We adopt absorbing boundaries at z = 0 and z = z_max = 500 pc.

The results are presented in the appendix in Figs A.1-1, A.2-2, A.3-3 for the WNM, A.4-1, A.5-2, A.6-3 for the CNM, and A.7-1, A.8-2, A.9-3 for the DiM and for the three different escape models (model 𝒫: 1, model ℱ: 2 and model 𝒮: 3). The color-code refers to the different ISM phases: yellow for WNM, green for CNM, and orange for DiM. The figures are separated in three parts. The upper part represents the spatial distribution of the CR pressure. The black dotted line represents the initial CRC pressure, and the blue, green, and red lines represent the simulation state at $t_{\frac{1}{2}}/4$, $t_{\frac{1}{2}}$, and $4t_{\frac{1}{2}}$, respectively, as in N16. We added when necessary two more solutions in brown and pink at $50t_{\frac{1}{2}}$ and ${10}^3{t}_{\frac{1}{2}}$, respectively, at high energies in the CNM and DiM phases to account for the slow dilution of the CRC when the damping rate drops. Dotted lines represent test-particle (TP) solutions while full lines represent numerical solutions. The middle part represents the spatial evolution of the CR pressure gradient along the mean magnetic field at the three above times. The bottom part represents the spatial evolution of the diffusion coefficient with time compared to D_ISM drawn with a black dotted line.

In the following sections, for each setup, CR leakage properties are described and associated qualitatively to a simplified propagation model as described in Table 2. The level of nonlinearity depends on the competition between the growth rate and damping rate. At t_1/2 we use “turbulence generation length” to refer to the width over which the turbulence growth rate exceeds the damping rate. This width is expected to decrease with time.

3.1. Warm neutral phase

Figure 2 shows our solutions in the WNM phase for model 𝒫. In this phase the CR escape begins at the start ST phase and ends at the start of the radiative phase. A general trend shows that the highest energies escape first, a result already obtained by N16. Figure A.1 shows the evolution of the CR pressure, the CR pressure gradient, and the diffusion coefficient evolution in space at three different times, t_1/2/4, t_1/2, and 4t_1/2.

At 10 GeV, CRs escape late ($t_{\frac{1}{2}}\approx 42$ kyr, see Fig. 2) when the shock radius is about 25 pc, at the beginning of the radiative phase. With the exception of early times, numerical solutions are close to TP solutions (see Fig. A.1). The maximum pressure gradient (10⁻³² erg/cm⁴) is reached at $t_{\frac{1}{2}}/4$ but is not strong enough to induce self-generated turbulence at rates much in excess of the background turbulence. This solution closely resembles regime 4, described above. The width of self-generated zone turbulence is about ∼50 pc.

At 100 GeV, CRs escape at ∼7 kyr when the shock radius is about 13 pc in the ST phase. This implies that the CR pressure gradient is higher, especially at earlier times (1.8 kyr) where it reaches ∼10⁻³⁰ erg/cm⁴. At these early times, CRs generate a noticeable amount of turbulence but as the SNR shock expands, the CR gradient drops and our solutions converge to the TP case. Between $t_{\frac{1}{2}}$ and $4t_{\frac{1}{2}}$ we observe that diffusion coefficients are close to the interstellar one. The width of generated turbulence is about ∼40 pc. Because the level of self-generated turbulence is at first high, it is best described by regime 2.

At 1 TeV, CRs escape at ∼1.5 kyr when the shock radius is about 6.4 pc still during the ST phase. At early times, CRs have a high pressure gradient but the generated turbulence is low because the most important part of them escape later than $t_{\frac{1}{2}}$. At this time, the value of the CR pressure gradient is about ∼5 × 10⁻³¹ erg cm⁻⁴ leading to a diffusion coefficient lower than 10²⁸ cm²/s. However, this turbulence is quickly damped because ions and neutrals are in the weakly coupled regime. The width of generated turbulence is ∼40 pc. It is again best described by regime 2.

In model ℱ, most CRs escape before the beginning of the radiative phase. Only 10 GeV CRs escape at the beginning of this phase at t_esc ≈ 24 kyr with a_esc ≈ 21 pc similarly to model 𝒫. As a consequence, the self-generated turbulence is stronger as showed in Fig. A.2. At t_esc/4, the diffusion coefficient is about ∼5 × 10²⁶ cm² s⁻¹ and increases gradually to 2 × 10²⁷ cm² s⁻¹ at t_esc to recover its interstellar value at 4t_esc. The generated turbulence width is similar to model 𝒫 but the ratio Γ_g/Γ_d is larger. This case is best described by regime 2.

In the model 𝒮, most CRs escape before the beginning of the radiative phase except CRs at 10 GeV which escape later during the transition t_PDS − t_MCS in the radiative phase, that is, at t_esc ≈ 60.7 kyr when the shock radius is about 28 pc. Our results are shown in Fig. A.3. In this case, the growth rate is weak and our solutions fall into regime 4. Numerical solutions follow thse of the TP case.

In summary, we find results similar to those obtained in N16. We show that at 10 TeV, CRs propagate following the TP case. We find, as expected, no strong differences among the three escape models with the exception of the propagation of 10 GeV CRs where nonlinear effects are stronger in model ℱ and weaker in model 𝒮.

3.2. Cold neutral phase

Figure 2 shows that in model 𝒫, in the CNM, CRs from 10 GeV to 10 TeV escape between the end of the ST and the beginning of the radiative phase. The solutions (a, t_1/2) do not follow a linear trend as in the case of the WNM. Indeed, around 1 TeV, CRs are more strongly confined as can be seen in Fig. A.4. Compared to the WNM case, CRs escape earlier in the CNM because of the SNR shock dynamics (see Eq. (26)) as well as the CRC expansion properties whose essential parameters are the initial pressure and the value of the nonlinear diffusion coefficient. Test-particle solutions are identical for all phases (see Eq. (24)).

At 10 GeV, CRs escape at t_1/2 ≈ 8 kyr when the cloud size is small (a = 5 pc). The SNR shock is in the radiative phase. Figure A.4 shows that the solutions are relatively close to the TP ones. At t_1/2/4, the CR pressure gradient is high and reaches ∼10⁻²⁸ erg cm⁻⁴ but then decreases rapidly to 10⁻³⁰ erg cm⁴ at t_1/2 and the self-generated turbulence level drops accordingly. At 4t_1/2, the self-generated Alfvén waves are completely damped and the linear propagation regime is recovered. This case is described by the regime 2.

At 100 GeV, CRs escape early (t_1/2 ≈ 1.2 kyr) during the ST phase when the cloud size is about 2.7 pc. As shown in Fig. A.4, the CR pressure gradient is ∼10⁻²⁸ erg cm⁻⁴ and then slightly decreases until t_1/2 to finally drop below 10⁻³⁰ erg cm⁻⁴ at 4t_1/2. Diffusion is substantially suppressed (D ∼ 10²⁶ − 10²⁷ cm² s⁻¹) until t_1/2. After this time, the solution converges to the TP case. This case is described by regime 2.

At 1 TeV, CRs escape at the same time as 100 GeV CRs (∼1.15 kyr). The generated turbulence rate is the same as for 100 GeV CRs and the diffusion is highly suppressed (D ∼ 3 × 10²⁶ cm² s⁻¹). The distance over which turbulence is generated however is different between the two energies (∼40 pc at 1 TeV versus ∼20 pc at 100 GeV). The associated turbulence regime is regime 1 at earlier times. Here the solutions are well approximated by the solutions of Malkov et al. (2013). At 4t_1/2, waves are damped and the linear propagation regime is recovered.

At 10 TeV, CRs escape during the ST phase (t_1/2 = 0.1 kyr) when the cloud size is about 1 pc. At t_1/2/4, the CR pressure gradient is high, about 10⁻²⁷ erg cm⁻⁴, and the diffusion coefficient is ∼10²⁷ cm² s⁻¹. The pressure gradient then decreases gradually during the simulation while the diffusion coefficient decreases gradually from ∼10²⁷ cm² s⁻¹ to 10²⁶ cm² s⁻¹ at 50t_1/2. This effect is a consequence of the weak damping rate of waves at 10 TeV (see Fig. 1). At this energy, the damping time of waves is about 100 kyr. Suppressing the damping term in the equation is a good approximation to describe the propagation which is described by regime 1 and the solutions derived in Malkov et al. (2013). Furthermore, the distance at which significant turbulence is generated is large, about ∼100 pc. We find a ratio of δB/B₀ ∼ 0.8. In order to prevent the level of self-generated turbulence from overtaking the quasi-linear regime we include the effect of nonlinear Landau damping and the effect of perpendicular diffusion by modulating the solution along the background magnetic field with a perpendicular dilution. At each time step, we calculate a new pressure given by

$$\begin{array}{c}\hfill P_{\mathrm{CR}}^{\mathrm{new}}(z,y=0,t)=\frac{1}{2}{P}_{\mathrm{CR}}^{\mathrm{old}}\times (\mathrm{erf}(-\mathrm{\Delta }{y}^2/4D_{\perp }\mathrm{\Delta }t)+\mathrm{erf}(\mathrm{\Delta }{y}^2/4D_{\perp }\mathrm{\Delta }t)),\end{array}$$(27)

where D_⊥ = I²D_∥ and the spatial step perpendicular to the mean magnetic field direction is $\mathrm{\Delta }y=\mathrm{\Delta }z=\sqrt{2C{D}_{\Vert }\mathrm{\Delta }t}$, where C = 0.2 is the CFL constant. In practice, the impact of the perpendicular diffusion is found to be limited. At late times, the CR distribution slowly converges to the TP solution. These solutions are however at the limit of validity of the quasi-linear theory. This point is discussed in Sect. 3.4.

Figure A.5 shows the solutions of model ℱ. Cosmic rays escape before the start of the radiative phase except at 10 GeV where the escape proceeds at the start of the radiative phase with t_esc ≈ 1.9 kyr when the shock radius is about ∼3 pc. In that case, numerical solutions are highly nonlinear over a large distance of ∼20 pc. Pressure gradients are high with ((∂P_CR/∂z)_max ≈ 10²⁷ erg cm⁻⁴, and D_min ≈ 5 × 10²⁴  cm² s⁻¹) testifying to a strong turbulence rate. In this model, the propagation regime is regime 2 but the growth rate exceeds the damping rate until 4t_esc.

Figure A.6 shows the solution of model 𝒮. Cosmic rays escape before the start of the radiative phase except at 10 GeV where CRs escape between the PDS and the MCS phase at t_esc ≈ 11 kyr when the cloud size is about 5.4 pc. Only the solution at t_esc/4 shows nonlinear behavior and reduced diffusion coefficients (∼10²⁶ cm² s) over a spatial range of 20 pc. At later times the turbulence is gradually relaxed and we tend to the propagation in regime 2 where Γ_d becomes higher than Γ_g. As expected however the solutions are less nonlinear with respect to model 𝒫.

In summary, model 𝒫 solutions at high energy (10 TeV) are in the nonlinear regime over a long time even if the growth rate decreases with time because the damping in the coupled regime has dropped. Comparing model 𝒫 with the two other models we find that the results are the same for 100 GeV, 1 TeV, and 10 TeV CRs for the three escape models. However, at 10 GeV, CRs are strongly confined in model ℱ. At this energy, diffusion is strongly suppressed by a factor exceeding 1000 with respect to standard values.

3.3. Diffuse molecular phase

Figure A.7 shows our results for model 𝒫 in the DiM phase. We find that here, CRs escape slightly earlier, and from a smaller cloud than in the CNM phase. This has important consequences on the level of self-generated turbulence and on the CRC expansion properties. The confinement time of CRs as a function of their energy does not evolve in a monotonic way. In particular we observe that at 1 TeV and 10 TeV CRs are well confined.

At 10 GeV, CRs escape at t_1/2 = 3 kyr during the radiative phase when the cloud size is about 2 pc. The solutions at t_1/2 and 4t_1/2 deviate from the TP ones. Diffusion is suppressed (D_min ≈ 5 × 10²⁵ − 2 × 10²⁶ cm² s⁻¹) but the self-generated turbulence drops at 4t_1/2 and this case falls into regime 2 with Γ_d >  Γ_g. The width over which turbulence is produced is about 17 pc.

At 1 TeV, CRs escape relatively late (t_1/2 ≈ 2.4 kyr) during the radiative phase when the cloud size is about ∼1.9 pc. High CR pressure gradients suggest that self-generated turbulence levels are high (∂P_CR/∂z_max ≈ 10⁻²⁸ erg cm⁻⁴ at t_1/2/4 and t_1/2) which is confirmed by the values of the diffusion coefficients (D_min ≈ 2 × 10²⁶ cm² s⁻¹ for t_1/2/4 and t_1/2). The turbulence level decreases very slowly at 4t_1/2 : (∂P_CR/∂z)_max ≈ 10⁻²⁹ erg cm⁻⁴ and D_min ≈ 10²⁷ cm² s⁻¹. This case is associated to regime 2 where Γ_g >  Γ_d.

At 10 TeV, CRs escape very early (t_1/2 = 0.07 kyr) during the ST stage when the cloud size is about 0.75 pc. The growth rates are high with (∂P_CR/∂z)_max ≈ 10⁻²⁷ erg cm⁻⁴ and D_min ≈ 3 × 10²⁶ − ×10²⁶ cm² s⁻¹. Turbulence production increases with time as in the case of 10 TeV CRs in the CNM and rapidly reaches a high turbulence regime with δB/B ∼ 0.8. Nonlinear Landau damping and dilution by perpendicular diffusion are both included. This case corresponds to regime 1 where Γ_g ≫ Γ_d. We note that the solution at t = 10³t_1/2 shows a slight increases at z ∼ 400 pc, this corresponds to a boundary effect which has no impact on the solution.

Figure A.8 shows the solutions for model ℱ. Cosmic rays at 100 GeV and 10 TeV escape during the ST phase while those at 10 GeV and 1 TeV escape at the beginning of the radiative stage at t_esc ≈ 0.5 kyr when the cloud size is about ∼1.2 pc. Turbulence levels are higher than in the model 𝒫: (∂P_CR/∂z)_max ≈ 10⁻²⁶ erg cm⁻⁴ and D_min ∼ 10²⁵ cm² s⁻¹ in the case of 10 GeV CRs. We observe that turbulence levels remain constant over a width of about 50 pc. These solutions correspond to propagation in regime 2 where Γ_g >  Γ_d. Solutions at 1 TeV are strongly nonlinear. The maximum value of the CR pressure gradient is about ∼10⁻²⁷ erg cm⁻⁴ and remains constant throughout the simulation while the turbulence rate increases with time (at t_esc/4, D_min ≈ 2 × 10²⁵ cm² s⁻¹ while at 4t_esc this value becomes lower than 10²⁵ cm² s⁻¹) over a width of about 40 pc. The CRC then relaxes towards the TP solution at later times. The waves growth rate is highly dominant and is in line with the propagation regime 1.

Figure A.9 shows the solutions for model 𝒮. At 100 GeV and 10 TeV, CRs escape during the ST phase while those at 10 GeV and 1 TeV escape during the radiative phase at t_esc = t_PDS = 4.7 kyr when the cloud size is about 2.3 pc. At 10 GeV, turbulence levels are relatively weak compared to model 𝒫. Only CRs escaping before t_esc/4 produce noticeable turbulence (D_min ∼ 10²⁶ cm² s⁻¹). The damping rate dominates at later times implying that the propagation regime is regime 4. At 1 TeV, CRs produce turbulence until relatively late times but finally relaxe to the TP solution.

In summary, we find the same trend as for the CNM at high energy: we believe the diffusion to be in the nonlinear regime because the damping rate drops. Comparisons between the models show that the results are the same for 100 GeV and 10 TeV CRs because they escape before the beginning of the SNR radiative stage and the shock velocity is higher than 150 km s⁻¹. However, the propagation properties of the 10 GeV and 1 TeV CRs depend on the escape model. In model ℱ the turbulence level is higher and increases with time compared to the solutions of model 𝒫. Model 𝒮 shows less nonlinear solutions, as expected.

3.4. Propagation at high energies

Here we comment over our solutions at 10 TeV in the CNM and DiM phases. At this energy diffusion is found to be highly nonlinear because of the drop of the damping term and because escape occurs with small CRC sizes. Relaxation towards the TP indeed occurs, although it does take a longer time with respect to smaller energies. However, our solutions are at the limit of validity of the quasi-linear theory and should be considered with some care. These cases require that the 2D diffusion be properly taken into account. This effect should contribute to diluting the local CR gradient and to decreasing the level of self-generated turbulence. In the meantime we also note that CRs at this energy have a large diffusing flux, meaning that they can trigger a large current that can further amplify magnetic perturbations through the triggering of the nonresonant streaming instability (Inoue 2019), producing further confinement. These two aspects however require a much more demanding numerical effort. Finally, our injection model is crude and more specific modeling is required find the phase of the SNR evolution at which CRs at high energy escape. All these aspects will be investigated in a forthcoming study.

4. Discussion: Residence time and grammage calculations

The grammage of CRs is defined by the density column of matter crossed by CRs along their pathway in the ISM. This quantity is highly dependent on the particle transport properties. The grammage is defined by X ≈ 1.4m_pn_Tcτ_res where τ_res is defined as the CR residence time. At galactic scales, this time is observationally inferred from secondary to primary ratios as is the case for boron (B) to carbon (C). However, this measurement is an average over the whole CR journey in the Galaxy and it includes the transit in the galactic disc and in the halo. In order to more appropriately calculate the CR pathway, it is important to evaluate the contribution of the propagation close to the CR sources where it is anticipated that the CR self-generated turbulence may control CR transport (D’Angelo et al. 2016).

We consider two ISM configurations. First, we calculate the residence time in the case of homogeneous phases as a function of the particle energy in Sect. 4.1, and subsequently we consider cases with multiphase ISM configurations in Sect. 4.3.

4.1. Single-phase studies

In our semi-cylindric 1D problem configuration the residence time is defined as the characteristic time for a CR to escape from the CRC to a fixed distance z_*. We derive τ_res using a method based on the fact that the CR trajectory outside the CRC is dominated by diffusion (an assessment which we have verified). The residence time is calculated from an average of the square of the diffusion distance (see N19):

$$\begin{array}{c}\hfill z_{\ast }^2=\frac{{\int }_0^{\infty }{z}^2{P}_{\mathrm{CR}}(E,z,{\tau }_{\mathrm{res}})\mathrm{d}z}{{\int }_0^{\infty }{P}_{\mathrm{CR}}(E,z,{\tau }_{\mathrm{res}})\mathrm{d}z}·\end{array}$$(28)

Once z_* is fixed, we deduce τ_res by balancing the two terms of Eq. (28).

We calculate the CR residence time (τ_res) and the associated grammage (X) in a given ISM configuration. The results are presented in Fig. 3 for the WNM (top), CNM (middle), and DiM (bottom). Solutions derived are represented by a continuous line. The positions z_* from the cloud center are indicated by different colors: black (10 pc), red (30 pc), blue (50 pc), and green (100 pc). All calculations are realized in the framework of model 𝒫 (see Fig. 2).

Fig. 3. Grammage (right vertical axis) and residence time (left vertical axis) of CRs as a function of their energy in three different phases of the ISM : WNM (top), CNM (middle), and DiM (bottom). Solid lines show the results obtained at different distances z* from the center of the source: green (100 pc), blue (50 pc), red (30 pc), and black (10 pc). The shaded region marks the no-wave propagation zone. Oscillations at high energy in the CNM and DiM are due to the limited energy resolution used to reconstruct the different curves.

4.1.1. Warm neutral medium

In the WNM at z_* = 100 pc (z_* = 50 pc and z_* = 30 pc respectively), the grammage evolves as E^−0.5 and is about X = 0.1 g cm⁻² (X = 4 × 10⁻² g cm⁻² and X = 2 × 10⁻² g cm⁻² , respectively) at 10 GeV and decreases linearly to reach X = 10⁻² g cm⁻² (X = 2 × 10⁻³ g cm⁻² and X = 1 × 10⁻³ g cm⁻² respectively) at 3 TeV.

In the WNM, we find in Sect. 3 (see also Fig. A.1) that nonlinear effects are weak at almost all energies because ion-neutral interactions proceed in the decoupled regime. This reflects the solutions obtained by the two methods as both show a trend with a grammage which can be parametrized as X ≃ X_WNM(E/10 GeV)^−0.5, with X_WNM ≃ 0.1 g cm⁻². This grammage evolution is characteristic of the transport in the background turbulence.

4.1.2. Cold neutral medium

In the CNM at z_* = 100 pc (z_* = 50 pc and z_* = 30 pc, respectively), the grammage evolves as E^−0.5 and is about X = 10 g cm⁻² (X = 3 g/cm² and X = 1.5 g cm⁻² respectively) at 10 GeV. The grammage in this model is not defined around 300 GeV due to the Alfvén waves nonpropagation band (see Fig. 1). We then see a slight softening of the slope around 3 TeV and a hardening up to 10 TeV where the grammage reaches X = 4 g cm⁻² (X = 2 g cm⁻² and X = 1 g cm⁻², respectively).

At low energy (E <  300 GeV), the grammage approximately follows X ∝ E^−0.5 which is characteristic of the propagation in background turbulence. The discussion in Sect. 3 and the results in Fig. A.4 show a trend very similar in the CNM and WNM at these energies. We do not have any solution around 300 GeV associated to the forbidden propagation zone. In that energy regime, the transport is completely controlled by the background turbulence. Above 300 GeV we can explain the change of behavior by the fact that diffusion is suppressed (D ∼ 5 × 10²⁶ cm² s⁻¹) until we reach the half-life time of the CRC. At high energies, the diffusion is even more suppressed over larger scales and longer times which results in an increase of the grammage. We emphasize however that the turbulence levels obtained in this configuration are at the limit the quasi-linear CR transport theory, and more realistic calculations (including 2D effects, improved escape model, and the triggering of other types of instabilities) should lead to a more accurate estimation.

4.1.3. Diffuse molecular

In the DiM at z_* = 50 pc (z_* = 30 pc and z_* = 10 pc, respectively), the grammage value is about X = 40 g cm⁻² (X = 20 g cm⁻² and X = 4 g/cm⁻², respectively) at 10 GeV. We do not have any values of the CR grammage between 20 and 100 GeV because of the absence of escaping solutions due to the nonpropagation of Alfvén waves. Beyond 3 GeV, we observe a softening up to 10 TeV for z_* = 30 pc and z_* = 50 pc. Close to the source (z_* = 10 pc) we observe that the grammage decreases more strongly up to 3 TeV and then hardens between 3 and 10 TeV, where it finally reaches X = 5 g cm⁻² (X = 1.5 g cm⁻² and X = 0.2 g cm⁻², respectively).

In the DiM, we observe that the grammage evolution with the CR energy is modified compared to the low-energy CNM values and those in the WNM. Softer slopes are induced by stronger CR confinement at high energy (see Appendix A.7). The confinement effect is also visible in the CR pressure distribution.

4.1.4. Single-phase study

In the WNM, the grammage does not exceed 0.1 g cm⁻² at 10 GeV and scales as E^−0.5. Still at 10 GeV, in the CNM and DiM the grammage is higher, namely ∼10 g cm⁻² for the CNM and ∼20 g cm⁻² for the DiM. These values are large, but it should be borne in mind that the CNM and DiM cover only a small fraction of the galactic disk volume, hence these media likely only contribute to a small fraction of the total grammage of the CR detected on Earth. Nonetheless, these results also show that the CR grammage may vary considerably depending on the location in the ISM. This could have implications for the production of light elements in the cold ISM located close to a CR source. The unexpected result we find is obtained at high energy where the grammage increases with the energy in the CNM and DiM. High values are found at 10 TeV of the order of 2 g cm⁻² and 5 g cm⁻² at 50 pc from the source in the CNM and DiM, respectively. Here again, due to the small volume-filling factor of these phases, we do not expect to have a strong contribution to the CR grammage at these energies. Furthermore, these values are likely overestimates since the results at these energies and in these phases are at the limit of the validity of quasi-linear theory. Finally, the results are obtained for a homogeneous ISM. More realistic multi-phase calculations are presented below for low CR energies and in a forthcoming work for all CR energies.

4.2. Hadronic collision losses

While propagating in dense media, CRs are subject to proton-proton (pp) collisions and energy losses. The typical loss time by pp interaction is

$$\begin{array}{c}\hfill t_{\mathrm{loss},\mathrm{pp}}={\left(K_{\mathrm{pp}}{\sigma }_{\mathrm{pp}}{n}_{\mathrm{H}}c\right)}^{-1},\end{array}$$(29)

where K_pp ≃ 0.2 and σ_pp ∼ 30 mbarn is the pp collision cross section at the interaction threshold and is only increasing with the logarithm of the particle kinetic energy (Kafexhiu et al. 2014). We have checked that even in the DiM at high energy we always find t_res/t_loss, pp ≪ 1. Cosmic rays then propagate in dense phases without suffering strong losses, which can then be a fortiori neglected in Eq. (14).

4.3. Multiple-phase studies

We extend our study to more realistic ISM environments by calculating the CR grammage in an inhomogeneous ISM. We restrict our analysis to CR energies at 10 GeV for three reasons: (1) it is in that energy domain that the grammage has been primarily derived from secondary to primary direct measurements, (2) in this energy domain the self-generated turbulence possibly controls the CR propagation in the ISM (Blasi et al. 2012), and therefore this is a good regime to test the impact of such turbulence over CR propagation in different media, and (3) in the decoupled ion-neutral regime the Alfvén speed is to be taken with respect to the ion density; it does not vary significantly from one phase to another. This means that CRs do not suffer from strong adiabatic losses while propagating from one medium to another, and therefore Eqs. (14) and (15) are still a good approximation to describe CR and wave-coupled evolution. At high energies, ion-neutral collisions occur in the coupled regime, and the Alfvén speed is to be taken with respect to the total gas density. The Alfvén speed varies considerably and then the kinetic equations for particles and waves have to be generalized. This requires a more complete modeling, which is postponed to a future work.

For the multiple-phase case, we propose two setups from the least-dense to the densest phase:

Setup 1: The CRC evolves first in a WNM phase up to a distance of 50 pc where it reaches a CNM phase up to 130 pc. Beyond this distance, the WNM is retrieved up to 200 pc which is the limit of our simulation box.

Setup 2 : This setup is structured in the same way as Setup 1 except that we set a DiM cloud in the middle of the CNM phase, between 80 and 100 pc.

In order to avoid strong pressure and diffusion coefficient discontinuities between phases, all space-dependent variables appearing in Eqs. (14)–(15) are smoothed with a hyperbolic tangent function of a width equal to 3 pc. This transition can be changed and we verified that it does not affect the overall results of the grammage calculation. Our results are presented in Figs. 4 and 5 for the space-time evolution of the CR pressure and diffusion coefficient and for the grammage, respectively. In both figures, for each setup, WNM, CNM, and DiM phases are represented by green, orange, and red colors, respectively. In Fig. 4, the plot in the top left hand corner shows the evolution of the CR pressure at times corresponding to those chosen for the simulation of the expansion of the CRC in the WNM at 10 GeV. The plot in the bottom left corner shows the evolution of the diffusion coefficient. Cosmic ray pressures and diffusion coefficients are normalized to the solution obtained in the WNM case. In the top right corner, CR pressures for Setup 1 (solid line) are compared to those of Setup 2 (dashed line). In the bottom right, the diffusion coefficients for the three setups (WNM: dotted line, WNM-CNM: solid line, WNM-CNM-DiM: dashed line) are compared to those of the TP solution. In Fig. 5 the CR residence time is represented in the left plot while the CR grammage is shown in the right one.

Fig. 4. Evolution of the CR pressure (top figures) and their associated diffusion coefficient (bottom panels) in space at different times of the simulation based on the CRC half-life time for 10 GeV CRs in the WNM: 10.52 kyr (blue), 42.06 kyr (green), and 168.25 kyr (red). The simulation setups are represented by the different line styles: solid line for the WNM-CNM, dashed line for the WNM-CNM-DiM, and dotted line for the homogeneous WNM. The vertical black lines represent the initial size of the CRC. The phase decomposition for each setup is represented between the figures by the colored lines: green for WNM, orange for CNM, and red for DiM. Top-left figure: relative evolution of the CR pressure for WNM-CNM and WNM-CNM-DiM setups compared to the CR pressure evolution in the case of propagation in a homogeneous WNM phase. Top-right panel: relative CR pressure evolution for the WNM-CNM-DiM (PCR, 2) and the WNM-CNM (PCR, 1) setups, respectively. Bottom-left panel: relative evolution of the CR diffusion coefficients compared to the evolution of the CR diffusion coefficient in the case of propagation in a homogeneous WNM phase. Bottom-right panel: normalized relative evolution of the CR diffusion coefficient around the CNM and CNM-DiM phases. We recall that $D_{\text{CR}}^0={10}^{28}{\text{cm}}^2{\text{s}}^{-1}$.

In order to simplify our analysis we use α_i  = P_CR, i(z, t) / P_CR, WNM(z, t) and β_i  =  D_CR, i(z, t)/D_CR, WNM(z, t) where i  =  1,  2: the ratio between CR pressure and diffusion coefficients of the setup i with the CR pressure and diffusion coefficients in the case of a propagation in an homogeneous WNM (see Fig. A.1).

4.3.1. Setup 1: WNM-CNM

For this setup, Fig. 4 shows that the CR pressure distribution is modified by ∼10% compared to the case of propagation in a homogeneous WNM. At 10.52 kyr, the curve (blue solid line) follows the same behavior as in the homogeneous WNM case up to a distance of ∼30 pc. Beyond 30 pc the ratio α₁ decreases and increases suddenly at the level of the transition WNM-CNM (50 pc). Here, α₁ reaches ∼1.15 around 70 pc before decreasing again and tends to ∼0.9 after the transition CNM-WNM. At 42.06 kyr, we observe that α₁ has a slightly lower value than 1 at the spatial origin of the simulation. The dip close to the WNM–CNM transition is then accentuated (α₁ ∼ 0.9) and α₁ increases again reaching its maximum (∼1.2) around 110–120 pc and finally decreases. At 168.25 kyr, α₁ is lower than 0.9 and is almost flat except at the CNM–WNM transition where it begins to increase.

The diffusion coefficient evolution is also affected by the structure of the ISM. At 10.52 kyr, the ratio β₁ is constant (and equals 1) up to the WNM–CNM transition (50 pc) and starts to increase to reach ∼1.3 before gradually decreasing to retrieve the value of 1 at the CNM–WNM transition. At 42.06 kyr, the ratio β₁ is lower than 1 (∼0.9) and starts to increase up to the WNM–CNM transition to reach ∼1.6 and decreases progressively to reach 1 near the CNM–WNM transition. At 168.25 kyr, β₁ is again close to 1 and starts to increase slowly from the WNM–CNM transition to reach 1.1 close to the CNM–WNM transition where it decreases abruptly to the value of 1 beyond.

The different behaviors between the case of a propagation in a homogeneous WNM and Setup 1 essentially come from the Alfvén waves damping timescales (see Fig. 1). In the WNM, the damping time is about 10 yr while in the CNM it is about 1 yr. In the master Eqs. (14) and (15), the right-hand term, 2(Γ_growth − Γ_d)I, decreases from the WNM to the CNM. As Π is a conserved quantity regardless of the phase, the figure can be explained easily: the ratio of the pressures starts at 1 in the WNM, then as the damping is more severe in the CNM, less waves are produced at the WNM–CNM transition. The CR pressure in the CNM goes first below and above the WNM solution because CRs diffuse faster in the CNM. Finally, the CNM solution catches the WNM ones beyond, but the ratio reaches 0.9 instead of 1 because of the conservation of Π. The peak in CR pressure in the CNM is slowly advected out with time as the solution relaxes to the TP case. The bottom right part of Fig. 4 shows that the difference between the self-generated and background diffusion coefficients is not large, unlike in homogeneous cases. This small difference at 10 GeV is due to the strong damping effect of ion-neutral collisions.

The CR residence time and grammage are presented by solid lines in Fig. 5 on the left and on the right, respectively. We observe that the residence time increases linearly from 22 kyrs at 30 pc to 200 kyr at 110 pc. The evolution curve of Setup 1 is slightly lower than that obtained in the case of a CR propagation in a homogeneous WNM. The associated grammage is relatively weak close to the source (X ∼ 0.02 g cm⁻² at 30 pc) and increases linearly by following the residence time trend up to the WNM–CNM transition where it increases by a factor of 100 ∼ n_T, CNM/n_T, WNM to reach the order of 10 g cm⁻² at 110 pc from the source.

Fig. 5. Evolution of the CR grammage (right part) and residence time (left part) in Setup 1: WNM-CNM and Setup 2: WNM-CNM-DiM represented by a solid curve and a dashed curve, respectively. These setups are compared with the case of the expansion of a CRC in a homogeneous WNM represented by the dotted lines. Each setup phase configuration is represented by a color code: WNM (green), CNM (orange), and DiM (red).

4.3.2. Setup 2: WNM-CNM-DiM

As presented by the left plots of Fig. 4, the observed trends in Setup 2 are almost equivalent to those observed in Setup 1. The modifications associated to the presence of the DiM, although small, are visible in the right-hand plots. The CR transport properties are essentially modified by the presence of the CNM. The effect of the DiM is negligible (∼0.2% of relative CR pressure variation). The only impact can be seen in Fig. 5 where the grammage increases due to the density contrast with respect to Setup 1.

It is important to note again that we did not take into account adiabatic losses associated to the spatial dependence of the Alfvén velocity. This aspect is left out for a future, more detailed work which may reveal a noticeable effect of the presence of the DiM phase on the CR propagation properties in the ISM.

5. Discussion and conclusion

In this study, we extend the work initiated in Nava et al. (2016) aimed at investigating the propagation of CRs in the ambient ISM surrounding SNRs. Here, we focus on cold, partially ionized atomic and molecular phases. To this aim, using the CRC model we simultaneously solve two 1D transport equations in 1D: an equation over the CR pressure P_CR (see Eq. (14)) and an equation over the Alfvén slab waves energy density I(k) (see Eq. (15)). We restrict our analysis to a resonant interaction where the particle Larmor radius and the wave number verify r_l = k⁻¹ in the quasi-linear approximation of CR transport; hence the lack of strong feed-back over the background turbulence. We determined the characteristic confinement time of CRs depending on their energy and the ISM phases: WNM, CNM, DiM. We discuss three different escape models once the SNR enters in the radiative phase, either considering a continuous escape (model 𝒫), an escape at the end of the free expansion phase (model ℱ), or an escape in the late radiative phase once the SNR forward shocks reach a speed of ∼150 km s (model 𝒮). We performed CRC evolution simulations for monoenergetic CR populations from 10 GeV to 10 TeV in both ISM phases in order to probe the evolution of the CR pressure distribution P_CR(E, z, t); the CR pressure gradient distribution ∂P_CR(E, z, t)/∂z, which is analogous to the force exerted by CR fluid on the plasma; and the CR diffusion coefficient D(E, z, t), which is intimately linked to the rate of excited Alfvén modes. Using the numerical solutions obtained we derived the grammage around a single cloud either in a homogeneous ISM or an ISM composed of multiple phases.

5.1. Results

Our results show that the Alfvén turbulence generated by the streaming of CRs can have an important effect on CR acceleration and propagation over distances of between a few tens and one hundred parsecs.

The CR escape solutions are presented in Fig. 2 for the three phases under study. We find a similar trend for the WNM as already discussed by N16 for the warm ionized phase. We show that in the CNM, at 1 TeV and 3 TeV, CRs escape at the same time as 100 GeV at ∼1 kyr while at 300 GeV CRs escape earlier at ∼0.4 kyr at the same time as 10 and 30 TeV. In the DiM, the radiative phase begins at ∼0.3 kyr and at 10 GeV, 300 GeV, and 1 TeV, CRs escape during the radiative stage at ∼2 kyr, while 30 TeV CRs escape at ∼0.05 kyr. The larger escape times at high energy are explained by the drop of the ion-neutral collisions in the coupled interaction regime.

We find that a non-negligible quantity of Alfvén modes are generated in the CNM (D ∼ 5 × 10²⁷ cm²/s at 10 TeV) and in the DiM (D ≈ ×10²⁶ at 10 TeV) after the CRC half-life time. These modes contribute to an enhanced CR confinement with respect to the transport by background turbulence. The model ℱ systematically shows stronger turbulence generation and longer confinement times with respect to the two others models.

In the WNM we show that the grammage behavior follows the observationally measured behavior and can be parametrized as X = X_WNM(z, 10 GeV <  E₀ <  10TeV)(E/E₀)^−0.5. The findings concerning propagation regimes in the CNM and the DiM imply that the CR grammage around SNRs can be enhanced at very high energy in these phases. In the CNM, the grammage has two regimes: the first one being between 10 GeV and 3 TeV with X = X_CNM(z, 10 GeV <  E₀ <  3 TeV)(E/E₀)^−0.5. Above 3 TeV, the grammage increases. In the DiM, the grammage decreases more slowly than in the WNM from 300 GeV to 10 TeV.

We perform a study at 10 GeV with multiphase setups: WNM-CNM and WNM-CNM-DiM. We find that the CR grammage (τ_res) is controlled by the densest phase and can reach ∼10 g cm⁻² in the WNM-CNM case and up to ∼10² g cm⁻² in the WNM-CNM-DiM case. The CR pressure evolution is however controlled by the CNM phase.

As stated above, our solutions at high energy are at the limit of validity of the quasi-linear theory. A more complete modeling including a 2D description is necessary for more accurate estimations of the grammage. We can also argue that if the grammage is expected to be smaller as extrapolated from [B/C] measurements, it reflects an average of the ensemble of CR pathways in the ISM. Partially ionized phases like the CNM and the DiM (and denser phases) do not cover a large fraction of the galactic volume, meaning that in total they do not significantly contribute to the average grammage.

5.2. Discussion and perspectives

We show that CRs accelerated at SNR shocks can have an important impact on the plasma turbulence properties around SNRs. Slightly unexpectedly, self-excited Alfvén modes can provide some confinement of CRs around the sources in the dense cold ISM environments beyond 1 TeV due to a drop of the ion-neutral collision damping in the coupled regime. At lower energies, ion-neutral collisions are in the decoupled regime and strongly damp self-generated waves. However, a substantial number of waves can still be generated, especially if CRs are released at the end of the free expansion phase as in model ℱ. In that case, diffusion can be suppressed by two to three orders of magnitude over distances of a few tens of parsecs over a rime period of a few thousand years. These waves are driven by strong CR gradients which in turn can have some impact on ISM dynamics (Commerçon et al. 2019).

However, several assumptions made in our model limit its implications. First, we only consider the resonant streaming instability and slab Alfvén modes. We do not investigate other types of instabilities like the nonresonant streaming instability uncovered by Bell (2004). Our model is only 1D and is therefore restricted to scales smaller than the coherence length of the background turbulence. One-dimensional geometry may also be problematic to investigate the inhomogeneous ISM, as is likely the case for atomic and molecular phases. The magnetic field lines can be mixed with turbulent motions in the super-Alfvénic regime which may occur for instance in the DiM phase, or possibly also for denser phases. A cylindrical 2D approach would be valuable to account for the geometry of the molecular clouds as well as the compression of the magnetic field lines (Chandran 2000). Two-dimensional calculations can also account for perpendicular CR diffusion. Another limit of our model is that it is restricted to medium where the speed of the scattering center remains approximately uniform. All these aspects, namely the geometry and adiabatic losses, will be accounted for in a future improved modeling.

Acknowledging these limitations, we argue that this simple model indicates that CR leakage from SNRs can induce a suppression of the diffusion with respect to what we could expect from the transport in the background turbulence. This effect should be taken into account in various aspects of the modeling of the ISM around CR sources.

Acknowledgments

The authors thank S.Gabici, L. Nava, S. Recchia for fruitful discussions and for their comments on the manuscript.

References

Appendix A: Cosmic ray propagation in the different ISM phases for the three models

A.1. Results for the WNM

Fig. A.1. Time evolution of a CRC of initial radius a in the WNM phase for model 𝒫. The results at 10 GeV, 100 GeV, 1 TeV, and 10 TeV are shown from the left to the right and the top to the bottom. Each panel presents three plots and shows CR pressure evolution as a function of space for three times (top), the CR pressure gradient evolution (middle), and diffusion coefficients (bottom). Blue, green, and red refer to t1/2/4, t1/2, and 4t1/2, respectively. The initial CR pressure distribution is represented by a thin black dashed line. The same typography is used for the background diffusion coefficient. Numerical solutions are represented with solid lines while test-particle solutions are presented by dashed lines.

Fig. A.2. Time evolution of a CRC of initial radius a in the WNM phase for model ℱ. Light gray plots show solutions differing from model 𝒫. See Fig. A.1 for more details.

Fig. A.3. Time evolution of a CR cloud of initial radius a in the WNM phase of the interstellar medium for the model 𝒮. Light gray plots show solutions differing from model 𝒫. See Fig. A.1 for more details.

A.2. Results for the CNM

Fig. A.4. Time evolution of a CR cloud of initial radius a in the CNM phase of the interstellar medium for model 𝒫. See Fig. A.1 for more details. In this medium we add two curves in brown and pink corresponding to solutions at 50 and 103t1/2.

Fig. A.5. Time evolution of a CR cloud of initial radius a in the CNM phase of the interstellar medium for model ℱ. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.

Fig. A.6. Time evolution of a CR cloud of initial radius a in the CNM phase of the interstellar medium for model 𝒮. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.

A.3. Results for the DiM medium

Fig. A.7. Time evolution of a CR cloud of initial radius a in the DiM phase of the interstellar medium for model 𝒫. See Fig. A.1 for more details. In this medium we add two curves in brown and pink corresponding to solutions at 50 and 103t1/2.

Fig. A.8. Time evolution of a CR cloud of initial radius a in the DiM phase of the interstellar medium for model ℱ. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.

Fig. A.9. Time evolution of a CR cloud of initial radius a in the DiM phase of the interstellar medium for model 𝒮. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.

All Tables

All Figures

Fig. 1. Ion-neutral damping rate (${\mathrm{\Gamma }}_{\mathrm{d}}^{\mathrm{in}}$, black curves), Farmer and Goldreich damping rate (${\mathrm{\Gamma }}_{\mathrm{d}}^{\mathrm{FG}}$, red curves), and Lazarian damping rate (${\mathrm{\Gamma }}_{\mathrm{d}}^{\mathrm{L}}$, blue curves) as a function of resonant CR energy. The shaded region refers to the range of energy where the real pulsation of Alfvén waves is not defined, which is why the turbulent damping, which depends on the Alfvén, is not defined in this scale range.

In the text

Fig. 2. Properties of the expansion of the CRC in the ISM. Left: evolution of the half-life time of the CRC as a function of its initial radius a in parsec units. Green, blue, red, and black solid lines refer to particle energies at 10, 102, 103, and 104 GeV. The black dashed line refers to the time evolution of the radius of the SNR front shock in the considered phase: WNM for the top figure, CNM for the middle one, and DiM for the bottom one. These plots have been obtained for Γ = 2.2. Right: evolution of the half-life time of the CRC derived from the left plots as a function of CR energy. The blue solid line refers to a CR energy spectral index of 2.2 while the red dashed line refers to a spectral index of 2.0. The gray zone corresponds to the no-wave propagation zone; it does not lead to any (a, t1/2) solution.

In the text

Fig. 3. Grammage (right vertical axis) and residence time (left vertical axis) of CRs as a function of their energy in three different phases of the ISM : WNM (top), CNM (middle), and DiM (bottom). Solid lines show the results obtained at different distances z* from the center of the source: green (100 pc), blue (50 pc), red (30 pc), and black (10 pc). The shaded region marks the no-wave propagation zone. Oscillations at high energy in the CNM and DiM are due to the limited energy resolution used to reconstruct the different curves.

In the text

Fig. 4. Evolution of the CR pressure (top figures) and their associated diffusion coefficient (bottom panels) in space at different times of the simulation based on the CRC half-life time for 10 GeV CRs in the WNM: 10.52 kyr (blue), 42.06 kyr (green), and 168.25 kyr (red). The simulation setups are represented by the different line styles: solid line for the WNM-CNM, dashed line for the WNM-CNM-DiM, and dotted line for the homogeneous WNM. The vertical black lines represent the initial size of the CRC. The phase decomposition for each setup is represented between the figures by the colored lines: green for WNM, orange for CNM, and red for DiM. Top-left figure: relative evolution of the CR pressure for WNM-CNM and WNM-CNM-DiM setups compared to the CR pressure evolution in the case of propagation in a homogeneous WNM phase. Top-right panel: relative CR pressure evolution for the WNM-CNM-DiM (PCR, 2) and the WNM-CNM (PCR, 1) setups, respectively. Bottom-left panel: relative evolution of the CR diffusion coefficients compared to the evolution of the CR diffusion coefficient in the case of propagation in a homogeneous WNM phase. Bottom-right panel: normalized relative evolution of the CR diffusion coefficient around the CNM and CNM-DiM phases. We recall that $D_{\text{CR}}^0={10}^{28}{\text{cm}}^2{\text{s}}^{-1}$.

In the text

Fig. 5. Evolution of the CR grammage (right part) and residence time (left part) in Setup 1: WNM-CNM and Setup 2: WNM-CNM-DiM represented by a solid curve and a dashed curve, respectively. These setups are compared with the case of the expansion of a CRC in a homogeneous WNM represented by the dotted lines. Each setup phase configuration is represented by a color code: WNM (green), CNM (orange), and DiM (red).

In the text

Fig. A.1. Time evolution of a CRC of initial radius a in the WNM phase for model 𝒫. The results at 10 GeV, 100 GeV, 1 TeV, and 10 TeV are shown from the left to the right and the top to the bottom. Each panel presents three plots and shows CR pressure evolution as a function of space for three times (top), the CR pressure gradient evolution (middle), and diffusion coefficients (bottom). Blue, green, and red refer to t1/2/4, t1/2, and 4t1/2, respectively. The initial CR pressure distribution is represented by a thin black dashed line. The same typography is used for the background diffusion coefficient. Numerical solutions are represented with solid lines while test-particle solutions are presented by dashed lines.

In the text

	Fig. A.2. Time evolution of a CRC of initial radius a in the WNM phase for model ℱ. Light gray plots show solutions differing from model 𝒫. See Fig. A.1 for more details.
In the text

	Fig. A.3. Time evolution of a CR cloud of initial radius a in the WNM phase of the interstellar medium for the model 𝒮. Light gray plots show solutions differing from model 𝒫. See Fig. A.1 for more details.
In the text

	Fig. A.4. Time evolution of a CR cloud of initial radius a in the CNM phase of the interstellar medium for model 𝒫. See Fig. A.1 for more details. In this medium we add two curves in brown and pink corresponding to solutions at 50 and 103t1/2.
In the text

	Fig. A.5. Time evolution of a CR cloud of initial radius a in the CNM phase of the interstellar medium for model ℱ. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.
In the text

	Fig. A.6. Time evolution of a CR cloud of initial radius a in the CNM phase of the interstellar medium for model 𝒮. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.
In the text

	Fig. A.7. Time evolution of a CR cloud of initial radius a in the DiM phase of the interstellar medium for model 𝒫. See Fig. A.1 for more details. In this medium we add two curves in brown and pink corresponding to solutions at 50 and 103t1/2.
In the text

	Fig. A.8. Time evolution of a CR cloud of initial radius a in the DiM phase of the interstellar medium for model ℱ. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.
In the text

	Fig. A.9. Time evolution of a CR cloud of initial radius a in the DiM phase of the interstellar medium for model 𝒮. Light gray plots refer to behaviors differing from model 𝒫. See Fig. A.1 for more details.
In the text