© ESO 2019

1. Introduction

Cosmic rays (CRs) are an important component of the interstellar medium (ISM); their pressure is close to equipartition with magnetic, gas, and radiation pressures, all on the order of 1 eV cm⁻³ in the local ISM (Grenier et al. 2015). The details of CR propagation in the ISM are still not well understood. CRs propagate partly by a succession of random walks and their mean free path between two scatterings is poorly known. Therefore, the process is controlled by the properties of local magnetic turbulence, which is poorly modelled. Furthermore, while pervading the ISM CRs can generate their own turbulence. This is one of the main ways they back-react over ISM turbulence and its dynamics (Kulsrud & Pearce 1969).

CRs have a multi-scale impact over galactic ISM dynamics (Grenier et al. 2015): they drive large-scale winds (Hanasz et al. 2013; Booth et al. 2013; Salem & Bryan 2014; Girichidis et al. 2016; Wiener et al. 2017), generate large-scale magnetic fields through the Parker instability (Hanasz et al. 2009; Butsky et al. 2017), inject turbulent magnetic motions at intermediate scales (Rogachevskii et al. 2012), contribute to ionise, diffuse, and heat molecular clouds or protostellar environments (Padovani et al. 2009, 2018; Rodgers-Lee et al. 2017), and finally they produce radio elements by spallation reaction in young stellar systems (Ceccarelli et al. 2014). There, CRs as charged fast moving particles have an impact on dust charge (Ivlev et al. 2015).

ISM turbulence is a key mechanism for models of star formation rate and stellar initial mass function (e.g. Hennebelle & Chabrier 2008, 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012). ISM turbulence is characterised by its statistical properties; in particular, as turbulent motions in the ISM are often supersonic, they can be characterised by the density variance-Mach number relation (Vazquez-Semadeni 1994; Krumholz & McKee 2005; Kritsuk et al. 2007). In studies to date, several physical effects have been included to investigate this relation: effects due to different turbulent forcing geometries (compressive versus solenoidal modes, Federrath et al. 2010), the impact of magnetic fields (Lemaster & Stone 2008; Molina et al. 2012), and thermodynamic properties of the turbulence in non-isothermal turbulent models (e.g. Nolan et al. 2015). It has been found that the turbulent forcing geometries and the magnetic fields affect the width of the log-normal distribution that fits the probability density function. In addition to the above-mentioned effects, CRs exert a force over the plasma along their pressure gradients. Hence, this force can also modify ISM dynamics and statistics. These effects should be expected to be particularly strong close to CR sources like supernova remnants or massive star clusters, over spatial scales and timescales that are strongly dependent on the ambient ISM phases surrounding the sources (Nava et al. 2016).

One of the goals of this study is to test the impact of anisotropic diffusion of CRs over the gas probability density function (PDF). It should be noticed that at smaller galactic scales corresponding to molecular clouds, little is known about the potential effect of CR gradients on the collapse of molecular clouds (see however Everett & Zweibel 2011), especially again if the cloud is near a CR source. This aspect will be treated in a forthcoming and associated study. In this study we first addressthe conditions in which the force due to a CR gradient can modify ISM turbulent statistics. We compute the response of a bi-stable fluid to different CR diffusivity laws. An associated study will investigate how CRs can modify this diffusivity through the generation of magnetic fluctuations and consider the effect of CR sources.

The lay-out of the paper is as follows. In Sect. 2, we discuss the necessary conditions on the CR spatial diffusion coefficient to produce a modification of ISM dynamics through the parallel CR gradient. In Sect. 3 we detail the simulation designed in this study. The main results of our parametric numerical study are presented in Sect. 4. We discuss our results and their limitations in Sect. 5, and we propose a list of different mechanisms that could lead to variations ofthe CR diffusion coefficient in the ISM. We conclude and present some perspectives on this study in Sect. 6.

2. Qualitative analysis of the properties of the turbulent ISM including cosmic rays

2.1. Turbulence and magnetic fields in the ISM

In this section, we look at the conditions under which CRs can be trapped within the different phases of the ISM. We focus in this short analytical study on the densest phases of the ISM: the warm neutral medium (WNM) and the cold neutral medium (CNM), essentially composed of neutral atomic hydrogen (HI). Dense ISM is also composed of molecular clouds, essentially composed of molecular hydrogen (H₂), which we also discuss qualitatively, but which we do not investigate using numerical experiments in the present study. Numerical models of molecular cloud formation and evolution will be the subject of a future study. This analysis is qualitative and is not aimed at covering the wide range of properties of the different gas phases in the ISM.

2.1.1. Turbulence

Turbulence is observed in the ISM at various scales and appears to be universal (Larson 1981; Heyer & Brunt 2004). Recent observations by Burkhart et al. (2015) report on transonic turbulence from HI emission PDF, which traces both the CNM and the WNM, while the HI absorption lines from the CNM show supersonic Mach numbers with a median value of 4.

Empiric relations between the gaseous velocity dispersion σ and the region size L are observed in the ISM, which are suggestive of a turbulent cascade (Larson 1979, 1981; Miville-Deschênes et al. 2017). The turbulent velocity dispersion increases with the size as $$\sigma \simeq {\sigma }_{\text{1} \text{pc}}{L}^q,$$(1)

where σ_1 pc is the velocity dispersion at a scale of 1 pc. For a Kolmogorov law, we have q = 1∕3.

In HI clouds, made of a mixture of WNM and CNM, observations report values of q ≃ 0.35 with σ_1 pc ≃ 1−1.5 km s⁻¹ (Larson 1979; Roy et al. 2008), which is close to the expected value for a Kolmogorov cascade. In the densest part of CNM, the molecular clouds, it is currently accepted that the gas follows scaling relations, named the Larson relations (Larson 1981), which relate the non-thermal velocity fluctuations to the size as well as to the mass of molecular clouds. Following the scaling given in Hennebelle & Falgarone (2012), the Larson relations read $$\begin{array}{lll}\sigma \hfill & \simeq \hfill & 1 \text{km} {\text{s}}^{-1}{\left(\frac{L}{\text{1} \text{pc}}\right)}^{0.5},\hfill \\ n\hfill & \simeq \hfill & 3000 {\text{cm}}^{-3}{\left(\frac{L}{\text{1} \text{pc}}\right)}^{-0.7},\hfill \end{array}$$

where n is the gas number density, and L is the size of the region in which the non-thermal velocity and the density are measured. It is generally assumed that this length scale represents the molecular cloud size. In this case, it follows that $$\sigma \simeq 1 \text{km} {\text{s}}^{-1}{\left(\frac{n}{3000 {\text{cm}}^{-3}}\right)}^{-5/7}.$$(4)

These scaling relations indicate that velocity fluctuations are reduced as the gas becomes denser, until self-gravity takes over to initiate the gravitational protostellar collapse and to re-accelerate the gas on much smaller scales. We do not consider molecular cloud and collapsing dense core scales in our numerical work. Hereafter, we assume that non-thermal motions are due to the turbulence only and consider that the non-thermal velocity is the rms velocity dispersion v_rms.

It is worth noting that the exponent values q ≃ 0.35−0.5 in the velocity-size relation are commonly retrieved in numerical works (e.g. Audit & Hennebelle 2010; Saury et al. 2014). In the following, we assume that the velocity dispersion and the size scale as $$v_{\text{rms}}\simeq 1 \text{km} {\text{s}}^{-1}{\left(\frac{L}{\text{1} \text{pc}}\right)}^q,$$(5)

with q varying in the different phases of the ISM.

2.1.2. Magnetic fields

Magnetic fields are observed at all scales in the ISM and show typical variations as a function of the density. Similarly to the turbulence, magnetic fields induce processes that are highly anisotropic. Magnetic fields are not amplified if the motion occurs along the field lines. In contrast, if a cloud is contracting, the anisotropic magnetic force will channel the gas to collapse first along the mean magnetic field direction, and then perpendicular to it. For a spherical contraction, the conservation of magnetic flux BR² for an enclosed mass M ∝ nR³ implies that B ∝ n^2∕3. The mass-to-flux conservation gives B ∝ nh, where h is the thickness of the cloud along the field lines. Equilibrium along the field lines implies that $c_{\text{s}}^2~\varphi$, with c_s the sound speed and ϕ the gravitational potential. Finally, the Poisson equation leads to ϕ ∝ c_sh². We thus get B ∝ n^1∕2 for an isothermal medium.

Measurements of the magnetic intensity via the Zeeman effect report on two different regimes for the magnetic intensity-density scaling (e.g. Crutcher 2012). At low density n < 300 cm⁻³, the magnetic intensity is insensitive to the density variations, with B_0,max ≃ 10 μG. While for densities n > 300 cm⁻³, that is within molecular clouds, the magnetic intensity scales with the density with an upper value given by B_max ∝ n^0.65 (which would correspond to a spherical contraction).

We can thus infer a simple relation that gives the magnetic intensity as a function of the density $$B\le 10 \mu \text{G} {\left[1+\left(\frac{n}{300 {\text{cm}}^{-3}}\right)\right]}^{0.65}.$$(6)

This relation is purely empirical and shows the most probable maximum values of the magnetic field amplitude (Crutcher 2012). It does not reflect the variations of the magnetic amplitude within the ISM, but it is useful to illustrate the different regimes that are expected in the ISM magnetohydrodynamics (MHD) turbulence.

2.1.3. MHD turbulence

We can derive two regimes for the MHD turbulence, which depends on the Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$, defined as ${\mathcal{M}}_{\text{A}}=v_{\text{rms}}/v_{\text{A}}$, with $v_{\text{A}}=B/{(4\pi {\mu }_{\text{gas}}{m}_{\text{H}}n)}^{0.5}$ the Alfvén speed and μ_gas the gas mean molecular weight. By combining the above-mentioned scaling relations for B and v_rms at the densities n < 100 cm⁻³ we interest us in this study, the Alfvénic Mach number can be estimated as $${\mathcal{M}}_{\text{A}}\simeq 0.9{\left(\frac{L}{\text{1} \text{pc}}\right)}^{0.35}{\left(\frac{n}{100 {\text{cm}}^{-3}}\right)}^{0.5},$$(7)

where we have taken a constant magnetic field amplitude of 6 μG according to the results of Heiles & Troland (2005) for the CNM. We further assume that this value is also valid for the WNM. The above relation only applies for the CNM and WNM scaling, that is, for q = 0.35.

Figure 1 illustrates the value of the Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$ as function of the gas density and the length scale from Eq. (7). This figure is only shown as an illustrative example of the range of Alfvénic Mach numbers that are observed (see Fig. 4 for instance). The MHD turbulence within the ISM switches between the regime of super- and sub-Alfvénic turbulence if the flow satisfies the above-mentioned scaling relations. At a gas density of ≤ 10 cm⁻³, the flow is always sub-Alfvénic, while at larger densities it may be super-Alfvénic. CRs propagate along magnetic field lines, and the topology of the lines thus affects the global propagation speed. In the super-Alfvénic regime, the magnetic field lines are tangled and the CRs take some time to escape from turbulent eddies. This reduces the effectivediffusion coefficient. In addition, the turbulence tends to isotropize the diffusion process so that the effective parallel and perpendicular coefficients are equal. On the contrary, CRs can propagate rapidly along smooth magnetic field lines in the sub-Alfvénic regime. The parallel diffusion coefficient is large, but the perpendicular diffusion coefficient is typically orders of magnitude smaller (e.g. Yan & Lazarian 2004).

Fig. 1. Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$ as a function of the gas density n and of the length scale L following Eq. (7). This plot is valid only if one assumes the scaling relation and magnetic field strength observed in the CNM and WNM.

2.2. The critical CR diffusion coefficient

We consider a blob of turbulent interstellar gas of length L, with a turbulence characterised by the root-mean-square velocity v_rms. The turbulent time is defined as the time at which a turbulent fluctuation travels through the blob, $$t_{\text{turb}}\equiv \frac{L}{v_{\text{rms}}}.$$(8)

We define the diffusion time t_diff as $$t_{\text{diff}}\equiv \frac{L^2}{D},$$(9)

where D is the CR diffusion coefficient.

It is straightforward to estimate a critical diffusion coefficient D_crit that corresponds to the regime where CR diffusion is not sufficient to prevent energy fluctuations. CRs then get trapped, which can lead to a significant gradient of CRs and, thus, a dynamical effect of the CR pressure onto the flow.

The critical diffusion is given at t_turb = t_diff $$D_{\text{crit}}\equiv L{v}_{\text{rms}}.$$(10)

Using the scaling relation in Eq. (5), the critical diffusion coefficient is $$D_{\text{crit}}\simeq 3.1\times {10}^{23} {\text{cm}}^2 {\text{s}}^{-1}{\left(\frac{L}{\text{1} \text{pc}}\right)}^{q+1}.$$(11)

The above calculation of D_crit is an order of estimate as the scaling relation for v_rms shows some dispersion and as we neglect the dimensionality of the diffusion process. Nevertheless, the value of q does not vary a lot between the WNM and CNM and the molecular clouds so that our prediction is robust and does not significantly change in the different phases of the ISM. For a cloud 10 pc in size, which is typical of the length used in studies of molecular clouds turbulence (e.g. Federrath & Klessen 2012), the critical diffusion coefficient is D_crit ≃ 10²⁵ cm² s⁻¹. A usual parametrisation of the CR parallel diffusion coefficient derived from CR data fitting is $$D_{\text{ISM}}~{10}^{28-29} {\text{cm}}^2 {\text{s}}^{-1}{\left(\frac{E}{1 \text{Gev}}\right)}^s,$$(12)

with s ≃ 0.3–0.6 (Strong et al. 2007). The value of D_crit we derived is small compared to the canonical value of D_ISM ≃ 10²⁸ cm² s⁻¹ at a kinetic energy of ≃1 GeV (e.g. Ptuskin et al. 2006) in which we are interested in this study. CR energy spectra peak around ≃1 GeV in Galactic environments (Grenier et al. 2015), which thus corresponds to the CR particles that contribute most to the CR pressure. We discuss in Sect. 5 some possible mechanisms that can modify the CR diffusion coefficient in the different phases of the ISM.

3. Numerical setup

3.1. Magnetohydrodynamics with cosmic-ray diffusion

Our numerical model integrates the equation of ideal MHD for a fluid mixture made of gas and CRs, including anisotropic CR diffusion. We describe in the following the details of our numerical tool.

We use the adaptive-mesh-refinement code RAMSES (Teyssier 2002), which is based on a finite-volume, second-order Godunov scheme and a constrained transport algorithm for ideal MHD (Fromang et al. 2006; Teyssier et al. 2006). We use the HLLD (Harten-Lax-Van Leer-Discontinuities) Riemann solver to compute the MHD flux, as well as the electric field to integrate the induction equation (Miyoshi & Kusano 2005). CRs are treated as a fluid following the implementation of anisotropic diffusion described by Dubois & Commerçon (2016). It consists in an advection-diffusion equation of the CR energy density, in which the CRs are advected with the fluid, contribute to the total energy, and propagate along the magnetic field lines. The full CR-MHD equations read $$\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot \left[\rho u\right]=0,\\ \frac{\partial \rho u}{\partial t}+\nabla \cdot \left[\rho uu+P_{\text{tot}}\mathbb{I}-\frac{BB}{4\pi }\right]=\rho f,\\ \frac{\partial E_{\text{tot}}}{\partial t}+\nabla \cdot \left[u\left(E_{\text{tot}}+P_{\text{tot}}\right)-\frac{\left(u\cdot B\right)B}{4\pi }\right]\\ =-P_{\text{CR}}\nabla \cdot u+\nabla \cdot \mathbb{D}\nabla E_{\text{CR}}\\ +\rho f\cdot u-\mathcal{L}(\rho ,T),\\ \frac{\partial E_{\text{CR}}}{\partial t}+\nabla \cdot \left[u{E}_{\text{CR}}\right]=-P_{\text{CR}}\nabla \cdot u+\nabla \cdot \mathbb{D}\nabla E_{\text{CR}},\\ \frac{\partial B}{\partial t}-\nabla \times \left[u\times B\right]=0,\end{array}$$

where ρ is the gas density, u is the fluid velocity, $\mathbb{I}$ is the identity matrix, P_tot the total pressure, that is, the sum of the gas thermal pressure P = (γ − 1)E and the CR pressure P_CR = (γ_CR − 1)E_CR where γ and γ_CR are, respectively, the gas and CR adiabatic index (and E and E_CR the gas and CR internal energy densities), and the magnetic pressure P_mag = B²∕(8π). The magnetic field vectors is B, E_tot is the total energy E_tot = E + 0.5ρu² + B²∕(8π) + E_CR. The anisotropic diffusion tensor of cosmic rays is $\mathbb{D}$, given by $$\mathbb{D}=D_{\perp }{\delta }_{ij}-(D_{\perp }-D_{\parallel })bb,$$(18)

where b = B∕||B||, D_∥ and D_⊥ are the diffusion coefficients, respectively, along and across the magnetic field. Their amplitude remains poorly constrained in the turbulent ISM. For numerical convenience, we set D_⊥ = 0.01D_∥ in order to avoid poor convergence of our implicit solver used for the anisotropic diffusion (Dubois & Commerçon 2016). For the thermal budget, we include heating and cooling processes $\mathcal{L}(\rho ,T)$ following Audit & Hennebelle (2005) to be able to describe the CNM and the WNM phases of the ISM, as well as the thermal instability (Field 1965). The impact of gas heating due to CR streaming remains to be quantified with respect to the classical heating and cooling processes included in this work. We also neglect the role of CR streaming and CR loss processes (see for instance Ruszkowski et al. 2017; Pfrommer et al. 2017, for alternative implementations of CR that account for some of these effects). The CR streaming velocity magnitude corresponds to the velocity of the Alfvén waves v_A. We thus expect that CR streaming may have an impact in the sub-Alfvénic regimes. In addition, CR streaming leads to a loss term in the CR energy equation, and thus to a heating term for the gas thermal energy. These two effects will be investigated in future works. Lastly, CR cooling rates in the ionised ISM are often accounted for (e.g. Guo & Oh 2008; Booth et al. 2013), but the cooling rates in the neutral ISM remain to be quantified.

The specific force f represents the force that is applied to drive the turbulence. This force field is generated in the Fourier space and its mode is given by a stochastic process based on the Ornstein-Uhlenbeck process (Eswaran & Pope 1988; Schmidt et al. 2006). For more details on the turbulence forcing module, we refer readers to the work of Schmidt et al. (2009) and Federrath et al. (2010) on which our implementation is based. We set the forcing to a mixture of compressive and solenoidal modes, with a compressive force power equal to 1∕3 of the total forcing power. We use a unique autocorrelation timescale T ≃ 0.6 Myr. We neglect self-gravity in this study.

The system is closed using the perfect gas relation, P = ρk_BT∕(μ_gasm_H), with μ_gas = 1.4. We set the gas adiabatic index to γ = 5∕3 and the CR fluid adiabatic index to γ_CR = 4∕3. Lastly, we define the modified sonic Mach number for the gas and CR mixture as $${\mathcal{M}}_{\text{s},\text{CR}}\equiv \frac{v}{c_{\text{s},\text{CR}}},$$(19)

where $$c_{\text{s},\text{CR}}\equiv \sqrt{\gamma \frac{P}{\rho }+{\gamma }_{\text{CR}}\frac{P_{\text{CR}}}{\rho }}$$(20)

is the modified sound speed which accounts for the CR pressure.

3.2. Initial conditions

Our initial set-up is very similar to that of Seifried et al. (2011) and Saury et al. (2014). We set a box of size L_box = 40 pc filled with gas that follows the observed scaling relations. The initial state of the gas corresponds to a thermally unstable state. The initial temperature is T₀ ≃ 4460 K, density n₀ = 2 cm⁻³, gas thermal pressure P ≃ 1.23 × 10⁻¹² dyne cm⁻². The initial sound speed is c_s ≃ 5.2 km s⁻¹. The initial cosmic-ray pressure is set by its ratio with the gas thermal pressure, which we define as ζ ≡ P_CR∕P. If not stated otherwise, we set ζ = 0.1. The initial magnetic field is set to B₀ = 1 μG, which gives an initial Alfvénic speed v_A ≃ 1.3 km s⁻¹ and a plasma beta β = P∕P_mag ≃ 30. Magnetic field amplification by turbulent dynamo is unlikely to take place in our models since we start with a magnetic field amplitude that corresponds to the values expected at saturation (e.g. Ferrière & Schmitt 2000; Rieder & Teyssier 2017).

The turbulence is driven with various mode wave numbers k that correspond to a fraction of the box size. For k_turb = 2, the turbulence is driven at a scale corresponding to half the box size. We explore four different driving scales, k_turb = 2, k_turb = 4, k_turb = 8, and k_turb = 16. This setup mimics the effect of different driving mechanisms such as supernovae at large driving scales and protostellar jets at the smallest scales. We use a parabolic function between k ∈ [k_turb − 1, k_turb + 1] to choose the wave numbers for the forcing module. The amplitude of the driving is kept constant between all the different scales, which means that the resulting v_rms varies as a function of k_turb.

3.3. Parametric study

We performed a set of 50 simulations to explore the initial conditions of the parameter space. We used a uniform fiducial resolution of 128³. Saury et al. (2014) performed a resolution study (from 128³ to 512³) and have shown that even though the scales to describe the cold CNM and the molecular clouds are barely resolved with 128³, this resolution is sufficient for a box of 40 pc to perform a parameter study. Going towards higher resolution will change the mass fraction in the cold gas, as well as its structure, but it will not affect the qualitative picture. In this paper, we are interested on the global effect of CRs on the two phases of ISM dynamics, and we do not detail the structure of the cold gas quantitatively. We present a resolution study in Appendix A that validates our choice of a 128³ resolution for a qualitative study.

We explored the parameter space with the five following parameters: D₀ the amplitude of the parallel diffusion coefficient, k_turb the wave number for turbulence driving, v_rms the amplitude of the turbulence, ζ the initial ratio between the CR and the gas pressure, and ${\mathcal{M}}_{\text{A}}$ the initial Alfvénic Mach number.

Table 1 summarises the different runs adopted for our study (the simulations performed in Appendix A are not reported here). The turbulence injection scale varies from 20 pc for k_turb = 2 to 2.5 pc for k_turb = 16. The amplitude of the rms velocity also varies with the driving scale (high driving wave numbers have weaker rms velocity). Most models are initially in a super-Alfvénic regime, while the initial turbulent Mach numbers vary from slightly supersonic with k_turb = 2 to subsonic with k_turb = 16 ($\mathcal{M}\simeq 0.4$).

If not stated otherwise, we analysed our simulation results at times beyond two turbulent crossing times t_cross = L_box∕v_rms. We determined the value of the turbulent velocity dispersion v_rms by running MHD models without CRsand averaging over ≃40 Myr. Changing the time of analysis or taking a unique absolute time does not change the qualitative results provided that the simulations are analysed once the turbulent motions have settled in a statistically stationary regime.

Similar to Sect. 2.2, we can infer a critical diffusion coefficient by assuming that the diffusion timescale is comparable to the eddy turnover timescale, t_dyn = L_inj∕v_rms, at the injection scale L_inj ≡ L_box∕k_turb. The critical diffusion coefficient in our setup reads $$D_{\text{crit}}=L_{\text{inj}}{v}_{\text{rms}}.$$(21)

Since we use a unique value of 40 pc for the size of the computational domain L_box, the critical diffusion coefficient is $$D_{\text{crit}}\simeq 1.2\times {10}^{25}{k}_{\text{turb}}^{-1}\left(\frac{v_{\text{rms}}}{1 \text{km} {\text{s}}^{\text{-1}}}\right) {\text{cm}}^2 {\text{s}}^{-1}.$$(22)

We give in Table 1 the values of the critical diffusion coefficient for each setup, which ranges from 1.4 × 10²⁴ to 6.3 × 10²⁵ cm² s⁻¹. We varied D₀ from 10²² to 10²⁸ cm² s⁻¹ in order to cover the transition regime.

We can also estimate the value of the critical length below which CR are trapped as a function of D₀ and v_rms, that is, L_crit = D₀∕v_rms, $$L_{\text{crit}}\simeq 3.2\times {10}^5\left(\frac{D_0}{{\text{10}}^{\text{28}} {\text{cm}}^{\text{2}} {\text{s}}^{\text{-1}}}\right){\left(\frac{v_{\text{rms}}}{1 \text{km} {\text{s}}^{\text{-1}}}\right)}^{-1} \text{pc}.$$(23)

In the case k_turb = 2 and v_rms = 5.7 km s⁻¹, the critical length varies from about 5.7 kpc for D₀ = 10²⁸ cm² s⁻¹ to 5.6 × 10⁻³ pc for D₀ = 10²² cm² s⁻¹. The condition for CR trapping is thus simply L_crit < L_inj.

4. Results

4.1. A fiducial case, k_turb = 2

We first discuss the fiducial set of simulations with a driving scale corresponding to k_turb = 2 and an amplitude of the turbulence v_rms ≃ 5.7 km s⁻¹, which corresponds to a turbulent crossing time t_cross = 6.9 Myr. We explore five values for the diffusion coefficient: 0, 10²², 10²⁴, 10²⁶, 10²⁸ cm² s⁻¹. We also compare the results with a model without CR (labelled No CR). From Eq. (11), we expect the critical value for the diffusion coefficient to be D_crit ≃ 8 × 10²⁵ cm² s⁻¹.

Figure 2 shows the time evolution of the clumping factor of the gas density, CR energy density, andAlfvénic Mach number for the five above-mentioned simulations with various diffusion coefficients D₀ and the model without CRs. The clumping factor of a quantity X is defined as C(X) = ⟨X²⟩∕⟨X^⟩2. A clumping factor C(X) > 1 indicates that the field X has significant variations. The evolution of the clumping factors quickly settles to a quasi-stationary regime within about one turbulent crossing time (≃7 Myr) once the turbulent motions are established. First, we verified that all simulations show large variations in the density, which indicates that the gas has undergone the thermal instability to settle in the two-phases CNM-WNM medium. The clumping factor of the gas density shows two distinct regimes. In the case where the CR diffusion coefficient is large, D₀ ≥ 10²⁶ cm² s⁻¹, the clumping factor is large (>4) and the results are very similar to that of the case without CRs, with an amplitude of a factor > 10³ between the minimum and maximum values. On the contrary, with D₀ ≤ 10²⁴ cm² s⁻¹ the clumping factor is close to 1, which indicates that the density field is smoother. The CR energy density shows the opposite behaviour. In the case with D₀ ≥ 10²⁶ cm² s⁻¹, the CR energy density does not vary much, with less than a factor of four to five between the minimum and maximum values, and the clumping factor is very close to unity. For D₀ ≤ 10²⁴ cm² s⁻¹, the CR energy density shows a clumping factor larger than unity corresponding to large amplitude variations of a factor > 100 at a given snapshot, which corresponds roughly to the ones found in the density field. Last, the Alfvénic Mach number variations do not depend on the CR diffusion coefficient, or on the presence of CRs. In addition, we observe that the Alfvén speed is independentof the CR pressure, implying that the magnetic fields-density relation is also independent of the CR physics. Since our numerical setup uses the ideal MHD framework, it is natural that the density and the magnetic fields remain well coupled.

The time evolution of the CR energy density shows consistent results, with large variations occurring for D₀ ≤ 10²⁴ cm² s⁻¹. Consequently, the gas density variations are tempered by the extra pressure support provided by the gradients in the cosmic-ray pressure, and the maximum gas density is reduced by about one order of magnitude. In the following, we will refer to this regime as the “trapped CR regime”.

Figure 3 shows maps of the gas density, CR energy density, and Alfvénic Mach number at a time corresponding to 2t_cross for the fiducial runs with k_turb = 2 with CR diffusion. We clearly recover the two different behaviours for D₀ ≤ 10²⁴ cm² s⁻¹ and D₀ ≥ 10²⁶ cm² s⁻¹ in the density and CR energy density maps. The CR energy and the gas density fields share the same morphology when CRs are trapped, while the density behaves independently of the CR pressure when D₀ ≥ 10²⁶ cm² s⁻¹ (CR pressure gradients are shallow). As already mentioned, the Alfvénic Mach number shows variations of the same amplitude in all cases.

Figure 4 shows 2D-histograms of the CR energy density, the Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$, the ratio ${\mathcal{M}}_{\text{s},\text{CR}}/{\mathcal{M}}_{\text{s}}$, and ζ = P_CR∕P as a function of the density for all the cells in the computational domain at the same time as in Fig. 3. The colour coding corresponds to the temperature, the magnetic field magnitude, and the sonic Mach number ${\mathcal{M}}_{\text{s}}$ (for the two bottom rows), respectively. The correlation between the CR energy density and the gas density is obvious for D₀ ≤ 10²⁴ cm² s⁻¹ and the CR fluid behaves as an adiabatic fluid with P_CR ∝ n^γ_CR, while the CR energy density is insensitive to the density variations for a larger CR diffusion coefficient. The temperature varies from ≃ 10⁴ K in the low density regions to T < 100 K in the high density regions, which indicates that the thermal instability has operated in all models to settle in the expected regime where the cold neutral medium and the warm neutral medium coexist. We note that we neglect CR radiative losses and CR streaming heating, which may affect the thermal balance. These effects will be investigated in future work. As mentioned before, the Alfvénic Mach number and magnetic field magnitude do not depend on the value of the CR diffusion coefficient. The ratio ${\mathcal{M}}_{\text{s},\text{CR}}/{\mathcal{M}}_{\text{s}}$ shows that when CRs are trapped, the CR pressure can substantially modify the effective sound speed of the gas and CR fluid mixture. In the densest regions, the Mach number can be reduced by a factor of two. In the opposite case, the CR pressure has a moderate effect for D₀ ≥ 10²⁶ cm² s⁻¹ in the low density regions but no influence at high densities. The behaviour of the pressures ratio ζ shows also opposite trends. For D₀ ≤ 10²⁴ cm² s⁻¹, the CR pressure outmatches the thermal one with increasing density, being the greatest in the highly supersonic regions (${\mathcal{M}}_{\text{s}}>1$). The CNM and WNM are roughly in pressure equilibrium and then the temperature decreases with increasing density in the CNM, so that the gas effective polytropic index is γ_eff < 1. The CR energy density is strongly correlated with the density. The CR pressure thus grows faster with the density than the thermal pressure. For D₀ ≥ 10²⁶ cm² s⁻¹, ζ decreases with the density. Since the CR pressure is roughly uniform, ζ varies with the density as ζ ∝ 1∕n^γ_eff. At high density, the gas is almost isothermal and thus we get γ_eff ≃ 1.

The volume-weighted PDF of the gas density, mass-weighted PDF of the gas temperature, and mass-weighted PDF of the CR energy density are shown in Fig. 5. In the remainder of the paper, all the PDFs are time-averaged, computed once the calculations have reached 2t_cross. The time average is done over a timescale varying from one to two dynamical times. For comparison, we also show the two extreme runs with D₀ = 0 (no diffusion) and without CRs. There is again an obvious distinction depending on whether the CR diffusion coefficient ishigher or lower than the critical value. The gas density and temperature PDFs are wider in models with D₀ ≥ 10²⁶ cm² s⁻¹, suggesting that the gas is more efficiently compressed (or expanded) compared to the lower diffusion coefficient values. From first principles, one would have expected that the gas gets more compressed since the mixture of gas and CRs has a global adiabatic index 4∕3 < γ < 5∕3 (see Eq. (20)). Nevertheless, the gas compression is primarily dependent on the thermal budget resulting from the heating and cooling processes included. As we just discussed, the effective polytropic index of the gas is even lower than 4∕3 (e.g. Wolfire et al. 2003). In the trapped CR regime, the CR pressure dominates so that the effectivepolytropic index of the mixture is close to the CR one and prevents strong compressionwith respect to the case where thermal pressure dominates (D₀≥ 10²⁶ cm² s⁻¹). The temperature variations are reduced in the models with D₀ ≤ 10²⁴ cm² s⁻¹. In the case without CRs, the temperature PDF barely shows the two expected peaks for the CNM and the WNM mixture at T ≃60 K (CNM phase) and T ≃ 5000 K (WNM phase). The maximum temperature increases by adiabatic heating to reach a higher value than the initial one, while the various thermal processes included in the heating and cooling function lead to an efficient cooling in the high density regions. The CR energy density PDF is more sensitive to the diffusion coefficient than the gas density and temperature ones. The CR energy variation amplitudes are anti-correlated with those of the gas density and temperature. For D₀ = 10²⁸ cm² s⁻¹, the CR energy shows no variation and remains uniform at an energy density corresponding to the initial state. The CR and gas fluids are decoupled and CRs propagate freely in the gas. The CR energy density begins to show variations of less than a decade for D₀ = 10²⁶ cm² s⁻¹. CR pressure gradients thus develop, but have no effect on the gas flow as indicated in the density PDF. For D₀ = 10²⁴ cm² s⁻¹, the CR energy density PDF is much wider (≃2.5 decades at PDF = 10⁻³), so that CR pressure gradients develop. The peak of the PDF shifts towards higher values compared to the initial value. The trapped CRs are compressed and then diffuse slowly so that the peak moves towards higher CR energy densities. Interestingly, the CR energy density PDFs vary substantially between D₀ = 10²² cm² s⁻¹ and D₀ = 10²⁴ cm² s⁻¹ while this is not the case for the gas density and the temperature. For D₀ = 10²² cm² s⁻¹, the three PDFs match closely the ones of the model without CR diffusion, that is where the CRs are perfectly trapped with the fluid. In the opposite case, the gas density and temperature of the case where D₀ = 10²⁸ cm² s⁻¹ correspond to those of the run without CRs.

The picture we can draw from these fiducial runs is the following: for diffusion coefficients D₀ ≤ 10²⁴ cm² s⁻¹, CRs are trapped and can provide extra pressure support to counterbalance the increase in gas density due to adiabatic compression and/or shocks, and the cooling becomes less efficient. In the model with D₀ ≤ 10²⁴ cm² s⁻¹, the thermal instability is thus less efficient in segregating the gas in the mixture of WNM and CNM. As already mentioned, we neglect in this work CR streaming and radiative losses. Further work will investigate the impact of these effects, in particular in the sub-Alfvénic regions, on the qualitative picture we present here. We observe that the width of the gas density PDF increases with the CR diffusion coefficient. As mentioned in Sect. 1, the width of the gas density PDF can be related to the fluid sonic Mach number. In our model, the fluid is made of a mixture of gas and CRs and the extra CR pressure support increases the effective sound speed of the mixture. As a consequence, the effective Mach number decreases as shown in Fig. 4 and the PDF gets tighter.

We finally note that given our choice of parameters, the critical length for CR trapping was ≃ 57 pc (larger than the box size of 40 pc) for D₀ = 10²⁶ cm² s⁻¹, which is indeed verified by our numerical results since the results show no evidence of CR trapping with this diffusion coefficient.

Fig. 2. Time evolution for six simulations with kturb = 2 of the clumping factor of the gas density n (top panel), the cosmic-ray energy density ECR (middle panel), and the Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$ (bottom panel). The clumping factor C(X) of a quantity X is defined as C(X) = ⟨X2⟩∕⟨X⟩2.

Fig. 3. Fiducial runs with kturb = 2: maps of the gas density (top panels), cosmic-ray energy density (middle panels), and Alfvénic Mach number (bottom panels) in the plane z = 20 pc at a time corresponding to 2tcross. From left to right panels: CR diffusion coefficient is D0 = 1022, 1024, 1026, and 1028 cm2 s−1. The arrows on the density maps represent the velocity vectors (normalised to a value of 14 km s−1, see upper right panel), and the streamlines on the CR energy density maps represent the magnetic field lines, both in the xy-plane.

Fig. 4. Two-dimensional histograms of, from top to bottom panels, the CR energy density, the Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$, the ratio ${\mathcal{M}}_{\text{s},\text{CR}}/{\mathcal{M}}_{\text{s}}$, and ζ = PCR∕P as a function of the density for the fiducial runs kturb = 2 with CR diffusion at the same time as in Fig. 3. The colour coding indicates the temperature, the magnetic field strength, and the sonic Mach number ${\mathcal{M}}_{\text{s}}$ (two bottom rows), respectively.

4.2. Parameter study

In this section, we explore the parameter space by varying the CR initial pressure, the driving scale, the driving amplitude, and the initial magnetisation level.

4.2.1. Effect of the CR initial pressure

The first parameter we test is the initial ratio ζ between the CR pressure and the gas pressure. In the fiducial case, we used ζ = 0.1. We now explore three additional cases, with ζ = 1, 10, and 100. We used the same initial conditions and turbulence driving properties as in the fiducial runs (only the CR pressure increases). For each value of ζ, we performed two simulations with D₀ = 10²⁴ and D₀ = 10²⁶ cm² s⁻¹, which surround the critical value D_crit.

Figure 6 shows the PDFs of the gas density, the temperature, and the CR energy density for the fiducial cases and the six runs with ζ = 1, 10, and 100. For a diffusion coefficient D₀ = 10²⁶ cm² s⁻¹ and ζ ≤ 1, the gas density and temperature PDFs are similar in width. We note some differences in the maximum density (and thus minimum temperature), which is lower (respectively, higher) with ζ = 1. The CR energy density is shifted by about one decade towards a higher value with ζ = 1, which corresponds to the initial difference between the CR pressure. The width of the PDF is narrower with ζ = 1, which indicates that the CRs are less sensitive to the gas flow fluctuations if their pressure is larger. Nevertheless, the CR pressure has a slight effect on the gas flow as indicated by the maximum density, which is a little bit lower. The CR pressure is strong enough to provide an extra support. For ζ > 1 and D₀ = 10²⁶ cm² s⁻¹, the three PDFs become tighter as ζ increases, up to the point where the CR pressure prevents any compression of the gas and thus the development of the thermal instability. As the gas density PDF can be well approximated by a log-normal distribution, which width is determined by the sonic Mach number, increasing the CR pressure leads to a decrease of the modified sound speed c_s,CR and, hence, to a decrease in the width of the PDF. In the D₀ = 10²⁴ cm² s⁻¹ case, the difference between ζ = 0.1 and ζ ≥ 1 is striking. The larger CR pressure inhibits the development of the thermal instability and the gas remains in a state that is thermally unstable. CNM gas formation, where n > 10 cm⁻³, is strongly altered if CRs are in pressure equilibrium with the gas (roughly energy equipartition) and if they get trapped. In this case, the transition between the two regimes is much more abrupt. If the initial CR pressure exceeds the thermal one, the gas does not evolve much and the thermal instability does not develop even though it is initially in a state that is thermally unstable. The classical criterion for the thermal instability can roughly be can be summarised as (d P∕dρ) < 0 (Field 1965). Our results indicate that this criterion must be revised if one accounts for the CR pressure. Values of ζ ≥ 1 are to be expected in the ISM surrounding CR sources while CRs are freshly injected in. Our results show that CR sources can strongly modify the ambient ISM dynamics even if the CR diffusion coefficient is larger than D_crit. This aspect has to be accounted for in simulations including supernova (SN) feedback.

Fig. 5. PDF of the gas density (top panel), temperature (middle panel), and CR energy density (bottom panel) for the fiducial runs with kturb = 2 at the same time as in Fig. 3, for D0 = 1022 (purple), 1024 (blue), 1026 (red), and 1028 cm2 s−1 (black). Two models with D0 = 0 (dashed cyan) and without CRs (dashed grey) are also represented for comparison.

Fig. 6. Same as Fig. 5 for the models with kturb = 2 and D0 = 1024 (blue), D0 = 1026 cm2 s−1 (red). The solid lines represent the fiducial model with ζ = 0.1, the dashed line the models with ζ = 1, the dotted line ζ = 10, and the dashed-dotted line ζ = 100.

Fig. 7. Same as Fig. 5 for kturb = 4.

4.2.2. Effect of the driving scale

We present in this section the results obtained by running models with the same initial conditions as in the fiducial case, but where we vary the scale for turbulence driving. We explore three smaller driving scales, corresponding to k_turb = 4, 8, and 16. Theaim is to probe the effect of a smaller injection scale for the turbulence. We keep the driving amplitude of the acceleration constant between the different models so that the rms Mach number decreases with increasing k_turb. We thus indirectly test at the same time the effect of varying the turbulence amplitude.

From the measured value of v_rms reported in Table 1, the critical length for CR trapping is 8.2, 11, and 17.8 pc for k_turb = 4, 8, 16 respectively (and for D₀ = 10²⁴ cm² s⁻¹), while the turbulence injection scale varies from 10 to 2.5 pc. This means that for the k_turb = 8 and 16 models, the driving scale L_inj is less than the expected critical length L_crit for D₀ = 10²⁴ cm² s⁻¹. In addition, since the turbulent energy decreases with increasing k_turb, D₀ = 10²⁴ cm² s⁻¹ is closer to thecritical diffusion coefficient reported in Table 1 than in the fiducial runs.

Figures 7–9 show the gas density, temperature, and CR energy density PDFs for the models with varying driving scalek_turb (4, 8, and 16, respectively). In all cases, we recover the two regimes, but the transition shifts towards lower values of D₀ as k_turb increases, as expected in the qualitative analysis above. For D₀ ≥ 10²⁶ cm² s⁻¹, the gas always settles in the two phases mixture of CNM and WNM. The width and global shape of the gas density and temperature remain similar for k_turb = 4 and 8, though with gas that is denser and warmer in the WNM as k_turb increases (compression is simply less efficient). In the k_turb = 16 case, the density PDF is shallower at high density. Most of the gas remains at a low density because the turbulence is not strong enough to push the gas along the Hugoniot curve towards a higher density. Nevertheless, we observe that the PDF of the temperature does not reflect that of the density in the k_turb = 16 models, with an excess of low temperature gas. The spread in density at a given temperature is indeed larger in the cases where the turbulence is stronger as the density fluctuations may get destroyed faster, so that the thermal instability cannot develop quickly and cooling becomes less efficient. In the case of low turbulence, the gas remains predominantly in thermal equilibrium and the fraction of cold gas increases (e.g. Audit & Hennebelle 2005; Seifried et al. 2011; Gazol & Kim 2013). Overall, the PDF of the CR energy density is more sensitive to the driving scale. For D₀ ≥ 10²⁶ cm² s⁻¹, the width gets narrower as the turbulence level decreases. In the trapped CR regime, the width of the CR energy density PDFs scales with the gas density and thus gets also narrower at higher k_turb. Lastly, we do not observe any particular effect when the driving scale is less than the critical length for D₀ = 10²⁴ cm² s⁻¹. The transition to the trapped CR regime remains on the same parameter range.

Fig. 8. Same as Fig. 5 for kturb = 8.

4.2.3. Effect of the driving amplitude or sonic Mach number

The next experiment consists in increasing by a factor ≃ 2 the amplitude of the turbulence driving in the fiducial case (k_turb = 2, i.e. the driving scale remains the same contrary to the preceding section) to reach v_rms = 10.5 km s⁻¹. We performed four simulations with a CR diffusion coefficient ranging from D₀ = 10²² to D₀ = 10²⁸ cm² s⁻¹.

Figure 10 shows the gas density, temperature, and CR energy density PDFs with v_rms = 10.5 km s⁻¹. The PDFs of the fiducial case with v_rms = 5.7 km s⁻¹ are also shown for comparison. From the gas density PDF, the main effect of a stronger turbulence is to produce lower density regions for all values of D₀, as a result of larger expansion. The high gas densities are not very sensitive to the amplitude of the turbulence driving. As a consequence, the temperature is also not very sensitive to the amplitude of the turbulence. Finally, the CR energy density shows the same qualitative behaviour as in the fiducial case, with trapped CRs for D₀ ≤ 10²⁴ cm² s⁻¹. The PDFs are nevertheless different, with larger width and higher peak values with stronger turbulence. The peak is shifted toward higher energy density as a result of stronger compression by turbulent motions. The minimum value of the CR energy density decreases with increasing v_rms as a result of local expansions, which create low energy regions density that cannot be refilled with surrounding CRs since diffusion is less efficient.

Fig. 9. Same as Fig. 5 for kturb = 16.

Fig. 10. Same as Fig. 5 for kturb = 2, but different amplitude in the turbulence driving. The solid lines represent the fiducial case results with vrms = 5.7 km s−1 and the dashed lines the ones obtained with vrms = 10.5 km s−1.

4.2.4. Effect of the magnetic field amplitude or Alfvénic Mach number

In this last experiment, the dependency on the initial Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$ is explored by changing the value of the initial magnetic field amplitude. We compared the fiducial case results, with ${\mathcal{M}}_{\text{A}}=4.5$ (super-Alfvénic), with the ones obtained with ${\mathcal{M}}_{\text{A}}=20$ (super-Alfvénic, weak magnetic fields), ${\mathcal{M}}_{\text{A}}=1.2$ (trans-Alfvénic), and ${\mathcal{M}}_{\text{A}}=0.6$ (sub-Alfvénic). We recall that magnetic field amplitude increases with decreasing ${\mathcal{M}}_{\text{A}}$. We use two values for the CR diffusioncoefficient, D₀ = 10²⁴ and D₀ = 10²⁶ cm² s⁻¹.

Figure 11 shows the PDFs of the eight different runs. Again, the qualitative behaviour is not modified, with a transition to the regime of trapped CRs for D₀ ≤ 10²⁶ cm² s⁻¹. The gas density and temperature PDFs do not depend much on ${\mathcal{M}}_{\text{A}}$ for D₀ = 10²⁶ cm² s⁻¹. A noticeable feature is the width of the gas density PDF that is the narrowest with ${\mathcal{M}}_{\text{A}}=0.6$. This is the expected behaviour since stronger magnetic fields provide support that softens the turbulent motions (e.g. Molina et al. 2012). The CR energy density PDF gets wider as ${\mathcal{M}}_{\text{A}}$ increases. The magnetic field lines are indeed more tangled in the super-Alfvénic case, which reduces the effective diffusion coefficient. For D₀ = 10²⁴ cm² s⁻¹, we see clearly two different regimes as a function of ${\mathcal{M}}_{\text{A}}$, and the results are qualitatively different from the case with D₀ = 10²⁶ cm² s⁻¹. In the super-Alfvénic regime, the gas density and temperature PDFs are narrower, meaning that as the magnetisation increases the gas gets more compressed. As we just mentioned, the effective CR diffusion coefficient decreases in the super-Alfvénic regime. CRs are thus more efficiently trapped and CR pressure gradients get stronger. As a consequence, the CR pressure support increases and counterbalances the effect of the increase in magnetisation level. Clearly, the coupled effect of the magnetisation level and of CR propagation is non-linear and should be investigated in future works with more details partly based on some specific CR transport modelling (see Sect. 5.2.1).

Fig. 11. Same as Fig. 5 for kturb = 2, but different amplitudes of the initial magnetic field strength. The solid lines represent the fiducial case with ${\mathcal{M}}_{\text{A}}=4.5$, the dashed lines ${\mathcal{M}}_{\text{A}}=1.2$, the dotted lines ${\mathcal{M}}_{\text{A}}=0.6$, and the dashed-dotted lines ${\mathcal{M}}_{\text{A}}=20$. The blue colour represents results with D0 = 1024 cm2 s−1, the red D0 = 1026 cm2 s−1.

5. Discussion

5.1. Dependence on the parameters

We have seen in the previous section that the critical diffusion coefficient for CR trapping D_crit is on the order of 10²⁴ –10²⁵ cm² s⁻¹ and is in good agreement with the theoretical expectations. In addition, we have observed a strong effect of the initial pressure ratio ζ. For large initial CR pressure, ζ = 10, the gas and CR mixture behaves differently even if CRs are not trapped. The non-thermal support provided by the CR pressure may be very efficient at preventing the gas compression by turbulent motions so that the thermal instability is inhibited in CR-pressure-dominated flow. We show in the following discussion that CR sources and CR acceleration at small scales are possible mechanisms for local enhancements of the CR pressure, which might then change the gas dynamics in the environment of these sources.

In addition, the initial Alfvénic Mach number has an effect on the gas compression. In the case of initial super-Alfvénic Mach numbers, the magnetic field lines are tangled and the CR diffusion tends to be isotropic with a reduced effective coefficient. The maximum CR energy density is thus lower and the CRs’ pressure gradients are weaker than in the case of sub-Alfvénic turbulence. We note that we retrieve this behaviour when an isotropic diffusion coefficient is used (see Appendix B). We have seen in Sect. 2 and Fig. 1 that the ISM MHD turbulence can be either super- or sub-Alfvénic, and we expect to recover both the isotropic and anisotropic diffusion regimes. We also discuss in Sect. 5.2.1 some models for the diffusion coefficient depending on the large-scale Alfvénic Mach number.

5.2. Possible mechanisms for CR diffusion coefficient variation

At a first glance it might be hopeless to investigate ISM dynamics under the effect of CRs because the usual ISM diffusion coefficient is D_ISM ≫ D_crit. However, several physical processes can contribute to decrease the effective CR diffusion coefficient below D_ISM and possibly below D_crit in some specific ISM environments. We provide a qualitative discussion of these processes below but as this first study aims only at a parametric investigation of the effect of CR pressure over ISM dynamics, we postpone to an associated paper (Paper II) a detailed study of these specific diffusion regimes. We also caution the reader that some processes we invoke in the following may be at scales we do not consider in our previous numerical experiments.

5.2.1. Large-scale turbulence models

Large-scale motions in the ISM are sources of turbulence (Mac Low & Klessen 2004). Among the potential processes that are sources of turbulence we find magneto-rotational instabilities produced by differential galactic rotation motions, supernova explosions either isolated or through a cumulative effect in superbubbles, and gravitational instabilities induced by galactic arm motions. At smaller galactic scales other sources of turbulence injection can contribute, especially protostellar jets, stellar wind feedback, or ionisation front motions associated to H_II regions. All these processes force turbulent motions with different geometries, either solenoidal or compressible or a combination of both, at some preferred scale L_inj. Once fluid motions are forced at L_inj, energy cascades down to smaller scales. This cascade can be described by different models. As the ISM is magnetised, large-scale motions can be described approximately using a MHD turbulent model. At the scales under investigation, we have seen in Sect. 2.1.1 that turbulent motions are transonic in the WNM and supersonic in the CNM, so compressible. MHD compressible turbulence has been mainly studied by means of numerical simulations (Lithwick & Goldreich 2001; Cho et al. 2002).

In particular, Cho et al. (2002) isolate different regimes of MHD turbulence depending on the regime of Alfvénic Mach number. Based on this turbulence model, Yan & Lazarian (2008) propose a derivation of CR diffusion coefficients induced by these large-scale-injected turbulent motions. Let us recall their main results.

Super-Alfvénic regime

If $M_{\text{A}}>1$, then large-scale motions are hydrodynamical and follow a Kolmogorov law. At a scale ${\mathcal{l}}_{\text{A}}=L_{\text{inj}}/M_{\text{A}}^3$ the turbulence becomes trans-Alfvénic. At scales below ℓ_A the turbulence is MHD and can be decomposed into three MHD mode cascades: Alfvén and slow magnetosonic cascades are anisotropic and follow a Goldreich–Sridhar model (Goldreich & Sridhar 1995), while the fast magnetosonic cascade is isotropic and follows a Kraichnan model (Cho et al. 2002). CR diffusion in super-Alfvénic turbulence is isotropic, hence D = D_∥ = D_⊥ but depends on λ_∥, the CR mean free path parallel to the mean magnetic field. If λ_∥ > L_inj, the effective diffusion coefficient is D = ℓ_Av∕3 while if λ_∥ < L_inj, D = λ_∥v∕3, where v = β_cc is the CR velocity. The calculation of λ_∥ in partially ionised phases is complex and involves a calculation of the turbulence damping scale by ion-neutral collisions. Xu et al. (2016) find a characteristic U shape for λ_∥ with respect to the CR energy. It appears that in the WNM and CNM phases, GeV CRs have Larmor radii smaller than MHD turbulence damping scales and are in the regime λ_∥≫ L_inj. Alfvénic Mach numbers in these partially ionised phases are not well known but it is possible to define a critical Alfvénic Mach number using Eq. (11), $$M_{\text{A},\text{crit}}\equiv {\left(\frac{v{L}_{\text{inj}}}{3D_{\text{crit}}}\right)}^{1/3}\simeq 46 {\left(\frac{{\beta }_{\text{c}}{L}_{\text{inj}}}{1 \text{pc}}\right)}^{1/3}{\left(\frac{L}{1 \text{pc}}\right)}^{-\frac{(q+1)}{3}}.$$(24)

If $M_{\text{A}}\ge {\mathcal{M}}_{\text{A},\text{crit}}$, CRs can have dynamical effects. Using L_inj = L∕2 = 20 pc, q = 1/3 or 1/2 as in our fiducial simulation we find for 1 GeV CRs $M_{\text{A},\text{crit}}\simeq 23$ for q = 1/3 and $M_{\text{A},\text{crit}}\simeq 19$ for q = 1/2. From Fig. 1 we can deduce that this regime is marginally obtained for the gas density and scale length considered in this work. However, smaller critical Alfvénic Mach numbers are expected if the turbulence is injected at smaller scales (at higher k_turb).

Sub-Alfvénic regime

If $M_{\text{A}}<1$, magnetic fields dominate gas dynamics at the injection scale. The turbulence is weak and cascades in the direction perpendicular to the mean magnetic field (Galtier et al. 2000). Below a scale ${\mathcal{l}}_{\text{w}}=L_{\text{inj}}{M}_{\text{a}}^2$ the Goldreich–Sridhar scaling prevails again. Sub-Alfvénic CR diffusion is no longer isotropic. The perpendicular diffusion coefficient depends again on λ_∥ (Yan & Lazarian 2008). If λ_∥ > L_inj then $D_{\perp }=D_{\parallel }{M}_{\text{A}}^4{L}_{\text{inj}}/{\lambda }_{\parallel }$, while if λ_∥ < L_inj, $D_{\perp }=D_{\parallel }{M}_{\text{A}}^4$; in each case, D_⊥ ≪ D_∥. The parallel diffusion coefficient is D_∥ = λ_∥v∕3. In that case one can define a critical parallel mean free path λ_∥,crit $${\lambda }_{\parallel ,\text{crit}}\simeq 3.1\times {10}^{13} \text{cm} {\beta }_{\text{c}}^{\text{-1}}{\left(\frac{L}{\text{1} \text{pc}}\right)}^{\text{q+1}},$$(25)

which in our fiducial case using Eq. (11) is λ_∥,crit ≃ 4.9 × 10¹⁵ cm for q = 1/3 or λ_∥,crit ≃ 9.0 × 10¹⁵ cm for q = 1/2. Calculations by Xu et al. (2016) indicate that at GeV energies λ_∥ > L_inj so likely larger than λ_∥,crit. We do not expect any strong CR effect in this regime.

5.2.2. CR streaming and turbulence generation

Low energy CRs are injected at the end of a supernova remnant’s (SNR) lifetime likely when the forward shock becomes trans-sonic but possibly before when the SNR enters the radiative phase when shock acceleration stops. As most SNRs explode behind spiral arms and usually in groups (thus forming superbubbles), they produce a mean CR density (or pressure) gradient in the galaxy both in vertical and radial directions. CRs while propagating at the speed of light along the mean magnetic field can trigger MHD perturbations preferentially along the magnetic field lines through the so-called streaming instability (Skilling 1971). Perturbations generated by CRs are resonant and have typical wave numbers $k_{\text{res}}~r_{\text{g}}^{-1}$, where r_g = E∕eB is the CR gyroradius for a CR kinetic energy E and a charge e in a magnetic field of strength B. The streaming instability growth rate Γ_g is controlled by the CR density or pressure gradient along the field lines. In the ISM, we can expect typical CR pressure gradients on the order of $P_{{\text{CR}}^{\prime }}~P_{\text{CR}}/L_{\text{CR}}$, where L_CR is the CR gradient scale. A reference value is $P_{\text{CR},{\text{ref}}^{\prime }}\simeq 5.2\times {10}^{-33} \text{dyne} {\text{cm}}^{\text{-3}} P_{\text{CR},\text{eVcc}}/L_{\text{CR},\text{100} \text{pc}}$, where P_{CR, eVcc} is the CR pressure in units of 1 eV  cm⁻³ and L_CR is the CR gradient scale in units of 100 pc. In the Milky Way, one expects $P_{{\text{CR}}^{\prime }}\le P_{\text{CR},{\text{ref}}^{\prime }}$ as the typical CR variation scale lengths are larger than the galactic disc height.

CRs can inject magnetic perturbations at small scale lengths through the streaming instability. The level of the perturbations saturate due to a balance of the growth rate with a damping rate. Depending on the ISM phase under consideration the dominant damping process changes (Wiener et al. 2013; Nava et al. 2016). In partially ionised phases, the dominant damping process is due to collisions between ions and neutrals. At GeV energies, the CR self-generated waves are at too high frequencies to prevent ions from being decoupled from neutrals. The damping rate is Γ_in = ν_in∕2, where again ν_in is the ion-neutral collision frequency. We take ${\nu }_{\text{in}}\simeq 1.4 {10}^{-9}{n}_{\text{n}} {(T/100\text{K})}^{0.5} {\text{cm}}^{\text{3}} {\text{s}}^{\text{-1}}$ and ${\nu }_{\text{in}}\simeq 1.6 {10}^{-9}{n}_{\text{n}} {\text{cm}}^{\text{3}} {\text{s}}^{\text{-1}}$ with n_n the number density of neutrals, in the WNM and the CNM, respectively (Jean et al. 2009). In the decoupled ion-neutral regime the streaming instability growth rate (in s⁻¹) can be written as ${\Gamma }_{\text{g}}(E)\simeq V_{\text{A},\text{i}}{P}_{{\text{CR}}^{\prime }}(E)/(2I(k)W_{\text{B}})$ (Nava et al. 2016), where $V_{\text{A},\text{i}}=B/\sqrt{4\pi {\rho }_{\text{i}}}$ is the ion Alfvén velocity, W_B = B²∕8π is the ambient mean magnetic energy density, and $I(k)={(\delta B(k)/B)}^2$ is the relative amplitude of perturbations at a wave number $k=k_{\text{res}}~r_{\text{g}}^{-1}$ with respect to the ambient mean magnetic field. The condition Γ_g = Γ_in sets the value of I(k). If I(k) ≪ 1 then parallel and perpendicular CR diffusion coefficients are given by the CR transport quasi-linear theory (Schlickeiser 2002) by $D_{\parallel }=4/3\pi r_{\text{g}}v/I(k=r_{\text{g}}^{-1})$ and D_⊥ ≃ $D_{\parallel }I{(k)}^2$, and again v = β_cc is the particle speed.

The parallel diffusion coefficient is then $$D_{\parallel }=\frac{4}{3\pi } \frac{r_{\text{g}}{\beta }_{\text{c}}c{\nu }_{\text{in}}{W}_{\text{B}}}{V_{\text{A},\text{i}}{P}_{\text{CR}}{(E)}^{\prime }},$$(26)

or comparing it with D_crit given by Eq. (11), $$\frac{D_{\parallel }}{D_{\text{crit}}}\simeq 4.2\times {10}^3\frac{E_{\text{GeV}}{\beta }_{\text{c}}{\nu }_{\text{in},\text{-8}}{({\mu }_{\text{i}}{n}_{\text{i},\text{-2}})}^{1/2}}{P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}}{L}_{\text{pc}}^{-(1+q)}.$$(27)

We have expressed the CR energy E_GeV in GeV units, the ion density n_i in units of 10⁻² cm⁻³, the ion-neutral collision frequency in units of 10⁻⁸ s⁻¹, and the ion mean mass is μ_i.

For numerical application we choose the following parameters for WNM and CNM. In the WNM we use T ≃ 6000 K, μ_i = 1, n_i,-2 ≃ 0.3, and ν_{in, -8} ≃ 0.3. In the CNM we use μ_i = 12 but provide some estimates for μ_i = 1 as well, n_i,-2 ≃ 2 and ν_{in, -8} ≃ 5. In order to have a substantial effect of CR self-generated waves over ISM dynamics, that is to have D_∥ ~ D_crit, we need a ratio $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}>1$. If we use the conditions adopted in our fiducial case, then L = 40 pc and q = 1/3 or q = 1/2. In the WNM we find $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}\ge 4$ at a CR kinetic energy of 1 GeV in order to get D_∥ = D_crit for q = 1/3 and $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}\ge 2$ for q = 1/2. In the CNM we find $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}\ge 390$ still at 1 GeV in order to get D_∥ ≃ D_crit for q = 1/3 and $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}\ge 210$ for q = 1/2. If we assume μ_i = 1 we find $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}\ge 110$ and ≥ 60 for q = 1/3 and q = 1/2, respectively.

A CR source can perturb the ambient ISM because it can inject some CR overdensity (or overpressure) over scales that can be smaller than typically 100 pc. Hence, the associated CR pressure gradient can exceed $P_{\text{CR},{\text{ref}}^{\prime }}$ by a factor 10–30 in SNR or even more in star-forming or in starburst galaxies. For instance, Nava et al. (2016) consider thepropagation of an SNR in warm ISM phases (warm ionised and warm neutral). They investigate the propagation regime of CR release from the forward shock. In that case, the release operates upstream of the shock because high energy CRs have a mean free path larger than a fraction of the shock radius. While escaping ahead of the shock, CRs trigger resonant MHD waves by the streaming instability. The diffusion coefficient is again set bybalancing the wave growth rate and the dominant damping process, which is again due to ion-neutral collisions. The authors find that at 10 GeV, so for CRs released late on in the SNR lifetime, diffusion is suppressed¹ because of scattering off the self-generated waves down to a factor 100 at time ~ 10 000 yr after the explosion, while higher CR energies released earlier have a suppression of their diffusion coefficient by a factor 10–30. Calculationsin the CNM show similar behaviours although over timescales shorter by a factor of a few kilo-years (Brahimi et al., in prep.). However, as far as we know, no calculations have investigated the case of GeV CRs released at the end of shock acceleration or SNR lifetime. To be conservative, at 1 GeV, considering an overpressure of P_CR ≃ 10–30 eV  cm⁻³ in the ISM over a typical scale-length R ≃ 10 pc, typical of a fraction of the radius of a supernova remnant at the end of its lifetime in these phases, we find $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}~100\text{--}300,$ enough with the above requirements to obtain D_∥ ≃ D_crit.

However, it is mandatory to check that I(k) ≪ 1 as the streaming instability growth rate has been obtained in the quasi-linear limit. From the condition Γ_g = Γ_in we deduce a turbulence level of $$I(k)\simeq 2.9\times {10}^{-5}\left(\frac{(P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}){({\mu }_{\text{i}}{n}_{\text{i},\text{-2}})}^{-1/2}}{B_{\mu \text{G}}{\nu }_{\text{in},\text{-8}}}\right).$$(28)

In the conditions that prevail in partially ionised phases, using B ≃ 6 μG, the ratio $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}$ cannot exceed ~3 and ~ 300 in WNM and CNM, respectively, for the quasi-linear condition (I(k) ≤ 10⁻⁴) to apply. However, this simple calculation made through a rate equilibration likely overestimates the level of magnetic field perturbations produced so the previous limits are likely larger. Non-linear calculations as performed in Nava et al. (2016) tend to confirm this assertion. Another caveat of these calculations is that low energy sub-relativistic CRs also released at the end of the acceleration phase or SNR lifetime will ionise the surrounding ISM (Padovani et al. 2009). This effect has two implications. Firstly, as the ionisation fraction increases, the ion density also increases and D_∥ ∕D_ISM increases as well. This first effect has then a negative feedback on CR ISM dynamical impact. Secondly, as the ion density increases, the ion-neutral collision rate decreases and so the damping of self-generated wave also decreases. This second effect has then a positive feedback on the CR ISM dynamical impact. As all quantities above depend on $n_{\text{i}}^{1/2}$, whereas ν_in ∝ n_n, it is tempting to say that the second effect will take over the first and that the CR confinement and wave generation around SNR will be enhanced. However, an accurate calculation of the propagation of each populations (i.e. low energy sub-relativistic CRs responsible for the ionisation and a mildly relativistic population responsible for the pressure) is needed before stating any conclusion. The first situation, where ISM ionisation properties are not modified, can still be valid if the mean free path of the sub-relativistic CRs is smaller than the mean free path of mildly relativistic CRs. How these mean free paths compare with each other strongly depends on the nature of the turbulence in the SNR and in its close environment. The latter in practice is strongly dependent on the star formation history that occurs at the location of CR release.

A detailed and dynamical calculation of I(k) and D_∥ with respectto the different types of partially ionised phases will be presented in Paper II. A major difficulty, also discussed in Paper II, in implementing this estimate of CR diffusion in any MHD code is that the scales of the self-generated waves cannot be resolved by the codes. This model of turbulence will be implemented as a sub-grid model.

5.2.3. CR impact over small star formation scales

Young stellar objects (YSO) are also potential sites of particle acceleration and can be sources of turbulence injection at ISM scales ≃ 0.1−1 pc. Padovani et al. (2015) have shown that background CRs cannot be responsible for the observed high ionisation rates observed in some YSOs (see references therein), which can exceed usual rates by one to four orders of magnitude. They propose that a process of strong in situ particle acceleration should occur in these objects. Padovani et al. (2016) have identified several potential acceleration sites in YSOs: the jet termination shocks in class 0–I YSOs, and accretion shocks in class 0–I or more evolved objects. Theregion of interaction between the accretion disc and stellar magnetic fields can also produce episodic reconnection events (de Gouveia Dal Pino et al. 2010). Several groups have investigated the impact of freshly accelerated energetic particles in isolated YSO (Rodgers-Lee et al. 2017; Rab et al. 2017) or in groups (Gaches & Offner 2018) of YSOs over their surrounding accretion discs or envelopes. These models find that YSOs are, in all these acceleration zones, able to accelerate these CRs in the relativistic regime. These CRs can provide both ionisation and pressure. If ionisation by CRs is an obvious feedback to take into account in models, the feedback produced by CR pressure gradients is a possible supplementary issue to consider in these environments. Of course, in accretion discs, envelopes, or in jets, particle propagation and self-generated wave production are strongly limited by ion-neutral collisions (Padovani et al. 2016). On the other hand, CR gradients can develop over scales between 1 and 1000 AU depending on the site, so leading to higher $P_{{\text{ref}}^{\prime }}$ values. Among the different acceleration sites invoked above, jets could contribute to drive strong CR gradients. For instance, taking a reference gradient produced by a CR pressure of 1 eV   cm⁻³ over 1000 AU leads to $P_{\text{CR},{\text{ref}}^{\prime }}\simeq 1.1\times {10}^{-28} \text{dyne} {\text{cm}}^{\text{-3}}$. Considering typical molecular cloud core properties as n_n ≃ 10⁴ cm⁻³, n_i ≃ 10⁻² cm⁻³, μ_i ≃ 29 (Xu et al. 2016) and a simulation box L ≃ 0.1 pc with q = 1/2, we find D_crit ≃ 10²² cm² s⁻¹² and D_∥ ≃ 6.7 × 10²⁵ cm² s⁻¹. A dynamical effect of CRs over jet environments seems then to require a very large $P_{{\text{CR}}^{\prime }}/P_{\text{CR},{\text{ref}}^{\prime }}$ ratio. This speculative estimate, however, requires a more detailed analysis postponed to a future investigation where the interior of molecular clouds will be considered.

5.3. Gamma-ray diffuse emission

Gamma-ray diffuse emission can provide some constraints over CR propagation in different parts of the ISM. Gamma rays are either produced by hadrons through pion production, or by leptons through Inverse Compton emission, or through Bremsstrahlung if the gas density is high enough. The Fermi-LAT detector has observed local molecular clouds (Ackermann et al. 2012b) and giant molecular clouds (Ackermann et al. 2012a), and the Fermi Collaboration has obtained a model of diffuse galactic emission (Ackermann et al. 2011; Acero et al. 2016). Among the main results extracted from these data we list: local molecular clouds have gamma-ray spectra compatible with the emission produced by the local CR spectrum observed on Earth but show a variation in normalisation by a factor of 20%; diffuse gamma-ray emission shows a trend with more intense and harder CR spectra towards small galactic radii to less intense and softer CR spectra in the outer galaxy (Acero et al. 2016). However, the gamma-ray signal in the outer galaxy is larger than expected from the standard CR diffusion model from known CR source distributions; this is the so-called CR gradient problem. The GeV gamma-ray signal from the Large Magellanic Cloud (LMC) requires a 1–100 GeV CR population with a density of ~30% larger than the local Galactic value.

It is difficult to compare our results with these data. Fermi observed molecular clouds usually in active star formation regions whereas our simulations only consider an ISM free of CR sources. It seems that very active star-forming regions such as the galactic centre or the LMC produce more intense CR components with respect to the local Galactic value. Observations of active CR sources may show that propagation of CRs is different in the environment of these sources. This is for instance the case for the W28 SNR detected at GeV by Fermi (Abdo et al. 2010) and at TeV by H.E.S.S. (Aharonian et al. 2008). Gamma-ray emission from clouds located near the SNR may be explained if CRs gradually escape from the accelerator, the highest energies being released first. In this framework, diffusion coefficients for TeV CRs required to explain the gamma rays need to be reduced by two orders of magnitude with respect to the standard values (Nava & Gabici 2013; Yan et al. 2012). Uchiyama et al. (2012) report on the Fermi detection of molecular gas around the W44 SNR. The observations may be interpreted by invoking anisotropic diffusion from the SNR and in that case they would also require some suppression of CR diffusion. These findings advocate in favour of a modified CR propagation in the vicinity of the CR accelerator. The suppression of diffusion is likely a time-dependent process, which depends on the CR energy and the properties of the surrounding ISM (see Nava et al. 2016). This transitory aspect makes a comparison between gamma-ray observations and properties of CR-modified ISM turbulence difficult and, at the stage of the present work, premature. The possible impact of CRs over ISM turbulence, in particular around CR sources requires us to include the effect of sources in several ISM numerical realisations and to develop a statistical analysis of the gamma-ray signal from synthetic ISM maps. This modelling is beyond the scope of this study and deserves a future analysis.

5.4. Model limitation and future works

In this work, we neglect the effect of CR advection due to the streaming, CR heating due to streaming instability, and CR radiative losses. It has been shown that these processes can have significant effects on large galactic scales for the launching of galactic winds (Ruszkowski et al. 2017; Wiener et al. 2017). It is nevertheless not clear whether they also affect CR propagation at the scales considered in this work. Figures 1, 3, and 4 indicate that the flow can be sub-Alfvénic, in particular within the WNM. As a consequence, we may underestimate the cosmic-ray propagation in the low density regions. These aspects will be investigated in a future work.

Another limitation comes from the rather idealised set-up used in this study. Although it is more realistic to study ISM dynamics with heating and cooling processes able to handle the thermal instability, the effects of large (galactic) scales and of self-consistent turbulence driving by supernovae or galactic shears are not considered here. It is currently admitted that CR pressure can modify the star formation rate (SFR) by providing extra support, but the details of CR propagation through the galactic scale to the molecular cloud scale are not yet understood. A step forward will be to consider the zooming technique as well as self-gravity, as proposed in Hennebelle (2018) for instance. However, this step is far beyond our current work.

Lastly, we show in Appendix A that our results may not have converged in the trapped CR regime. The CNM gas fraction is underestimated at low resolution and the effect of CR pressure gradients is overestimated. A quantitative analysis has to be done, in particular with respect to the impact the CR can have on the SFR. Nevertheless, we note that CR trapping in the high density regions is robust at high resolution. Local enhancements of the CR energy density are thus still expected, which may have an influence on the ionisation budget of star-forming clouds. In particular, CR ionisation remains the most efficient ionisation process within collapsing dense cores, and appears to control the coupling between the gas and the magnetic fields within the centre of the collapsing region (Marchand et al. 2016; Wurster et al. 2018). Future work is needed in order to quantify the variations of the CR energy density within molecular clouds.

6. Conclusion

We have studied the conditions for (GeV) CR trapping within a turbulent and magnetised ISM. First, we have derived an analytical prediction in Eq. (22) for the critical diffusion coefficient D_crit below whichCRs are trapped. In the case of a turbulent ISM following the observed scaling relations, we get typical values of 10²⁴ −10²⁵ cm² s⁻¹ for D_crit. Then, we presented numerical experiments of the evolution of a mixture of CRs and gas, which is subject to the thermal instability. We explored the parameter space and varied the initial magnetisation, the initial CR energy density, as well as the turbulence amplitude and driving scale. Our results indicate that our analytical prediction is very robust over the entire explored parameter space. The two regimes of trapped versus free streaming CRs are always recovered, with a transition occurring around 10²⁴ −10²⁵ cm² s⁻¹. The initial CR pressure is found to have a dramatic effect if it exceeds the gas thermal pressure by a factor of ten. In this case, the non-thermalsupport provided by the CRs is sufficient to inhibit the development of the thermal instability even if the CR diffusion coefficient is larger than the critical value. In the trapped CR regime, the width of the gas density PDF gets narrower, as a result of the larger modified sound speed that accounts for CR pressure. The width of the gas density PDF can be related to the SFR (e.g. Hennebelle & Chabrier 2011), and any physical effect, such as turbulence and magnetic fields, can then have an impact on star formation. We have shown that in the trapped CR regime, as well as in CR-pressure-dominated ISM, the width of the PDF is reduced, so that CRs can have a dynamical impact on star formation. Lastly, we discussed possible mechanisms that can lower the diffusion coefficient to values that are smaller than the canonical admitted ones of a few 10²⁸ cm² s⁻¹. The turbulence generation by CR streaming especially around CR sources, as well as models for large-scale turbulence, appear to be potential mechanisms able to modify diffusion coefficient values.

To summarise, we have shown that CR pressure gradients can develop within the ISM forming molecular clouds, if the condition for CR trapping is satisfied or if the CR pressure largely exceeds the gas thermal pressure. CR pressure gradients affect the gas dynamics, as well as the propagation of CRs. This feedback loop needs to be better understood in future works, which will account for local variations of the CR diffusion coefficient as a function of large-scale turbulent models, as well as sub-grid models for turbulence generated by CR streaming.

Acknowledgements

We thank the anonymous referee for his or her constructive comments that helped to improve the quality of the manuscript. B.C. thanks Andrew McLeod for his contribution to the driven turbulence module. This work was supported by the CNRS programmes“Programme National de Cosmologie et Galaxies” (PNCG), “Physique et Chimie du Milieu Interstellaire” (PCMI), and “Programme national des hautes énergies” (PNHE). This project was supported by the IDEXLyon project (contract No. ANR-16-IDEX-0005) under the auspices of the University of Lyon. This work was granted access to the HPC resources of CINES (Occigen) under the allocation DARI A0020407247 made by GENCI. Some of the figures were created using the OSIRIS³ visualisation package for RAMSES.

Appendix A: Convergence with the resolution

In this appendix, we perform a resolution study for the fiducial case, with k_turb = 2, and for D₀ = 10²⁴ and D₀ = 10²⁶ cm² s⁻¹. We perform four calculations with D₀ = 10²⁴ cm² s⁻¹ and a resolution ranging from 64³ to 512³, and three calculations with D₀ = 10²⁶ cm² s⁻¹ and a resolution ranging from 64³ to 256³. We do not present results at a resolution of 512³ with D₀ = 10²⁶ cm² s⁻¹ because the CPU time needed to reach only one dynamical time is too expensive. This comes from the bad efficiency of our implicit solver when dealing with rather smooth CR energy density and highly anisotropic diffusion.

Figure A.1 shows gas density and CR energy density maps for the four calculations done with D₀ = 10²⁴ cm² s⁻¹. The denser structures get more compact as resolution increases. We attribute this behaviour to the poor sampling of the CNM structures with a resolution less than 256³ for a 40 pc box. Audit & Hennebelle (2005) argued that the typical size of the fragments of CNM is ≃ 0.1 pc. We are clearly sub-sampling this scale with our fiducial resolution, which indicates that we cannot properly describe the structure of the densest gas, that is, the molecular clouds. Nevertheless, the cooling length in the WNM, ≃ 1−3 pc at a density of 2 cm⁻³ (Hennebelle & Pérault 1999), is always resolved. The high CR energy density regions also exhibit smaller structures as resolution increases, but not as much as the gas density ones.

Fig. A.1. Gas density (top panels) and CR energy density (bottom panels) maps for the fiducial runs with D0 = 1024 cm2 s−1 in the plane z = 20 pc at a time corresponding to 1.2tcross and a resolution of 643, 1283, 2563, and 5123 (from left to right).

Figure A.2 shows the PDFs of the gas density, temperature, and CR energy density of the seven different calculations at times greater than1t_cross. The gas density and temperature PDFs broaden as resolution increases. In particular, the PDFs of the case with D₀ = 10²⁴ cm² s⁻¹ do not show convergence as resolution increases. The gas fraction in the CNM increases significantly in the 512³ run. This feature was also reported in other works that report more pronounced phase transition at higher resolution (Seifried et al. 2011). Nevertheless, the CR energy density PDFs are not very sensitive to the resolution, except for the high energy tail. Importantly, the qualitative behaviour of the CR energy density does not change if resolution increases, with a transition to the trapped CRs regime with a low diffusion coefficient. We note, however, that given the fiducial resolution we use of 128³, our study cannot be quantitative for the CNM gas fraction and structures.

Fig. A.2. Same as Fig. 5 for different numerical resolutions in the fiducial case (kturb = 2) with D0 = 1024 cm2 s−1 (blue) and D0 = 1026 cm2 s−1 (red), but at times greater than 1tcross.

Appendix B: Isotrop versus anisotrop CR diffusion

Cosmic rays are charged particles and their propagation is therefore highly dependent on the magnetic field configuration, which is thus anisotropic in nature. We show in this appendix to what extent the anisotropic case differs from a simple isotropic diffusion case. We perform four runs of the fiducial case, with a diffusion coefficient ranging from 10²² to 10²⁸ cm² s⁻¹, but we set the perpendicular diffusion coefficient equal to the parallel one.

Figure B.1 compares the PDFs of the gas density, gas temperature, and CR energy density for the cases with anisotropic diffusion and with isotropic diffusion. For the two extreme cases, D₀ = 10²² and D₀ = 10²⁸ cm² s⁻¹, the PDFs are not very sensitive to the nature of the diffusion process. In the first case, the diffusion coefficient is too low and CRs are trapped anyway, while in the second case diffusion is too fast. We find differences in the CR energy density PDF with D₀ = 10²⁶ cm² s⁻¹. The PDF is indeed narrower with isotropic diffusion because the diffusion process is more efficient. With D₀ = 10²⁴ cm² s⁻¹, the three PDFs are different if isotropic diffusion is used. The CRs can escape more efficiently and cannot build up pressure gradients that act on the gas motions. Since this value is close to the critical value, a slight change of roughly afactor of three in the overall diffusion coefficient makes important differences. Overall, these tests indicate that anisotropic diffusion has to be taken into account if one aims at describing properly the transition to CR trapping.

Fig. B.1. Same as Fig. 5. Effect of anisotropic diffusion versus isotropic for the fiducial runs.

References

All Tables

All Figures

	Fig. 1. Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$ as a function of the gas density n and of the length scale L following Eq. (7). This plot is valid only if one assumes the scaling relation and magnetic field strength observed in the CNM and WNM.
In the text

	Fig. 2. Time evolution for six simulations with kturb = 2 of the clumping factor of the gas density n (top panel), the cosmic-ray energy density ECR (middle panel), and the Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$ (bottom panel). The clumping factor C(X) of a quantity X is defined as C(X) = ⟨X2⟩∕⟨X⟩2.
In the text

Fig. 3. Fiducial runs with kturb = 2: maps of the gas density (top panels), cosmic-ray energy density (middle panels), and Alfvénic Mach number (bottom panels) in the plane z = 20 pc at a time corresponding to 2tcross. From left to right panels: CR diffusion coefficient is D0 = 1022, 1024, 1026, and 1028 cm2 s−1. The arrows on the density maps represent the velocity vectors (normalised to a value of 14 km s−1, see upper right panel), and the streamlines on the CR energy density maps represent the magnetic field lines, both in the xy-plane.

In the text

Fig. 4. Two-dimensional histograms of, from top to bottom panels, the CR energy density, the Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$, the ratio ${\mathcal{M}}_{\text{s},\text{CR}}/{\mathcal{M}}_{\text{s}}$, and ζ = PCR∕P as a function of the density for the fiducial runs kturb = 2 with CR diffusion at the same time as in Fig. 3. The colour coding indicates the temperature, the magnetic field strength, and the sonic Mach number ${\mathcal{M}}_{\text{s}}$ (two bottom rows), respectively.

In the text

	Fig. 5. PDF of the gas density (top panel), temperature (middle panel), and CR energy density (bottom panel) for the fiducial runs with kturb = 2 at the same time as in Fig. 3, for D0 = 1022 (purple), 1024 (blue), 1026 (red), and 1028 cm2 s−1 (black). Two models with D0 = 0 (dashed cyan) and without CRs (dashed grey) are also represented for comparison.
In the text

	Fig. 6. Same as Fig. 5 for the models with kturb = 2 and D0 = 1024 (blue), D0 = 1026 cm2 s−1 (red). The solid lines represent the fiducial model with ζ = 0.1, the dashed line the models with ζ = 1, the dotted line ζ = 10, and the dashed-dotted line ζ = 100.
In the text

	Fig. 7. Same as Fig. 5 for kturb = 4.
In the text

	Fig. 8. Same as Fig. 5 for kturb = 8.
In the text

	Fig. 9. Same as Fig. 5 for kturb = 16.
In the text

	Fig. 10. Same as Fig. 5 for kturb = 2, but different amplitude in the turbulence driving. The solid lines represent the fiducial case results with vrms = 5.7 km s−1 and the dashed lines the ones obtained with vrms = 10.5 km s−1.
In the text

Fig. 11. Same as Fig. 5 for kturb = 2, but different amplitudes of the initial magnetic field strength. The solid lines represent the fiducial case with ${\mathcal{M}}_{\text{A}}=4.5$, the dashed lines ${\mathcal{M}}_{\text{A}}=1.2$, the dotted lines ${\mathcal{M}}_{\text{A}}=0.6$, and the dashed-dotted lines ${\mathcal{M}}_{\text{A}}=20$. The blue colour represents results with D0 = 1024 cm2 s−1, the red D0 = 1026 cm2 s−1.

In the text

	Fig. A.1. Gas density (top panels) and CR energy density (bottom panels) maps for the fiducial runs with D0 = 1024 cm2 s−1 in the plane z = 20 pc at a time corresponding to 1.2tcross and a resolution of 643, 1283, 2563, and 5123 (from left to right).
In the text

	Fig. A.2. Same as Fig. 5 for different numerical resolutions in the fiducial case (kturb = 2) with D0 = 1024 cm2 s−1 (blue) and D0 = 1026 cm2 s−1 (red), but at times greater than 1tcross.
In the text

	Fig. B.1. Same as Fig. 5. Effect of anisotropic diffusion versus isotropic for the fiducial runs.
In the text