Issue 
A&A
Volume 585, January 2016



Article Number  A138  
Number of page(s)  10  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201527126  
Published online  12 January 2016 
An implicit scheme for solving the anisotropic diffusion of heat and cosmic rays in the RAMSES code
^{1}
Institut d’Astrophysique de Paris, UMR 7095, CNRS, UPMC Univ. Paris VI,
98bis boulevard
Arago,
75014
Paris,
France
email:
dubois@iap.fr
^{2}
École Normale Supérieure de Lyon, CRAL, UMR 5574 du CNRS,
Université de Lyon I, 46 allée
d’Italie, 69364
Lyon Cedex 07,
France
email:
benoit.commercon@enslyon.fr
Received: 5 August 2015
Accepted: 31 October 2015
Astrophysical plasmas are subject to a tight connection between magnetic fields and the diffusion of particles, which leads to an anisotropic transport of energy. Under the fluid assumption, this effect can be reduced to an advectiondiffusion equation, thereby augmenting the equations of magnetohydrodynamics. We introduce a new method for solving the anisotropic diffusion equation using an implicit finitevolume method with adaptive mesh refinement and adaptive timestepping in the ramses code. We apply this numerical solver to the diffusion of cosmic ray energy and diffusion of heat carried by electrons, which couple to the ion temperature. We test this new implementation against several numerical experiments and apply it to a simple supernova explosion with a uniform magnetic field.
Key words: conduction / diffusion / magnetohydrodynamics (MHD) / cosmic rays / plasmas / methods: numerical
© ESO, 2016
1. Introduction
Diffusive processes are relevant at various scales in astrophysical systems from the diffusion of cosmic rays (CRs) in the interstellar medium (e.g. Strong et al. 2007) to the the diffusion of heat and CRs in cluster of galaxies (e.g. Rosner & Tucker 1989; Fabian 1994; Narayan & Medvedev 2001). As the diffusion of these charged particles operate in a magnetised plasma, the diffusion process becomes anisotropic since it is conducted along magnetic field lines. Furthermore, for a fully ionised plasma, like the lighter population of particles, electrons conduct heat much faster than protons, which are virtually in the limit of nondiffusion. In some situations, it leads to a nonthermal equilibrium between ions and electrons if the collisional coupling time of the two temperatures is longer than the diffusion timescale of electrons. Thus, there could be fluctuations in electron temperature without the corresponding fluctuation in the ion temperature.
Those processes have mostly been ignored from the numerical astrophysics perspective with a few exceptions. In the intracluster medium, it has been shown that ions and electrons are not perfectly at thermal equilibrium (Chièze et al. 1998; Courty & Alimi 2004; Rudd & Nagai 2009; Wong & Sarazin 2009; Gaspari & Churazov 2013) and that the transport of heat produces an electron precursor at the virial radius of massive dark matter halos (Teyssier et al. 1998). The transport of heat eventually reduces the cooling flow at the core of galaxy clusters, transporting the cosmic shockdriven internal energy of the outskirts towards the centre (Ruszkowski & Begelman 2002; Jubelgas et al. 2004). Thermal and magnetohydrodynamical instabilities form because of the anisotropic nature of the conduction, the socalled magnetothermal instability (Balbus 2000; Parrish et al. 2008) and the heatingbuoyancy instability (Parrish & Quataert 2008; Bogdanović et al. 2009; Parrish et al. 2009), and their variations with the adjunction of the diffusion of CRs (Rasera & Chandran 2008). Conduction may also play an important role in redistributing the mechanical energy deposited by active galactic nuclei jets (Ruszkowski & Begelman 2002; Brighenti & Mathews 2003; Brüggen et al. 2005; McNamara & Nulsen 2007), helping to solve the cooling catastrophe in galaxy clusters. A similar heating source in the intracluster medium is that of comic rays deposited in active galactic nuclei or at shocks and redistributed by diffusion (Miniati et al. 2001; Enßlin et al. 2011).
On galactic scales, since CR energy in the interstellar medium (ISM) is at equipartition with kinetic, thermal, and magnetic energies (Beck & Krause 2005), they may play a crucial role in the selfregulation of the ISM dynamics and, thus, of star formation. For instance, the injection of CRs within remnants of supernovae lead to stronger galacticscale winds (Jubelgas et al. 2008; Uhlig et al. 2012; Hanasz et al. 2013; Booth et al. 2013; Salem & Bryan 2014) and to amplified galactic dynamos with anisotropic diffusion (Hanasz et al. 2004). Anisotropic conduction in the ISM affects the shape and size of cold clouds (Koyama & Inutsuka 2004; Piontek & Ostriker 2004; Choi & Stone 2012), and the expansion of supernova remnants (Tilley et al. 2006; Balsara et al. 2008a). Also, CRs can penetrate deep inside cold dense cores and, as a consequence, can regulate the ionisation rates of the gas (e.g. Spitzer & Tomasko 1968; Padovani et al. 2009).
One problem with the implementation of diffusion in numerical simulation is that the stability criterion is Δt_{diff} = Δx^{2}/ (2D_{diff}), where Δx is the cell size and D_{diff} the diffusion coefficient. For comparison, the hydrodynamical Courant Friedrich Levy (CFL) condition is Δt_{h} = CΔx/ (u + c_{s}), where u is the gas velocity, c_{s} the fluid sound speed, and C ≤ 0.8 is the Courant factor. Since the diffusion stability time step does not scale linearly with resolution, unlike the hydrodynamical time step, it could be a bottleneck to employ an explicit diffusion scheme that must verify this time step condition at any time, especially for multiscale problems where gravity shapes strong contrasts in gas densities and triggers refinements in resolution.
Implicit numerical solvers in these situations could be favoured. They do not need to fulfil the time step constraint on diffusion, and any time step can be chosen (such as the CFL time step), at the expense of numerical complexity and intensity compared to an explicit solver. There is a variety of numerical implementations of anisotropic diffusion from centred symmetric to centred asymmetric schemes (Günter et al. 2005) using various slope limiters to preserve the monotonicity of the system (Sharma & Hammett 2007). Amongst all these implementations, implicit and explicit solvers have been developed for astrophysical codes: Hanasz & Lesch (2003), Parrish & Stone (2005), and Rasera & Chandran (2008) implemented an explicit method for anisotropic conduction and diffusion, Balsara et al. (2008b) developed a semiimplicit solver, Yokoyama & Shibata (2001) and Meyer et al. (2012) a fully implicit solver.
In this paper, which uses the ramses code (Teyssier 2002), we present the first implementation of an implicit solver for anisotropic diffusion or conduction on adaptive mesh refinement (AMR) with adaptive timestepping. This new numerical implementation is augmented by modelling a multitemperature and energy component with temperature coupling of ions and electrons. In Sect. 2, we present the new numerical solver for diffusion or conduction and temperature coupling. In Sect. 3, we test our method with several numerical experiments. We discuss future developments and applications of this method in Sect. 4.
2. Numerical setup
2.1. Magnetohydrodynamics with heat and cosmic ray diffusion
The set of differential equations to be solved for magnetohydrodynamics of a fluid made of ions and electrons with heat conduction and CRs diffusion is where ρ is the gas mass density, u the gas velocity, B the magnetic field, e = 0.5ρu^{2} + e_{th} + e_{cr} + B^{2}/ 8π is the total energy density, e_{th} = e_{I} + e_{E} is the thermal energy density, e_{E} and e_{I} are the electron and ion energy densities, e_{cr} is the energy density of CRs, and total pressure p_{tot} = (γ − 1)e_{th} + (γ_{cr} − 1)e_{cr} + 0.5B^{2}/ 4π, where γ and γ_{cr} are the adiabatic indexes of the gas and CRs, respectively. All energy components e_{i} are energies per unit volume e_{i} = E_{i}/ Δx^{3}, where Δx is the cell size. The heat flux F_{cond} is carried by electrons alone, and the temperature of ions and electrons couple through the energy exchange term ℋ_{EI}. Since ions are heavier than electrons, their thermal velocity is lower by a factor ∝. Therefore, for a given increase in thermal velocity at a shock, it corresponds to a higher temperature increase for ions than for electrons. It justifies the fact that electron energy is treated in a separate energy equation that does not fulfil the jump condition. The same arguments apply to CRs compared to the thermal gas since velocities of CRs are much higher than thermal velocities of the ions (a.k.a. CRs are a nonthermal gas component).
We assume that the thermal components of the gas are made only of ons and electrons, meaning that the gas is fully ionised. In some situations, this approximation will break down (starforming regions), but it is a good approximation for the warm or hot phases of the ISM and for the intergalactic gas. CRs also diffuse through the diffusion flux term F_{CR}. Equation (6) can be repeated for as many CR energy components required to decompose the CR energy power spectrum (with different diffusion coefficients). Introducing multiple CR energy bands requires including adiabatic and radiative losses for a proper treatment of the CR power spectrum, so we defer this to future work and simply assume that CRs are represented by one fluid with one effective diffusion coefficient and one effective adiabatic index.
We use the AMR ramses code detailed in Teyssier (2002). Equations (5) and (6) for electrons and CRs, respectively, are added to the set of MHD equations. The full set of equations is solved with the standard MHD solver of ramses described in Fromang et al. (2006), where the righthand terms of Eq. (3) are treated separately as source terms. The ion internal energy is obtained for free by subtracting the electron internal energy from the total thermal energy. The induction (Eq. (4)) is solved using constrained transport (Teyssier et al. 2006). Godunov fluxes are modified to account for the extra energy components and total pressure made of ions, electrons, and CRs. Therefore, the effective sound speed used for the CFL time step condition accounts for the extra pressure components (i.e. total pressure of the fluid).
2.2. Anisotropic diffusion
We explain the method for the conduction of heat through electrons, but this can be applied to a single temperature model (ions and electrons at thermal equilibrium) or for the diffusion of CRs. The difference for CRs is that the diffusion coefficient D_{CR} is constant throughout the simulation volume and that it applies to the gradient of CR energy instead of the gradient of temperature, e.g. F_{CR} = −D_{CR}b(b.∇)e_{CR}, where b is the magnetic field unit vector.
In the presence of magnetic fields, the conduction of heat in a plasma is (7)where κ_{iso} and κ_{∥} are the isotropic conduction and parallel along the magnetic field lines coefficients respectively, with κ_{∥} = κ_{Sp} − κ_{iso}, and κ_{iso} = f_{iso}κ_{Sp} with f_{iso} an arbitrary small coefficient for anisotropic diffusion (f_{iso} = 1 is for isotropic diffusion). In many astrophysical cases, κ_{iso}/κ_{∥} ≪ 1 since the Larmor radius is much smaller than the mean free path of electrons. For instance, in the hot gas of galaxy clusters with temperature T_{E} = 3 keV, electron density n_{E} = 10^{2} cm^{3} and B = 1 μG, the Larmor radius is λ_{L} = 10^{8} cm and the mean free path λ_{mfp} = 10^{21} cm. However, to ensure numerical stability, we set up the isotropic conduction term to be around one per cent of the parallel conduction coefficient.
The conduction coefficient for electrons is the Spitzer (1956) value (8)where k_{B} is the Boltzmann constant and D_{C} the thermal diffusivity (9)
2.2.1. Implicit scheme
This diffusion step is performed after the magnetohydrodynamics step. Here, we explain how the anisotropic conduction and diffusion equation is solved for a uniform grid, and the AMR part of the solver is explained in Sect. 2.2.2. Given the restrictive stability condition for an explicit scheme on the Eq. (7), Δt_{diff} = Δx^{2}/ (2D), compared to the behaviour of the CFL condition for the hydrodynamics time step Δt_{h} = Δx/ (u + c_{s}), using an implicit solver on the equation diffusion will alleviate this numerical stability constraint. Thus, the diffusion term is solved over one hydrodynamical time step using implicit integration. Discretising Eq. (7) in 1D leads to (we omit subscripts E and cond for clarity) (10)where the subscript i is for the cell position, and the superscript n for the time. Each of the fluxes F^{n + 1} is evaluated at the cell interface (2 in 1D, 4 in 2D, and 6 in 3D) at time n + 1. This can be rewritten as (11)where e_{i} = n_{i}k_{B}T/ (γ − 1) (n_{i} here is the gas number density, not time index n). Equation (11) implicitly assumes that the cells are cubic, which is the case in ramses. We use the value of κ at time n since we employ an implicit solver, which requires the constancy of the coefficient during the integration step. Thus, defining and C_{v,i} = n_{i}k_{B}/ (γ − 1), we obtain (12)We discretise Eq. (7) in 2D (3D can be obtained from 2D with little effort) assuming the cells have the same extent in x, y, z directions (13)for cell position i,j. These quantities are evaluated with the centred symmetric scheme proposed by Günter et al. (2005) for the anisotropic part of the flux. The anisotropic flux at cell interfaces and are evaluated from their cell corner fluxes , thus The anisotropic cell corner flux is (14)where barred quantities are arithmetic averages over the cells connected to the corner; i.e., For the isotropic part of the flux, we use a classical discretisation, where the fluxes of Eq. (13) are simply obtained from the left and right states of the cell interface; i.e., (15)With this numerical scheme, the matrix system Ax = c formed by Eq. (12), with A the matrix including the C^{n} and coefficients, x the vector of temperature T^{n + 1}, and c the vector of energy densities can be generalised in multidimensions by Eq. (13). The sparse matrix A is positive definite and symmetric, so we can use the conjugate gradient algorithm to solve this system of linearised equations as in Commerçon et al. (2011; for radiation hydrodynamics in their case).
2.2.2. Diffusion with adaptive mesh refinement
We follow the procedure introduced in Commerçon et al. (2014) for solving the diffusion equation with AMR and adaptive timestepping on a levelbylevel basis. Each level ℓ is evolved with a time step Δt^{ℓ} that is twice less than the coarser level ℓ − 1, Δt^{ℓ} = Δt^{ℓ − 1}/2, meaning that level ℓ evolves with two consecutive time steps before doing one time step of level ℓ − 1. The adaptive timestepping is performed for all solvers in ramses including the diffusion solver presented in this paper, and finer levels are updated first.
At a given level of refinement ℓ, there are two types of nonuniform interfaces: the finetocoarse interface (interface between a cell at level ℓ and a cell at ℓ − 1) and the coarsetofine interface (interface between a cell at level ℓ and a cell at ℓ + 1). In both cases, we use Dirichlet boundary conditions: cell values at level boundaries are imposed. One could choose Neumann boundary conditions, which impose the fluxes and guarantee energy conservation as what is done for the hydrodynamical solver. However, as shown in Commerçon et al. (2014), this sort of imposed flux conditions could lead to negative values for lang time steps.
For the finetocoarse interface, we use values of ℓ − 1 at time n as imposed boundary conditions for level ℓ. With the minmod scheme (van Leer 1979; monotonised central is also a valid choice), we interpolate the values of level ℓ − 1 on a finer virtual grid at level ℓ to impose the finetocoarse boundary. For the coarsetofine interface, we use values of ℓ + 1 at time n + 1 as imposed boundary conditions for level ℓ. We restrict the value of the boundary coarse cell at level ℓ to the average value of the 2^{dim} cells of level ℓ + 1 to impose the coarsetofine boundary. In any case, Eq. (13) remains correct since the levelbylevel diffusion solver only deals with data estimated at the same level of refinement. The combination of the use of Dirichlet boundary conditions at level interfaces and of interpolation or restriction operations does not break the symmetry of matrix A. The imposed values of the neighbouring cells at different refinement levels go into the righthand side (vector c) of the matrix system Ax = c, because it is the case for the computational domain imposed boundary conditions.
2.2.3. Limitations of the method
There are two main limitations to our current anisotropic diffusion solver: the nonenergy conservation due to Dirichlet boundary at level interfaces and the nonmonotonicity preserving nature of the solver. To circumvent the first limitation, one could adopt a unique time step strategy, where diffusion is solved for all levels at once. This guarantees energy conservation but slows down the calculation. A good compromise would be to do the diffusion step rather than every fine time step of the simulation, but only over coarse time steps (or even every 10, 100, etc. coarse time steps). Nevertheless, as shown in Commerçon et al. (2014), Dirichlet boundary conditions are robust in most cases, and in practice, energy conservation is acceptable (see tests below).
Concerning the nonmonotonicity preserving nature of the diffusion solver, since the flux at a cell interface (Eq. (14)) includes a transverse gradient of temperature, it is not guaranteed that the flux flows from the high temperature cell to the low temperature cell. Thus, the temperature can become negative, and the monotonicity of the solution is not preserved (see for instance the test problem illustrated in Fig. 5 of Sharma & Hammett 2007). Sharma & Hammett (2007) propose a slope limiter for the centred symmetric scheme with explicit time integration that preserves the monotonicity of the solution. However, employing a slope limiter would not allow us to use the conjugate gradient algorithm any longer since matrix A becomes nonsymmetric. In our case, we prefer to conserve the simple and fast solution offered by the conjugate gradient algorithm and to fix the negative temperatures created by the anisotropic conduction solver with tiny positive values. Future developments will include slope limiters, and the matrix system will be solved using a biconjugate gradient algorithm instead (as in González et al. 2015).
2.3. Ion and electron temperature coupling
We now explain how we model the coupling of the temperature of ions and electrons in ramses. This coupling step is performed after the diffusion step. The energy exchange rate ℋ_{EI} introduced in Eq. (5) is written for a gas that we assume to be fully ionised and composed of hydrogen and helium (16)with the equilibration timescales (17)and (18)where m_{E} and m_{I} stand for the electron and ion mass, respectively; n_{I} = ρ/ (μ_{I}m_{p}) is the ion number density with m_{p} the proton mass and μ_{I} being the ion mean molecular weight; Z_{I} is the ion charge number; q_{E} the electron charge; and lnΛ = 40 is the Coulomb logarithm. Since m_{E} ≪ m_{I}, we neglect the T_{I}/m_{I} term in the energy exchange rates. For a mixture of hydrogen and helium, the timescale of the temperature coupling between hydrogen and helium would be proportional to . Thus, the equilibration between hydrogen and helium operates 10^{5} times faster than the equilibration timescale between hydrogen and electrons or helium and electrons. For hydrogen or helium, the ratio of also simplifies to m_{p} (not the case for heavier elements). Therefore, protons and HeIII can be considered as one single population sharing the same temperature.
Writing , we see that (19)The variation in energy is symmetric between the transfer of ion energy and electron energy, but this is not the case for the variation in temperature, which is . Therefore, the difference in temperature change between ions and electrons will be a factor μ_{I}/μ_{E}, where μ_{E} is the “mean molecular weight of electrons” μ_{E} = ρ/ (n_{E}m_{p}). For astrophysical plasmas essentially composed of hydrogen and helium, the change in temperature between ions and electrons is close to symmetric.
To solve for the system of two nonlinear coupled equations, we rewrite (16) as (20)where . We introduce the vector , and rewrite the previous two equations as The goal is now to iterate on ΔX to find the value of G_{1}(X^{k} + ΔX) = 0 and G_{2}(X^{k} + ΔX) = 0 as required by equations (20). The derivatives and are obtained by differentiating equations (20) with respect to and , where the superscript k replaces the superscript n + 1 in (20). The values of ΔX are obtained with Cramer’s rule, and ΔX is added to the value of X^{k} until the variation in ΔX/X^{k}< 10^{4}.
3. Numerical tests
If not indicated, we use a courant factor of C = 0.8 with adaptive timestepping between the different levels of refinement.
Fig. 1 Top: diffusion of a pressure step function as a function of position at t = 9.3 × 10^{4} (black), t = 1.9 × 10^{3} (red), t = 5.6 × 10^{3} (blue) for the simulation (symbols) and the analytic solution (solid lines). The extra levels of refinement are indicated as dashed lines. Bottom: relative errors of the simulation with respect to the analytic solution. 
3.1. Diffusion of a step function
We performed a simple 1D test of the diffusion of a step function with a constant diffusion coefficient of D = 1 without hydrodynamics. The initial setup is a pressure of P_{1} = 0.8 between x = [ 0.25,0.75 ] and P_{0} = 0.4 outside with C_{V} = ρk_{B}/ (μm_{p}) = 1. We start with a minimum level of refinement of ℓ_{min} = 3 and refine up to level ℓ_{max} = 7 wherever the relative pressure variation is larger than 10%. The time step is chosen to be Δt_{7} = 7 × Δt_{diff} (for level 7), where Δt_{diff} is the stability time step for diffusion (only useful for explicit schemes).
The time evolution of the step function has the following analytical solution for a constant diffusion coefficient (21)where x_{0} = 0.25. Figure 1 shows the result of the diffusion of the step function with the analytical solution at three different times. The analytical solution is reproduced well, below a 10% relative error.
Fig. 2 Electron temperature (red or diamonds) and ion temperature (blue or pluses) as a function of time for the one cell problem. The solution for a time step of τ_{eq,EI,0}/ 50 is represented in solid lines and for the time step of 20τ_{eq,EI,0} with symbols. 
3.2. Temperature coupling
Figure 2 shows the result of the one cell problem with two different initial temperatures for electrons T_{E,0} = 10^{10} K and ions T_{I,0} = 10^{8} K. The solid lines indicate the solution with a time step of τ_{eq,EI,0}/ 50, and points are the solution with a time step of 20τ_{eq,EI,0}. They are in excellent agreement, even for the run where the time step is 20 times the timescale of temperature coupling.
Fig. 3 Solution of the Sod (1978) shock tube problem with two temperatures. There is no diffusion and coupling of the temperatures. We plot the velocity (top left), pressure (bottom left), and density (bottom right) at t = 0.245. The squares correspond to the result of the simulation and the red lines to the analytical solution. The light blue and dark blue squares in the pressure plot stand for the first and second species, respectively, while the black squares are for the total pressure. The level of refinement is plotted in the top right panel. 
Fig. 4 Stationary cosmological shock experience including conduction and twotemperature coupling. The result of the twotemperature run (blue) is represented with dotted lines for the ion temperature (pressure), in dashed for the electron temperature (pressure), and solid for the average temperature (total pressure) of the gas at t = 10 Gyr. For comparison, the results of the single temperature runs are plotted in red with conduction and in black without conduction. An electron precursor forms in the upstream part of the shock front. 
3.3. Sod test with two energy components
We set up the conditions for a Sod (1978) shock tube with two temperatures and without conduction or temperature coupling. The initial left and right gas density is ρ_{L} = 1, ρ_{R} = 0.125, the first pressure P_{1,L} = 0.34, P_{1,R} = 0.066 (gas or ion), the second pressure P_{2,L} = 0.66, P_{2,R} = 0.034 (respectively CR or electron), thus P_{L} = 1 and P_{R} = 0.1, and zero velocity. The adiabatic index for the first and second temperatures is γ_{1,2} = 1.4. The minimum level of refinement is ℓ_{min} = 3, and we refine up to level ℓ_{max} = 10 wherever the relative gradient of any hydro variable is greater than 5%. Figure 3 shows the result of the simulation, which exhibits extremely good agreement with the analytical solution. At the contact discontinuity x = 0.7, the total pressure is continuous, but not the individual pressure terms. We reproduce the structure of a twocomponent shock well (see e.g. Pfrommer et al. 2006, for comparison).
3.4. Shock with two temperatures and heat conduction
We model the formation of a 1D stationary shock typically arising at the virial radius of a cluster of galaxies. The gas is of primordial cosmological composition with a 0.76 fraction of hydrogen and 0.24 of helium and is fully ionised. Thus, the mean molecular weight of ions, electrons, and gas are μ_{I} = 1.22, μ_{E} = 1.13, and μ = 0.59. The right state is with gas density n_{R} = 10^{4} cm^{3}, gas velocity u_{R} = 2000 km s^{1} oriented in the negative xdirection, and ion and electron temperatures at equilibrium T_{I,R} = T_{E,R} = T_{R} = 4 × 10^{7} K. Using the RankineHugoniot jump conditions, we find that the corresponding left state is n_{L} ≃ 2.3 × 10^{4} cm^{3}, u_{L} = 875 km s^{1}, and T_{I,R} = T_{E,R} = T_{R} ≃ 8.4 × 10^{7} K for an adiabatic index of γ = 5/3. The box length is 8 Mpc, with a minimum and a maximum level of refinement of ℓ_{min} = 4 and ℓ_{max} = 9, respectively, and with refinement triggered wherever the relative gradient of any hydro variable is larger than 10%. The left and right boundaries are imposed with the initial values of left and right states.
Figure 4 shows the result of the simulation with the two temperatures including conduction of the electrons and the temperature coupling between ions and electrons, as well as the solution for a single temperature with or without conduction. For the single temperature experiment without conduction, the initial left and right states are preserved at time t = 10 Gyr. With the addition of conduction, a thermal precursor forms ~1 Mpc ahead of the shock, which is present for both the single and the twotemperature experiments.
In the case of the two temperatures, we see that the electrons have the higher temperature value in the thermal precursor compared to the ions with a slow rise similar to the thermal precursor in the singletemperature experiment. Ion temperature in the precursor lies belowwith a sharp rise at the shock interfacebecause the coupling of the two temperatures is not instantaneous. The increase in ion temperature right after the shock interface compared to the average gas temperature is because the ions absorb the upstream kinetic energy at the shock, and electrons only lag behind and thermally recouple with the ions. Far away from the shock front, both temperatures are at equilibrium. The reader can refer to Zel’dovich & Raizer (1967) for the detailed description of a twotemperature shock with conduction.
3.5. Conduction in a circular loop
Anisotropy can be tested by the propagation of heat through a magnetic loop, here in 2D (see for comparison e.g. Sharma & Hammett 2007; Rasera & Chandran 2008). We initialised a patch of temperature of K wherever 2 <r< 3 kpc and 165 <θ< 195° in a background gas with temperature K and density n = 1 cm^{3}. The size of the box is 10×10 kpc, and we performed two runs either with a uniform grid 128^{2} or with AMR with ℓ_{min} = 4 and ℓ_{max} = 8 and refinement triggered if the relative gradient of pressure is larger than 20%. We used a Spitzer conductivity with κ_{∥} = 0.99κ_{Sp} and κ_{iso} = 0.01κ_{Sp}.
The result of this numerical experiment is shown in Fig. 5 for two runs with and without AMR. The initial patch of temperature has diffused anisotropically along the circularloop shape of the magnetic field lines. The results of the AMR and of the uniform grid runs are in very good agreement, as shown by the white and red contours.
We also explored the effect of varying the perpendicular conductivity coefficient with κ_{iso} = [ 0.001,0.01,0.1 ] κ_{Sp}. The resulting temperature maps are shown in Fig. 6. As expected, the spread of temperature perpendicular to the magnetic field lines increases with the increase in the isotropic component of the diffusion coefficient.
Fig. 5 Temperature map at t = 2 Gyr for the heat conduction in a 2D magnetic loop without AMR (grey levels). The solution with AMR is plotted with red contours for T = 5 × 10^{6} K and T = 10^{7} K (white contours without AMR). 
Fig. 6 Temperature map at t = 2 Gyr for the heat conduction in a 2D magnetic loop with AMR and different fractions of perpendicular conductivity coefficients κ_{iso} = 0.1κ_{Sp} (red contours), κ_{iso} = 0.01κ_{Sp} (grey levels and white contours), and κ_{iso} = 0.001κ_{Sp} (blue contours). Two contours are plotted for each simulation result corresponding to temperatures T = 5 × 10^{6} K and T = 10^{7} K (except for κ_{iso} = 0.1κ_{Sp} where only T = 5 × 10^{6} K is reached). 
3.6. Sovinec test
Sovinec et al. (2004) designed a test to measure the numerical perpendicular diffusion κ_{num} for the anisotropic diffusion in 2D. The energy evolution follows (22)where is a heating source term of the form . The magnetic field is chosen such that B.∇T_{E} = 0, i.e. the parallel heat flux component is suppressed, and only the perpendicular component remains. The analytical stationary solution at the centre is . For this test κ_{iso} = 1 (κ_{Sp} = 100), ρ = 1, , B_{x} = cos(πx)sin(πy), and B_{y} = −sin(πx)cos(πy). The box size ranges over [−0.5,0.5 ] × [ −0.5,0.5 ]. As in Rasera & Chandran (2008), we imposed the temperature at the boundaries of the simulation box with T_{bound} = 0.
We ran the simulation for various uniform grid resolutions from 16^{2} to 256^{2} up to the steadystate solution, and we measured the value of T_{mea}(0,0). We followed Sharma & Hammett (2007) to obtain the value of κ_{num} as follows κ_{num} = T_{iso}(0,0) /T_{mea}(0,0)−1, where T_{iso}(0,0) is the the central temperature calculated in the isotropic limit (κ_{iso} = κ_{Sp}). This measurement is more accurate than assuming that the isotropic diffusion gives T_{iso}(0,0) = 1 exactly. The result is shown in Fig. 7, and we see that the perpendicular numerical diffusion coefficient scales quadratically with resolution ∝ Δx^{2} and that, even at low resolution Δx = 1/16, the numerical diffusion is a factor 100 below our explicit perpendicular diffusion coefficient κ_{iso}.
Fig. 7 Results of the perpendicular numerical diffusion coefficient κ_{num} as a function of resolution Δx for the Sovinec et al. (2004) test problem in a solid line with square symbols. The numerical diffusion scales with Δx^{2} (dotted line). 
Fig. 8 Spherically averaged values of density (left), pressure (middle), and velocity (right) as a function of radius at t = 4 kyr for the 3D Sedov explosion without diffusion (top), and with isotropic diffusion (bottom). The solid lines are the result of the simulation, and the dashed line is the analytic prediction for the standard Sedov solution. The pressure is decomposed into total pressure (black), ion pressure (blue), electron pressure (red), and CR pressure (magenta). Without diffusion, the analytic solution is reproduced well. With diffusion, there are both a CR precursor and an electron precursor. 
Fig. 9 Same as Fig. 8 at t = 100 kyr. The shock front takes over the conduction front of electrons. A CR pressure precursor is still visible in the conduction and diffusion case with a subtle effect on the Sedov shock front position, which slightly lags behind the solution without conduction and diffusion. 
3.7. Supernova explosion in a 3phase medium: ions, electrons, and cosmic rays
To fully validate the implementation of anisotropic diffusion of CRs and conduction of heat through electrons and their coupling with ions in an astrophysical context, we model the explosion of a supernova in 3D. The explosion is triggered in a homogeneous medium of density n = 1 cm^{3} and with temperatures of electrons and ions at equilibrium T = 10^{4} K, and the background energy density of CRs is onethird of the total energy. We assume that the gas is totally ionised and composed of a cosmic mixture of hydrogen and helium, thus, μ_{I} = 1.22, μ_{E} = 1.13, and μ = 0.59. The magnetic field is uniform along the xaxis of the box and with a magnitude of 1 μG, thus the ratio of internal to magnetic pressure is initially of β = 45. We inject in the eight central cells a total of E_{SN} = 10^{51} erg of energy with 1/3 in CR energy, 1/3 in ion energy, and 1/3 in electron energy.
The temperature of electrons is conducted at the Spitzer rate and couples with the ion temperature. Both ion and electron adiabatic index are equal to γ = 5/3. CRs are diffused with a diffusion coefficient D_{CR} = 3 × 10^{27} cm^{2} s^{1} and have an adiabatic index of γ_{CR} = 4/3 corresponding to ultrarelativistic particles. As currently implemented, our solver could allow for several CR energy components with various diffusion coefficients and CR adiabatic indexes. We use a box size of 500 pc with a coarse grid of level ℓ_{min} = 5 and refine up to ℓ_{max} = 9 wherever the relative gradient of any hydro variable is more than 20%.
Fig. 10 Map cuts in the xz plane for the 3D Sedov explosion including anisotropic electron heat conduction and anisotropic CR diffusion at t = 5 kyr (first two rows) and t = 100 kyr (last two rows). From left to right and top to bottom: the ratio of electron to ion temperature, the ion temperature, the electron temperature, the gas density, the magnetic field amplitude, the gas velocity amplitude, the ratio of CR to internal energy, and the CR energy density. The black dotted circles are the shock front positions of the standard Sedov solution R(t) ≃ 1.1527(E_{SN}/ρ)^{1/5}t^{2/5}. The explosion exhibits a nozzlelike heat/CR precursor in the horizontal direction due to the anisotropic nature of the diffusion (along the horizontal magnetic field). 
For the first test, we did electron and ion coupling without diffusion or conduction, and the second test was with isotropic diffusion and conduction. Figure 8 shows the result of the two simulations after time t = 4 kyr and is compared to the analytic prediction from Sedov (1959). In the absence of conduction and diffusion, the analytical solution is reproduced well, except that the shock interface is smeared by the numerical diffusion. The electrons and ions are not yet at thermal equilibrium within the bubble. The solution is slightly modified with isotropic diffusion of CR energy and conduction of electron temperature. The gas density exhibits two peaks: one ahead of the shock front position and another behind. The ahead overdensity is the product of the conduction of electrons over a few parsec, while the second overdensity lags behind. The total pressure within the bubble is decreased since the CR and electron energies have leaked out of the bubble. Figure 9 shows the result at time t = 100 kyr. Again, the experiment without conduction and diffusion reproduces the analytical prediction well, and ions and electrons are at thermal equilibrium. At this time, the Sedov shock has taken over the electron precursor, which has disappeared. It turns out that electrons and ions are at thermal equilibrium. However, a large amount of the CR energy has escaped the bubble, so that there is less pressure within the bubble, explaining why the shock front slightly lags behind.
The final experiment was run with anisotropic conduction of electrons and anisotropic diffusion of CRs. Figure 10 shows 2D slices in the middle xz plane of several quantities: temperatures, density, CR energy, velocity, and magnetic field amplitude at two different times t = 5 kyr and t = 100 kyr. We see that the electron temperature and CR energy are diffused along the xaxis. The ion temperature catches up with the electron temperature in the heat precursor, but it takes some time for the two species to reach thermal equilibrium. A jetlike nozzle forms in the horizontal direction due to the anisotropic overpressurised region, although it propagates more slowly than the genuine shock front. Nevertheless, at later times, the explosion approaches spherical with still some degree of anisotropy in the bubble expansion. In particular, the shock front moves faster along magnetic field lines (xaxis) than perpendicular to them (note the amplification of the magnetic field due to the compression along the zaxis). This results in an aspherical bubble shape, which is in advance in the xdirection but behind in the zdirection (as for the ydirection).
4. Conclusion
We have presented the implementation of a new solver for the anisotropic diffusion along magnetic field lines of CRs and of heat with twotemperature components, electrons, and ions in the AMR code ramses. This solver is implicit in time and supports adaptivetimestepping on an AMR levelbylevel basis. The implicit anisotropic diffusion scheme uses a conjugate gradient method adapted from Commerçon et al. (2011, 2014). We tested this new implementation against several basic numerical experiments with, for some cases, typical astrophysical conditions, as well as a first simple application to a supernova explosion in a uniform medium. This preliminary work indicates that anisotropic conduction or diffusion, together with temperature coupling, is relevant to several astrophysical problems, as already highlighted in the literature, and, in particular, may affect the propagation of remnants of supernovae.
Future applications of this new part of the ramses code will include the study of active galactic nuclei jets in the intracluster medium, virial shocks in cosmological simulations of galaxy clusters, CRdriven galactic winds, supernovae explosions, and thermal instabilities in the ISM, amongst others.
Acknowledgments
We thank J. Masson for enlightening discussion of the temperature coupling. We thank R. Teyssier, Y. Rasera, I. Parrish, T. Bogdanovich, and P. Hennebelle for useful informal discussions. We warmly thank S. Rouberol for smoothly running the Horizon cluster on which several simulations were run. We thank the anonymous referee for help clarifying several aspects of the paper. This work was supported by the CNRS programmes “Programme National de Cosmologie et Galaxies” (PNCG) and “Physique et Chimie du Milieu Interstellaire” (PCMI). BC gratefully acknowledges support from the French ANR Retour Postdoc programme (ANR11PDOC0031).
References
 Balbus, S. A. 2000, ApJ, 534, 420 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., Bendinelli, A. J., Tilley, D. A., Massari, A. R., & Howk, J. C. 2008a, MNRAS, 386, 642 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., Tilley, D. A., & Howk, J. C. 2008b, MNRAS, 386, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, R., & Krause, M. 2005, Astron. Nachr., 326, 414 [NASA ADS] [CrossRef] [Google Scholar]
 Bogdanović, T., Reynolds, C. S., Balbus, S. A., & Parrish, I. J. 2009, ApJ, 704, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJ, 777, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Brighenti, F., & Mathews, W. G. 2003, ApJ, 587, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Brüggen, M., Ruszkowski, M., & Hallman, E. 2005, ApJ, 630, 740 [NASA ADS] [CrossRef] [Google Scholar]
 Chièze, J.P., Alimi, J.M., & Teyssier, R. 1998, ApJ, 495, 630 [NASA ADS] [CrossRef] [Google Scholar]
 Choi, E., & Stone, J. M. 2012, ApJ, 747, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Courty, S., & Alimi, J. M. 2004, A&A, 416, 875 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Enßlin, T., Pfrommer, C., Miniati, F., & Subramanian, K. 2011, A&A, 527, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabian, A. C. 1994, ARA&A, 32, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaspari, M., & Churazov, E. 2013, A&A, 559, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Günter, S., Yu, Q., Krüger, J., & Lackner, K. 2005, J. Comput. Phys., 209, 354 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasz, M., & Lesch, H. 2003, A&A, 412, 331 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hanasz, M., Kowal, G., OtmianowskaMazur, K., & Lesch, H. 2004, ApJ, 605, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasz, M., Lesch, H., Naab, T., et al. 2013, ApJ, 777, L38 [NASA ADS] [CrossRef] [Google Scholar]
 Jubelgas, M., Springel, V., & Dolag, K. 2004, MNRAS, 351, 423 [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]
 Koyama, H., & Inutsuka, S.I. 2004, ApJ, 602, L25 [NASA ADS] [CrossRef] [Google Scholar]
 McNamara, B. R., & Nulsen, P. E. J. 2007, ARA&A, 45, 117 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2012, MNRAS, 422, 2102 [NASA ADS] [CrossRef] [Google Scholar]
 Miniati, F., Ryu, D., Kang, H., & Jones, T. W. 2001, ApJ, 559, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., & Medvedev, M. V. 2001, ApJ, 562, L129 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parrish, I. J., & Quataert, E. 2008, ApJ, 677, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Parrish, I. J., & Stone, J. M. 2005, ApJ, 633, 334 [NASA ADS] [CrossRef] [Google Scholar]
 Parrish, I. J., Stone, J. M., & Lemaster, N. 2008, ApJ, 688, 905 [NASA ADS] [CrossRef] [Google Scholar]
 Parrish, I. J., Quataert, E., & Sharma, P. 2009, ApJ, 703, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Pfrommer, C., Springel, V., Enßlin, T. A., & Jubelgas, M. 2006, MNRAS, 367, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Piontek, R. A., & Ostriker, E. C. 2004, ApJ, 601, 905 [NASA ADS] [CrossRef] [Google Scholar]
 Rasera, Y., & Chandran, B. 2008, ApJ, 685, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Rosner, R., & Tucker, W. H. 1989, ApJ, 338, 761 [NASA ADS] [CrossRef] [Google Scholar]
 Rudd, D. H., & Nagai, D. 2009, ApJ, 701, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Ruszkowski, M., & Begelman, M. C. 2002, ApJ, 581, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312 [NASA ADS] [CrossRef] [Google Scholar]
 Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York: Academic Press) [Google Scholar]
 Sharma, P., & Hammett, G. W. 2007, J. Comput. Phys., 227, 123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sod, G. A. 1978, J. Comput. Phys., 27, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Sovinec, C. R., Glasser, A. H., Gianakon, T. A., et al. 2004, J. Comput. Phys., 195, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, L. 1956, Physics of Fully Ionized Gases (New York: Interscience Publishers) [Google Scholar]
 Spitzer, Jr., L., & Tomasko, M. G. 1968, ApJ, 152, 971 [NASA ADS] [CrossRef] [Google Scholar]
 Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Ann. Rev. Nucl. Part. Sci., 57, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teyssier, R., Chièze, J.P., & Alimi, J.M. 1998, ApJ, 509, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Tilley, D. A., Balsara, D. S., & Howk, J. C. 2006, MNRAS, 371, 1106 [NASA ADS] [CrossRef] [Google Scholar]
 Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1979, J. Comput. Phys., 32, 101 [Google Scholar]
 Wong, K.W., & Sarazin, C. L. 2009, ApJ, 707, 1141 [NASA ADS] [CrossRef] [Google Scholar]
 Yokoyama, T., & Shibata, K. 2001, ApJ, 549, 1160 [NASA ADS] [CrossRef] [Google Scholar]
 Zel’dovich, Y. B., & Raizer, Y. P. 1967, Physics of shock waves and hightemperature hydrodynamic phenomena (New York: Academic Press) [Google Scholar]
All Figures
Fig. 1 Top: diffusion of a pressure step function as a function of position at t = 9.3 × 10^{4} (black), t = 1.9 × 10^{3} (red), t = 5.6 × 10^{3} (blue) for the simulation (symbols) and the analytic solution (solid lines). The extra levels of refinement are indicated as dashed lines. Bottom: relative errors of the simulation with respect to the analytic solution. 

In the text 
Fig. 2 Electron temperature (red or diamonds) and ion temperature (blue or pluses) as a function of time for the one cell problem. The solution for a time step of τ_{eq,EI,0}/ 50 is represented in solid lines and for the time step of 20τ_{eq,EI,0} with symbols. 

In the text 
Fig. 3 Solution of the Sod (1978) shock tube problem with two temperatures. There is no diffusion and coupling of the temperatures. We plot the velocity (top left), pressure (bottom left), and density (bottom right) at t = 0.245. The squares correspond to the result of the simulation and the red lines to the analytical solution. The light blue and dark blue squares in the pressure plot stand for the first and second species, respectively, while the black squares are for the total pressure. The level of refinement is plotted in the top right panel. 

In the text 
Fig. 4 Stationary cosmological shock experience including conduction and twotemperature coupling. The result of the twotemperature run (blue) is represented with dotted lines for the ion temperature (pressure), in dashed for the electron temperature (pressure), and solid for the average temperature (total pressure) of the gas at t = 10 Gyr. For comparison, the results of the single temperature runs are plotted in red with conduction and in black without conduction. An electron precursor forms in the upstream part of the shock front. 

In the text 
Fig. 5 Temperature map at t = 2 Gyr for the heat conduction in a 2D magnetic loop without AMR (grey levels). The solution with AMR is plotted with red contours for T = 5 × 10^{6} K and T = 10^{7} K (white contours without AMR). 

In the text 
Fig. 6 Temperature map at t = 2 Gyr for the heat conduction in a 2D magnetic loop with AMR and different fractions of perpendicular conductivity coefficients κ_{iso} = 0.1κ_{Sp} (red contours), κ_{iso} = 0.01κ_{Sp} (grey levels and white contours), and κ_{iso} = 0.001κ_{Sp} (blue contours). Two contours are plotted for each simulation result corresponding to temperatures T = 5 × 10^{6} K and T = 10^{7} K (except for κ_{iso} = 0.1κ_{Sp} where only T = 5 × 10^{6} K is reached). 

In the text 
Fig. 7 Results of the perpendicular numerical diffusion coefficient κ_{num} as a function of resolution Δx for the Sovinec et al. (2004) test problem in a solid line with square symbols. The numerical diffusion scales with Δx^{2} (dotted line). 

In the text 
Fig. 8 Spherically averaged values of density (left), pressure (middle), and velocity (right) as a function of radius at t = 4 kyr for the 3D Sedov explosion without diffusion (top), and with isotropic diffusion (bottom). The solid lines are the result of the simulation, and the dashed line is the analytic prediction for the standard Sedov solution. The pressure is decomposed into total pressure (black), ion pressure (blue), electron pressure (red), and CR pressure (magenta). Without diffusion, the analytic solution is reproduced well. With diffusion, there are both a CR precursor and an electron precursor. 

In the text 
Fig. 9 Same as Fig. 8 at t = 100 kyr. The shock front takes over the conduction front of electrons. A CR pressure precursor is still visible in the conduction and diffusion case with a subtle effect on the Sedov shock front position, which slightly lags behind the solution without conduction and diffusion. 

In the text 
Fig. 10 Map cuts in the xz plane for the 3D Sedov explosion including anisotropic electron heat conduction and anisotropic CR diffusion at t = 5 kyr (first two rows) and t = 100 kyr (last two rows). From left to right and top to bottom: the ratio of electron to ion temperature, the ion temperature, the electron temperature, the gas density, the magnetic field amplitude, the gas velocity amplitude, the ratio of CR to internal energy, and the CR energy density. The black dotted circles are the shock front positions of the standard Sedov solution R(t) ≃ 1.1527(E_{SN}/ρ)^{1/5}t^{2/5}. The explosion exhibits a nozzlelike heat/CR precursor in the horizontal direction due to the anisotropic nature of the diffusion (along the horizontal magnetic field). 

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.