© 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 non-diffusion. In some situations, it leads to a non-thermal 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 intra-cluster 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 shock-driven 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 so-called magneto-thermal instability (Balbus 2000; Parrish et al. 2008) and the heating-buoyancy 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 intra-cluster 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 self-regulation of the ISM dynamics and, thus, of star formation. For instance, the injection of CRs within remnants of supernovae lead to stronger galactic-scale 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²/ (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 multi-scale 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 semi-implicit 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 time-stepping. This new numerical implementation is augmented by modelling a multi-temperature 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 set-up

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² + e_th + e_cr + B²/ 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²/ 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³, 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 non-thermal 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 (star-forming 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 right-hand 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_CRb(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⁸ cm and the mean free path λ_mfp = 10²¹ 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²/ (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_ik_BT/ (γ − 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_ik_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ⁿ and coefficients, x the vector of temperature T^{n + 1}, and c the vector of energy densities can be generalised in multi-dimensions 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 time-stepping on a level-by-level 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 time-stepping 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 non-uniform interfaces: the fine-to-coarse interface (interface between a cell at level ℓ and a cell at ℓ − 1) and the coarse-to-fine 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 fine-to-coarse 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 fine-to-coarse boundary. For the coarse-to-fine 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 coarse-to-fine boundary. In any case, Eq. (13) remains correct since the level-by-level 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 right-hand 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 non-energy conservation due to Dirichlet boundary at level interfaces and the non-monotonicity 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 non-monotonicity 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 non-symmetric. 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 bi-conjugate 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 = ρ/ (μ_Im_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⁵ 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_Em_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 non-linear 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₁(X^k + ΔX) = 0 and G₂(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 time-stepping 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₁ = 0.8 between x = [ 0.25,0.75 ] and P₀ = 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 × Δ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.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¹⁰ K and ions T_I,0 = 10⁸ 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 two-temperature coupling. The result of the two-temperature 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 two-component 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 x-direction, and ion and electron temperatures at equilibrium T_I,R = T_E,R = T_R = 4 × 10⁷ K. Using the Rankine-Hugoniot 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⁷ 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 two-temperature 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 single-temperature experiment. Ion temperature in the precursor lies below-with a sharp rise at the shock interface-because 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 two-temperature 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² 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 circular-loop 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 × 106 K and T = 107 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 × 106 K and T = 107 K (except for κiso = 0.1κSp where only T = 5 × 106 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² to 256² up to the steady-state 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² 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 Δx2 (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 3-phase 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⁴ K, and the background energy density of CRs is one-third 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 x-axis 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⁵¹ 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²⁷ cm² s^-1 and have an adiabatic index of γ_CR = 4/3 corresponding to ultra-relativistic 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(ESN/ρ)1/5t2/5. The explosion exhibits a nozzle-like 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 over-density is the product of the conduction of electrons over a few parsec, while the second over-density 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 x-axis. 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 jet-like nozzle forms in the horizontal direction due to the anisotropic over-pressurised 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 (x-axis) than perpendicular to them (note the amplification of the magnetic field due to the compression along the z-axis). This results in an aspherical bubble shape, which is in advance in the x-direction but behind in the z-direction (as for the y-direction).

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 two-temperature components, electrons, and ions in the AMR code ramses. This solver is implicit in time and supports adaptive-time-stepping on an AMR level-by-level 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 intra-cluster medium, virial shocks in cosmological simulations of galaxy clusters, CR-driven 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 (ANR-11-PDOC-0031).

References

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 two-temperature coupling. The result of the two-temperature 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 × 106 K and T = 107 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 × 106 K and T = 107 K (except for κiso = 0.1κSp where only T = 5 × 106 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 Δx2 (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(ESN/ρ)1/5t2/5. The explosion exhibits a nozzle-like heat/CR precursor in the horizontal direction due to the anisotropic nature of the diffusion (along the horizontal magnetic field).

In the text