Issue 
A&A
Volume 622, February 2019



Article Number  A143  
Number of page(s)  18  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201833809  
Published online  12 February 2019 
Cosmicray propagation in the bistable interstellar medium
I. Conditions for cosmicray trapping
^{1}
Centre de Recherche Astrophysique de Lyon UMR5574, ENS de Lyon, Univ. Lyon1, CNRS, Université de Lyon,
69007
Lyon,
France
email: benoit.commercon@enslyon.fr
^{2}
Laboratoire Univers et Particules de Montpellier (LUPM), Université Montpellier, CNRS/IN2P3, CC72, place Eugène Bataillon,
34095
Montpellier Cedex 5,
France
email: Alexandre.Marcowith@umontpellier.fr
^{3}
Institut d’Astrophysique de Paris, UMR 7095, CNRS, UPMC Université Paris VI,
98 bis boulevard Arago,
75014
Paris,
France
email: dubois@iap.fr
Received:
9
July
2018
Accepted:
23
November
2018
Context. Cosmic rays propagate through the galactic scales down to the smaller scales at which stars form. Cosmic rays are close to energy equipartition with the other components of the interstellar medium and can provide a support against gravity if pressure gradients develop.
Aims. We study the propagation of cosmic rays within the turbulent and magnetised bistable interstellar gas. The conditions necessary for cosmicray trapping and cosmicray pressure gradient development are investigated.
Methods. We derived an analytical value of the critical diffusion coefficient for cosmicray trapping within a turbulent medium, which follows the observed scaling relations. We then presented a numerical study using 3D simulations of the evolution of a mixture of interstellar gas and cosmic rays, in which turbulence is driven at varying scales by stochastic forcing within a box of 40 pc. We explored a large parameter space in which the cosmicray diffusion coefficient, the magnetisation, the driving scale, and the amplitude of the turbulence forcing, as well as the initial cosmicray energy density, vary.
Results. We identify a clear transition in the interstellar dynamics for cosmicray diffusion coefficients below a critical value deduced from observed scaling relations. This critical diffusion depends on the characteristic length scale L of D_{crit} ≃ 3.1 × 10^{23} cm^{2} s^{−1}(L/1 pc)^{q+1}, where the exponent q relates the turbulent velocity dispersion σ to the length scale as σ ~ L^{q}. Hence, in our simulations this transition occurs around D_{crit} ≃ 10^{24}–10^{25} cm^{2} s^{−1}. The transition is recovered in all cases of our parameter study and is in very good agreement with our simple analytical estimate. In the trapped cosmicray regime, the induced cosmicray pressure gradients can modify the gas flow and provide a support against the thermal instability development. We discuss possible mechanisms that can significantly reduce the cosmicray diffusion coefficients within the interstellar medium.
Conclusions. Cosmicray pressure gradients can develop and modify the evolution of thermally bistable gas for diffusion coefficients D_{0} ≤ 10^{25} cm^{2} s^{−1} or in regions where the cosmicray pressure exceeds the thermal one by more than a factor of ten. This study provides the basis for further works including more realistic cosmicray diffusion coefficients, as well as local cosmicray sources.
Key words: magnetohydrodynamics (MHD) / methods: numerical / ISM: structure / diffusion / cosmic rays / ISM: individual objects: molecular clouds
© 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^{−3} 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 backreact over ISM turbulence and its dynamics (Kulsrud & Pearce 1969).
CRs have a multiscale impact over galactic ISM dynamics (Grenier et al. 2015): they drive largescale winds (Hanasz et al. 2013; Booth et al. 2013; Salem & Bryan 2014; Girichidis et al. 2016; Wiener et al. 2017), generate largescale 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; RodgersLee 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 varianceMach number relation (VazquezSemadeni 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 nonisothermal 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 lognormal distribution that fits the probability density function. In addition to the abovementioned 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 bistable 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 layout 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_{2}), 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; MivilleDeschênes et al. 2017). The turbulent velocity dispersion increases with the size as (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^{−1} (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 nonthermal 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
where n is the gas number density, and L is the size of the region in which the nonthermal 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 (4)
These scaling relations indicate that velocity fluctuations are reduced as the gas becomes denser, until selfgravity takes over to initiate the gravitational protostellar collapse and to reaccelerate 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 nonthermal motions are due to the turbulence only and consider that the nonthermal velocity is the rms velocity dispersion v_{rms}.
It is worth noting that the exponent values q ≃ 0.35−0.5 in the velocitysize 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 (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^{2} for an enclosed mass M ∝ nR^{3} implies that B ∝ n^{2∕3}. The masstoflux conservation gives B ∝ nh, where h is the thickness of the cloud along the field lines. Equilibrium along the field lines implies that , with c_{s} the sound speed and ϕ the gravitational potential. Finally, the Poisson equation leads to ϕ ∝ c_{s}h^{2}. 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 intensitydensity scaling (e.g. Crutcher 2012). At low density n < 300 cm^{−3}, the magnetic intensity is insensitive to the density variations, with B_{0,max} ≃ 10 μG. While for densities n > 300 cm^{−3}, 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 (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 , defined as , with the Alfvén speed and μ_{gas} the gas mean molecular weight. By combining the abovementioned scaling relations for B and v_{rms} at the densities n < 100 cm^{−3} we interest us in this study, the Alfvénic Mach number can be estimated as (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 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 subAlfvénic turbulence if the flow satisfies the abovementioned scaling relations. At a gas density of ≤ 10 cm^{−3}, the flow is always subAlfvénic, while at larger densities it may be superAlfvénic. CRs propagate along magnetic field lines, and the topology of the lines thus affects the global propagation speed. In the superAlfvé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 subAlfvé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 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 rootmeansquare velocity v_{rms}. The turbulent time is defined as the time at which a turbulent fluctuation travels through the blob, (8)
We define the diffusion time t_{diff} as (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} (10)
Using the scaling relation in Eq. (5), the critical diffusion coefficient is (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^{25} cm^{2} s^{−1}. A usual parametrisation of the CR parallel diffusion coefficient derived from CR data fitting is (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^{28} cm^{2} s^{−1} 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 cosmicray 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 adaptivemeshrefinement code RAMSES (Teyssier 2002), which is based on a finitevolume, secondorder Godunov scheme and a constrained transport algorithm for ideal MHD (Fromang et al. 2006; Teyssier et al. 2006). We use the HLLD (HartenLaxVan LeerDiscontinuities) 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 advectiondiffusion 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 CRMHD equations read
where ρ is the gas density, u is the fluid velocity, 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^{2}∕(8π). The magnetic field vectors is B, E_{tot} is the total energy E_{tot} = E + 0.5ρu^{2} + B^{2}∕(8π) + E_{CR}. The anisotropic diffusion tensor of cosmic rays is , given by (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 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 subAlfvé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 OrnsteinUhlenbeck 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 selfgravity in this study.
The system is closed using the perfect gas relation, P = ρk_{B}T∕(μ_{gas}m_{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 (19)
is the modified sound speed which accounts for the CR pressure.
3.2. Initial conditions
Our initial setup 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_{0} ≃ 4460 K, density n_{0} = 2 cm^{−3}, gas thermal pressure P ≃ 1.23 × 10^{−12} dyne cm^{−2}. The initial sound speed is c_{s} ≃ 5.2 km s^{−1}. The initial cosmicray 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_{0} = 1 μG, which gives an initial Alfvénic speed v_{A} ≃ 1.3 km s^{−1} 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^{3}. Saury et al. (2014) performed a resolution study (from 128^{3} to 512^{3}) and have shown that even though the scales to describe the cold CNM and the molecular clouds are barely resolved with 128^{3}, 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^{3} resolution for a qualitative study.
We explored the parameter space with the five following parameters: D_{0} 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 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 superAlfvénic regime, while the initial turbulent Mach numbers vary from slightly supersonic with k_{turb} = 2 to subsonic with k_{turb} = 16 ().
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 (21)
Since we use a unique value of 40 pc for the size of the computational domain L_{box}, the critical diffusion coefficient is (22)
We give in Table 1 the values of the critical diffusion coefficient for each setup, which ranges from 1.4 × 10^{24} to 6.3 × 10^{25} cm^{2} s^{−1}. We varied D_{0} from 10^{22} to 10^{28} cm^{2} s^{−1} 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_{0} and v_{rms}, that is, L_{crit} = D_{0}∕v_{rms,} (23)
In the case k_{turb} = 2 and v_{rms} = 5.7 km s^{−1}, the critical length varies from about 5.7 kpc for D_{0} = 10^{28} cm^{2} s^{−1} to 5.6 × 10^{−3} pc for D_{0} = 10^{22} cm^{2} s^{−1}. The condition for CR trapping is thus simply L_{crit} < L_{inj}.
Summary of the initial conditions of all the simulations.
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^{−1}, which corresponds to a turbulent crossing time t_{cross} = 6.9 Myr. We explore five values for the diffusion coefficient: 0, 10^{22}, 10^{24}, 10^{26}, 10^{28} cm^{2} s^{−1}. 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^{25} cm^{2} s^{−1}.
Figure 2 shows the time evolution of the clumping factor of the gas density, CR energy density, andAlfvénic Mach number for the five abovementioned simulations with various diffusion coefficients D_{0} and the model without CRs. The clumping factor of a quantity X is defined as C(X) = ⟨X^{2}⟩∕⟨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 quasistationary 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 twophases CNMWNM medium. The clumping factor of the gas density shows two distinct regimes. In the case where the CR diffusion coefficient is large, D_{0} ≥ 10^{26} cm^{2} s^{−1}, 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^{3} between the minimum and maximum values. On the contrary, with D_{0} ≤ 10^{24} cm^{2} s^{−1} 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_{0} ≥ 10^{26} cm^{2} s^{−1}, 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_{0} ≤ 10^{24} cm^{2} s^{−1}, 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 fieldsdensity 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_{0} ≤ 10^{24} cm^{2} s^{−1}. Consequently, the gas density variations are tempered by the extra pressure support provided by the gradients in the cosmicray 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_{0} ≤ 10^{24} cm^{2} s^{−1} and D_{0} ≥ 10^{26} cm^{2} s^{−1} 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_{0} ≥ 10^{26} cm^{2} s^{−1} (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 2Dhistograms of the CR energy density, the Alfvénic Mach number , the ratio , 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 (for the two bottom rows), respectively. The correlation between the CR energy density and the gas density is obvious for D_{0} ≤ 10^{24} cm^{2} s^{−1} 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^{4} 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 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_{0} ≥ 10^{26} cm^{2} s^{−1} in the low density regions but no influence at high densities. The behaviour of the pressures ratio ζ shows also opposite trends. For D_{0} ≤ 10^{24} cm^{2} s^{−1}, the CR pressure outmatches the thermal one with increasing density, being the greatest in the highly supersonic regions (). 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_{0} ≥ 10^{26} cm^{2} s^{−1}, ζ 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 volumeweighted PDF of the gas density, massweighted PDF of the gas temperature, and massweighted PDF of the CR energy density are shown in Fig. 5. In the remainder of the paper, all the PDFs are timeaveraged, 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} = 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_{0} ≥ 10^{26} cm^{2} s^{−1}, 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_{0}≥ 10^{26} cm^{2} s^{−1}). The temperature variations are reduced in the models with D_{0} ≤ 10^{24} cm^{2} s^{−1}. 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 anticorrelated with those of the gas density and temperature. For D_{0} = 10^{28} cm^{2} s^{−1}, 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_{0} = 10^{26} cm^{2} s^{−1}. CR pressure gradients thus develop, but have no effect on the gas flow as indicated in the density PDF. For D_{0} = 10^{24} cm^{2} s^{−1}, the CR energy density PDF is much wider (≃2.5 decades at PDF = 10^{−3}), 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_{0} = 10^{22} cm^{2} s^{−1} and D_{0} = 10^{24} cm^{2} s^{−1} while this is not the case for the gas density and the temperature. For D_{0} = 10^{22} cm^{2} s^{−1}, 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_{0} = 10^{28} cm^{2} s^{−1} correspond to those of the run without CRs.
The picture we can draw from these fiducial runs is the following: for diffusion coefficients D_{0} ≤ 10^{24} cm^{2} s^{−1}, 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_{0} ≤ 10^{24} cm^{2} s^{−1}, 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 subAlfvé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_{0} = 10^{26} cm^{2} s^{−1}, 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 k_{turb} = 2 of the clumping factor of the gas density n (top panel), the cosmicray energy density E_{CR} (middle panel), and the Alfvénic Mach number (bottom panel). The clumping factor C(X) of a quantity X is defined as C(X) = ⟨X^{2}⟩∕⟨X^{⟩2}. 
Fig. 3. Fiducial runs with k_{turb} = 2: maps of the gas density (top panels), cosmicray energy density (middle panels), and Alfvénic Mach number (bottom panels) in the plane z = 20 pc at a time corresponding to 2t_{cross}. From left to right panels: CR diffusion coefficient is D_{0} = 10^{22}, 10^{24}, 10^{26}, and 10^{28} cm^{2} 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 xyplane. 
Fig. 4. Twodimensional histograms of, from top to bottom panels, the CR energy density, the Alfvénic Mach number , the ratio , and ζ = P_{CR}∕P as a function of the density for the fiducial runs k_{turb} = 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 (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_{0} = 10^{24} and D_{0} = 10^{26} cm^{2} s^{−1}, 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_{0} = 10^{26} cm^{2} s^{−1} 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_{0} = 10^{26} cm^{2} s^{−1}, 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 lognormal 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_{0} = 10^{24} cm^{2} s^{−1} 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^{−3}, 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 k_{turb} = 2 at the same time as in Fig. 3, for D_{0} = 10^{22} (purple), 10^{24} (blue), 10^{26} (red), and 10^{28} cm^{2} s^{−1} (black). Two models with D_{0} = 0 (dashed cyan) and without CRs (dashed grey) are also represented for comparison. 
Fig. 6. Same as Fig. 5 for the models with k_{turb} = 2 and D_{0} = 10^{24} (blue), D_{0} = 10^{26} cm^{2} 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 dasheddotted line ζ = 100. 
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_{0} = 10^{24} cm^{2} s^{−1}), 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_{0} = 10^{24} cm^{2} s^{−1}. In addition, since the turbulent energy decreases with increasing k_{turb}, D_{0} = 10^{24} cm^{2} s^{−1} 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_{0} as k_{turb} increases, as expected in the qualitative analysis above. For D_{0} ≥ 10^{26} cm^{2} s^{−1}, 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_{0} ≥ 10^{26} cm^{2} s^{−1}, 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_{0} = 10^{24} cm^{2} s^{−1}. The transition to the trapped CR regime remains on the same parameter range.
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^{−1}. We performed four simulations with a CR diffusion coefficient ranging from D_{0} = 10^{22} to D_{0} = 10^{28} cm^{2} s^{−1}.
Figure 10 shows the gas density, temperature, and CR energy density PDFs with v_{rms} = 10.5 km s^{−1}. The PDFs of the fiducial case with v_{rms} = 5.7 km s^{−1} 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_{0}, 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_{0} ≤ 10^{24} cm^{2} s^{−1}. 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. 10. Same as Fig. 5 for k_{turb} = 2, but different amplitude in the turbulence driving. The solid lines represent the fiducial case results with v_{rms} = 5.7 km s^{−1} and the dashed lines the ones obtained with v_{rms} = 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 is explored by changing the value of the initial magnetic field amplitude. We compared the fiducial case results, with (superAlfvénic), with the ones obtained with (superAlfvénic, weak magnetic fields), (transAlfvénic), and (subAlfvénic). We recall that magnetic field amplitude increases with decreasing . We use two values for the CR diffusioncoefficient, D_{0} = 10^{24} and D_{0} = 10^{26} cm^{2} s^{−1}.
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_{0} ≤ 10^{26} cm^{2} s^{−1}. The gas density and temperature PDFs do not depend much on for D_{0} = 10^{26} cm^{2} s^{−1}. A noticeable feature is the width of the gas density PDF that is the narrowest with . 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 increases. The magnetic field lines are indeed more tangled in the superAlfvénic case, which reduces the effective diffusion coefficient. For D_{0} = 10^{24} cm^{2} s^{−1}, we see clearly two different regimes as a function of , and the results are qualitatively different from the case with D_{0} = 10^{26} cm^{2} s^{−1}. In the superAlfvé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 superAlfvé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 nonlinear 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 k_{turb} = 2, but different amplitudes of the initial magnetic field strength. The solid lines represent the fiducial case with , the dashed lines , the dotted lines , and the dasheddotted lines . The blue colour represents results with D_{0} = 10^{24} cm^{2} s^{−1}, the red D_{0} = 10^{26} cm^{2} 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^{24} –10^{25} cm^{2} s^{−1} 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 nonthermal 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 CRpressuredominated 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 superAlfvé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 subAlfvé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 subAlfvé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 largescale 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. Largescale turbulence models
Largescale motions in the ISM are sources of turbulence (Mac Low & Klessen 2004). Among the potential processes that are sources of turbulence we find magnetorotational 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, largescale 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 largescaleinjected turbulent motions. Let us recall their main results.
SuperAlfvénic regime
If , then largescale motions are hydrodynamical and follow a Kolmogorov law. At a scale the turbulence becomes transAlfvé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 superAlfvé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 = ℓ_{A}v∕3 while if λ_{∥} < L_{inj}, D = λ_{∥}v∕3, where v = β_{c}c is the CR velocity. The calculation of λ_{∥} in partially ionised phases is complex and involves a calculation of the turbulence damping scale by ionneutral 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), (24)
If , 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 for q = 1/3 and 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}).
SubAlfvénic regime
If , 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 the Goldreich–Sridhar scaling prevails again. SubAlfvénic CR diffusion is no longer isotropic. The perpendicular diffusion coefficient depends again on λ_{∥} (Yan & Lazarian 2008). If λ_{∥} > L_{inj} then , while if λ_{∥} < L_{inj}, ; 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} (25)
which in our fiducial case using Eq. (11) is λ_{∥,crit} ≃ 4.9 × 10^{15} cm for q = 1/3 or λ_{∥,crit} ≃ 9.0 × 10^{15} 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 transsonic 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 socalled streaming instability (Skilling 1971). Perturbations generated by CRs are resonant and have typical wave numbers , 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 , where L_{CR} is the CR gradient scale. A reference value is , where P_{CR, eVcc} is the CR pressure in units of 1 eV cm^{−3} and L_{CR} is the CR gradient scale in units of 100 pc. In the Milky Way, one expects 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 selfgenerated 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 ionneutral collision frequency. We take and with n_{n} the number density of neutrals, in the WNM and the CNM, respectively (Jean et al. 2009). In the decoupled ionneutral regime the streaming instability growth rate (in s^{−1}) can be written as (Nava et al. 2016), where is the ion Alfvén velocity, W_{B} = B^{2}∕8π is the ambient mean magnetic energy density, and is the relative amplitude of perturbations at a wave number 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 quasilinear theory (Schlickeiser 2002) by and D_{⊥} ≃ , and again v = β_{c}c is the particle speed.
The parallel diffusion coefficient is then (26)
or comparing it with D_{crit} given by Eq. (11), (27)
We have expressed the CR energy E_{GeV} in GeV units, the ion density n_{i} in units of 10^{−2} cm^{−3}, the ionneutral collision frequency in units of 10^{−8} s^{−1}, 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 selfgenerated waves over ISM dynamics, that is to have D_{∥} ~ D_{crit}, we need a ratio . 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 at a CR kinetic energy of 1 GeV in order to get D_{∥} = D_{crit} for q = 1/3 and for q = 1/2. In the CNM we find still at 1 GeV in order to get D_{∥} ≃ D_{crit} for q = 1/3 and for q = 1/2. If we assume μ_{i} = 1 we find 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 by a factor 10–30 in SNR or even more in starforming 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 ionneutral collisions. The authors find that at 10 GeV, so for CRs released late on in the SNR lifetime, diffusion is suppressed^{1} because of scattering off the selfgenerated 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 kiloyears (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^{−3} in the ISM over a typical scalelength 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 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 quasilinear limit. From the condition Γ_{g} = Γ_{in} we deduce a turbulence level of (28)
In the conditions that prevail in partially ionised phases, using B ≃ 6 μG, the ratio cannot exceed ~3 and ~ 300 in WNM and CNM, respectively, for the quasilinear condition (I(k) ≤ 10^{−4}) 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. Nonlinear calculations as performed in Nava et al. (2016) tend to confirm this assertion. Another caveat of these calculations is that low energy subrelativistic 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 ionneutral collision rate decreases and so the damping of selfgenerated wave also decreases. This second effect has then a positive feedback on the CR ISM dynamical impact. As all quantities above depend on , 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 subrelativistic 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 subrelativistic 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 selfgenerated waves cannot be resolved by the codes. This model of turbulence will be implemented as a subgrid 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 (RodgersLee 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 selfgenerated wave production are strongly limited by ionneutral 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 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^{−3} over 1000 AU leads to . Considering typical molecular cloud core properties as n_{n} ≃ 10^{4} cm^{−3}, n_{i} ≃ 10^{−2} cm^{−3}, μ_{i} ≃ 29 (Xu et al. 2016) and a simulation box L ≃ 0.1 pc with q = 1/2, we find D_{crit} ≃ 10^{22} cm^{2} s^{−1}^{2} and D_{∥} ≃ 6.7 × 10^{25} cm^{2} s^{−1}. A dynamical effect of CRs over jet environments seems then to require a very large 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. Gammaray diffuse emission
Gammaray 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 FermiLAT 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 gammaray 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 gammaray 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 gammaray signal in the outer galaxy is larger than expected from the standard CR diffusion model from known CR source distributions; this is the socalled CR gradient problem. The GeV gammaray 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 starforming 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). Gammaray 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 timedependent 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 gammaray observations and properties of CRmodified 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 gammaray 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 subAlfvénic, in particular within the WNM. As a consequence, we may underestimate the cosmicray propagation in the low density regions. These aspects will be investigated in a future work.
Another limitation comes from the rather idealised setup 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 selfconsistent 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 selfgravity, 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 starforming 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^{24} −10^{25} cm^{2} s^{−1} 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^{24} −10^{25} cm^{2} s^{−1}. 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 nonthermalsupport 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 CRpressuredominated 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^{28} cm^{2} s^{−1}. The turbulence generation by CR streaming especially around CR sources, as well as models for largescale 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 largescale turbulent models, as well as subgrid 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. ANR16IDEX0005) 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^{3} 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_{0} = 10^{24} and D_{0} = 10^{26} cm^{2} s^{−1}. We perform four calculations with D_{0} = 10^{24} cm^{2} s^{−1} and a resolution ranging from 64^{3} to 512^{3}, and three calculations with D_{0} = 10^{26} cm^{2} s^{−1} and a resolution ranging from 64^{3} to 256^{3}. We do not present results at a resolution of 512^{3} with D_{0} = 10^{26} cm^{2} s^{−1} 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_{0} = 10^{24} cm^{2} s^{−1}. 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^{3} 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 subsampling 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^{−3} (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 D_{0} = 10^{24} cm^{2} s^{−1} in the plane z = 20 pc at a time corresponding to 1.2t_{cross} and a resolution of 64^{3}, 128^{3}, 256^{3}, and 512^{3} (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_{0} = 10^{24} cm^{2} s^{−1} do not show convergence as resolution increases. The gas fraction in the CNM increases significantly in the 512^{3} 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^{3}, 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 (k_{turb} = 2) with D_{0} = 10^{24} cm^{2} s^{−1} (blue) and D_{0} = 10^{26} cm^{2} s^{−1} (red), but at times greater than 1t_{cross}. 
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^{22} to 10^{28} cm^{2} s^{−1}, 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_{0} = 10^{22} and D_{0} = 10^{28} cm^{2} s^{−1}, 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_{0} = 10^{26} cm^{2} s^{−1}. The PDF is indeed narrower with isotropic diffusion because the diffusion process is more efficient. With D_{0} = 10^{24} cm^{2} s^{−1}, 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.
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 718, 348 [NASA ADS] [CrossRef] [Google Scholar]
 Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 223, 26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ackermann, M., Ajello, M., Baldini, L., et al. 2011, ApJ, 726, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Allafort, A., et al. 2012a, ApJ, 756, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Allafort, A., et al. 2012b, ApJ, 755, 22 [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]
 Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Audit, E., & Hennebelle, P. 2010, A&A, 511, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJ, 777, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Burkhart, B., Lee, M.Y., Murray, C. E., & Stanimirović, S. 2015, ApJ, 811, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Butsky, I., Zrake, J., Kim, J.h., Yang, H.I., & Abel, T. 2017, ApJ, 843, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Ceccarelli, C., Dominik, C., LópezSepulcre, A., et al. 2014, ApJ, 790, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, J., Lazarian, A., & Vishniac, E. T. 2002, ApJ, 564, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Crutcher, R. M. 2012, ARA&A, 50, 29 [NASA ADS] [CrossRef] [Google Scholar]
 de Gouveia Dal, E. M., Piovezan, P. P., & Kadowaki, L. H. S. 2010, A&A, 518, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubois, Y., & Commerçon B. 2016, A&A, 585, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eswaran, V., & Pope, S. B. 1988, Comput. Fluids, 16, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Everett, J. E., & Zweibel, E. G. 2011, ApJ, 739, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., RomanDuval, J., Klessen, R. S., Schmidt, W., & Mac Low M.M. 2010, A&A, 512, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferrière, K., & Schmitt, D. 2000, A&A, 358, 125 [NASA ADS] [Google Scholar]
 Field, G. B. 1965, ApJ, 142, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaches, B. A. L., & Offner, S. S. R. 2018, ApJ, 861, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Galtier, S., Nazarenko, S. V., Newell, A. C., & Pouquet, A. 2000, J. Plasma Phys., 63, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Gazol, A., & Kim, J. 2013, ApJ, 765, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJ, 816, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Grenier, I. A., Black, J. H., & Strong, A. W. 2015, ARA&A, 53, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasz, M., Wóltański, D., & Kowalik, K. 2009, ApJ, 706, L155 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasz, M., Lesch, H., Naab, T., et al. 2013, ApJ, 777, L38 [NASA ADS] [CrossRef] [Google Scholar]
 Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P. 2018, A&A, 611, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Pérault, M. 1999, A&A, 351, 309 [NASA ADS] [Google Scholar]
 Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Ivlev, A. V., Padovani, M., Galli, D., & Caselli, P. 2015, ApJ, 812, 135 [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]
 Kritsuk, A. G., Norman, M. L. P., & Wagner, R. 2007, ApJ, 665, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1979, MNRAS, 186, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1981, MNRAS, 194, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Lemaster, M. N., & Stone, J. M. 2008, ApJ, 682, L97 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., & Goldreich, P. 2001, ApJ, 562, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M.M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MivilleDeschênes, M.A., Murray, N., & Lee, E. J. 2017, ApJ, 834, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Molina, F. Z., Glover, S. C. O., Federrath, C., & Klessen, R. S. 2012, MNRAS, 423, 2680 [NASA ADS] [CrossRef] [Google Scholar]
 Nava, L., & Gabici, S. 2013, MNRAS, 429, 1643 [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]
 Nolan, C. A., Federrath, C., & Sutherland, R. S. 2015, MNRAS, 451, 1380 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Hennebelle, P., Marcowith, A., & Ferrière, K. 2015, A&A, 582, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., & Springel, V. 2017, MNRAS, 465, 4500 [NASA ADS] [CrossRef] [Google Scholar]
 Ptuskin, V. S., Moskalenko, I. V., Jones, F. C., Strong, A. W., & Zirakashvili, V. N. 2006, ApJ, 642, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Rab, C., Güdel, M., Padovani, M., et al. 2017, A&A, 603, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieder, M., & Teyssier, R. 2017, MNRAS, 471, 2674 [NASA ADS] [CrossRef] [Google Scholar]
 RodgersLee, D., Taylor, A. M., Ray, T. P., & Downes, T. P. 2017, MNRAS, 472, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Rogachevskii, I., Kleeorin, N., Brandenburg, A., & Eichler, D. 2012, ApJ, 753, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Roy, N., Peedikakkandy, L., & Chengalur, J. N. 2008, MNRAS, 387, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Ruszkowski, M., Yang, H. Y. K., & Zweibel, E. 2017, ApJ, 834, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312 [NASA ADS] [CrossRef] [Google Scholar]
 Saury, E., MivilleDeschênes, M.A., Hennebelle, P., Audit, E., & Schmidt, W. 2014, A&A, 567, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Berlin: Springer) [CrossRef] [Google Scholar]
 Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2006, Comput. Fluids, 35, 353 [CrossRef] [Google Scholar]
 Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seifried, D., Schmidt, W., & Niemeyer, J. C. 2011, A&A, 526, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Skilling, J. 1971, ApJ, 170, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annu. Rev. Nucl. Sci., 57, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Uchiyama, Y., Funk, S., Katagiri, H., et al. 2012, ApJ, 749, L35 [NASA ADS] [CrossRef] [Google Scholar]
 VazquezSemadeni, E. 1994, ApJ, 423, 681 [NASA ADS] [CrossRef] [Google Scholar]
 Wiener, J., Zweibel, E. G., & Oh, S. P. 2013, ApJ, 767, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467, 906 [NASA ADS] [Google Scholar]
 Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 [NASA ADS] [CrossRef] [Google Scholar]
 Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 476, 2063 [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]
 Yan, H., & Lazarian, A. 2008, ApJ, 673, 942 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., Lazarian, A., & Schlickeiser, R. 2012, ApJ, 745, 140 [NASA ADS] [CrossRef] [Google Scholar]
We still assume that Larson laws apply at a scale of 0.1–0.01 pc and we use Eq. (2).
All Tables
All Figures
Fig. 1. Alfvénic Mach number 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 k_{turb} = 2 of the clumping factor of the gas density n (top panel), the cosmicray energy density E_{CR} (middle panel), and the Alfvénic Mach number (bottom panel). The clumping factor C(X) of a quantity X is defined as C(X) = ⟨X^{2}⟩∕⟨X^{⟩2}. 

In the text 
Fig. 3. Fiducial runs with k_{turb} = 2: maps of the gas density (top panels), cosmicray energy density (middle panels), and Alfvénic Mach number (bottom panels) in the plane z = 20 pc at a time corresponding to 2t_{cross}. From left to right panels: CR diffusion coefficient is D_{0} = 10^{22}, 10^{24}, 10^{26}, and 10^{28} cm^{2} 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 xyplane. 

In the text 
Fig. 4. Twodimensional histograms of, from top to bottom panels, the CR energy density, the Alfvénic Mach number , the ratio , and ζ = P_{CR}∕P as a function of the density for the fiducial runs k_{turb} = 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 (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 k_{turb} = 2 at the same time as in Fig. 3, for D_{0} = 10^{22} (purple), 10^{24} (blue), 10^{26} (red), and 10^{28} cm^{2} s^{−1} (black). Two models with D_{0} = 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 k_{turb} = 2 and D_{0} = 10^{24} (blue), D_{0} = 10^{26} cm^{2} 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 dasheddotted line ζ = 100. 

In the text 
Fig. 7. Same as Fig. 5 for k_{turb} = 4. 

In the text 
Fig. 8. Same as Fig. 5 for k_{turb} = 8. 

In the text 
Fig. 9. Same as Fig. 5 for k_{turb} = 16. 

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

In the text 
Fig. 11. Same as Fig. 5 for k_{turb} = 2, but different amplitudes of the initial magnetic field strength. The solid lines represent the fiducial case with , the dashed lines , the dotted lines , and the dasheddotted lines . The blue colour represents results with D_{0} = 10^{24} cm^{2} s^{−1}, the red D_{0} = 10^{26} cm^{2} s^{−1}. 

In the text 
Fig. A.1. Gas density (top panels) and CR energy density (bottom panels) maps for the fiducial runs with D_{0} = 10^{24} cm^{2} s^{−1} in the plane z = 20 pc at a time corresponding to 1.2t_{cross} and a resolution of 64^{3}, 128^{3}, 256^{3}, and 512^{3} (from left to right). 

In the text 
Fig. A.2. Same as Fig. 5 for different numerical resolutions in the fiducial case (k_{turb} = 2) with D_{0} = 10^{24} cm^{2} s^{−1} (blue) and D_{0} = 10^{26} cm^{2} s^{−1} (red), but at times greater than 1t_{cross}. 

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

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.