© Y. Dubois et al. 2019

1 Introduction

Cosmic rays (CR) are understood to play an important role in astrophysical plasmas due to their capacity to ionise the interstellar matter (Padovani et al. 2009) and their non-negligible pressure support to gas dynamics according to evolutionary processes that differ substantially from the thermal component since they diffuse efficiently and have different dissipation timescales. CRs are likely produced at shocks through the process of diffuse shock acceleration (DSA; see Bell 1978; Drury 1983; Blandford & Eichler 1987; Jones & Ellison 1991; Berezhko & Ellison 1999 and Marcowith et al. 2016 for a recent review). Recent advances in the numerical modelling of DSA through hybrid particle-in-cell codes (Caprioli & Spitkovsky 2014) have provided accurate predictions about the amount of CRs injected at shocks as a function of various properties of the shock including the Mach number, the obliquity of the magnetic field, or the pre-existing amount of CRs (Caprioli et al. 2018). There is a large body of evidence for CRs accelerated in the shocked-shell material of supernova (SN) explosions (e.g. Koyama et al. 1995; Decourchelle et al. 2000; Aharonian et al. 2004; Warren et al. 2005; Helder et al. 2009; Ackermann et al. 2013) and it has been shown that they have a significant impact on the shell structure and dynamics (Chevalier 1983; Dorfi 1990; Zank et al. 1993; Wagner et al. 2009; Ferrand et al. 2010; Castro et al. 2011; Pfrommer et al. 2017; Pais et al. 2018; Diesing & Caprioli 2018). Supernova remnants (SNRs) are expected to be the main source of CRs permeating the entire interstellar medium (ISM) of galaxies (Aguilar et al. 2015), though the consistency of the accelerated CR spectrum in a SNR with that of entire galaxies is still intensely debated (see Blasi 2013 for a review).

Cosmic rays likely have an important dynamical impact over the ISM on all galactic scales. On small scales, while released by a SNR, CRs possess enough pressure to overcome the background magnetic and gas pressures and trigger different types of plasma instabilities which result in the production of waves and turbulence (Ptuskin et al. 2008; Malkov et al. 2013). This self-generated turbulence can confine CRs over distances and amounts of time that depend on the conditions prevailing in the ISM, especially the ionisation degree (Nava et al. 2016, 2019). The generation of waves contribute to locally heating the warm ionised medium (Wiener et al. 2013b). On larger galactic scales, comparable to the disc height, CR gradients can modify the dynamics of Jeans unstable regions in the atomic phase (Commerçon et al. 2019), and they can propel cold galactic-wide outflows (Jubelgas et al. 2008; Wadepuhl & Springel 2011; Uhlig et al. 2012; Hanasz et al. 2013; Salem & Bryan 2014; Salem et al. 2014; Girichidis et al. 2016, 2018; Simpson et al. 2016; Recchia et al. 2017; Fujita & Mac Low 2018; Mao & Ostriker 2018) with a preferential impact in low-mass galaxies (Booth et al. 2013; Jacob et al. 2018, Dashyan & Dubois 2019). However, the capability of winds to carry mass and momentum depends on the detailed CR physics such as streaming (Ruszkowski et al. 2017b; Wiener et al. 2017; Holguin et al. 2019; Butsky & Quinn 2018), or taking into account the unresolved multi-phase nature of the gas and its impact on CR transport (Farber et al. 2018). CRs also boost the dynamo amplification of the magnetic field in disc galaxies (Hanasz et al. 2004, 2009a,b; Pakmor et al. 2016).

On very large cosmological scales CRs are released in shocks (Miniati et al. 2000, 2001; Ryu et al. 2003; Skillman et al. 2008; Pfrommer et al. 2007, 2008, 2017; Vazza et al. 2009, 2012) with external cosmological infall of gas producing the strongest shocks, while pre-processed internal shocks in halos drive the bulk of the shock distribution in the more moderate strength regime.

Similarly, strong shocks are produced in jets from active galactic nuclei; they release large amounts of CRs as observed in radio emission (Fanaroff & Riley 1974; Pierre Auger Collaboration 2007; Croston et al. 2009) and help to release the feedback back to the hot gas from galaxy clusters (Croston et al. 2008; Guo & Oh 2008; Sijacki et al. 2008; Guo & Mathews 2011; Fujita & Ohira 2011; Jacob & Pfrommer 2017; Ruszkowski et al. 2017a; Ehlert et al. 2018). However, again, their impact might significantly differ depending on which CR dynamical processes are modelled and which ignored.

In a previous work (Dubois & Commerçon 2016), we introduced a numerical model for anisotropic CR diffusion. Here, we extend it byincluding a model of the CR streaming instability and CR injection at shocks through DSA in the adaptive mesh refinement code RAMSES (Teyssier 2002). In another work (Brahimi et al. in prep.) we introduce new diffusive transport for CRs accounting for the generation of turbulence produced by the streaming. This ensemble of work aims to provide a consistent description of CR dynamical effect on the interstellar or intergalactic media. In the same view, a recent model has been proposed by Thomas & Pfrommer (2019).

In Sect. 2, we introduce the full set of CR magneto-hydrodynamics including the streaming and acceleration terms, whose numerical modelling and tests are respectively tackled in Sects. 3 and 4. We finally test CR acceleration and streaming combined in turbulent ISM experiments in Sect. 5.

2 Magneto-hydrodynamics with cosmic rays

By taking the energy moment of the Fokker–Planck CR transport equation (Drury & Voelk 1981), the following set of differential equations tobe solved for cosmic-ray magneto-hydrodynamics (CRMHD) of a fluid mixture made of thermal particles and CRs can be obtained: $$\begin{array}{lll}\hfill & \hfill & \frac{\partial \rho }{\partial t}+\nabla .(\rho u)=0,\hfill \\ \hfill & \hfill & \frac{\partial \rho u}{\partial t}+\nabla .\left(\rho uu+P_{\text{tot}}-\frac{BB}{4\pi }\right)=0,\hfill \\ \hfill & \hfill & \frac{\partial e}{\partial t}+\nabla .\left((e+P_{\text{tot}})u-\frac{B(B.u)}{4\pi }\right)\hfill \\ \hfill & \hfill & =-P_{\text{CR}}\nabla .u-\nabla .F_{\text{CR},\text{d}}+{\mathcal{L}}_{\text{rad}},\hfill \\ \hfill & \hfill & \frac{\partial B}{\partial t}-\nabla \times (u\times B)=0,\hfill \\ \hfill & \hfill & \frac{\partial e_{\text{CR}}}{\partial t}+\nabla .\left(e_{\text{CR}}u+(e_{\text{CR}}+P_{\text{CR}})u_{\text{st}}\right)\hfill \\ \hfill & \hfill & =-P_{\text{CR}}\nabla .u-\nabla .F_{\text{CR},\text{d}}+{\mathcal{L}}_{\text{st}}+{\mathcal{H}}_{\text{acc}}+{\mathcal{L}}_{\text{rad},\text{CR}}.\hfill \end{array}$$

Here ρ is the gas mass density, u is the gas velocity, u_st is the streaming velocity, B is the magnetic field, e = 0.5ρu² + e_th + e_CR + B²∕8π is the total energy density, e_th is the thermal energy density, and e_CR is the CR energy density; P_tot = P_th + P_CR + P_mag is the sum of thermal P_th = (γ − 1)e_th, CR P_CR = (γ_CR − 1)e_cr, and magnetic P_mag = 0.5B²∕(4π) pressures, where γ and γ_CR are the adiabatic indexes of the thermal and CR components, respectively. We note that all energy components e_i are energies per unit volume e_i = E_i∕Δx³, where Δx is the cell size. The terms on the right-hand side of the equations are treated as source terms with P_CR ∇.u the CR pressure work term, F_CR,d = −D₀b(b.∇e_CR) the anisotropic diffusion flux term, D₀ the diffusion coefficient (usually taken as a constant value for simplicity, but it can also be a function of local MHD quantities), b = B∕||B|| the magnetic unity vector, and a total radiative loss term ${\mathcal{L}}_{\text{rad}}={\mathcal{L}}_{\text{rad},\text{th}}+{\mathcal{L}}_{\text{rad},\text{CR-}>\text{th}}$ composed of the thermal ${\mathcal{L}}_{\text{rad},\text{th}}$ and CR ${\mathcal{L}}_{\text{rad},\text{CR-}>\text{th}}$ radiative loss terms, where the CR loss term (${\mathcal{L}}_{\text{rad},\text{CR-}>\text{th}}={\mathcal{L}}_{\text{rad},\text{CR}}+{\mathcal{H}}_{\text{rad},\text{CR-}>\text{th}}$) is the non-conserving sum of radiative losses from cosmic rays ${\mathcal{L}}_{\text{rad},\text{CR}}$ turning as a heating rate ${\mathcal{H}}_{\text{rad},\text{CR-}>\text{th}}$ for the thermal component. Finally, and this is the core of this paper, we detail how the streaming instability terms $\nabla .\left((e_{\text{CR}}+P_{\text{CR}})u_{\text{st}}\right)$ (advection-diffusion term) and ${\mathcal{L}}_{\text{st}}$ (heating term), and the CR acceleration at shocks ${\mathcal{H}}_{\text{acc}}$ are modelled.

We use the RAMSES code detailed in Teyssier (2002) to solve these equations with adaptive mesh refinement (AMR). The full set of equations is solved with the standard MHD solver of RAMSES described in Fromang et al. (2006), where the right-hand side terms of Eq. (3) are treated separately as source terms. The induction equation (Eq. (4)) is solved using constrained transport (Teyssier et al. 2006), which by construction guarantees at all times that ∇.B ≃ 0 at machine precision. Godunov fluxes are solved with the approximate Harten–Lax–van Leer Discontinuities (HLLD) Riemann solver (Miyoshi & Kusano 2005) and the minmod total variation diminishing slope limiter are modified to account for the extra energy components and total pressure made of the thermal and CR component. Accordingly, the effective sound speed used for the Courant–Friedrichs–Lewy time-step condition accounts for the extra pressure components (i.e. total pressure of the fluid). The implementation of the anisotropic CR diffusion in RAMSES, which our new implementation of CR streaming relies on, is described in Dubois & Commerçon (2016).

It should be noted that Eq. (5) can be expanded to as many CR energy bins as required to sample a full spectrum of CRs in energy-momentum space with source terms communicating the energy fluxes between the various energy bins (see Miniati 2001; Girichidis et al. 2014; Winner et al. 2019, for such efforts in those directions). We ignore this extra level of complexity to represent the entire spectrum of CR energy by a single bin of energy. For sake of completeness, we introduced the anisotropic diffusion term as well as the CR radiative loss terms (trivially modelled as a simple density and CR energy-dependent term; see e.g. Enßlin et al. 2007; Guo & Oh 2008) in the equations; we do not make use of them in the various tests of this paper, i.e. D₀ = 0 and ${\mathcal{L}}_{\text{rad},\text{CR}}=0$.

3 Cosmic-ray streaming

3.1 Numerical implementation

Cosmic rays propagating faster than the Alfvén velocity $u_{\text{A}}=B/\text{}\sqrt{4\pi \rho }$ excite Alfvén waves, which in turn drive the scattering of the CR pitch angle with magnetic field lines. This coupling leads to a reduced CR bulk velocity at the Alfvén velocity and confines the CR streaming transport along the field lines and their own gradient of pressure (Wentzel 1968; Kulsrud & Pearce 1969; Skilling 1975). Several damping mechanisms, such as ion-neutral damping, non-linear Landau damping, or turbulence damping (Kulsrud & Pearce 1969; Yan & Lazarian 2002; Farmer & Goldreich 2004; Lazarian & Beresnyak 2006; Wiener et al. 2013a), can lead to a significant suppression of these self-excited Aflvén waves and can increase the effective value at which CRs are allowed to stream down their own gradient at super-Alfvénic velocities u_st = −f_SAu_Asign(b.∇e_CR), where f_SA ≥ 1 is the super-Aflvénic boost factor of the streaming velocity.

In addition, while CRs scatter onto the Aflvén waves, they experience a drag force, whose work is transferred to the thermal pool at the following rate: $${\mathcal{L}}_{\text{st}}=-\text{sign}(b.\nabla e_{\text{CR}})u_{\text{A}}.\nabla P_{\text{CR}}.$$(6)

We note that this heating term has f_SA = 1 since only the Alfvén waves mediate the energy exchange between CRs and the thermal component (see e.g. Ruszkowski et al. 2017b). This term, which is by construction always a heating (resp. loss) term for the thermal (resp. CR) component, is obtained by simply differentiating the values of the CR energy density with neighbouring cells.

For simplicity, in the rest of this work, whose aim is to test the implementation of CR streaming, we systematically assume f_SA = 1. The advection or diffusion term of streaming ∇.((e_CR + P_CR)u_st) can be solved via two distinct approaches. One is to update the CR energy density using an explicit upwind method; however, since the streaming velocity can become discontinuous at extrema of e_CR, it modifies the condition of stability of the solution to Δt ∝ Δx³ (Sharma et al. 2009). Sharma et al. (2009) proposed regularising the streaming velocity by replacing sign (b.∇e_CR) by tanh (hb.∇e_CR∕e_CR) in order to obtain a less constraining time-step condition of Δt = hΔx²∕(2e_CRu_A), and where h should be the size of a few cells. Nonetheless, this time-step condition is still too constraining due to the quadratic dependency on cell size, and it is necessary to rely on a different strategy to make such a numerical implementation practicable in all possible situations. Sharma et al. (2009) suggested using an implicit solver for the regularised upwind method. Here we decided to take a different route that relies on the modelling of the anisotropic diffusion with an implicit solver, as done in Dubois & Commerçon (2016).

We can rewrite the streaming velocity as $$u_{\text{st}}=-\frac{b.\nabla e_{\text{CR}}}{|b.\nabla e_{\text{CR}}|}{u}_{\text{A}},$$(7)

which, whenrecast into ∇.((e_CR + P_CR)u_st), can be rewritten as a diffusion term (see also Uhlig et al. 2012, where the same diffusion approach for the isotropic version of CR streaming is used): $$\begin{array}{lll}\nabla .F_{\text{CR},\text{s}}\hfill & =\hfill & \nabla .(-D_{\text{st}}b(b.\nabla e_{\text{CR}}))\hfill \\ \hfill & =\hfill & \nabla .\left(-\frac{(e_{\text{CR}}+P_{\text{CR}})\left|B\right|}{|b.\nabla e_{\text{CR}}|\sqrt{4\pi \rho }}b(b.\nabla e_{\text{CR}})\right).\hfill \end{array}$$(8)

Therefore, this advection-diffusion part of the streaming instability can be treated as an addition to the standard F_CR,d CR diffusion term (F_{CR, ds} = F_{CR, d} + F_{CR, s}), for clarity hereafter written as follows: $$\nabla .F_{\text{CR},\text{ds}}=\nabla .\left(-Db(b.\nabla e_{\text{CR}})\right),$$(9)

where D = D₀ + D_st. The F_CR,ds diffusion flux can be arbitrarily decomposed into an anisotropic and isotropic part $$\nabla .F_{\text{CR},\text{ds}}=\nabla .\left(-D_{\parallel }b(b.\nabla e_{\text{CR}})-D_{\text{iso}}\nabla e_{\text{CR}}\right),$$(10)

where D_∥ = (1 − f_iso)D, D_iso = f_isoD, and f_iso ≤ 1. We briefly recallthe framework of the implicit solver developed in Dubois & Commerçon (2016). For the 2D case, the time update of the CR energy by the anisotropic part (the isotropic part is trivially obtained) of the diffusion flux is $$e_{i,j}^{n+1}+\Delta t\frac{F_{i+\frac{1}{2},j}^{n+1}+F_{i,j+\frac{1}{2}}^{n+1}-F_{i-\frac{1}{2},j}^{n+1}-F_{i,j-\frac{1}{2}}^{n+1}}{\Delta x}=e_{i,j}^n,$$(11)

where the cell-centred fluxes are computed with cell-cornered values using the symmetric scheme from Günter et al. (2005): $$\begin{array}{lll}{F}_{i+\frac{1}{2},j}^{\text{ani}}\hfill & =\hfill & \frac{F_{i+\frac{1}{2},j-\frac{1}{2}}^{\text{ani}}+F_{i+\frac{1}{2},j+\frac{1}{2}}^{\text{ani}}}{2},\hfill \\ F_{i,j+\frac{1}{2}}^{\text{ani}}\hfill & =\hfill & \frac{F_{i-\frac{1}{2},j+\frac{1}{2}}^{\text{ani}}+F_{i+\frac{1}{2},j+\frac{1}{2}}^{\text{ani}}}{2}.\hfill \end{array}$$

The anisotropic cell corner flux is $$F_{i+\frac{1}{2},j+\frac{1}{2}}^{\text{ani}}=-\overline{D}{\overline{b}}_x\left({\overline{b}}_x\overline{\frac{\partial e}{\partial x}}+{\overline{b}}_y\overline{\frac{\partial e}{\partial y}}\right),$$(12)

where barred quantities are arithmetic averages over the cells connected to the corner, i.e. $$\begin{array}{lll}\hfill & \hfill & {\overline{b}}_x=\frac{b_{x,i+\frac{1}{2},j}^n+b_{x,i+\frac{1}{2},j+1}^n}{2},\hfill \\ \hfill & \hfill & {\overline{b}}_y=\frac{b_{y,i,j+\frac{1}{2}}^n+b_{y,i+1,j+\frac{1}{2}}^n}{2},\hfill \\ \hfill & \hfill & \overline{\frac{\partial e}{\partial x}}=\frac{e_{i+1,j+1}^{n+1}+e_{i+1,j}^{n+1}-e_{i,j+1}^{n+1}-e_{i,j}^{n+1}}{2\Delta x},\hfill \\ \hfill & \hfill & \overline{\frac{\partial e}{\partial y}}=\frac{e_{i+1,j+1}^{n+1}+e_{i,j+1}^{n+1}-e_{i+1,j}^{n+1}-e_{i,j}^{n+1}}{2\Delta x},\hfill \\ \hfill & \hfill & \overline{D}=\frac{D_{i,j}^n+D_{i+1,j}^n+D_{i,j+1}^n+D_{i+1,j+1}^n}{4}.\hfill \end{array}$$

We notethat all hydrodynamical variables in RAMSES are cell-centred except for the magnetic field which is face-centred. The streaming diffusion coefficient is computed as $$D_{i,j}^n=\frac{(e_{i,j}^n+P_{i,j}^n)}{\sqrt{4\pi {\rho }_{i,j}}}\frac{|\tilde{B}{|}_{i,j}^n}{{\widetilde{|b.\nabla e|}}_{i,j}^n},$$(13)

where upper tilde quantities stand for cell-centred quantities reconstructed from a combination of cell-centred and face-centred quantities: $$\begin{array}{lll}\hfill & \hfill & |\tilde{B}{|}_{i,j}^n=\frac{1}{2}\sqrt{{\left(B_{x,i-\frac{1}{2},j}^n+B_{x,i+\frac{1}{2},j}^n\right)}^2+{\left(B_{y,i,j-\frac{1}{2}}^n+B_{y,i,j+\frac{1}{2}}^n\right)}^2},\hfill \\ \hfill & \hfill & {\widetilde{|b.\nabla e|}}_{i,j}^n=\frac{1}{4\Delta x|\tilde{B}{|}_{i,j}^n}|\left(B_{x,i-\frac{1}{2},j}^n+B_{x,i+\frac{1}{2},j}^n\right)\left(e_{x,i+1,j}^n-e_{x,i-1,j}^n\right)\hfill \\ \hfill & \hfill & +\left(B_{y,i,j-\frac{1}{2}}^n+B_{y,i,j+\frac{1}{2}}^n\right)\left(e_{y,i,j+1}^n-e_{y,i,j-1}^n\right)|.\hfill \end{array}$$(14)

It should be noted that, in principle, the solver can deal with any arbitrary large values of the diffusion coefficient; however, the number of iterative steps of the implicit solver to converge towards the solution can be large for a large diffusion coefficient, typically at extrema of |b.∇e_CR| where this value can become close to zero. In practice, we cap the value of the streaming diffusion coefficient to 10²⁸ cm² s⁻¹ in all practical astrophysical applications to reduce the spectral condition number of the matrix involved in the implicit solver in order to save computational iterations. From the 2D case, the method is trivially expanded into three dimensions.

Fig. 1 Evolution of a 1D sinusoid of CR energy density with streaming advection only as a function of position with 512 cells and imposing a constant Alfvén velocity of 1. The solution is made of two plateaus as the maxima are capped over time, due to the infinite streaming diffusion coefficient, while the two regions between the two plateaus move at a velocity of ± γCR = ±1.4.

3.2 Tests of CR streaming

3.2.1 One-dimensional sinusoid

In order to test the implementation of the CR advection-diffusion streaming term, a 1D sinusoid experiment is set up where the rest of the physics is deactivated, and with γ_CR = 1.4 similar to the test proposed by Sharma et al. (2009). Unfortunately, there is no known analytical solution to that experiment, but we can test the numerical convergence of the implementation to test its self-consistency. The initial condition for CR energy density is e_CR = 1 + 0.5sin(2πx), and we assume that the Alfvén velocity equals 1 oriented along the x-axis. In this 1D test we set the maximum streaming diffusion coefficient to be no larger than 100. As shown in Fig. 1 for this 1D test problem using 512 cells (level 9), the evolved solution is a sinusoid where the extrema are cropped and where the regions of maximum slope are advected at γ_CR u_st (i.e. − 1.4 if ∂E_CR∕∂x > 0 and + 1.4 if ∂E_CR∕∂x < 0). A more evolved time shows a higher cropped fraction of the high and low part of the sinusoid. We perform a consistency test by varying the resolution of the simulation from 16 cells to 1024 cells, where the highest resolution simulation is used as a reference for comparison. Figure 2 shows the solution at time t = 0.02 for 16, 32, 64, 128, 256, and 1024 cells, and their relative variation to the reference run. The solution shows very good numerical convergence towards the high-resolution reference solution, which never exceeds a few percentage points relative variation even when using only 16 cells to resolve the wavelength of the sinusoid. Finally, the L2 norm (again using the 1024-cell run as a reference) is computed and has a convergence with a scaling of Δx^1.87±0.08, as shown in Fig. 3.

3.2.2 Two-dimensional sinusoid in a looped magnetic field

In this test case we try to mimic the 1D sinusoid problem embedded in a non-uniform magnetic configuration. We initialise a 2D looped magnetic field centred on the middle of the box, hence in the circular coordinate system the magnetic field is purely tangential. We also initialise the CR energy density in the same way as the previous 1D test case with a θ angle dependency e_CR = 1 + 0.5sin(θ) for a radius 0.15 < r < 0.35 and e_CR = 10⁻⁵ for r ≤ 0.15 and r ≥ 0.35. In this 2D test we set the maximum streaming diffusion coefficient to be no larger than 1 and an isotropic component of f_iso = 10⁻²; we discuss the effect of changing these values on the solution in Appendix A. We choose an Alfvén velocity of 1, and again we deactivate the rest of the hydrodynamics. Figure 4 shows the result at times t = 0 and t = 0.02. The solution shows a similar angle-dependent pattern for the evolved solution at t = 0.02 to that of the 1D case at the same time (i.e. the value of energy density is close to uniform around regions of initial extrema). We note that the capping of extrema is slightly late in this 2D configuration with respect to the 1D test: compared with Fig. 1, where the maximum and minimum are respectively 1.2 and 0.7 at time t = 0.02, here in 2D we obtain 1.28 and 0.6, respectively. We also tested the 2D streaming for a ten times wider range of initial CR energy density. Theresult, not shown here, is qualitatively similar to that of our reference test.

Fig. 2 Solution at t = 0.02 of the sinusoid experiment with different resolution from uniform level 4 to 8 (from light red to dark red) and level 10 (in black). The relative errors are compared to the reference numerical solution of level 10. Even for very low resolution the relative error is never larger than a few percentage points.

Fig. 3 Convergence of the L2 norm for the sinusoid experiment using the solution at t = 0.02. The norm is compared to the reference numerical solution of level 10. The L2 norm scales with Δx1.87±0.08 as indicated by the dashed line.

Fig. 4 Cosmic-ray energy density maps at t = 0 (top left) and t = 0.02 (top right) for an initial angle-dependent sinusoid within a purely circular magnetic field with an Alfvén velocity of 1 and a resolution of 1282 cells. The energy is evolved with the streaming advection-diffusion term only. The bottom panel shows the radially averaged energy in theradius interval r = [0.15, 0.35] as a functionof the polar angle θ.

4 Shock-accelerated CRs

4.1 Shock finder algorithm

Our shock finder algorithm relies on several criteria. A shock cell is identified as such when all of the following conditions are met: (i) ∇T.∇S > 0 (Ryu et al. 2003, where S = T∕n^2∕3 is the pseudo-entropy) and ∇T.∇ρ > 0 (which filters out tangential discontinuities, Schaal & Springel 2015); (ii) ∇.u is negative (compression region); (iii) ∇.u is a local minimum along the normal to n_s = −∇T∕|∇T| (where the local value of ∇.u is compared to the cloud-in-cell interpolated value of ∇.u at one Δx local cell distance in the upstream and downstream of the local cell); and (iv) the Mach number is larger $\mathcal{M}>{\mathcal{M}}_{\text{min}}$, with ${\mathcal{M}}_{\text{min}}\simeq 1.5$. Keeping in mind these conditions, the Mach number of eligible cells is computed according to the criteria using upstream (pre-shock) and downstream (post-shock) fluid variables. Using the Rankine-Hugoniot shock jump relations, the Mach number can be computed from density, temperature, or pressure values. For instance, the Mach number for a single thermal component can be obtained from the ratio ${\mathcal{R}}_{\text{P}}=P_{\text{2}}/P_{\text{1}}$ of the downstream to upstream pressures (here and in the following we keep the 1 and 2 subscripts for the upstream and downstream quantities), leading to $${\mathcal{M}}^2=\frac{1}{2\gamma }\left[(\gamma -1)+(\gamma +1){\mathcal{R}}_{\text{P}}\right]\text{}.$$(15)

We note that it is also possible to employ the jump relations for density or velocity; however, they quickly saturate at high Mach numbers, while pressure jumps offer better leverage for probing the values of the Mach number.

Since our aim is to apply this shock finder to a thermal–CR mixture, the following relation (Pfrommer et al. 2017) should be used instead: $${\mathcal{M}}^2=\frac{1}{{\gamma }_{\text{e}}}\frac{{\mathcal{R}}_{\text{P}}\mathcal{C}}{\mathcal{C}-\left[({\gamma }_1+1)+({\gamma }_1-1){\mathcal{R}}_{\text{P}}\right]({\gamma }_2-1)}.$$(16)

Here $\mathcal{C}=[({\gamma }_2+1){\mathcal{R}}_{\text{P}}+({\gamma }_2-1)]({\gamma }_1-1)$, γ_i = P_i∕ϵ_i + 1 for i = {1, 2} (respectively upstream and downstream) and γ_e = (γP_th,2 + γ_CRP_CR,2)∕P₂ for the downstream region. In the limit where the weighted adiabatic indexes are equal γ_e = γ₁ = γ₂ this formula for the Mach number is equal to the classical formulation of Eq. (15).

The normal to the shock is provided by the gradient of temperature n_s. A first guess of the upstream and downstream values of pressure are obtained by cloud-in-cell interpolating the values of the 2^D cell pressure (where D is the dimensionality of the system to simulate), one cell and two cells away from the shocked cell candidates along n_s and −n_s for the upstream and downstream quantities, respectively. The upstream and downstream pressures are respectively the minimum and maximum of pressures obtained from the one cell and two cell distances away from the shocked cell. This first guess of the Mach number is kept for cells with moderate Mach numbers $\mathcal{M}<5$, while cells with higher Mach numbers require probing regions further than two cells away from the shocked cell to properly evaluate their Mach numbers. As we see in the tests, the stronger the shock, the larger the number of cells to sample the discontinuity, and we thus need to probe more distant cells to accurately capture the true upstream and downstream values of the shock. This first guess is limited to two cells to fully exploit the code structure of RAMSES that tracks at each time the 3^D − 1 neighbouring octs of each cell (anoct contains 2^D cells), including virtual octs that belong to another domain (hence, going further away requires communication between CPU domains and can be prohibitive, which is why we limit this search to the strongest shocked cells).

The second guess of the Mach number, and other related quantities (see next section), is obtained by moving forward along the normal to the shock by steps of Δx up to four cells distance, thus probing both 3Δx and 4Δx in the upstream and the downstream regions. For the new value of upstream and downstream pressures (and other related quantities) to be accepted for the calculation of the new Mach number, we check that the slope of the thermal energy is getting shallower (the profile must flatten as we are moving outwards) by computing the new gradient of thermal energy and comparing to its value from the previous distance step, and that the total pressure and the density both have a new extremum (either an upstream minimum or a downstream maximum). Our experiments with Mach numbers as strong as 1000 has lead us to use up to four cells distance to probe the estimated Mach number of strong shocks, hence we always use this maximum value in the following, but our implementation can work with arbitrarily larger distances.

4.2 Cosmic-ray acceleration at shocks

At shocks the kinetic energy flux of the upstream flow ${\varphi }_{\text{K},\text{1}}=0.5{\rho }_{\text{1}}{u}_{\text{1}}^3$ (where the velocities are measured in the moving shock frame) is dissipated by the shock interface into a thermal energy flux ϕ_th,2 = e_th,dissu₂, CR energy ϕ_CR,2 = e_CR,dissu₂ and the remaining into kinetic and magnetic energy. For classical strong shocks without CR acceleration, the ratio of post-shock thermal (dissipated) energy to the pre-shock kinetic energy $e_{\text{th},\text{diss}}/(0.5\rho u_1^2)$ can be obtained from the Rankine–Hugoniot jump relations, and tends towards 0.56 for γ = 5∕3. Once shocked cells are identified, the amount of accelerated CRs is obtained with the CR flux following $${\varphi }_{\text{CR}}=\eta (\mathcal{M},X_{\text{CR}},{\theta }_{\text{B}})e_{\text{diss}}{u}_{\text{2}},$$(17)

where e_diss = e_th,diss + e_CR,diss is the dissipated internal energy of the gas, u₂ is the downstream velocity in the frame of the moving shock, and $\eta (\mathcal{M},X_{\text{CR}},{\theta }_{\text{B}})$ is the acceleration efficiency of CRs at shocks, which is a function of the Mach number, the upstream CR-to-thermal ratio X_CR = P_CR,1∕P_th,1, and the magnetic obliquity to the normal of the shock θ_B. Instead of measuring the downstream velocity in the shock frame (which requires knowing both the upstream and downstream velocities in the lab frame, as well as the jump density ratio ${\mathcal{R}}_{\rho }$), we replace u₂ by $\mathcal{M}{c}_{\text{s},\text{1}}/{\mathcal{R}}_{\rho }$, where c_s,1 is the upstream sound speed. The dissipated energy can be directly measured from the upstream and downstream thermal and CR energy densities $$e_{\text{diss}}=e_{\text{th},\text{2}}+e_{\text{CR},\text{2}}-e_{\text{th},\text{1}}{\mathcal{R}}_{\rho }^{\gamma }-e_{\text{CR},\text{1}}{\mathcal{R}}_{\rho }^{{\gamma }_{\text{CR}}},$$(18)

where e_{th, 2} and e_{th, 1} are respectively the downstream and upstream thermal energy densities, e_{CR, 2} and e_{CR, 1} the downstream and upstream CR energy densities, and ${\mathcal{R}}_{\rho }$ the jump density ratio. The jump density ratio is obtained from the direct evaluation of the upstream and downstream densities $${\mathcal{R}}_{\rho }=\frac{{\rho }_2}{{\rho }_1}.$$(19)

The ${\mathcal{R}}_{\rho }^{\gamma }$ and ${\mathcal{R}}_{\rho }^{{\gamma }_{\text{CR}}}$ terms account for the fact that the upstream thermal and CR energies are also adiabatically compressed at the shock. Finally, the new CR energy is updated using Δe_CR = ϕ_CRΔt∕Δx.

According to detailed simulations of accelerated CRs at shocks (Caprioli & Spitkovsky 2014), their acceleration efficiency depends on both the Mach number of the shock and the upstream magnetic field orientation with respect to the normal to the shock θ_B = arccos(b₁.n_s). The dependency of the efficiency of CR acceleration with this so-called “magnetic obliquity” can be factorised out, $\eta (\mathcal{M},X_{\text{CR}},{\theta }_{\text{B}})={\eta }_0\xi (\mathcal{M},X_{\text{CR}})\zeta ({\theta }_{\text{B}}),$ and approximated by the following functional form (Pais et al. 2018): $$\zeta ({\theta }_{\text{B}})=\frac{1}{2}\left[\tanh \left(\frac{{\theta }_{\text{crit}}-{\theta }_{\text{B}}}{{\delta }_{\theta }}\right)+1\right],$$(20)

where θ_crit = π∕4 and δ_θ = π∕18. Therefore, we probe the angle θ_B by evaluating the orientation of the magnetic vector in the upstream region using the cell that defines the value of the upstream pressure as defined in the previous section.

The dependency of the acceleration $\xi (\mathcal{M},X_{\text{CR}})$ is obtained from the results of Kang & Ryu (2013), and is an increasing function of both $\mathcal{M}$ and X_CR. They provide values of the acceleration efficiency fortwo values of X_CR, namely 0 and 0.05, and ten values of the Mach number (from 1.5 to 100). Since, to the best of our knowledge, no work has explored the cases with X_CR > 0.05, in order to explore the full range of admissible values of X_CR we simply interpolate and extrapolate the values of $\xi (\mathcal{M},X_{\text{CR}})$ from X_CR = 0 and 0.05, sampling values of X_CR = 0.025, 0.1, 0.2, 0.5, and 1. In addition, we fix those sampling values so that ξ is a monotonic increasing function of $\mathcal{M}$ and X_CR. We note that their obtained values of the acceleration efficiency saturates at η₀ = 0.225, a factor of ~ 2 larger than the maximum values obtained by Caprioli & Spitkovsky (2014) for parallel shocks (θ_B = 0). We thus renormalise $\xi (\mathcal{M},X_{\text{CR}})$ by 0.225 so that the maximum allowed efficiency is explicitly controlled by η₀. The values of ξ are shown in Fig. 5 and are available as tabulated values in Appendix B. We note that obliquity-dependent CR acceleration simulations conducted by Caprioli et al. (2018) with a pre-existing population of CRs in the upstream region suggest that the transition of the obliquity-dependent part of the efficiency ζ(θ_B) from the efficient to the inefficient regime is displaced from θ_crit = π∕4 to θ_crit = π∕3. We neglect this effect at the moment.

Finally, we decided to inject the CR energy accelerated at shocks a few cells away from the shock cell. We were guided by the fact that numerical shocks are not pure discontinuities and are in fact numerically broadened; therefore, any CR pressure deposited in the numerically broadened shock layer experiences a work P_CR∇.u of pressure forces. For this reason, the CR energy is deposited in the cell of minimum |∇.u| in the post-shock direction up to four cells away from the shock cell. We emphasize that this choice is crucial to obtaining the correct amount of CR energy density in the post-shock region, and our experiments have taught us that the direct injection in the shock systematically overestimates the resulting CR energy density in the post-shock region by a large factor even in the simplest 1D test case (e.g. by a factor of ~ 2 for the Sod test).

Fig. 5 Acceleration efficiency $\xi (\mathcal{M},X_{\text{CR}})$ as a functionof the Mach number $\mathcal{M}$ for different values of the upstream CR-to-thermal pressure ratio XCR. The values are obtained from the XCR = 0 and 0.025 values of Kang & Ryu (2013) and renormalised to a maximum value of 1.

4.3 One-dimensional Sod shock tube

4.3.1 Convergence of the shock Mach number

In this first test for the convergence of the evaluated shock Mach number, we used the standard Sod shock tube initial conditions for a Mach of 10; in other words, we started with initial left and right states separated by a virtual interface at x = 5 in a box of size of 10 with thermal pressure P_th,L = 63.499 and P_th,R = 0.1, density ρ_L = 1 and ρ_R = 0.125, velocity u_L = u_R = 0. This test was run without any initial or accelerated CR component (i.e. free of CR pressure), and we adopted an adiabatic index of the gas of 5∕3. In addition we also explored more aggressive shock tube initial conditions to probe Mach of 100 (P_{th, L} = 6349.9), and Mach of 1000 (P_{th, L} = 634 990). We employed a base grid of level 5 with up to three additional levels of refinements triggered in regions where the relative cell-to-cell variation of either the density, velocity, or pressure is larger than 10%.

Figure 6 shows the quality of the Mach number evaluation with the statistics of its value relative to the exact analytical value for various shock tube tests, changing the strength of the shock by two orders of magnitude. We tested two maximum values of the extent of the pre-shock and post-shock quantities, either probing up to ncell_max = 2 cells or ncell_max = 4 away from the shock cell. We note that we removed the estimates of the Mach number for the first 15 time steps of the simulations (over the 263 available time steps, reaching final times t = 0.35, t = 0.035, and t = 0.0035 for Mach numbers of 10, 100, and 1000, respectively), where the shock, contact, and rarefaction waves are not yet sufficiently separated to correctly capture the Mach number of the shock. It shows that ncell_max = 2 cells can be sufficient to obtain Mach numbers accurate to a level of a few percentage points up to Mach numbers of the order of ~ 100, even though it is systematically underevaluated; however, Mach numbers of 1000 are almost never correctly captured. On the contrary, going up to ncell_max = 4 cells distance to measure hydrodynamical quantities involved in the reconstruction of the Mach number allows a precision of better than 0.1% in this simple 1D shock tube test. This behaviour is the natural outcome of the larger numerical broadening of shock discontinuities for stronger shocks (see Appendix C): strong shocks require more cells to resolve the entire shock layer. We note that increasing the level of refinement does not cure the problem; the shocks are narrower in physical extent, but the number of cells required to describe the shock jump remains the sam.

Fig. 6 Statistics of the numerical Sod shock Mach number relative to its expected value for different shock Mach numbers, 10 (black), 100 (blue), and 1000 (red), using either a maximum of ncellmax = 2 cells (left panel) or ncellmax = 4 cells (right panel) to probe hydrodynamical values in the post-shock and pre-shock regions. These shock tube tests do not model CRs. Pre-shock and post-shock regions need to be probed up to four cells away from the shock cell location for the strongest Mach numbers to be captured accurately.

4.3.2 Cosmic-ray acceleration with constant efficiency

In this test we set up the previous 1D Sod shock tube test with Mach number $\mathcal{M}=10$, and allowed for CR acceleration with a constant efficiency of η = 0.5 (the exact Mach number accounting for CRs added at the shock is $\mathcal{M}=9.56$ for this particular efficiency). We used an adiabatic index for the thermal and CR components of respectively γ = 5∕3 and γ_CR = 4∕3. All the Sod experiments were run without streaming and without radiative thermal or CR losses. The analytical solution with accelerated CRs was provided by Pfrommer et al. 2017 (see their Appendix B).

Figure 7 shows the result of the numerical calculation where the analytical solution is nicely reproduced with the correct Mach number of $\mathcal{M}\simeq 9.56$ positioned at the shock front in one of the cell sampling the numerically broadened discontinuity. Right after the shock discontinuity, in the post-shock region, the thermal pressure shows a few cells that overshoot the expected value. This effect is due to our choice of depositing the accelerated CR energy density a few cells beyond the exact shock location (a strategy we employ to avoid the Pd V compression). Apart from this expected effect, pressures, velocity, density, and the effective adiabatic index of the gas are accurately reproduced.

Fig. 7 Sod shock tube experiment with CR acceleration efficiency of η = 0.5, zero initial CR pressure and γCR = 4∕3 at t = 0.35. Left panels: solution over the full box. Right panels: zoomed-in region over the shock and contact discontinuities for better clarity of the CR shock-accelerated region. From top to bottom: pressures (black: total, blue: thermal, red: CR), the density, the velocity, the Mach number, the effective adiabatic index, and the level of refinement. The symbols stand for the numerical solution, while the solid lines are forthe analytical solution. The exact Sod solution with accelerated CRs is reproduced well by our numerical implementation.

Fig. 8 Sod shock tube experiment with zero initial CR pressure and γCR = 4∕3 with obliquity-dependent CR acceleration efficiency (θB is the so-called obliquity: angle of the pre-shock B field with the normal to the shock) η0 = 0.5ζ(θB) for θB = 30, 45, 60° from topto bottom. The panels show the pressures (black: total, blue: thermal, red: CR) at t = 0.35 over a zoomed-in region over the shock and contact discontinuities for better clarity of the CR shock-accelerated region. The symbols stand for the numerical solution, while the solid lines are for the analytical solution. As expected, the amount of CRs produced at the shock decreases with obliquity, and reproduces well the exact solution.

4.3.3 Cosmic-ray acceleration with magnetic obliquity dependency

In this Sod test, we let the acceleration efficiency η(θ_B) vary with the pre-shock magnetic obliquity angle θ_B and imposedη = 0.5ζ(θ_B) (the previous Sod test wasrun with θ_B = 0°, i.e. the efficiency was η = η₀ = 0.5). We ran three experiments with θ_B = 30, 45, and 60° (i.e. ζ ≃ 0.95, 0.5, and 0.05 respectively), starting with an initial magnetic field with components (B_x, B_y, B_z) = (10⁻¹⁰, 0, 5.77 × 10⁻¹¹), (10⁻¹⁰, 0, 10⁻¹⁰), (5.77 × 10⁻¹¹, 0, 10⁻¹⁰), respectively. Magnetic field magnitudes were chosen to be arbitrarily small so that the magnetic field had no dynamical impact on the gas (i.e. B² ≪ P). The results are shown in Fig. 8, where we see that the expected values of the CR pressure in the shock are reproduced well for any of the adopted magnetic obliquity. We note that the exact location of the shock jump is modified, due to the modified shock velocity, which is governed by the effective adiabatic index in the shock that depends on the amount of accelerated CRs.

4.4 Three-dimensional Sedov explosion

We set up a3D Sedov explosion with the following unitless values: a background at rest with gas density of ρ = 1, P_th = 10⁻⁴, and a point-like explosion of energy E_th = 1 spread over the eight central cells in a box of size unity¹. There are no CRs initially, and only those accelerated into the shock with a constant acceleration efficiency of η = 0.5 will necessarily contribute to the CR distribution. The adiabatic index of the thermal component is γ = 5∕3, and γ_CR = 4∕3 for CRs. In a box of size unity, we start with a base grid of level 6 and allow for 2 extra levels of refinement wherever the cell-to-cell density and pressure variations are larger than 20 and 50%, respectively. The criterion for density is used only where the gas density is higher than that of the background in order to avoid excessive refinement into the hot interior, and instead we focus on the shocked swept-up shell material. For this particular test it is customary to employ a more diffusive solver than HLLD (or Harten–Lax–van Leer–Contact for a pure hydro run) to avoid the formation of the carbuncle phenomenon in shocked cells around the x-, y-, or z-axis of the box, hence, we use, here, the Lax–Friedrich approximate Riemann solver. All Sedov experiments are run without streaming and without radiative thermal or CR losses.

Figure 9 (left panels) shows the density and CR pressure in a thin slice through the centre of the explosion at time t = 0.05. The swept-up material accumulates in a thin shocked layer of gas where CRs are accelerated and they propagate backward through a reverse shock in the bubble interior. We can see finger-like features in the shocked material, which are produced by the discretised nature of the grid; amongst the post-shock cells receiving the accelerated CR energy, some of them can indeed receive energy from several shock cells, while some others receive it only once. We note that Pfrommer et al. (2017) also noticed this effect in their unstructured mesh code, the difference is that their features are randomly located in angle, while here, due to the structured Cartesian nature of our grid, these features follow some π∕2 periodic pattern.

As expected, due to the high adopted value of acceleration efficiency η = 0.5, there is a very significant amount of CRs produced into the dissipation layer of the shock as seen in the spherically averaged radial profiles from Fig. 10. The pressure in the shock layer is a mixture of CRs and thermal particles, while the CR pressure completely dominates the total pressure in the diffuse bubble interior.

It leads to a sharp transition of the effective adiabatic index of the gas from purely thermal outside of the explosion γ_e = γ to purely CR-like in the diffuse bubble γ_e = γ_CR. What matters for the shock dynamics is the effective adiabatic index in the swept-up shock layer that can be inferred from the exact Sedov shock dynamics given a value of γ_e. For analytical guidance, with enthalpy arguments Chevalier (1983) provides the solution for the effective adiabatic index as a function of the fraction of CR pressure w = P_CR∕P_tot in the shocked shell (not to be confused with the acceleration efficiency)

Fig. 9 Sedov explosion with accelerated CRs with η = 0.5 (left panels),and obliquity-dependent acceleration efficiency η = 0.5ζ(θB) with either a uniform magnetic field (middle panels) or a random magnetic field (right panels). Top and bottom panels: respectively slices of density and CR pressure at time t = 0.05, with the solid circle line indicating the position of the Sedov shock front for the exact solution with γe = 7∕5, which are reproduced in all panels to guide the eye throughout (the random magnetic field configuration is better fitted with γe = 1.55), and with magnetic unit vectors overplotted as black segments (the length scale of the random magnetic field corresponds to the size of two large arrows). In the simulation without obliquity dependent acceleration, CR production is close to uniform in the shell except for small numerical grid artefacts. With obliquity dependency, CRs accumulate in polar caps for a uniform magnetic field, and in small patches for the random magnetic field corresponding to the length scale of the field. The position and shape of the shell are also affected by the presence and the configuration of the magnetic field with respectto the obliquity-independent case.

$${\gamma }_{\text{e}}=\frac{5+3w}{3(1+w)}$$(21)

for γ_CR = 4∕3. In agreement with Pfrommer et al. (2017), we find that for the same set-up, an effective adiabatic index in the shock of γ_e = 7∕5 for the exact solution leads to a good recovery of the numerical solution in both total pressure and density, though the maximum values are less pronounced at the shock because of the limited resolution. Increasing the resolution naturally captures theshock profile more faithfully.

We ran two extra simulations with the acceleration efficiency depending on magnetic obliquity η = 0.5ζ(θ_B) and changing from an initial initially uniform magnetic field with (B_x, B_y, B_z) = (10⁻¹⁰, 0, 0) or a randommagnetic field configuration (see Appendix D for details) with a typical coherence length of λ_B = 1∕16 and a similar magnitude of 10⁻¹⁰. For the uniform magnetic field configuration, CRs are accelerated around polar caps along the x-axis of the box with maximum efficiency, and go to zero along the y-axis (or z-axis) as a result of magnetic obliquity (see middle panels of Fig. 9). It results in an ellipsoid shape of the explosion: the position of the shell where CR acceleration is close to zero (y- and z-axes) is further away than where CRs are produced (x-axis) as a result of the higher (resp. lower) effective adiabatic index of the gas mixture in the shell. We note that the exact shape of the ellipsoid is a function of the obliquity-independent part of the acceleration efficiency: the larger ξ is, the more stretched the explosion is (see Pais et al. 2018, for a thorough analysis of this effect). As expected, the density is also higher along the x-direction than along the y-direction (z-direction) as a result of the dependency of the density jump to the adiabatic index of the gas (for strong shocks, ${\mathcal{R}}_{\rho }=4$ for γ_e = 5∕3 and ${\mathcal{R}}_{\rho }=6$ for γ_e = 7∕5).

Finally, the random magnetic field set-up shows a shell mass distribution close to spherical with significant fluctuations with angle (right panels of Fig. 9). It reflects the underlying patchy acceleration and distribution of CR pressure in the swept-up shock layer. On average, the acceleration efficiency is reduced by a factor $\langle \zeta \rangle ={\displaystyle {\int }_0^{\pi /2}\zeta }({\theta }_{\text{B}})\sin {\theta }_{\text{B}}d{\theta }_{\text{B}}\simeq 0.302$ for a purelyrandom upstream magnetic field orientation (see Fig. 11) compared to the simulation without obliquity dependency, and thus to an effective acceleration parameter of η_e ≃ 0.15. Therefore, there is a smaller amount of CRs produced in the shock, and as expected from Chevalier (1983) (see also Castro et al. 2011; Bell 2015), the exact solution is now better reproduced for a lower effective adiabatic index of γ_e = 1.55 (see Fig. 12) and leads to a shock front in advance compared to the obliquity-independent simulation.

Fig. 10 Spherically averaged radial profiles for the 3D Sedov explosion with CR acceleration with constant acceleration efficiency of η = 0.5 of the pressure (blue: thermal pressure, red: CR pressure), density, and effective adiabatic index of the thermal–CR mixture from top to bottom at time t = 0.1. Solid lines stand for the result of the numerical simulation, while the dashed lines of the pressure and density plots are the exact solution of the self-similar profile for an effective adiabatic index of 7∕5 in black (the exact density profile for γ = 5∕3 is also shown as a dashed blue line). The blue and red dashed lines in γe stand for the adiabatic index used for the thermal and CR component, respectively. The thermal–CR mixture produces an explosion similar to a Sedov solution with effective adiabatic index of γe = 7∕5, which delays the position of the shock due to the lower pressure work exerted by the shocked shell.

5 Turbulent box of the ISM

We ran turbulent ISM boxes in the same spirit of Commerçon et al. (2019) except that here we started with negligible CR pressure (10⁻¹⁰ that of the thermal pressure) and let it build through the turbulence-generated shocks. The simulations have a uniform 128³ Cartesian resolution in a box of 50 pc, leading to a spatial resolution of 0.4 pc. The initial gas density and temperature are 2 cm⁻³ and 4460 K, respectively, with a mean molecular weight of μ = 1.4 assumed throughout. We started with an initial thermal pressure of P_th,0 = 1.2 × 10⁻¹² erg cm⁻³. The initial magnetic field was uniform and was set up in the x-direction of the box with a magnitude of 0.1 μG, leading to a plasma beta parameter of β = P_th,0∕P_mag,0 ≃ 3 × 10³. We did not allow for self-gravity of the gas or for any refinement. Cooling proceeded on the thermal component following Audit & Hennebelle (2005), while we neglected the role of Coulomb and hadronic losses of CR protons (Enßlin et al. 2007; Guo & Oh 2008).

The turbulence is forced at all times with an injection scale of k_turb = 2 (i.e. corresponding to half the size of the box) and with a parabolic shape in the Fourier space $\tilde{f}(k)\propto 1-{(k-k_{\text{turb}})}^2$ with k sampled in the range k = [1, 3]. The turbulence is applied intermittently with an auto-correlation time of 0.5 Myr and with a compression-to-solenoidal ratio of 1 (see Commerçon et al. 2019, for more details).

Fig. 11 Stacked PDF of the magnetic obliquity in the Sedov experiment between t = 0.05 − 0.1 for the random magnetic field configuration (solid histogram), compared to the random distribution (black dashed line). The distribution of magnetic obliquity is compatible with a purely random field as expected, thus leading to a reduced efficiency of ⟨ζ⟩ = 0.302.

Fig. 12 Similar to Fig. 10, but for the random magnetic field configuration and with obliquity dependency of the CR acceleration efficiency η = 0.5ζ(θB). Here the Sedov profile is better fitted with an effective adiabatic index of γe = 1.55.

5.1 $\mathcal{M}$- and X_CR-independent acceleration efficiency

We start with a batch of simulations where the acceleration efficiency does not depend on $\mathcal{M}$ and X_CR (i.e. ξ = 1). We set up three different simulations: (i) without CR acceleration (i.e. η₀ = 0, NoShock); (ii) with CR acceleration and η = η₀ = 0.1 (i.e. where CR acceleration does not depend on magnetic, NoThetaB); (iii) with CR acceleration and η₀ = 0.1 (i.e. where CR acceleration depends on magnetic obliquity, ThetaB); and with η₀ = 0.1 and CR streaming (Streaming). We note that we use rather large values of CR acceleration efficiencies given the moderate Mach numbers of only 2–4 (e.g. Kang & Jones 2005; Kang & Ryu 2013) obtained in that experiment. This somewhat reflects the more typicalSN-generated CR acceleration efficiencies corresponding to much larger values of the shock Mach number than we can capture here with this simplified set-up. For the sake of a testable set-up for our new implemented algorithm, these values allow us to reach an appreciable amount of CR energy density in the simulated volume over a few turbulent crossing times t_cross = 6.7 Myr, where it is the box length divided by the rms velocity u_rms = 7.3 km s⁻¹ (here measured at t = 20 Myr for the Streaming run).

Shocks are driven in sheets with moderate Mach numbers of $\mathcal{M}\simeq 3$–4, as can be seen in Fig. 13 for the Streaming run (other simulations show similar features) at time t = 10 Myr, which dissipates the energy of shocks with a typical range of flux values of e_dissu₂ ≃ 10⁴⁴–10⁴⁵ erg Myr⁻¹ pc⁻². Figure 14 shows maps of the CR pressure at two different times t = 10 and 20 Myr for the simulation NoThetaB, ThetaB, andStreaming. At t = 10 Myr the CR pressure has already built up to appreciable levels thanks to turbulence-generated shocks in the box, with clustered regions of pressure at levels similar to or above the initial thermal pressure (P_th,0 ≃ 10⁻¹² erg cm⁻³). The NoThetaB simulation has, as expected, the highest values of CR pressure since CR acceleration efficiency is always equal to η = 0.1, while in the two other runs it can only reach this value for a perfectly aligned pre-shock magnetic field with the normal to the shock. At this early stage of the simulation, the effect of streaming is still very moderate on the CR pressure distribution. It reduces the range of the lowest and highest values of pressure mimicking the effect of a diffusion process; nonetheless, the geometrical features are easily recognizable between the ThetaB and Streaming runs (and NoThetaB as well).

Figure 15 (top panel) shows the thermal and CR energies in the simulated volumes as a function of time. The total thermal energy in the box is quickly reduced in 3 Myr by nearly a factor of 3 with very negligible differences by the end of the simulation between the four simulations. The total CR energy builds up almost linearly with time as a result of nearly constant dissipated energy and acceleration efficiency over time, once passed the first 5 Myr. This CR pressure provides a support to the total pressure close to the thermal pressure, if not above (NoThetaB case at t = 20 Myr). The magnetic energy quickly increases early on and saturates at a plasma beta β ≃10 similar for the four different simulations. We note that this level of magnetic field is crucial for the CR streaming to have an appreciable effect on the CR pressure distribution as the streaming velocity scales with the Alfvén velocity.

As we discussed in Sect. 4.4, the average obliquity-dependent part of the CR acceleration efficiency must be ⟨ζ⟩ ≃ 0.302 for a purely random field, which seems supported by the apparent randomness of magnetic vectors (white arrows in Fig. 14), but we show that this is not the case. Figure 15 (bottom panel) shows the dissipated energy per unit time in the form of thermal or CR energy. Dissipated thermal energies are very similar for the three simulations, although there is a slight deviation at late times for the Streaming run. However, the dissipated CR energy shows a larger than a factor 3 difference between the non-θ_B and the θ_B dependencies, closer to a factor 6–8 difference between the NoThetaB and ThetaB runs. This is indirect evidence that pre-shock magnetic fields are not randomly oriented, but show preferentially within-shock-plane orientations. To clarify further, we measure the probability density function (PDF) of the obliquity for the ThetaB and Streaming runs at time t = 20 Myr in Fig. 16, which shows that the PDF is skewed towards larger angles: upstream magnetic fields are more likely to be perpendicular to the normal of shocks than for a random field, in agreement with the estimated reduced efficiency of CR acceleration.

We also note that at time t = 10 Myr, the CR energy density is a factor 2 lower with streaming, while the CR dissipated energy before t ≤ 10 Myr is similar to that of the simulation without streaming. Therefore, this difference in CR energy density is directly due to streaming (as opposed to streaming reducing shock strengths) putting CRs away from compressed regions (shocks or not) where the adiabatic compression can further enhance the overall CR pressure.

At time t = 20 Myr, the distributions of CR pressure (Fig. 14) in the three simulations differ very significantly. While the NoThetaB and ThetaB runs look like a renormalised versions of one another, albeit with different specific locations of voids and plume-like features, the Streaming run has lost most of its CR structure with a closer to uniform distribution of CR pressure in the box.

These distinct CR pressure evolutions and distributions lead to very important differences in the way the matter is compressed into overdense regions of the flow. Figure 17 shows the time evolution of the mass fraction of dense gas, which is arbitrarily chosen at five times the initial gas density (i.e. for n > 10 cm⁻³, but the results are qualitatively independent of this choice). Since only the thermal pressure is affected by radiative losses, which are larger at high gas densities, it is the CR pressure that accumulates in regions of high gas densities that can provide the support against compression. Therefore, it shows that the simulations with the largest total CR energy are the simulations with the lowest amount of dense gas. However, the streaming introduces a subtle but significant difference to this overall picture. Since streaming smooths theCR pressure in the ISM, the high gas density is much less clustered for a given total energy in the box. At t = 20 Myr in the Streaming run, the total CR energy is indeed equal to that at t = 18 Myr in the ThetaB run; nonetheless, the mass fraction of dense gas is respectively 40% higher in the Streaming run. Recast into an “effective” diffusion framework, we can deduce that streaming behaves like anisotropic diffusion with an effective diffusion coefficient to be determined through comparison with the corresponding simulations, which we defer to a future work.

Fig. 13 Projection of the Mach number $\mathcal{M}$ (top) and dissipated energy flux edissu2 (bottom) for the Streaming turbulent box, with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$, at time t = 10 Myr over a box thickness of half the size of the box centred on the middle of the box. Shocks are driven in sheets with a bulk of the Mach number of moderate values $\mathcal{M}\simeq 3$–4.

Fig. 14 Cosmic-ray pressure maps of the turbulent box simulation, with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$, in a thin plane within the x-plane of the middle of the box at time t = 10 Myr (top panels) and t = 20 Myr (bottom panels) for the simulation without CR streaming and without (left panels) or with (right panels) obliquity dependency for CR acceleration, and with obliquity and CR streaming (right panels). The black segments depict the orientation of the unitary magnetic vectors. The simulation without obliquity builds the CR pressure faster. The presence of the streaming instability allows for a more uniform distribution of CRs in the simulated volume.

5.2 $\mathcal{M}$- and X_CR-dependent acceleration efficiency

We show here the results of the turbulent box experiments, where this time the efficiency dependency $\xi (\mathcal{M},X_{\text{CR}})$ is not assumed to be equal to 1, but varies according to the scaled values of Kang & Ryu (2013). We ran two numerical experiments, free of CR streaming, with and without the magnetic obliquity dependency ζ(θ_B), called ThetaB_KR13 and NoThetaB_KR13, respectively. We recall that we started with an initial CR pressure of almost zero so that X_CR = 10⁻¹⁰ everywhere in the box at time t = 0, and that a normalisation (maximum) acceleration efficiency η₀ = 0.1 was used throughout.

Figure 18 shows the evolution of the CR flux-weighted mean value of $\mathcal{M}$ (top panel) and X_CR (middle panel), and the evolution of the energy flux-weighted mean acceleration efficiencies (bottom panel) $\eta ={\eta }_0\xi (\mathcal{M},X_{\text{CR}})\zeta ({\theta }_{\text{B}})$, $\xi (\mathcal{M},X_{\text{CR}})$, and ζ(θ_B) as a function of time. The bulk of the CR energy is produced in shocks of $\mathcal{M}\simeq 3-4$ with a slight decrease over time. As CRs are produced, the upstream CR-to-thermal pressure ratio rises to values close to X_CR ≃ 0.1−0.2 at time t = 20 Myr. The corresponding CR acceleration efficiencies also evolve with time since ξ varies significantly for this range of moderate Mach number as a function of X_CR reaching ξ ≃ 0.03 and 0.1 at t = 20 Myr for the ThetaB_KR13 and NoThetaB_KR13 runs respectively. In particular, there is an increase between 10 and 20 Myr of the acceleration efficiency by one order of magnitude in both simulations. The difference between the two simulations is that the obliquity dependent run has a lower overall acceleration efficiency η since nearlyrandom magnetic fields (see Fig. 16) reduce the ζ component to ≃ 0.2. We note that the choice of starting with X_CR = 0 for educative purposes makes these simulations extremely unrepresentative of the ISM of normal galaxies (though it might apply for proto-galaxies), and delay the build-up of the CR pressure. Nonetheless, we show that our implementation of the $\mathcal{M}$, X_CR (and θ_B) dependency of η leads to interesting results in the build-up of the CR pressure through shocks, and might be useful for a broad range of applications.

Fig. 15 Top panel: time evolution of total thermal (solid lines), CR energy (dashed lines), and magnetic energy (dot-dashed lines) in the simulated turbulent ISM boxes for the simulations without shock-acceleration (black), without CR streaming, and without (red) or with (green) obliquity dependency for CR acceleration, and with obliquity and CR streaming (blue) with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$. Bottom panel: evolution of the dissipated thermal (solid) and CR (dashed) energy rates at shocks for the same simulations.

Fig. 16 Probability density function of the magnetic obliquity in the ISM boxes, and with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$, at time t = 20 Myr with CR streaming (blue) or without (green), compared to the random distribution in black dashed. Those simulations are more likely to have magnetic field perpendicular to the normal of shocks than for a random distribution, therefore, lowering the CR acceleration efficiency compared to the averaged random distribution, i.e. ⟨ζ⟩ = 0.302.

Fig. 17 Time evolution of the mass fraction of dense gas in the simulated turbulent ISM boxes for the simulations without shock-acceleration (black), without CR streaming, and without (red) or with (green) obliquity dependency for CR acceleration, and with obliquity and CR streaming (blue), and with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$.

6 Conclusion

We have introduced a new modelling of anisotropic CR streaming and dynamical CR shock-acceleration for the AMR code RAMSES (Teyssier 2002). Streaming is solved with a diffusion approach where the diffusion step is performed with a time implicit scheme (Dubois & Commerçon 2016), and can handle complex multi-dimensional problems with non-trivial magnetic field geometries. CR acceleration at shocks through the DSA mechanism is obtained by accurately detecting shocks, and measuring their Mach number and magnetic obliquity. We have shown that our numerically CR accelerated solutions faithfully reproduces exact 1D Sod shock tube solutions. CR-modified 3D Sedov solutions with accelerated CRs have been tested with various background magnetic field configurations (hence, obliquities). They show very good agreement with previous numerical experiments (Pfrommer et al. 2017) with CRs reducing the effective adiabatic index and slowing down the motion of the shell. Obliquity dependency of the acceleration leads to a significant modification of the CR distribution in the shell of the Sedov explosion with either a polar or patchy distribution when the coherence length of the background magnetic field is respectively larger or smaller than the bubble size. This also has consequences on the final shape of the bubble, with a significant elongation of the bubble when the magnetic field has a large field coherence with respect to the bubble size (Pais et al. 2018).

Finally, the effect of CR streaming and CR acceleration has been tested in a turbulent box mimicking the motions within the ISM on scales of tens of pc (Commerçon et al. 2019). CRs are produced at shock surfaces and are spread throughout the entire volume by convection and streaming. CRs have important consequences on the reservoir of cold gas available as they provide a long-term pressure support against compressed material, and streaming substantially modifies the small-scale distribution of CRs, and in turn the clustering of gas. The obliquity of the field produces a strong suppression of the effective acceleration efficiency, a factor of ~ 2 beyond the pure random case as a result of the preferential alignment of magnetic fields with shock surfaces.

These new CR physics modules embedded in the RAMSES code make it useful for the study of the impact of CRs in a wide variety of situations, such as the acceleration of CRs by cosmic shocks, galactic-wide outflows driven by CRs (Dashyan & Dubois 2019), the release of CRs in galaxy clusters by active galactic nuclei, studies of supernova remnants, and the release of CRs in the supernova-driven turbulence of the ISM, which we defer to future work.

Fig. 18 Top to bottom: cosmic-ray flux-weighted mean Mach number $\mathcal{M}$, CR flux-weighted mean CR-to-thermal pressure ratio XCR, and dissipated energy flux-weighted efficiencies $\eta ={\eta }_0\xi (\mathcal{M},X_{\text{CR}})\zeta ({\theta }_{\text{B}})$, $\xi (\mathcal{M},X_{\text{CR}}),$ and ζ(θB) as a function of time for the simulation with (red) and without (green) obliquity dependency. We recall that η0 = 0.1 is used in those simulations and that ξ is a function of $\mathcal{M}$ and XCR as extrapolated from the values of Kang & Ryu (2013).

Acknowledgements

We thank G. Dashyan, M. Lemoine, C. Pfrommer, R. Teyssier, and A. Wagner for enlightening discussions. We warmly thank S. Rouberol for smoothly running the Horizon cluster on which several simulations were run. This work was supported by the ANR grant LYRICS (ANR-16-CE31-0011) and the CNRS programs “Programme National de Cosmologie et Galaxies” (PNCG) and “Physique et Chimie du Milieu Interstellaire” (PCMI).

Appendix A Effect of perpendicular diffusion on streaming

Fig. A.1 Energy density maps at t = 0.02 for the same set-up as for the 2D sinusoid loop described in Sect. 3.2.2 with fiso = 10−1 (left) or fiso = 10−3 (right).

Appendix B Tabulated values of ξ($\mathcal{M}$, XCR)

Table B.1 shows the tabulated values of Kang & Ryu (2013) renormalised to 1 (see Sect. 4.2 for details).

Here we vary the value of the isotropic component of the streaming diffusion term from f_iso = 10⁻³ to f_iso = 10⁻¹ (to be compared with the value of f_iso = 10⁻² used by default in Sect. 3.2.2) with respect to pure anisotropy. Figure A.1 shows that increasing the value of f_iso to 10⁻¹ leads to more diffusion outside of the loop, which decreases the values of the maximum, while f_iso = 10⁻³ produces numerically driven finger-like features but allows a more contained CR distribution in the loop.

Appendix C Shock numerical broadening

We show in Fig. C.1 a zoomed-in view of the shock discontinuity for the Sod shock tube experimentsdescribed in Sect. 4.3.1 (i.e. without CRs) and for the three different Mach numbers $\mathcal{M}=10$, 100, and 1000. Instead of a pure discontinuity (the exact solution is shown as a solid line) the numerical shock is broadened by numerical diffusion with typically 4–5 cells; the number of cells in the discontinuity to match the exact pre- and post-shock pressures increases with the value of the Mach number, and given the quadratic increase in pressure jump with Mach number, any error is strongly amplified. In the strongest shock example shown in the bottom panel, using only two cells away from the shock would lead to underestimating the Mach number by a factor of 10 (Mach number scales with $R_{\text{P}}^{1/2}$ and the upstream value two cells away from the shock is ≃100 times that of the true value).

Fig. C.1 Pressure profiles at time $t\simeq 3.5/\mathcal{M}$ around the shock discontinuity for the Mach $\mathcal{M}=$ 10 (top), 100 (middle), and 1000 (bottom) experiments. The result of the numerical solution is shown as diamonds; the red symbol highlights the position of the shock cell given by a shock finder algorithm. The solid line is the exact numerical solution. We see that the numerical shock tends to broaden with increasing Mach number, and given the largest error made on the post- and pre-shock regions, the error on the evaluated Mach number becomes larger for a small kernel (ncell,max = 2).

Appendix D Random magnetic fields

In order to set a random magnetic field fulfilling the ∇.B = 0 constraint, we first set up a random potential vector on the nodes of a Cartesian grid of arbitrary resolution $n_{\text{pot}}^3$ cells (there are actually ${(n_{\text{pot}}+1)}^3$ values of potential vectors drawn at nodes of the $n_{\text{pot}}^3$ sampling cells), with the right-, top-, and back-most boundaries being replicates of the left-, bottom-, and front-most boundaries to ensure the correct periodicity of the (staggered) magnetic field. In the cases simulated in this paper the AMR cell size is smaller than or equal to 1∕n_pot, which means that the vector potential is the trilinear interpolation of the surrounding node vector potentials projected along the AMR cell edge. Once these reconstructed vector potentials are obtained along AMR cell edges, the staggered magnetic field (one B-field perpendicular to each face of AMR cells) is obtained by taking the rotational of the potential vector of the face-surrounding edges. This procedure guarantees that the magnetic field is random, ∇.B = 0, and the consistency of the coarse-to-fine values of the B-field. We note that we took the initial random potential vector as a white noise vector, but this can be modified to account for any given spectrum of the vector potential (or magnetic field), and to obtain any desired shape of the magnetic power spectrum, as the power spectrum of B scales as k (i.e. the wave number) times the power spectrum of A.

References

All Tables

All Figures

	Fig. 1 Evolution of a 1D sinusoid of CR energy density with streaming advection only as a function of position with 512 cells and imposing a constant Alfvén velocity of 1. The solution is made of two plateaus as the maxima are capped over time, due to the infinite streaming diffusion coefficient, while the two regions between the two plateaus move at a velocity of ± γCR = ±1.4.
In the text

	Fig. 2 Solution at t = 0.02 of the sinusoid experiment with different resolution from uniform level 4 to 8 (from light red to dark red) and level 10 (in black). The relative errors are compared to the reference numerical solution of level 10. Even for very low resolution the relative error is never larger than a few percentage points.
In the text

	Fig. 3 Convergence of the L2 norm for the sinusoid experiment using the solution at t = 0.02. The norm is compared to the reference numerical solution of level 10. The L2 norm scales with Δx1.87±0.08 as indicated by the dashed line.
In the text

Fig. 4 Cosmic-ray energy density maps at t = 0 (top left) and t = 0.02 (top right) for an initial angle-dependent sinusoid within a purely circular magnetic field with an Alfvén velocity of 1 and a resolution of 1282 cells. The energy is evolved with the streaming advection-diffusion term only. The bottom panel shows the radially averaged energy in theradius interval r = [0.15, 0.35] as a functionof the polar angle θ.

In the text

	Fig. 5 Acceleration efficiency $\xi (\mathcal{M},X_{\text{CR}})$ as a functionof the Mach number $\mathcal{M}$ for different values of the upstream CR-to-thermal pressure ratio XCR. The values are obtained from the XCR = 0 and 0.025 values of Kang & Ryu (2013) and renormalised to a maximum value of 1.
In the text

Fig. 6 Statistics of the numerical Sod shock Mach number relative to its expected value for different shock Mach numbers, 10 (black), 100 (blue), and 1000 (red), using either a maximum of ncellmax = 2 cells (left panel) or ncellmax = 4 cells (right panel) to probe hydrodynamical values in the post-shock and pre-shock regions. These shock tube tests do not model CRs. Pre-shock and post-shock regions need to be probed up to four cells away from the shock cell location for the strongest Mach numbers to be captured accurately.

In the text

Fig. 7 Sod shock tube experiment with CR acceleration efficiency of η = 0.5, zero initial CR pressure and γCR = 4∕3 at t = 0.35. Left panels: solution over the full box. Right panels: zoomed-in region over the shock and contact discontinuities for better clarity of the CR shock-accelerated region. From top to bottom: pressures (black: total, blue: thermal, red: CR), the density, the velocity, the Mach number, the effective adiabatic index, and the level of refinement. The symbols stand for the numerical solution, while the solid lines are forthe analytical solution. The exact Sod solution with accelerated CRs is reproduced well by our numerical implementation.

In the text

Fig. 8 Sod shock tube experiment with zero initial CR pressure and γCR = 4∕3 with obliquity-dependent CR acceleration efficiency (θB is the so-called obliquity: angle of the pre-shock B field with the normal to the shock) η0 = 0.5ζ(θB) for θB = 30, 45, 60° from topto bottom. The panels show the pressures (black: total, blue: thermal, red: CR) at t = 0.35 over a zoomed-in region over the shock and contact discontinuities for better clarity of the CR shock-accelerated region. The symbols stand for the numerical solution, while the solid lines are for the analytical solution. As expected, the amount of CRs produced at the shock decreases with obliquity, and reproduces well the exact solution.

In the text

Fig. 9 Sedov explosion with accelerated CRs with η = 0.5 (left panels),and obliquity-dependent acceleration efficiency η = 0.5ζ(θB) with either a uniform magnetic field (middle panels) or a random magnetic field (right panels). Top and bottom panels: respectively slices of density and CR pressure at time t = 0.05, with the solid circle line indicating the position of the Sedov shock front for the exact solution with γe = 7∕5, which are reproduced in all panels to guide the eye throughout (the random magnetic field configuration is better fitted with γe = 1.55), and with magnetic unit vectors overplotted as black segments (the length scale of the random magnetic field corresponds to the size of two large arrows). In the simulation without obliquity dependent acceleration, CR production is close to uniform in the shell except for small numerical grid artefacts. With obliquity dependency, CRs accumulate in polar caps for a uniform magnetic field, and in small patches for the random magnetic field corresponding to the length scale of the field. The position and shape of the shell are also affected by the presence and the configuration of the magnetic field with respectto the obliquity-independent case.

In the text

Fig. 10 Spherically averaged radial profiles for the 3D Sedov explosion with CR acceleration with constant acceleration efficiency of η = 0.5 of the pressure (blue: thermal pressure, red: CR pressure), density, and effective adiabatic index of the thermal–CR mixture from top to bottom at time t = 0.1. Solid lines stand for the result of the numerical simulation, while the dashed lines of the pressure and density plots are the exact solution of the self-similar profile for an effective adiabatic index of 7∕5 in black (the exact density profile for γ = 5∕3 is also shown as a dashed blue line). The blue and red dashed lines in γe stand for the adiabatic index used for the thermal and CR component, respectively. The thermal–CR mixture produces an explosion similar to a Sedov solution with effective adiabatic index of γe = 7∕5, which delays the position of the shock due to the lower pressure work exerted by the shocked shell.

In the text

	Fig. 11 Stacked PDF of the magnetic obliquity in the Sedov experiment between t = 0.05 − 0.1 for the random magnetic field configuration (solid histogram), compared to the random distribution (black dashed line). The distribution of magnetic obliquity is compatible with a purely random field as expected, thus leading to a reduced efficiency of ⟨ζ⟩ = 0.302.
In the text

	Fig. 12 Similar to Fig. 10, but for the random magnetic field configuration and with obliquity dependency of the CR acceleration efficiency η = 0.5ζ(θB). Here the Sedov profile is better fitted with an effective adiabatic index of γe = 1.55.
In the text

	Fig. 13 Projection of the Mach number $\mathcal{M}$ (top) and dissipated energy flux edissu2 (bottom) for the Streaming turbulent box, with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$, at time t = 10 Myr over a box thickness of half the size of the box centred on the middle of the box. Shocks are driven in sheets with a bulk of the Mach number of moderate values $\mathcal{M}\simeq 3$–4.
In the text

Fig. 14 Cosmic-ray pressure maps of the turbulent box simulation, with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$, in a thin plane within the x-plane of the middle of the box at time t = 10 Myr (top panels) and t = 20 Myr (bottom panels) for the simulation without CR streaming and without (left panels) or with (right panels) obliquity dependency for CR acceleration, and with obliquity and CR streaming (right panels). The black segments depict the orientation of the unitary magnetic vectors. The simulation without obliquity builds the CR pressure faster. The presence of the streaming instability allows for a more uniform distribution of CRs in the simulated volume.

In the text

Fig. 15 Top panel: time evolution of total thermal (solid lines), CR energy (dashed lines), and magnetic energy (dot-dashed lines) in the simulated turbulent ISM boxes for the simulations without shock-acceleration (black), without CR streaming, and without (red) or with (green) obliquity dependency for CR acceleration, and with obliquity and CR streaming (blue) with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$. Bottom panel: evolution of the dissipated thermal (solid) and CR (dashed) energy rates at shocks for the same simulations.

In the text

Fig. 16 Probability density function of the magnetic obliquity in the ISM boxes, and with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$, at time t = 20 Myr with CR streaming (blue) or without (green), compared to the random distribution in black dashed. Those simulations are more likely to have magnetic field perpendicular to the normal of shocks than for a random distribution, therefore, lowering the CR acceleration efficiency compared to the averaged random distribution, i.e. ⟨ζ⟩ = 0.302.

In the text

	Fig. 17 Time evolution of the mass fraction of dense gas in the simulated turbulent ISM boxes for the simulations without shock-acceleration (black), without CR streaming, and without (red) or with (green) obliquity dependency for CR acceleration, and with obliquity and CR streaming (blue), and with η0 = 0.1 and $\xi (\mathcal{M},X_{\text{CR}})=1$.
In the text

Fig. 18 Top to bottom: cosmic-ray flux-weighted mean Mach number $\mathcal{M}$, CR flux-weighted mean CR-to-thermal pressure ratio XCR, and dissipated energy flux-weighted efficiencies $\eta ={\eta }_0\xi (\mathcal{M},X_{\text{CR}})\zeta ({\theta }_{\text{B}})$, $\xi (\mathcal{M},X_{\text{CR}}),$ and ζ(θB) as a function of time for the simulation with (red) and without (green) obliquity dependency. We recall that η0 = 0.1 is used in those simulations and that ξ is a function of $\mathcal{M}$ and XCR as extrapolated from the values of Kang & Ryu (2013).

In the text

	Fig. A.1 Energy density maps at t = 0.02 for the same set-up as for the 2D sinusoid loop described in Sect. 3.2.2 with fiso = 10−1 (left) or fiso = 10−3 (right).
In the text

Fig. C.1 Pressure profiles at time $t\simeq 3.5/\mathcal{M}$ around the shock discontinuity for the Mach $\mathcal{M}=$ 10 (top), 100 (middle), and 1000 (bottom) experiments. The result of the numerical solution is shown as diamonds; the red symbol highlights the position of the shock cell given by a shock finder algorithm. The solid line is the exact numerical solution. We see that the numerical shock tends to broaden with increasing Mach number, and given the largest error made on the post- and pre-shock regions, the error on the evaluated Mach number becomes larger for a small kernel (ncell,max = 2).

In the text