Shock-accelerated cosmic rays and streaming instability in the adaptive mesh refinement code Ramses

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 non-thermal 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 obliquity-dependent 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.


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 andMarcowith et al. 2016 for a recent review). Recent advances in the numerical modelling of DSA through hybrid particle-incell 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 ). 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 galax-ies (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(Nava et al. , 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;Girichidis et al. 2016Girichidis et al. , 2018Simpson 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, sub.). 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. 2018;Butsky & Quinn 2018), or taking into account the unresolved multi-phase nature of the gas and its impact on CR transport . CRs also boost the dynamo amplification of the magnetic field in disc galaxies (Hanasz et al. 2004(Hanasz et al. , 2009aPakmor et al. 2016).
On very large cosmological scales CRs are released in shocks (Miniati et al. 2000Ryu et al. 2003;Skillman et al. 2008;Pfrommer et al. 2007Pfrommer et al. , 2008Pfrommer et al. , 2017Vazza et al. 2009Vazza et al. , 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 et al. 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 by including 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 section 2, we introduce the full set of CR magnetohydrodynamics including the streaming and acceleration terms, whose numerical modelling and tests are respectively tackled in Sections 3 and 4. We finally test CR acceleration and streaming combined in turbulent interstellar medium experiments in section 5.

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 to be solved for cosmic-ray magneto-hydrodynamics (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 right-hand 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 L rad = L rad,th + L rad,CR−>th composed of the thermal L rad,th and CR L rad,CR−>th radiative loss terms, where the CR loss term (L rad,CR−>th = L rad,CR + H rad,CR−>th ) is the non-conserving sum of radiative losses from cosmic rays L rad,CR turning as a heating rate H rad,CR−>th for the thermal component. Finally, and this is the core of this paper, we detail how the streaming instability terms ∇. ((e CR + P CR )u st ) (advection-diffusion term) and L st (heating term), and the CR acceleration at shocks H 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 equation (3) are treated separately as source terms. The induction equation (equation 4) is solved using constrained transport , 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 equation (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 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 L rad,CR = 0.

Numerical implementation
Cosmic rays propagating faster than the Alfvén velocity u A = B/ 4πρ 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 SA u A sign(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: 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 time-step condition of ∆t = h∆x 2 /(2e CR u 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 which, when recast 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): 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: where D = D 0 + D st . The F CR,ds diffusion flux can be arbitrarily decomposed into an anisotropic and isotropic part where D = (1 − f iso )D, D iso = f iso D, and f iso ≤ 1. We briefly recall the 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 where the cell-centred fluxes are computed with cell-cornered values using the symmetric scheme from Günter et al. (2005): The anisotropic cell corner flux is F ani where barred quantities are arithmetic averages over the cells connected to the corner, i.e.
We note that all hydrodynamical variables in ramses are cellcentred except for the magnetic field which is face-centred. The streaming diffusion coefficient is computed as where upper tilde quantities stand for cell-centred quantities reconstructed from a combination of cell-centred and face-centred quantities: 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.

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.5 sin(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.

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.5 sin(θ) 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   tern 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. The result, not shown here, is qualitatively similar to that of our reference test.

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 M > M min , with M min 1.5. Keeping in mind these conditions, the Mach number of eligible cells is computed according to the criteria using upstream (preshock) 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 R P = P 2 /P 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 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 ) should be used instead: Here , 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 equation (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 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 (an oct 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.

Cosmic-ray acceleration at shocks
At shocks the kinetic energy flux of the upstream flow φ K,1 = 0.5ρ 1 u 3 1 (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 pre-shock kinetic energy e th,diss /(0.5ρu 2 1 ) 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 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 η(M, X CR , θ 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 R ρ ), we replace u 2 by Mc s,1 /R ρ , 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 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 R ρ the jump density ratio. The jump density ratio is obtained from the direct evaluation of the upstream and downstream densities The R γ ρ and R γ 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 1 .n s ). The dependency of the efficiency of CR acceleration with this so-called 'magnetic obliquity' can be factorised out, η(M, X CR , θ B ) = η 0 ξ(M, X CR )ζ(θ B ), and approximated by the following functional form (Pais et al. 2018): 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 ξ(M, X CR ) is obtained from the results of Kang & Ryu (2013), and is an increasing function of both M and X CR . They provide values of the acceleration efficiency for two 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 ξ(M, X 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 M 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 ξ(M, X CR ) 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 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).

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 = 634990). 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.

Cosmic-ray acceleration with constant efficiency
In this test we set up the previous 1D Sod shock tube test with Mach number 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 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   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 top to 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.
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 PdV compression). Apart from this expected effect, pressures, velocity, density, and the effective adiabatic index of the gas are accurately reproduced.

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 was run 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.

Three-dimensional Sedov explosion
We set up a 3D Sedov explosion with the following unitless values: a background at rest with gas density of ρ = 1, P th = 10 −4 , and a point-like 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 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-Laxvan 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 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 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). The top and bottom panels show 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 respect to the obliquity-independent case. the acceleration efficiency) 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 the shock 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 random magnetic 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 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 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 x-direction than along the y-direction (zdirection) as a result of the dependency of the density jump to the adiabatic index of the gas (for strong shocks, R ρ = 4 for γ e = 5/3 and R ρ = 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 < ζ >= π/2 0 ζ(θ B ) sin θ B dθ B 0.302 for a purely random 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 Article number, page 9 of 17 A&A proofs: manuscript no. article  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. 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.

Turbulent box of the interstellar medium
We ran turbulent interstellar medium (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 turbulence-generated 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 ×  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. 10 −12 erg cm −3 . 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 3 . 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 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 spacef (k) ∝ 1 − (k − k 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).

Mand X CR -independent acceleration efficiency
We start with a batch of simulations where the acceleration efficiency does not depend on M 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, No-ThetaB); 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 typical SN-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 −1 (here measured at t = 20 Myr for the Streaming run).
Shocks are driven in sheets with moderate Mach numbers of M 3-4, as can be seen in Fig. 13 for the Streaming run (other  Figure 14 shows maps of the CR pressure at two different times t = 10 and 20 Myr for the simulation NoThetaB, ThetaB, and Streaming. 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 −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 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 section 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 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  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 No-ThetaB 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 −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 the CR 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.

Mand X CR -dependent acceleration efficiency
We show here the results of the turbulent box experiments, where this time the efficiency dependency ξ(M, X 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  Fig. 16. Probability density function of the magnetic obliquity in the ISM boxes, and with η 0 = 0.1 and ξ(M, X 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 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 ξ(M, X CR ) = 1. 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 flux-weighted mean value of M (top panel) and X CR (middle panel), and the evolution of the energy flux-weighted mean acceleration efficiencies (bottom panel) η = η 0 ξ(M, X CR )ζ(θ B ), ξ(M, X CR ), and ζ(θ B ) as a function of time. The bulk of the CR energy is produced in shocks of M 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 nearly random 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 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.

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 mecha-  Fig. 18. Top to bottom: Cosmic-ray flux-weighted mean Mach number M, CR flux-weighted mean CR-to-thermal pressure ratio X CR , and dissipated energy flux-weighted efficiencies η = η 0 ξ(M, X CR )ζ(θ B ), ξ(M, X 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 M and X CR as extrapolated from the values of Kang & Ryu (2013). nism 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  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 interstellar medium 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 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, galactic-wide outflows driven by CRs (Dashyan & Dubois, sub.), 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.
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 3 pot cells (there are actually (n pot + 1) 3 values of potential vectors drawn at nodes of the n 3 pot 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 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 (n cell,max = 2). 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.  X CR = 0 X CR = 0.025 X CR = 0.05 X CR = 0.1 X CR = 0.2 X CR = 0.5 X CR = 1 M = 2 4.44 × 10 −4 3.80 × 10 −2 7.55 × 10 −2 1.51 × 10 −1 3.01 × 10 −1 7.51 × 10 −1 1.00 M = 3 2.66 × 10 −2 1.47 × 10 −1 2.66 × 10 −1 5.06 × 10 −1 9.86 × 10 −1 1.00 1.00 M = 4 2.00 × 10 −1 4.11 × 10 −1 6.22 × 10 −1 9.08 × 10 −1 1.00 1.00 1.00 M = 5 4.44 × 10 −1 5.60 × 10 −1 6.76 × 10 −1 9.08 × 10 −1 1.00 1.00 1.00 M = 7 6.66 × 10 −1 7.33 × 10 −1 8.00 × 10 −1 9.33 × 10 −1 1.00 1.00 1.00 M = 10 8.66 × 10 −1 8.89 × 10 −1 9.10 × 10 −1 9.55 × 10 −1 1.00 1.00 1.00 M = 20 1.00 1.00 1.00 1.00 1.00 1.00 1.00 M = 30 1.00 1.00 1.00 1.00 1.00 1.00 1.00 M = 50 1.00 1.00 1.00 1.00 1.00 1.00 1.00 M = 100 1.00 1.00 1.00 1.00 1.00 1.00 1.00