Issue 
A&A
Volume 631, November 2019



Article Number  A121  
Number of page(s)  17  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201936275  
Published online  05 November 2019 
Shockaccelerated cosmic rays and streaming instability in the adaptive mesh refinement code Ramses
^{1}
Institut d’Astrophysique de Paris,
UMR 7095,
CNRS, UPMC University of Paris VI,
98 bis boulevard Arago,
75014 Paris, France
email: dubois@iap.fr
^{2}
Centre de Recherche Astrophysique de Lyon UMR5574, ENS de Lyon, Universite Lyon 1, CNRS, Université de Lyon,
69007 Lyon, France
^{3}
Laboratoire Univers et Particules de Montpellier (LUPM), Université Montpellier, CNRS/IN2P3, CC72, place Eugène Bataillon,
34095 Montpellier Cedex 5, France
Received:
9
July
2019
Accepted:
19
September
2019
Cosmic rays (CRs) are thought to play a dynamically important role in several key aspects of galaxy evolution, including the structure of the interstellar medium, the formation of galactic winds, and the nonthermal pressure support of halos. We introduce a numerical model solving for the CR streaming instability and acceleration of CRs at shocks with a fluid approach in the adaptive mesh refinement code RAMSES. CR streaming is solved with a diffusion approach and its anisotropic nature is naturally captured. We introduce a shock finder for the RAMSES code that automatically detects shock discontinuities in the flow. Shocks are the loci for CR injection, and their efficiency of CR acceleration is made dependent on the upstream magnetic obliquity according to the diffuse shock acceleration mechanism. We show that the shock finder accurately captures shock locations and estimates the shock Mach number for several problems. The obliquitydependent injection of CRs in the Sedov solution leads to situations where the supernova bubble exhibits large polar caps (homogeneous background magnetic field), or a patchy structure of the CR distribution (inhomogeneous background magnetic field). Finally, we combine both accelerated CRs with streaming in a simple turbulent interstellar medium box, and show that the presence of CRs significantly modifies the structure of the gas.
Key words: magnetohydrodynamics / methods: numerical / cosmic rays / shock waves / ISM: supernova remnants / ISM: structure
© Y. Dubois et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Cosmic rays (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 nonnegligible 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 particleincell 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 preexisting amount of CRs (Caprioli et al. 2018). There is a large body of evidence for CRs accelerated in the shockedshell 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 selfgenerated 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 galacticwide 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 lowmass 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 multiphase 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 preprocessed 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 magnetohydrodynamics 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 Magnetohydrodynamics 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 cosmicray magnetohydrodynamics (CRMHD) of a fluid mixture made of thermal particles and CRs can be obtained:
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^{2} + e_{th} + e_{CR} + B^{2}∕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^{2}∕(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^{3}, where Δx is the cell size. The terms on the righthand side of the equations are treated as source terms with P_{CR} ∇.u the CR pressure work term, F_{CR,d} = −D_{0}b(b.∇e_{CR}) the anisotropic diffusion flux term, D_{0} 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 composed of the thermal and CR radiative loss terms, where the CR loss term () is the nonconserving sum of radiative losses from cosmic rays turning as a heating rate for the thermal component. Finally, and this is the core of this paper, we detail how the streaming instability terms (advectiondiffusion term) and (heating term), and the CR acceleration at shocks 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 righthand 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 timestep 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 energymomentum 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 energydependent 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} = 0 and .
3 Cosmicray streaming
3.1 Numerical implementation
Cosmic rays propagating faster than the Alfvén velocity 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 ionneutral damping, nonlinear 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 selfexcited Aflvén waves and can increase the effective value at which CRs are allowed to stream down their own gradient at superAlfvénic velocities u_{st} = −f_{SA}u_{A}sign(b.∇e_{CR}), where f_{SA} ≥ 1 is the superAflvé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: (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^{3} (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 timestep condition of Δt = hΔx^{2}∕(2e_{CR}u_{A}), and where h should be the size of a few cells. Nonetheless, this timestep 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 (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): (8)
Therefore, this advectiondiffusion 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: (9)
where D = D_{0} + D_{st}. The F_{CR,ds} diffusion flux can be arbitrarily decomposed into an anisotropic and isotropic part (10)
where D_{∥} = (1 − f_{iso})D, D_{iso} = f_{iso}D, 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 (11)
where the cellcentred fluxes are computed with cellcornered values using the symmetric scheme from Günter et al. (2005):
The anisotropic cell corner flux is (12)
where barred quantities are arithmetic averages over the cells connected to the corner, i.e.
We notethat all hydrodynamical variables in RAMSES are cellcentred except for the magnetic field which is facecentred. The streaming diffusion coefficient is computed as (13)
where upper tilde quantities stand for cellcentred quantities reconstructed from a combination of cellcentred and facecentred quantities: (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^{28} cm^{2} s^{−1} 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 Onedimensional sinusoid
In order to test the implementation of the CR advectiondiffusion 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 selfconsistency. 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 xaxis. 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 highresolution 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 1024cell 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 Twodimensional sinusoid in a looped magnetic field
In this test case we try to mimic the 1D sinusoid problem embedded in a nonuniform 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^{−5} 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^{−2}; 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 angledependent 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 Δx^{1.87±0.08} as indicated by the dashed line. 
Fig. 4 Cosmicray energy density maps at t = 0 (top left) and t = 0.02 (top right) for an initial angledependent sinusoid within a purely circular magnetic field with an Alfvén velocity of 1 and a resolution of 128^{2} cells. The energy is evolved with the streaming advectiondiffusion 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 Shockaccelerated 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 pseudoentropy) 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 cloudincell 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 , with . Keeping in mind these conditions, the Mach number of eligible cells is computed according to the criteria using upstream (preshock) and downstream (postshock) fluid variables. Using the RankineHugoniot 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 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 (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: (16)
Here , γ_{i} = P_{i}∕ϵ_{i} + 1 for i = {1, 2} (respectively upstream and downstream) and γ_{e} = (γP_{th,2} + γ_{CR}P_{CR,2})∕P_{2} for the downstream region. In the limit where the weighted adiabatic indexes are equal γ_{e} = γ_{1} = γ_{2} 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 cloudincell 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 , 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 Cosmicray acceleration at shocks
At shocks the kinetic energy flux of the upstream flow (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,diss}u_{2}, CR energy ϕ_{CR,2} = e_{CR,diss}u_{2} and the remaining into kinetic and magnetic energy. For classical strong shocks without CR acceleration, the ratio of postshock thermal (dissipated) energy to the preshock kinetic energy 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 (17)
where e_{diss} = e_{th,diss} + e_{CR,diss} is the dissipated internal energy of the gas, u_{2} is the downstream velocity in the frame of the moving shock, and is the acceleration efficiency of CRs at shocks, which is a function of the Mach number, the upstream CRtothermal 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 ), we replace u_{2} by , 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 (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 the jump density ratio. The jump density ratio is obtained from the direct evaluation of the upstream and downstream densities (19)
The and 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_{1}.n_{s}). The dependency of the efficiency of CR acceleration with this socalled “magnetic obliquity” can be factorised out, and approximated by the following functional form (Pais et al. 2018): (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 is obtained from the results of Kang & Ryu (2013), and is an increasing function of both 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 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 and X_{CR}. We note that their obtained values of the acceleration efficiency saturates at η_{0} = 0.225, a factor of ~ 2 larger than the maximum values obtained by Caprioli & Spitkovsky (2014) for parallel shocks (θ_{B} = 0). We thus renormalise by 0.225 so that the maximum allowed efficiency is explicitly controlled by η_{0}. The values of ξ are shown in Fig. 5 and are available as tabulated values in Appendix B. We note that obliquitydependent CR acceleration simulations conducted by Caprioli et al. (2018) with a preexisting population of CRs in the upstream region suggest that the transition of the obliquitydependent 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 postshock 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 postshock region, and our experiments have taught us that the direct injection in the shock systematically overestimates the resulting CR energy density in the postshock 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 as a functionof the Mach number for different values of the upstream CRtothermal pressure ratio X_{CR}. The values are obtained from the X_{CR} = 0 and 0.025 values of Kang & Ryu (2013) and renormalised to a maximum value of 1. 
4.3 Onedimensional 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 celltocell 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 preshock and postshock 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 ncell_{max} = 2 cells (left panel) or ncell_{max} = 4 cells (right panel) to probe hydrodynamical values in the postshock and preshock regions. These shock tube tests do not model CRs. Preshock and postshock 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 Cosmicray acceleration with constant efficiency
In this test we set up the previous 1D Sod shock tube test with Mach number , and allowed for CR acceleration with a constant efficiency of η = 0.5 (the exact Mach number accounting for CRs added at the shock is 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 positioned at the shock front in one of the cell sampling the numerically broadened discontinuity. Right after the shock discontinuity, in the postshock 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: zoomedin region over the shock and contact discontinuities for better clarity of the CR shockaccelerated 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 obliquitydependent CR acceleration efficiency (θ_{B} is the socalled obliquity: angle of the preshock 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 zoomedin region over the shock and contact discontinuities for better clarity of the CR shockaccelerated 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 Cosmicray acceleration with magnetic obliquity dependency
In this Sod test, we let the acceleration efficiency η(θ_{B}) vary with the preshock magnetic obliquity angle θ_{B} and imposedη = 0.5ζ(θ_{B}) (the previous Sod test wasrun with θ_{B} = 0°, i.e. the efficiency was η = η_{0} = 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^{−10}, 0, 5.77 × 10^{−11}), (10^{−10}, 0, 10^{−10}), (5.77 × 10^{−11}, 0, 10^{−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^{2} ≪ 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 Threedimensional 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^{−4}, and a pointlike explosion of energy E_{th} = 1 spread over the eight central cells in a box of size unity^{1}. 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 celltocell 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 sweptup 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 zaxis 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 sweptup 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 fingerlike features in the shocked material, which are produced by the discretised nature of the grid; amongst the postshock 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 CRlike in the diffuse bubble γ_{e} = γ_{CR}. What matters for the shock dynamics is the effective adiabatic index in the sweptup 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 obliquitydependent 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 obliquityindependent case. 
for γ_{CR} = 4∕3. In agreement with Pfrommer et al. (2017), we find that for the same setup, 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^{−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^{−10}. For the uniform magnetic field configuration, CRs are accelerated around polar caps along the xaxis of the box with maximum efficiency, and go to zero along the yaxis (or zaxis) 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 zaxes) is further away than where CRs are produced (xaxis) 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 obliquityindependent 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 xdirection than along the ydirection (zdirection) as a result of the dependency of the density jump to the adiabatic index of the gas (for strong shocks, for γ_{e} = 5∕3 and for γ_{e} = 7∕5).
Finally, the random magnetic field setup 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 sweptup shock layer. On average, the acceleration efficiency is reduced by a factor 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 obliquityindependent 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 selfsimilar 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^{−10} that of the thermal pressure) and let it build through the turbulencegenerated shocks. The simulations have a uniform 128^{3} 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^{−3} 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^{−12} erg cm^{−3}. The initial magnetic field was uniform and was set up in the xdirection 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^{3}. We did not allow for selfgravity 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 with k sampled in the range k = [1, 3]. The turbulence is applied intermittently with an autocorrelation time of 0.5 Myr and with a compressiontosolenoidal 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  and X_{CR}independent acceleration efficiency
We start with a batch of simulations where the acceleration efficiency does not depend on and X_{CR} (i.e. ξ = 1). We set up three different simulations: (i) without CR acceleration (i.e. η_{0} = 0, NoShock); (ii) with CR acceleration and η = η_{0} = 0.1 (i.e. where CR acceleration does not depend on magnetic, NoThetaB); (iii) with CR acceleration and η_{0} = 0.1 (i.e. where CR acceleration depends on magnetic obliquity, ThetaB); and with η_{0} = 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 typicalSNgenerated CR acceleration efficiencies corresponding to much larger values of the shock Mach number than we can capture here with this simplified setup. For the sake of a testable setup 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^{−1} (here measured at t = 20 Myr for the Streaming run).
Shocks are driven in sheets with moderate Mach numbers of –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_{diss}u_{2} ≃ 10^{44}–10^{45} erg Myr^{−1} pc^{−2}. 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 turbulencegenerated shocks in the box, with clustered regions of pressure at levels similar to or above the initial thermal pressure (P_{th,0} ≃ 10^{−12} erg cm^{−3}). 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 preshock 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 obliquitydependent 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 preshock magnetic fields are not randomly oriented, but show preferentially withinshockplane 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 plumelike 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^{−3}, 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 (top) and dissipated energy flux e_{diss}u_{2} (bottom) for the Streaming turbulent box, with η_{0} = 0.1 and , 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 –4. 
Fig. 14 Cosmicray pressure maps of the turbulent box simulation, with η_{0} = 0.1 and , in a thin plane within the xplane 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  and X_{CR}dependent acceleration efficiency
We show here the results of the turbulent box experiments, where this time the efficiency dependency 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^{−10} everywhere in the box at time t = 0, and that a normalisation (maximum) acceleration efficiency η_{0} = 0.1 was used throughout.
Figure 18 shows the evolution of the CR fluxweighted mean value of (top panel) and X_{CR} (middle panel), and the evolution of the energy fluxweighted mean acceleration efficiencies (bottom panel) , , and ζ(θ_{B}) as a function of time. The bulk of the CR energy is produced in shocks of with a slight decrease over time. As CRs are produced, the upstream CRtothermal 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 protogalaxies), and delay the buildup of the CR pressure. Nonetheless, we show that our implementation of the , X_{CR} (and θ_{B}) dependency of η leads to interesting results in the buildup 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 (dotdashed lines) in the simulated turbulent ISM boxes for the simulations without shockacceleration (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 . 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 , 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 shockacceleration (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 . 
6 Conclusion
We have introduced a new modelling of anisotropic CR streaming and dynamical CR shockacceleration 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 multidimensional problems with nontrivial 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. CRmodified 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 longterm pressure support against compressed material, and streaming substantially modifies the smallscale 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, galacticwide 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 supernovadriven turbulence of the ISM, which we defer to future work.
Fig. 18 Top to bottom: cosmicray fluxweighted mean Mach number , CR fluxweighted mean CRtothermal pressure ratio X_{CR}, and dissipated energy fluxweighted efficiencies , 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 and X_{CR} 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 (ANR16CE310011) 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 setup as for the 2D sinusoid loop described in Sect. 3.2.2 with f_{iso} = 10^{−1} (left) or f_{iso} = 10^{−3} (right). 
Appendix B Tabulated values of ξ(, 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^{−3} to f_{iso} = 10^{−1} (to be compared with the value of f_{iso} = 10^{−2} 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^{−1} leads to more diffusion outside of the loop, which decreases the values of the maximum, while f_{iso} = 10^{−3} produces numerically driven fingerlike features but allows a more contained CR distribution in the loop.
Acceleration efficiencies interpolated values of from Kang & Ryu (2013).
Appendix C Shock numerical broadening
We show in Fig. C.1 a zoomedin 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 , 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 postshock 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 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 around the shock discontinuity for the Mach 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 preshock regions, the error on the evaluated Mach number becomes larger for a small kernel (n_{cell,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 cells (there are actually values of potential vectors drawn at nodes of the sampling cells), with the right, top, and backmost boundaries being replicates of the left, bottom, and frontmost 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 Bfield perpendicular to each face of AMR cells) is obtained by taking the rotational of the potential vector of the facesurrounding edges. This procedure guarantees that the magnetic field is random, ∇.B = 0, and the consistency of the coarsetofine values of the Bfield. 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
 Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aguilar, M., Aisa, D., Alpat, B., et al. 2015, Phys. Rev. Lett., 114, 171103 [CrossRef] [PubMed] [Google Scholar]
 Aharonian, F. A., Akhperjanian, A. G., Aye, K. M., et al. 2004, Nature, 432, 75 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bell, A. R. 1978, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 2015, MNRAS, 447, 2224 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhko, E. G., & Ellison, D. C. 1999, ApJ, 526, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1 [Google Scholar]
 Blasi, P. 2013, A&ARv, 21, 70 [Google Scholar]
 Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJ, 777, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Butsky, I. S., & Quinn, T. R. 2018, ApJ, 868, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Caprioli, D., Zhang, H., & Spitkovsky, A. 2018, J. Plasma Phys., 84, 715840301 [CrossRef] [Google Scholar]
 Castro, D., Slane, P., Patnaude, D. J., & Ellison, D. C. 2011, ApJ, 734, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Chevalier, R. A. 1983, ApJ, 272, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Marcowith, A., & Dubois, Y. 2019, A&A, 622, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Croston, J. H., Hardcastle, M. J., Birkinshaw, M., Worrall, D. M., & Laing, R. A. 2008, MNRAS, 386, 1709 [NASA ADS] [CrossRef] [Google Scholar]
 Croston, J. H., Kraft, R. P., Hardcastle, M. J., et al. 2009, MNRAS, 395, 1999 [NASA ADS] [CrossRef] [Google Scholar]
 Dashyan, G., & Dubois, Y. 2019, A&A, submitted [Google Scholar]
 Decourchelle, A., Ellison, D. C., & Ballet, J. 2000, ApJ, 543, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Diesing, R., & Caprioli, D. 2018, Phys. Rev. Lett., 121, 091101 [NASA ADS] [CrossRef] [Google Scholar]
 Dorfi, E. A. 1990, A&A, 234, 419 [NASA ADS] [Google Scholar]
 Drury, L. O. 1983, Rep. Prog. Phys., 46, 973 [NASA ADS] [CrossRef] [Google Scholar]
 Drury, L. O., & Voelk, J. H. 1981, ApJ, 248, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Dubois, Y., & Commerçon B. 2016, A&A, 585, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ehlert, K., Weinberger, R., Pfrommer, C., Pakmor, R., & Springel, V. 2018, MNRAS, 481, 2878 [NASA ADS] [CrossRef] [Google Scholar]
 Enßlin, T. A., Pfrommer, C., Springel, V., & Jubelgas, M. 2007, A&A, 473, 41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P [NASA ADS] [CrossRef] [Google Scholar]
 Farber, R., Ruszkowski, M., Yang, H.Y. K., & Zweibel, E. G. 2018, ApJ, 856, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Farmer, A. J.,& Goldreich, P. 2004, ApJ, 604, 671 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrand, G., Decourchelle, A., Ballet, J., Teyssier, R., & Fraschetti, F. 2010, A&A, 509, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fujita, A., & Mac Low, M.M. 2018, MNRAS, 477, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Fujita, Y., & Ohira, Y. 2011, ApJ, 738, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Girichidis, P., Naab, T., Walch, S., & Hanasz, M. 2014, ArXiv eprints, unpublished [arXiv:1406.4861] [Google Scholar]
 Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJ, 816, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Girichidis, P., Naab, T., Hanasz, M., & Walch, S. 2018, MNRAS, 479, 3042 [NASA ADS] [CrossRef] [Google Scholar]
 Günter, S., Yu, Q., Krüger, J., & Lackner, K. 2005, J. Comput. Phys., 209, 354 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, F., & Mathews, W. G. 2011, ApJ, 728, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasz, M., Kowal, G., OtmianowskaMazur, K., & Lesch, H. 2004, ApJ, 605, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasz, M., OtmianowskaMazur, K., Kowal, G., & Lesch, H. 2009a, A&A, 498, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hanasz, M., Wóltański, D., & Kowalik, K. 2009b, 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]
 Helder, E. A., Vink, J., Bassa, C. G., et al. 2009, Science, 325, 719 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Holguin, F., Ruszkowski, M., Lazarian, A., Farber, R., & Yang, H.Y. K. 2019, MNRAS, 490, 1271 [NASA ADS] [CrossRef] [Google Scholar]
 Jacob, S., & Pfrommer, C. 2017, MNRAS, 467, 1478 [NASA ADS] [Google Scholar]
 Jacob, S., Pakmor, R., Simpson, C. M., Springel, V., & Pfrommer, C. 2018, MNRAS, 475, 570 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, F. C., & Ellison, D. C. 1991, Space Sci. Rev., 58, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Jubelgas, M., Springel, V., Enßlin, T., & Pfrommer, C. 2008, A&A, 481, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kang, H., & Jones, T. W. 2005, ApJ, 620, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Kang, H., & Ryu, D. 2013, ApJ, 764, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, K., Petre, R., Gotthelf, E. V., et al. 1995, Nature, 378, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Lazarian, A., & Beresnyak, A. 2006, MNRAS, 373, 1195 [NASA ADS] [CrossRef] [Google Scholar]
 Malkov, M. A., Diamond, P. H., Sagdeev, R. Z., Aharonian, F. A., & Moskalenko, I. V. 2013, ApJ, 768, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Mao, S. A., & Ostriker, E. C. 2018, ApJ, 854, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Marcowith, A., Bret, A., Bykov, A., et al. 2016, Rep. Prog. Phys., 79, 046901 [NASA ADS] [CrossRef] [Google Scholar]
 Miniati, F. 2001, Comput. Phys. Commun., 141, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Miniati, F., Ryu, D., Kang, H., et al. 2000, ApJ, 542, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Miniati, F., Ryu, D., Kang, H., & Jones, T. W. 2001, ApJ, 559, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Nava, L., Gabici, S., Marcowith, A., Morlino, G., & Ptuskin, V. S. 2016, MNRAS, 461, 3552 [NASA ADS] [CrossRef] [Google Scholar]
 Nava, L., Recchia, S., Gabici, S., et al. 2019, MNRAS, 484, 2684 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pais, M., Pfrommer, C., Ehlert, K., & Pakmor, R. 2018, MNRAS, 478, 5278 [NASA ADS] [CrossRef] [Google Scholar]
 Pakmor, R., Pfrommer, C., Simpson, C. M., & Springel, V. 2016, ApJ, 824, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Pfrommer, C., Enßlin, T. A., Springel, V., Jubelgas, M., & Dolag, K. 2007, MNRAS, 378, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Pfrommer, C., Enßlin, T. A., & Springel, V. 2008, MNRAS, 385, 1211 [NASA ADS] [CrossRef] [Google Scholar]
 Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., & Springel, V. 2017, MNRAS, 465, 4500 [NASA ADS] [CrossRef] [Google Scholar]
 Pierre Auger Collaboration (Abraham, J., et al.) 2007, Science, 318, 938 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Ptuskin, V. S., Zirakashvili, V. N., & Plesser, A. A. 2008, Adv. Space Res., 42, 486 [NASA ADS] [CrossRef] [Google Scholar]
 Recchia, S., Blasi, P., & Morlino, G. 2017, MNRAS, 470, 865 [NASA ADS] [CrossRef] [Google Scholar]
 Ruszkowski, M., Yang, H. Y. K., & Reynolds, C. S. 2017a, ApJ, 844, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Ruszkowski, M., Yang, H.Y. K., & Zweibel, E. 2017b, ApJ, 834, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Ryu, D., Kang, H., Hallman, E., & Jones, T. W. 2003, ApJ, 593, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312 [NASA ADS] [CrossRef] [Google Scholar]
 Salem, M., Bryan, G. L., & Hummels, C. 2014, ApJ, 797, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Schaal, K., & Springel, V. 2015, MNRAS, 446, 3992 [NASA ADS] [CrossRef] [Google Scholar]
 Sharma, P., Colella, P., & Martin, D. F. 2009, ArXiv eprints [arXiv:0909.5426] [Google Scholar]
 Sijacki, D., Pfrommer, C., Springel, V., & Enßlin, T. A. 2008, MNRAS, 387, 1403 [NASA ADS] [CrossRef] [Google Scholar]
 Simpson, C. M., Pakmor, R., Marinacci, F., et al. 2016, ApJ, 827, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Skillman, S. W., O’Shea, B. W., Hallman, E. J., Burns, J. O., & Norman, M. L. 2008, ApJ, 689, 1063 [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]
 Thomas, T., & Pfrommer, C. 2019, MNRAS, 485, 2977 [NASA ADS] [CrossRef] [Google Scholar]
 Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374 [NASA ADS] [CrossRef] [Google Scholar]
 Vazza, F., Brunetti, G., & Gheller, C. 2009, MNRAS, 395, 1333 [NASA ADS] [CrossRef] [Google Scholar]
 Vazza, F., Brüggen, M., Gheller, C., & Brunetti, G. 2012, MNRAS, 421, 3375 [NASA ADS] [CrossRef] [Google Scholar]
 Wadepuhl, M., & Springel, V. 2011, MNRAS, 410, 1975 [Google Scholar]
 Wagner, A. Y., Lee, J. J., Raymond, J. C., Hartquist, T. W., & Falle, S. A. E. G. 2009, ApJ, 690, 1412 [NASA ADS] [CrossRef] [Google Scholar]
 Warren, J. S., Hughes, J. P., Badenes, C., et al. 2005, ApJ, 634, 376 [NASA ADS] [CrossRef] [Google Scholar]
 Wentzel, D. G. 1968, ApJ, 152, 987 [NASA ADS] [CrossRef] [Google Scholar]
 Wiener, J., Oh, S. P., & Guo, F. 2013a, MNRAS, 434, 2209 [NASA ADS] [CrossRef] [Google Scholar]
 Wiener, J., Zweibel, E. G., & Oh, S. P. 2013b, ApJ, 767, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467, 906 [NASA ADS] [Google Scholar]
 Winner, G., Pfrommer, C., Girichidis, P., & Pakmor, R. 2019, MNRAS, 488, 2235 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., & Lazarian, A. 2002, Phys. Rev. Lett., 89, 281102 [NASA ADS] [CrossRef] [Google Scholar]
 Zank, G. P., Webb, G. M., & Donohue, D. J. 1993, ApJ, 406, 67 [NASA ADS] [CrossRef] [Google Scholar]
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 Δx^{1.87±0.08} as indicated by the dashed line. 

In the text 
Fig. 4 Cosmicray energy density maps at t = 0 (top left) and t = 0.02 (top right) for an initial angledependent sinusoid within a purely circular magnetic field with an Alfvén velocity of 1 and a resolution of 128^{2} cells. The energy is evolved with the streaming advectiondiffusion 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 as a functionof the Mach number for different values of the upstream CRtothermal pressure ratio X_{CR}. The values are obtained from the X_{CR} = 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 ncell_{max} = 2 cells (left panel) or ncell_{max} = 4 cells (right panel) to probe hydrodynamical values in the postshock and preshock regions. These shock tube tests do not model CRs. Preshock and postshock 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: zoomedin region over the shock and contact discontinuities for better clarity of the CR shockaccelerated 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 obliquitydependent CR acceleration efficiency (θ_{B} is the socalled obliquity: angle of the preshock 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 zoomedin region over the shock and contact discontinuities for better clarity of the CR shockaccelerated 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 obliquitydependent 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 obliquityindependent 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 selfsimilar 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 (top) and dissipated energy flux e_{diss}u_{2} (bottom) for the Streaming turbulent box, with η_{0} = 0.1 and , 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 –4. 

In the text 
Fig. 14 Cosmicray pressure maps of the turbulent box simulation, with η_{0} = 0.1 and , in a thin plane within the xplane 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 (dotdashed lines) in the simulated turbulent ISM boxes for the simulations without shockacceleration (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 . 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 , 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 shockacceleration (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 . 

In the text 
Fig. 18 Top to bottom: cosmicray fluxweighted mean Mach number , CR fluxweighted mean CRtothermal pressure ratio X_{CR}, and dissipated energy fluxweighted efficiencies , 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 and X_{CR} 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 setup as for the 2D sinusoid loop described in Sect. 3.2.2 with f_{iso} = 10^{−1} (left) or f_{iso} = 10^{−3} (right). 

In the text 
Fig. C.1 Pressure profiles at time around the shock discontinuity for the Mach 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 preshock regions, the error on the evaluated Mach number becomes larger for a small kernel (n_{cell,max} = 2). 

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.