Issue 
A&A
Volume 633, January 2020



Article Number  A72  
Number of page(s)  25  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201936166  
Published online  14 January 2020 
Nonlinear diffusion of cosmic rays escaping from supernova remnants: Cold partially neutral atomic and molecular phases
^{1}
Laboratoire Univers et Particules de Montpellier (LUPM) Université Montpellier, CNRS/IN2P3, CC72, place Eugène Bataillon, 34095 Montpellier Cedex 5, France
email: Loann.Brahimi@umontpellier.fr
^{2}
Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radiowave Propagation, 108840 Troitsk, Moscow, Russia
Received:
24
June
2019
Accepted:
6
September
2019
Aims. We aim to elucidate cosmic ray (CR) propagation in the weakly ionized environments of supernova remnants (SNRs) basing our analysis on the cosmic ray cloud (CRC) model.
Methods. We solved two transport equations simultaneously: one for the CR pressure and one for the Alfvén wave energy density where CRs are initially confined in the SNR shock. Cosmic rays trigger a streaming instability and produce slabtype resonant Alfvén modes. The selfgenerated turbulence is damped by ionneutral collisions and by noncorrelated interaction with Alfvén modes generated at large scales.
Results. We show that CRs leaking in cold dense phases such as those found in cold neutral medium (CNM) and diffuse molecular medium (DiM) can still be confined over distances of a few tens of parsecs from the CRC center for a few thousand years. At 10 TeV, CR diffusion can be suppressed by two or three orders of magnitude. This effect results from a reduced ionneutral collision damping in the decoupled regime. We calculate the grammage of CRs in these environments. We find that in both single and multiphase setups at 10 GeV, CNM and DiM media can produce grammage in the range 10–20 g cm^{−2} in the CNM and DiM phases. At 10 TeV, because of nonlinear propagation the grammage increases to values in the range 0.5–20 g cm^{−2} in these two phases. We also present preliminary calculations in inhomogeneous interstellar medium combining two or three different phases where we obtain the same trends.
Key words: instabilities / turbulence / cosmic rays
© L. Brahimi et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 nonnegligible 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 largescale plasma motions and/or by the waves that they selfgenerate 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 gammaray observations, such as for example the recent H.E.S.S. Tcherenkov telescope observations from the forward shock of the SNR RX J 1713.73946 (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, lowenergy CRs constitute an important source of ionization (Phan et al. 2018)^{1}, 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)^{2}. We assume in this study that the amplitude of the selfgenerated waves is weak and treat the problem according to the quasilinear theory of CR transport (Schlickeiser 2002). In partially ionized phases, the waves triggered by CRs are damped by ionneutral collisions and by nonlinear interactions with the largescale 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 overdensity in the neutral/molecular ISM? What is the spacetimedependent rate of selfgenerated waves around a SNR accelerator? How can we evaluate the impact of CR selfconfinement 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 wavedamping 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 twofluid 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.
Fiducial physical quantities for the WNM, CNM, and DiM phases adopted in this work.
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^{3}) with a high equilibrium kinetic temperature (6 000 − 10 000 K) whereas CNM is denser (n_{T} = 20 − 50 cm^{−3}) 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^{−2} 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 largescale 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_{2} 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_{2}). As enough interstellar radiation is present to photodissociate 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^{−3}); 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 multiphase ISM, which can require multidimensional simulations. This aspect is discussed in Sect. 5 and merits future dedicated study. However, in Sect. 4 we consider an example of multiphase ISM in the context of the calculation of the CR grammage.
2.1.2. Magnetohydrodynamic wave damping in partially ionized media
Ionneutral damping. Alfvén waves are supported by the motion of ions. In a twofluid model energy transfer rate operates between neutrals and ions which produces a damping of the Alfvén waves. In the WNM phase, the ionneutral collision rate ν_{in} expresses the energy transfer rate from ions to neutrals and is given in Eq. (19) of N16 by
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^{−3}, 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 ionneutral collision rate is expressed by (Jean et al. 2009)
The neutralion collision rate ν_{ni} expresses the energy transfer from neutrals to ions; it is linked with ν_{in} by the relation (Xu et al. 2015)
where χ = (m_{n}/m_{i})(X^{−1} − 1) = ρ_{n}/ρ_{i} is the ratio between the mass density of neutrals and ions. Ionion ν_{ii} and neutralneutral ν_{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 twofluid 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_{0} (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)^{2})^{−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_{0}. We observe two asymptotic behaviors. At low CR energy (E ≪ eB_{0}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
Fig. 1. Ionneutral damping rate (, black curves), Farmer and Goldreich damping rate (, red curves), and Lazarian damping rate (, 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 . At high CR energy (E ≫ eB_{0}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
where ξ_{n} = ρ_{n}/(ρ_{n} + ρ_{i}) is the neutral fraction in the gas, and the Alfvén speed is given by . Finally, in the energy range [eB_{0}V_{Ai}/ν_{in}, eB_{0}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 largescale 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 CRgenerated 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 threewave 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 subAlfvé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 subAlfvenic 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 the turbulence becomes strong and is described . 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 subAlfvenic, strong turbulence regime, (M_{A} < 1 and ) the damping term is given by
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 subAlfvenic, weak turbulence regime (M_{A} < 1 and ) the damping term is given by
for . Alternatively, if k^{−1} > LM_{A}, then it is given by
and
In the superAlfvenic turbulence regime (), the turbulence at scales larger than 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 . The damping rate is given by
if , and
if .
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 transAlfvé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, ionneutral 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 transAlfvé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 1100 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 transAlfvénic regime, the turbulent damping rate is given by
The difference between the FG04 model and the L16 model appears at
corresponding to the lowenergy 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 ionneutral 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). Highenergy 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 selfconsistent 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_{0}). 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_{0}). 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 slabtype Alfvén waves through the resonant streaming instability. The generation of selfgenerated turbulence produces a reduction of the CR diffusion coefficient. Finally, the model is developed in the quasilinear approximation, so it is valid only in the case of welldeveloped, 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
for the CR pressure P_{CR}, and
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 spacedependent 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 quasilinear theory. Any modelling of the backreaction of the selfgenerated 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_{0} = D_{ISM}(E/10 GeV)^{0.5} and has been fixed to D_{ISM} = 10^{28} cm^{2} s^{−1} 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:
where corresponds to the magnetic energy of the medium. The second RHS term in Eq. (15) describes the damping of the waves due to ionneutral 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Γ_{d}I_{0} where I_{0} is the background turbulence level. Hereafter, the latter is assumed to be identical whatever the ISM phase. A detailed investigation of the largescaleinjected 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
where f(E) is the CR energy space phase function defined by
where R_{s} is the shock radius and the CR energy spectrum dN/dE is given by
if Γ ≠ 2 and
if Γ = 2. Here, W_{CR} = 10^{50} erg is the amount of energy of the SN explosion imparted into the CR pressure, and E_{min} = 2 × 10^{−4} GeV and E_{max} = 2 × 10^{4} 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 . 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_{0}. Following Malkov et al. (2013), we introduce the partial pressure Π given by
where 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^{2}/D corresponds to the CR diffusion timescale and τ_{damp} ≈ [2Γ_{d}]^{−1} 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
where D_{ISM} = D_{B}/I_{ISM} and
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
In the above framework, Malkov et al. (2013) defined the CRC halflife time 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
where is the initial CR pressure in the cloud. To each CR energy corresponds a particular halflife time. In order to determine at which moment and at which SNR radius a given population of CR escapes, N16 proposes to intersect the halflife time evolution function resulting from the Eq. (25) 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 corecollapse 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,
SedovTaylor (ST),
pressuredriven snowplow (PDS),
momentumconserving snowplow (MCS),
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.
The SNR radius depends on the SN mechanical energy E_{SN} in units of 10^{51} ergs, the ejecta mass M_{ej,⊙} in units of solar mass, the ambient density n_{T} in units of cm^{−3}, 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 halflife 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, 10^{2}, 10^{3}, and 10^{4} 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 halflife 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 nowave propagation zone; it does not lead to any (a, t_{1/2}) solution. 
Warm neutral medium. In the WNM phase, we recover the results obtained in N16 (see Fig. 2 upper row). We observe that lowenergy CRs escape late, at 10 GeV R_{esc} ≈ 25 pc and , 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^{3}, 2 × 10^{3}, 0.5 × 10^{3} ≫ 1) at 10, 100, and 1000 GeV respectively meaning that the selfgenerated 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^{3}, 4 × 10^{3}, 25, 4 ≫ 1 or > 1 at 10, 100, 10^{3} and 10^{4} 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 ionneutral 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 selfgenerated 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 halflife 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^{3}, 16, and 9 at 10, 10^{3}, and 10^{4} GeV, respectively.
It is also possible to study the CRC expansion properties by observing the cloud halflife time evolution as a function of the CR energy (see left parts of Fig. 2). In the WNM, the halflife time decreases linearly with the energy of CRs. Lowenergy CRs are efficiently confined while highenergy 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 selfgenerated turbulence connected with the drop of the damping rate. At very high energies, the halflife time decreases again because CRs are injected with an energy spectrum with powerlaw 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 halflife time decreases at low energy because of the decreasing value of the pressure but this effect is not strong (see the righthand 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 selfgenerated turbulence can be seen to depend on the wave growth rate and the wavedamping rate imposed by the density and the ionization rate of the external medium. Let us briefly describe them.
Different selfgenerated turbulence regimes as functions of the CR streaming instability growth rate and the Alfven waves damping rate by the medium.
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 selfgenerated 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^{−1} 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, ionneutral 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 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 𝒫.
Values of the halflife time and their associated escaping radii for the three models in this study at low energy.
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^{−1}. 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 raypropagation 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 timedependent transition between the CRC and the CR background. We note that in order to avoid any nonphysical pressure gradients at small escape times () 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.11, A.22, A.33 for the WNM, A.41, A.52, A.63 for the CNM, and A.71, A.82, A.93 for the DiM and for the three different escape models (model 𝒫: 1, model ℱ: 2 and model 𝒮: 3). The colorcode 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 , , and , respectively, as in N16. We added when necessary two more solutions in brown and pink at and , 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 testparticle (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 ( 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^{−32} erg/cm^{4}) is reached at but is not strong enough to induce selfgenerated turbulence at rates much in excess of the background turbulence. This solution closely resembles regime 4, described above. The width of selfgenerated 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^{−30} erg/cm^{4}. 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 and we observe that diffusion coefficients are close to the interstellar one. The width of generated turbulence is about ∼40 pc. Because the level of selfgenerated 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 . At this time, the value of the CR pressure gradient is about ∼5 × 10^{−31} erg cm^{−4} leading to a diffusion coefficient lower than 10^{28} cm^{2}/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 selfgenerated turbulence is stronger as showed in Fig. A.2. At t_{esc}/4, the diffusion coefficient is about ∼5 × 10^{26} cm^{2} s^{−1} and increases gradually to 2 × 10^{27} cm^{2} s^{−1} 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. Testparticle 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^{−28} erg cm^{−4} but then decreases rapidly to 10^{−30} erg cm^{4} at t_{1/2} and the selfgenerated turbulence level drops accordingly. At 4t_{1/2}, the selfgenerated 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^{−28} erg cm^{−4} and then slightly decreases until t_{1/2} to finally drop below 10^{−30} erg cm^{−4} at 4t_{1/2}. Diffusion is substantially suppressed (D ∼ 10^{26} − 10^{27} cm^{2} s^{−1}) 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^{26} cm^{2} s^{−1}). 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^{−27} erg cm^{−4}, and the diffusion coefficient is ∼10^{27} cm^{2} s^{−1}. The pressure gradient then decreases gradually during the simulation while the diffusion coefficient decreases gradually from ∼10^{27} cm^{2} s^{−1} to 10^{26} cm^{2} s^{−1} 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} ∼ 0.8. In order to prevent the level of selfgenerated turbulence from overtaking the quasilinear 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
where D_{⊥} = I^{2}D_{∥} and the spatial step perpendicular to the mean magnetic field direction is , 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 quasilinear 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^{27} erg cm^{−4}, and D_{min} ≈ 5 × 10^{24} cm^{2} s^{−1}) 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^{26} cm^{2} 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 selfgenerated 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^{25} − 2 × 10^{26} cm^{2} s^{−1}) but the selfgenerated 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 selfgenerated turbulence levels are high (∂P_{CR}/∂z_{max} ≈ 10^{−28} erg cm^{−4} at t_{1/2}/4 and t_{1/2}) which is confirmed by the values of the diffusion coefficients (D_{min} ≈ 2 × 10^{26} cm^{2} s^{−1} for t_{1/2}/4 and t_{1/2}). The turbulence level decreases very slowly at 4t_{1/2} : (∂P_{CR}/∂z)_{max} ≈ 10^{−29} erg cm^{−4} and D_{min} ≈ 10^{27} cm^{2} s^{−1}. 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^{−27} erg cm^{−4} and D_{min} ≈ 3 × 10^{26} − ×10^{26} cm^{2} s^{−1}. 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^{3}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^{−26} erg cm^{−4} and D_{min} ∼ 10^{25} cm^{2} s^{−1} 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^{−27} erg cm^{−4} and remains constant throughout the simulation while the turbulence rate increases with time (at t_{esc}/4, D_{min} ≈ 2 × 10^{25} cm^{2} s^{−1} while at 4t_{esc} this value becomes lower than 10^{25} cm^{2} s^{−1}) 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^{26} cm^{2} s^{−1}). 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^{−1}. 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 quasilinear 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 selfgenerated 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_{p}n_{T}cτ_{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 selfgenerated 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. Singlephase studies
In our semicylindric 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):
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 nowave 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^{−2} (X = 4 × 10^{−2} g cm^{−2} and X = 2 × 10^{−2} g cm^{−2} , respectively) at 10 GeV and decreases linearly to reach X = 10^{−2} g cm^{−2} (X = 2 × 10^{−3} g cm^{−2} and X = 1 × 10^{−3} g cm^{−2} 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 ionneutral 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^{−2}. 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^{−2} (X = 3 g/cm^{2} and X = 1.5 g cm^{−2} 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^{−2} (X = 2 g cm^{−2} and X = 1 g cm^{−2}, 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^{26} cm^{2} s^{−1}) until we reach the halflife 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 quasilinear 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^{−2} (X = 20 g cm^{−2} and X = 4 g/cm^{−2}, 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^{−2} (X = 1.5 g cm^{−2} and X = 0.2 g cm^{−2}, respectively).
In the DiM, we observe that the grammage evolution with the CR energy is modified compared to the lowenergy 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. Singlephase study
In the WNM, the grammage does not exceed 0.1 g cm^{−2} 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^{−2} for the CNM and ∼20 g cm^{−2} 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^{−2} and 5 g cm^{−2} at 50 pc from the source in the CNM and DiM, respectively. Here again, due to the small volumefilling 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 quasilinear theory. Finally, the results are obtained for a homogeneous ISM. More realistic multiphase 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 protonproton (pp) collisions and energy losses. The typical loss time by pp interaction is
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. Multiplephase 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 selfgenerated 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 ionneutral 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 wavecoupled evolution. At high energies, ionneutral 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 multiplephase case, we propose two setups from the leastdense 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 spacedependent 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 spacetime 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, WNMCNM: solid line, WNMCNMDiM: 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 halflife 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 WNMCNM, dashed line for the WNMCNMDiM, 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. Topleft figure: relative evolution of the CR pressure for WNMCNM and WNMCNMDiM setups compared to the CR pressure evolution in the case of propagation in a homogeneous WNM phase. Topright panel: relative CR pressure evolution for the WNMCNMDiM (P_{CR, 2}) and the WNMCNM (P_{CR, 1}) setups, respectively. Bottomleft 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. Bottomright panel: normalized relative evolution of the CR diffusion coefficient around the CNM and CNMDiM phases. We recall that . 
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: WNMCNM
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 α_{1} decreases and increases suddenly at the level of the transition WNMCNM (50 pc). Here, α_{1} reaches ∼1.15 around 70 pc before decreasing again and tends to ∼0.9 after the transition CNMWNM. At 42.06 kyr, we observe that α_{1} 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 (α_{1} ∼ 0.9) and α_{1} increases again reaching its maximum (∼1.2) around 110–120 pc and finally decreases. At 168.25 kyr, α_{1} 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 β_{1} 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 β_{1} 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, β_{1} 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 righthand 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 selfgenerated 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 ionneutral 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^{−2} 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^{−2} at 110 pc from the source.
Fig. 5. Evolution of the CR grammage (right part) and residence time (left part) in Setup 1: WNMCNM and Setup 2: WNMCNMDiM 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: WNMCNMDiM
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 righthand 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^{−1} in the quasilinear approximation of CR transport; hence the lack of strong feedback 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 ionneutral collisions in the coupled interaction regime.
We find that a nonnegligible quantity of Alfvén modes are generated in the CNM (D ∼ 5 × 10^{27} cm^{2}/s at 10 TeV) and in the DiM (D ≈ ×10^{26} at 10 TeV) after the CRC halflife 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_{0} < 10TeV)(E/E_{0})^{−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_{0} < 3 TeV)(E/E_{0})^{−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: WNMCNM and WNMCNMDiM. We find that the CR grammage (τ_{res}) is controlled by the densest phase and can reach ∼10 g cm^{−2} in the WNMCNM case and up to ∼10^{2} g cm^{−2} in the WNMCNMDiM 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 quasilinear 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, selfexcited 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 ionneutral collision damping in the coupled regime. At lower energies, ionneutral collisions are in the decoupled regime and strongly damp selfgenerated 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. Onedimensional 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 superAlfvé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). Twodimensional 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.
Energetic particles with a kinetic energy in the range between MeV and 100 MeV for protons and 100eV and keV for electrons (Padovani et al. 2009).
We restrict our investigation to the case of the generation of resonant slab Alfvén waves, and defer to a future study the case of the generation of other types of kinetic or convective instabilities; see however Inoue (2019).
Acknowledgments
The authors thank S.Gabici, L. Nava, S. Recchia for fruitful discussions and for their comments on the manuscript.
References
 Adhikari, L., Zank, G. P., Bruno, R., et al. 2015, ApJ, 805, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2008, A&A, 481, 401 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20 [Google Scholar]
 Bell, A. R. 1978a, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 1978b, MNRAS, 182, 443 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 2004, MNRAS, 353, 550 [Google Scholar]
 Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhko, E. G. 1996, Astropart. Phys., 5, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhko, E. G., & Ellison, D. C. 1999, ApJ, 526, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Blasi, P., Amato, E., & Serpico, P. D. 2012, Phys. Rev. Lett., 109, 061101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Burlaga, L. F., Florinski, V., & Ness, N. F. 2015, ApJ, 804, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Chandran, B. D. G. 2000, ApJ, 529, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Marcowith, A., & Dubois, Y. 2019, A&A, 622, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 D’Angelo, M., Blasi, P., & Amato, E. 2016, Phys. Rev. D, 94, 083003 [NASA ADS] [CrossRef] [Google Scholar]
 Drury, L. O., Ellison, D. E., Aharonian, F. A., et al. 2001, Space Sci. Rev., 99, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Everett, J. E., & Zweibel, E. G. 2011, ApJ, 739, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Farmer, A. J., & Goldreich, P. 2004, ApJ, 604, 671 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrière, K. M. 2001, Rev. Mod. Phys., 73, 1031 [NASA ADS] [CrossRef] [Google Scholar]
 Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, BAAS, 1, 240 [NASA ADS] [Google Scholar]
 Gaensler, B. M., Madsen, G. J., Chatterjee, S., & Mao, S. A. 2008, PASA, 25, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Giacinti, G., & Kirk, J. G. 2017, ApJ, 835, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJ, 816, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Grenier, I. A., Black, J. H., & Strong, A. W. 2015, ARA&A, 53, 199 [NASA ADS] [CrossRef] [Google Scholar]
 H.E.S.S. Collaboration (Abdalla, H., et al.) 2018, A&A, 612, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haverkorn, M., Brown, J. C., Gaensler, B. M., & McClureGriffiths, N. M. 2008, ApJ, 680, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Heiles, C., & Troland, T. H. 2003, ApJ, 586, 1067 [NASA ADS] [CrossRef] [Google Scholar]
 Inoue, T. 2019, ApJ, 872, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Jean, P., Gillard, W., Marcowith, A., & Ferrière, K. 2009, A&A, 508, 1099 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014 [NASA ADS] [CrossRef] [Google Scholar]
 Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Lagage, P. O., & Cesarsky, C. J. 1983, A&A, 125, 249 [NASA ADS] [Google Scholar]
 Larson, R. B. 1979, MNRAS, 186, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Lazarian, A. 2016, ApJ, 833, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700 [NASA ADS] [CrossRef] [Google Scholar]
 Lazarian, A., Vishniac, E. T., & Cho, J. 2004, ApJ, 603, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J.J., Koo, B.C., Raymond, J., et al. 2007, ApJ, 659, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M.M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Malkov, M. A., Diamond, P. H., Sagdeev, R. Z., Aharonian, F. A., & Moskalenko, I. V. 2013, ApJ, 768, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Matthaeus, W. H., Oughton, S., Pontius, Jr., D. H., & Zhou, Y. 1994, J. Geophys. Res., 99, 19267 [NASA ADS] [CrossRef] [Google Scholar]
 McIvor, I. 1977, MNRAS, 178, 85 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Nava, L., Gabici, S., Marcowith, A., Morlino, G., & Ptuskin, V. S. 2016, MNRAS, 461, 3552 [NASA ADS] [CrossRef] [Google Scholar]
 Nava, L., Recchia, S., Gabici, S., et al. 2019, MNRAS, 484, 2684 [NASA ADS] [CrossRef] [Google Scholar]
 Neufeld, D. A., Wolfire, M. G., & Schilke, P. 2005, ApJ, 628, 260 [NASA ADS] [CrossRef] [Google Scholar]
 Ng, C. S., & Bhattacharjee, A. 1997, Phys. Plasmas, 4, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Phan, V. H. M., Morlino, G., & Gabici, S. 2018, MNRAS, 480, 5167 [NASA ADS] [Google Scholar]
 Ptuskin, V. S., Zirakashvili, V. N., & Plesser, A. A. 2008, Adv. Space Res., 42, 486 [NASA ADS] [CrossRef] [Google Scholar]
 Recchia, S., Blasi, P., & Morlino, G. 2017, MNRAS, 470, 865 [NASA ADS] [CrossRef] [Google Scholar]
 Reville, B., Kirk, J. G., & Duffy, P. 2009, ApJ, 694, 951 [NASA ADS] [CrossRef] [Google Scholar]
 Reynolds, R. J. 1989, ApJ, 345, 811 [NASA ADS] [CrossRef] [Google Scholar]
 Roy, N., Peedikakkandy, L., & Chengalur, J. N. 2008, MNRAS, 387, L18 [CrossRef] [Google Scholar]
 Schlickeiser, R. 2002, Cosmic Ray Astrophysics [CrossRef] [Google Scholar]
 Schure, K. M., & Bell, A. R. 2013, MNRAS, 435, 1174 [NASA ADS] [CrossRef] [Google Scholar]
 Shull, J. M., & McKee, C. F. 1979, ApJ, 227, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Skilling, J. 1971, ApJ, 170, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Snow, T. P., & McCall, B. J. 2006, ARA&A, 44, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Carbonell, M., & Ballester, J. L. 2013, ApJS, 209, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Sutherland, R. S., & Dopita, M. A. 2017, ApJS, 229, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2012, A&A, 541, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Truelove, J. K., & McKee, C. F. 1999, ApJS, 120, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Vaupré, S., HilyBlant, P., Ceccarelli, C., et al. 2014, A&A, 568, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wentzel, D. G. 1968, ApJ, 152, 987 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, S., Lazarian, A., & Yan, H. 2015, ApJ, 810, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, S., Yan, H., & Lazarian, A. 2016, ApJ, 826, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., & Lazarian, A. 2004, ApJ, 614, 757 [NASA ADS] [CrossRef] [Google Scholar]
 Zank, G. P., Dosch, A., Hunana, P., et al. 2012, ApJ, 745, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Zank, G. P., Matthaeus, W. H., & Smith, C. W. 1996, J. Geophys. Res., 101, 17093 [NASA ADS] [CrossRef] [Google Scholar]
 Zirakashvili, V. N., & Ptuskin, V. S. 2012, Astropart. Phys., 39, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Zirakashvili, V. N., & Ptuskin, V. S. 2018, Astron. Lett., 44, 769 [NASA ADS] [CrossRef] [Google Scholar]
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 t_{1/2}/4, t_{1/2}, and 4t_{1/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 testparticle 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 10^{3}t_{1/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 10^{3}t_{1/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
Fiducial physical quantities for the WNM, CNM, and DiM phases adopted in this work.
Different selfgenerated turbulence regimes as functions of the CR streaming instability growth rate and the Alfven waves damping rate by the medium.
Values of the halflife time and their associated escaping radii for the three models in this study at low energy.
All Figures
Fig. 1. Ionneutral damping rate (, black curves), Farmer and Goldreich damping rate (, red curves), and Lazarian damping rate (, 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 halflife 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, 10^{2}, 10^{3}, and 10^{4} 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 halflife 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 nowave propagation zone; it does not lead to any (a, t_{1/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 nowave 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 halflife 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 WNMCNM, dashed line for the WNMCNMDiM, 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. Topleft figure: relative evolution of the CR pressure for WNMCNM and WNMCNMDiM setups compared to the CR pressure evolution in the case of propagation in a homogeneous WNM phase. Topright panel: relative CR pressure evolution for the WNMCNMDiM (P_{CR, 2}) and the WNMCNM (P_{CR, 1}) setups, respectively. Bottomleft 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. Bottomright panel: normalized relative evolution of the CR diffusion coefficient around the CNM and CNMDiM phases. We recall that . 

In the text 
Fig. 5. Evolution of the CR grammage (right part) and residence time (left part) in Setup 1: WNMCNM and Setup 2: WNMCNMDiM 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 t_{1/2}/4, t_{1/2}, and 4t_{1/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 testparticle 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 10^{3}t_{1/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 10^{3}t_{1/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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.