© A. Navarro et al. 2022

1 Introduction

Thermal conductivity is an essential physical process by which pressure and temperature perturbations are spread out through a plasma. It corresponds to a macroscopic manifestation of a transport process (energy transfer by particle collisions) and its strength depends on the temperature and temperature gradient. It plays an important role in diverse range of conditions from those found in thermoelectric materials (Yu et al. 2022; Liu et al. 2022) to supernova remnants (Balsara et al. 2008). In the solar atmosphere, thermal conduction is one of the dominant processes that shapes the dynamics of the hot plasma (Ye et al. 2020). For instance, it contributes to the energy redistributions in plasma eruptions (Bradshaw et al. 2012; Liu et al. 2009; Ye et al. 2020; Navarro et al. 2021). Heat conduction is a key aspect in the numerical analysis of the coronal heating problem (Bingert & Peter 2011, 2013; Gudiksen et al. 2011; Bourdin et al. 2013; Chen et al. 2014).

The ratio of the collisional frequency to the cyclotron frequency, which depends on the magnetic field, determines the direction of the heat flux with respect to the magnetic field direction. In the corona, where the magnetic field is strong and the collision frequency is low, the thermal conduction is aligned with the magnetic field and is commonly modeled by Spitzer’s expression for fully ionized plasmas. Likewise, Braginskii (1965) developed a more complete description of the heat flux which depends on the ratio between the cyclotron frequency and the collisional frequency of the plasma. In his model, the heat flux is decomposed into three parts, a parallel component, a perpendicular component and a transverse component, each of them for ions and electrons. Therefore it allows a smooth transition between field-aligned conductivity and isotropic conductivity in regions with a low or null magnetic field or regions with high collisional frequency. Furthermore, the parallel component of the electron heat flux has roughly the same value as the Spitzer’s expression. Accordingly, Braginskii (1965) model is appropriate for the description of the complete solar atmosphere.

Modeling the thermal conductivity numerically is a challenging endeavor since explicit solvers for a parabolic problem lead to very small time-steps according to the Courant–Friedrichs–Lewy (CFL) condition in comparison to a hyperbolic problem defined by ideal magnetohydrodynamics (MHD) equations. Different schemes have been developed to overcome this restriction, for instance: Super Time-Stepping methods (Meyer et al. 2012) (used in the MPI-AMRVAC code, see Xia et al. 2018), implicit–explicit schemes (Balsara et al. 2008), semi-implicit directionally-split methods (Sharma & Hammett 2011; Ye et al. 2020), symmetric and asymmetric discretization schemes (Günter et al. 2005; Parrish & Stone 2005). With the same purpose, Rempel (2017) and Warnecke & Bingert (2020) introduced an efficient way to deal with the numerical time step constraints by solving the hyperbolic equation for the heat flux for the MURAM and PENCIL codes respectively. In the current work, we implemented into the code MANCHA3D (Felipe et al. 2010; González-Morales et al. 2018; Khomenko et al. 2018) a combination of both schemes, the hyperbolic equation to evolve the parallel to the magnetic field component of the heat flux, and the parabolic equation for the perpendicular and transverse directions.

The paper is organized as follows. In Sect. 3, we describe two different schemes for the implementation of thermal conductivity in the code. From Sects. 4.1 to 4.4 we validate the numerical schemes with different tests in one, two and three dimensions. The general Braginskii heat flux is presented in Sect. 5, we describe its implementation in the code and compare the thermal conductivity in parallel, perpendicular and transverse directions, from ions and electrons in a quiet Sun stratified atmosphere. Finally, in Sect. 5 we propose a two-dimensional test with a realistic solar stratified atmosphere and in Sect. 6 we present our conclusions.

2 Brief Description of the MANCHA3D Code

MANCHA3D is a three-dimensional numerical code designed to simulate complex dynamics of the solar atmosphere plasma. It solves the non-linear MHD equations taking into account non ideal effects as the ambipolar diffusion, Hall effect, the battery effect and radiative losses, although, in the simulations that are described in this work, such effects are not taken into account. A gravitationally stratified plasma in the ideal regime can be represented by the following set of equations $$\frac{\partial \rho }{\partial t}+\nabla .\left(\rho \text{v}\right)=0,$$(1) $$\frac{\partial \text{B}}{\partial t}=\nabla \times \left(\text{v}\times \text{B}\right),$$(2) $$\frac{\partial \left(\rho \text{v}\right)}{\partial t}+\nabla .\left[\rho \text{vv+}\left(p+\frac{{\text{B}}^2}{2{\mu }_0}\right)\text{I-}\frac{\text{BB}}{{\mu }_0}\right]=\rho \text{g,}$$(3) $$\frac{\partial e}{\partial t}+\nabla \cdot \left[\left(e+p+\frac{{\text{B}}^2}{2{\mu }_0}\right)\text{v-}\frac{\text{B}}{{\mu }_0}\left(\text{B}\cdot \text{v}\right)\right]=\rho v\cdot \text{g},$$(4)

where ρ is the density, v is the velocity, p is the gas pressure, B is the magnetic field, g is the gravitational acceleration, e is the total energy per unit volume and I is the identity tensor The dot “.” represents the scalar product of vectors, while the notation “BB” stands for the tensor product.

In the code, the equations are discretized with a centered finite differences scheme of sixth order and integrated in time with an explicit Runge–Kutta method. For numerical stability, the code uses hyperdiffusivity and filtering. A more complete description of the code can be found in the works of Khomenko & Collados (2006, 2008); Felipe et al. (2010); González-Morales et al. (2018); Khomenko et al. (2018).

3 Numerical Schemes for Thermal Conductivity

The effects of thermal conductivity can be represented by adding the divergence of the heat flux vector q into the energy equation (4). Accordingly, the new equation for the energy reads as $$\frac{\partial e}{\partial t}+\nabla \cdot \left[\left(e+p+\frac{{\text{B}}^2}{2{\mu }_0}\right)\text{v-}\frac{\text{B}}{{\mu }_0}\left(\text{B}\cdot \text{v}\right)\right]=\rho v\cdot \text{g-}\nabla \cdot \text{q.}$$(5)

In the strongly magnetized limit, for the fully ionized plasma, the heat flux is aligned with the magnetic field $$\text{q}=-{\kappa }_{\parallel }{\nabla }_{\parallel }T=-{\kappa }_0{T}^{5/2}\widehat{\text{b}}\left(\widehat{\text{b}}\cdot \nabla \right)T,$$(6)

where $\widehat{b}$ is the unit vector along the magnetic field and the term ($\widehat{b}\cdot \nabla T$) represents the temperature gradient along the magnetic field, κ_|| is the parallel conductivity and к₀ is a constant that depends on the properties of the plasma. This equation is referred as Spitzer’s thermal conduction vector (Spitzer 1956).

The integration time step of a fully explicit finite difference method is limited by the CFL condition. For the heat conduction problem, the value of the limiting time step can be obtained by comparing the time-derivative and the thermal conduction terms in the energy equation, assuming no variations in density and the ideal gas law $$\frac{\partial e}{\partial t}\approx \frac{\left(\gamma -1\right){\kappa }_0{T}^{7/2}}{p}{\nabla }^{2}e$$(7)

therefore, the time step is given by $$\text{d}t\le {\text{cfl}}_{\text{TC}}\frac{\min {\left(\nabla x,\nabla y,\nabla z\right)}^2}{\max \left(\left(\gamma -1\right){\kappa }_0{T}^{7/2}/p\right)},$$(8)

where ∆x, ∆y and ∆z are the grid sizes and cfl_TC is a stability constant that takes values of 0.5, 0.5 and 1/3 for 1D, 2D and 3D cases, respectively.

In some cases, the thermal conduction time step can become much smaller than the time step given by the MHD CFL condition. To overcome this difficulty, Rempel (2017) and Warnecke & Bingert (2020) proposed to solve the hyperbolic equation for the heat flux instead of the classical parabolic equation, i.e., $$\frac{\partial q}{\partial t}=\left(-f_{\text{sat}}{\kappa }_0{T}^{5/2}\left(\widehat{\text{b}}\cdot \nabla \right)T-q\right),$$(9) $$\frac{\partial q}{\partial t}=[\dots ]-\nabla \cdot \left(q\widehat{\text{b}}\right).$$(10)

The factor f_sat sets the saturation of the conductive heat flux, preventing it from taking highest values, and τ is the relaxation time. Following Fisher et al. (1985) and Meyer et al. (2012) the saturation factor can be written as $$f_{\text{sat}}=\left(1+\frac{\left|{\kappa }_0{T}^{5/2}\left(\text{b}\cdot \nabla \right)T\right|}{1.5\rho c_s^3}\right),$$(11)

where $c_s=\sqrt{\gamma p/\rho }$ denotes the speed of sound and γ is the adiabatic index.

The hyperbolic heat-conduction equation corresponds to the modification of the classical Fourier heat flux by the addition of an additional extra thermal inertia term. This modification was proposed independently by Vernotte (1958) and Cattaneo (1958) in order to suppress the paradox of the infinite speed of propagation of the heat conduction in the Fourier theory, and implying the existence of thermal waves. This approach has been widely studied, both theoretically and experimentally, see Abdel-Hamid (1999) and references therein.

In Eq. (9), the relaxation time depends on the mechanism of heat transport, and represents the time lag needed to establish steady-state heat conduction in an element of volume when a temperature gradient is suddenly applied to that element (Abdel-Hamid 1999). Different definitions for the relaxation time are proposed by Rempel (2017) and Warnecke & Bingert (2020), however they do not allow stable running of all our samples, therefore we use a fixed value for the relaxation time, τ = 4dt in Eq. (9).

Fig. 1 Left panel shows the evolution of temperature in the one-dimensional thermal conductivity test computed with the scheme #1 (explicit evolution of the parabolic term) and with the scheme #2 (hyperbolic equation for heat flux) with three different time steps. Right panel shows a comparison between simulations using scheme #2, one including the saturation factor of Eq. (11) and another one without it.

4 Tests for Validation of the Numerical Scheme

4.1 One-Dimensional Problem

For the first test, we have considered the following one-dimensional heat conduction problem proposed by Rempel (2017). The initial non-equilibrium temperature profile is given by $$T_0\left(x\right)=0.1+0.9x^5,$$(12)

in non-dimensional units. In these units, we use density ρ = 1.0, magnetic field B_x = 1.0, B_y = B_z = 0.0, and thermal conductivity coefficient к₀ = 1.0. The domain extends over 0 ≤ x ≤ 1 and is covered by 250 points. The boundary conditions set all the quantities to be fixed in time. The asymptotic stationary solution is given by $$T(x)={\left[{0.1}^{3.5}+\left(1-{0.1}^{3.5}\right)x\right]}^{\text{2/7}}.$$(13)

To analyze this problem, we perform four different simulations. On the one hand, we have performed a run with a standard parabolic scheme (we will refer to this scheme as scheme #1 over the paper), where we solve the Eqs. (1), (2), (3) and (5) using the field-aligned heat flux $$\text{q}=-f_{\text{sat}}{\kappa }_0{T}^{5/2}\widehat{\text{b}}\left(\widehat{\text{b}}\cdot \nabla \right)T\text{.}$$(14)

Here we have included the saturation factor f_sat given by Eq. (11) to the Spitzer heat flux (Eq. (6)). This is needed for consistency in comparison of the simulations with different schemes. On the other hand, we run simulations with the hyperbolic scheme using three different time steps, dt = dt_TC, dt = 50dt_TC and dt = 100dt_TC. We will refer to the latter one as scheme #2.

Figure 1 shows different stages of the time evolution, t = 0.25, t = 0.5, t = 0.75 and t = 1.0 for all the simulations. The solid black lines represent the initial state (t = 0) and the asymptotic solution (13). In all the cases, the reference solution is recovered around t = 1 and the intermediate steps of the evolution of the temperature are quite similar between different simulation cases. It should be emphasized that the temperature evolves with the same speed for both parabolic and hyperbolic schemes, because both schemes use the same saturation factor. It decreases the magnitude of the heat flux enhancing stability of the code. Without the saturation factor, the temperature evolves faster with respect to the simulation time, but it requires a smaller time step for a stable run. This faster evolution is clearly shown at left panel of Fig. 1. The saturation factor allows running simulation with much larger time step, however if the time step is too high, the simulation becomes unstable, as shown in the top left panel of Fig. 6.

To estimate the accuracy of our results, we compare them to the reference solution (13). Table 1 presents the errors L¹, L² and L^∞, computed according to the following definitions, $$L_N^1=\frac{1}{N}{\displaystyle \sum _{i=1,N}\left|T_i-T_i^{ref}\right|}\text{,}$$(15) $$L_N^2=\frac{1}{N}\sqrt{{{\displaystyle \sum _{i=1,N}\left|T_i-T_i^{ref}\right|}}^2}\text{,}$$(16) $$L_N^{\infty }=\max \left(\left|T_i-T_i^{\text{ref}}\right|\right),$$(17)

where N is the total number of points in the domain. We computed the values of L¹, L² and L^∞ at t = 1. The Table gives values for the simulations with the different schemes and, for the scheme #2, the values for large time steps are presented. The errors slightly increase between scheme #1 and #2, and also increase for larger time steps. However, they stay within the same order of magnitude. It suggests that we can speed up the simulations using the hyperbolic treatment and still get rather accurate modeling of thermal conduction problems.

Fig. 2 Two-dimensional static ring problem. Top panels show colormaps of the temperature and magnetic field lines at the initial and final states. The bottom left panel shows a comparison of the temperature at y = 0 between the different simulations at the final stage. Bottom right panel gives a graphical representation of the computational times.

4.2 Two-Dimensional Static Ring

The ring diffusion test proposed by Parrish & Stone (2005) is a standard test to evaluate the monotonicity properties of the anisotropic conduction. The setup consists of hot patch diffusing in a fixed circular-shaped magnetic field. This test is important to check the numerical schemes since the magnetic field lines align in all possible angles with respect to the Cartesian grid. The initial temperature is given by $$T=\{{}_{10\text{\hspace{0.17em}\hspace{0.17em}otherwish}}^{12\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}}0.5\text{\hspace{0.17em}\hspace{0.17em}<\hspace{0.17em}}r<0.7\text{\hspace{0.17em}and }\frac{11}{12}\pi <\theta <\frac{13}{12}\pi },$$(18)

in non-dimensional units, where $r=\sqrt{x^2+y^2}$ and tan θ = y/x. The density is constant and set to ρ = 1, and the circular magnetic field is given by $$B_x={10}^{-5}\cos \left(\theta +\pi /2\right)/r,$$(19) $$B_y={10}^{-5}\sin \left(\theta +\pi /2\right)/r.$$(20)

Initially, the system is at rest with υ_x = υ_y = υ_z = 0. To study this problem we use a numerical box of the size [−1, 1] × [−1, 1], covered by 200 × 200 points, we evolve only energy equation with a parallel thermal conductivity to κ_‖ = 0.01, and outflow boundary conditions. We run this problem using the scheme #1 with dt = dt_TC and, to evaluate the scheme #2, we performed simulations using the time steps dt = dt_TC and dt = 10dt_TC, dt = 20dt_TC, and dt = 30dt_TC. Top panels of Fig. 2 show the colormaps of the temperature and the magnetic field lines at the initial time moment (t = 0) on the left and at the final stage (t = 400) on the right.

The bottom left plot of Fig. 2 compares the temperature at y = 0 for the different runs, and shows nearly identical results. The panel at the bottom right of Fig. 2 displays the computational times for the different numerical schemes and it clearly demonstrates that the hyperbolic treatment is much more efficient then the parabolic one. When we set the time factor dt = 10dt_TC the simulation is 10.2 times faster compared with the scheme #1. Accordingly, when dt = 20dt_TC the efficiency increases 20.3 times and for dt = 30dt_TC, 30.3 times. An unstable example for this test is given at the top right panel of Fig. 6, which shows a colormap for a simulation computed with larger timestep dt = 33dt_TC. One can observe that an evident ringing is produced in this case, since the time step is beyond the stable limit.

According to the analytical solution, at advanced simulation times, the interior of the ring should have a uniform temperature of 10.1667 and outside the ring the temperature should be equal to 10. In Table 2 we provide the accuracy measures computed for the final snapshot of our simulations (t = 400). We calculate the maximum and the minimum temperatures in the domain, and the error L¹ for the maximum temperature by comparing it with the reference value. The magnitudes of the error are within the same order, which indicates that the results obtained with the different schemes and time steps are similar with the results produced by other codes, see Sharma & Hammett (2007); Xia et al. (2018) and Meyer et al. (2012).

The perpendicular numerical diffusion may cause non-negligible cross-field conduction. To estimate its value, we perform simulations with different explicit perpendicular conductivity K_⊥ and calculate the difference in the maximum temperature between the final and initial stages, ∆T. Similar to the calculations in Meier et al. (2010). Figure 3 shows ∆T for different values of κ_⊥. The numerical data are fit with a hyperbolic tangent function (solid line). We assign the perpendicular numerical conductivity to the value of κ_⊥ at location where the first derivative of the hyperbolic tangent function is 100 times below the maximum value of this derivative. The dashed line in Fig. 3 corresponds to the value of ∆T in a simulation where explicit κ_⊥. was set to zero (∆T_num = ∆T (κ_⊥ = 0)). It is shown in the plot to compare visually the precision of the calculation. The ratio between the numerical perpendicular conductivity and the physical parallel conductivity (κ_⊥,num/k_‖) for the different simulations is listed in the last column of Table 2. It is of the order ~10⁻⁵, which implies in one hand that we could not study the effects of anisotropic conduction in laboratory plasmas since in that case the ratio should be ~10⁻⁹, otherwise perpendicular numerical diffusion will swamp the true perpendicular diffusion. On the other hand, since the ratio it is smaller than ~10⁻³ it should be adequate to study qualitatively the effects of anisotropic conduction on dilute astrophysical plasmas (Sharma & Hammett 2007).

Fig. 3 Graphical representation of the method used to calculate the numerical perpendicular conductivity. The dots represent the final maximum temperature minus initial maximum temperature (∆T) in the simulations with different values of k⊥ shown at the horizontal axis. The dashed line is the difference of temperature (∆Tnum) in the simulation with κ⊥ = 0. We fit the data to a hyperbolic tangent function (solid line). The numerical perpendicular conduction coefficient к⊥,num is estimated by locating the position where the first derivative of the fit is 100 times smaller than its maximum value.

4.3 Two-Dimensional Hot Plate Test

A two-dimensional test designed to evaluate the anisotropic heat conduction along oblique magnetic field lines was proposed by Jiang et al. (2012) and reproduced by Navarro et al. (2017). The initial data consists of a small hot rectangular “plate” at the bottom of the domain. The initial pressure is set to a homogeneous value of p = 0.1, the system is at rest, υ_x = υ_y = υ_z = 0, and the magnetic field is set constant and inclined by a 45° angle, i.e., B_x = B_y = 1, B_z = 0. The initial density is given by $$\rho =\{{}_{0.1\text{\hspace{0.17em}elsewhere,}}^{0.01\text{\hspace{0.17em}for\hspace{0.17em}}\left|x\right|<\text{\hspace{0.17em}and\hspace{0.17em}}y=-0.5,}$$(21)

in a domain extended over [−0.5,0.5] × [−0.5,0.5] with 200 × 200 points. The boundary conditions are periodic in x, fixed to the initial values at the bottom surface at z = −0.5, and outflow at the top. As before, we perform simulations with both schemes. For the hyperbolic scheme, we use four times steps of dt = 10dt_tc, dt = 100dt_tc, dt = 200dt_tc and dt = 500dt_te. Figure 4 shows the colormaps of the temperature and the magnetic field lines at t = 1 for the different schemes. In all the cases the temperature is advected along the magnetic field lines at the same pace and with negligible transverse conduction. In order to compare the differences in more detail, the left-bottom plot shows a 1D cut of the temperature at y = −0.15. We find that all the numerical solutions are quite similar. The right-bottom plot compares the cpu times to the solution. It again confirms that the the hyperbolic treatment produces very similar results with a great advantage in computational time. For instance, for the run where the time step is dt = 10dt_tc times larger, the solution is 8.23 times faster. Similarly, the efficiency ratio for the factors 100, 200, and 500 are 76.57,156.7 and 405.9, respectively. In these simulations the advection time step is 1550 times larger then the conduction one. An example of an unstable simulation for this test is given in the bottom left panel of Fig. 6. It displays a colormap in a simulation using a time step beyond the stability limit (dt = 500dt_tc), a clear ringing is present from an early stage.

Fig. 4 Left top to right middle plot: colormaps of the temperature, together with magnetic field lines, at t = 1.0, obtained in the different simulations, using scheme #1 and scheme #2 with the time steps dt = 10dttc, dt = 100dttc, dt = 200dttc and dt = 400dttc. The bottom left plot shows a 1D cut of the temperature at z = −0.27, and the right bottom panel shows the simulation time for each run.

4.4 Three-Dimensional Static Ring

A 3D generalization of the static ring problem was proposed by Xia et al. (2018). Starting from the set up described in Sect. 4.2, we extend the hot patch over a width of 0.4 in z direction. The whole system, including the magnetic field, is rotated by π/4 around the x-axis and by π/4 around the z-axis. This rotation causes a misalignment between the field and the coordinate axes, so that the thermal flux evolves in the three dimensions. We run four different simulations in the numerical box of [−1, 1] × [−1, 1] × [−1, 1] with 200 × 200 × 200 points. The first one is with the scheme #1, and remaining three are with scheme #2 with time steps dt = 10dt_TC, dt = 20dt_TC, dt = 30dt_TC and dt = 40dt_TC.

Top panels of Fig. 5 present a contour of the temperature together with magnetic field lines at the initial and last stage of the simulation with scheme #1. Bottom left panel allows us to compare the results between the different runs showing the temperature at x = −0.5 and y = 0.21 at the final stage t = 400. One can observe that there is a slight difference in the temperature profiles. In particular, the runs with the scheme #2 with larger time step (dt = 20dt_tc and dt = 30dt_tc) exhibit less perpendicular dissipation since the maximum values of the temperature in those cases are bigger. Additionally, the bottom right panel of Fig. 5 displays the computational times for the simulations, showing the big advantage in computational times of the hyperbolic treatment. In these simulations the advection time step is 150 times larger than the conduction one. For the simulation where the time step is dt = 10dt_tc, the efficiency ratio is 4.9. Similarly, the efficiency ratio for the factors 20 and 30, are of 9.45, and 14.6 respectively. An example of a unstable simulation using scheme #2 is presented in bottom right panel of Fig. 6, which shows a ringing produced by the large time step, in this case, dt = 45dt_tc.

Table 3 presents the maximum temperature reached at the final stage and provides the L¹ error computed by comparing this maximum temperature with the theoretical value T = 10.1667. The values of the error are quite similar in the different simulations and with the results obtained in the 2D case, proving that the code resolves with sufficient accuracy three-dimensional thermal conduction problems with complex magnetic field configurations.

Fig. 5 Results for the three-dimensional static ring problem. The top panels correspond to the contour of temperature at the initial and final stages, together with the magnetic field lines, obtained with the simulation with scheme #1. Bottom left panel shows a comparison of one snapshot of the temperature at x = −0.5 and y = 0.21 at t = 10 for different simulations using scheme #2 and different time steps. In the bottom right panel, a graphical representation of the computational times of the different runs is given.

5 Thermal Conduction in the Solar Atmosphere

In a magnetized plasma, the efficiency of the heat conduction depends on the magnetic field direction. It is therefore convenient to decompose the conductivity vector using the projections into the parallel and perpendicular to the field directions, $$\text{q\hspace{0.17em}}=-{\kappa }_{\parallel }{\nabla }_{\parallel }T-{\kappa }_{\perp }T+{\kappa }_{\times }\text{b}\times {\nabla }_{\perp }T,$$(22)

where ${\nabla }_{\parallel }=\widehat{b}\left(\widehat{b}\cdot \nabla \right)$ gives the parallel projection to the field, ∇_⊥= ∇−∇_‖ gives the projection in the perpendicular direction, and the last term is the projection in the transverse direction (second perpendicular direction to the magnetic field). For example, similar decomposition was recently used in Hunana et al. (2022). Braginskii (1965) deduced the following general expressions for the conductivity coefficients used in the expressions of the heat conduction vector for electrons $${\kappa }_{\parallel }^e=3.1616\frac{k_B{p}_e}{v_{ei}{m}_e},$$(23) $${\kappa }_{\perp }^e=\frac{k_B{p}_e}{v_{ei}{m}_e}\frac{4.66x_e^2+11.92}{x_e^4+14.79x_e^2+3.77},$$(24) $${\kappa }_{\times }^e=\frac{k_B{p}_e}{v_{ei}{m}_e}{x}_e\frac{\frac{5}{2}{x}_e^2+21.67}{x_e^4+14.79x_e^2+3.77},$$(25)

and for ions, $${\kappa }_{\parallel }^i=3.906\frac{k_{\text{B}}pi}{v_{ii}{m}_i},$$(26) $${\kappa }_{\perp }^i=\frac{k_B{p}_i}{v_{ii}{m}_i}\frac{2x_i^2+2.645}{x_i^4+2.70x_i^2+0.677},$$(27) $${\kappa }_{\times }^i=\frac{k_B{p}_i}{v_{ii}{m}_i}{x}_i\frac{\frac{5}{2}{x}_i^2+4.65}{x_i^4+2.70x_i^2+0.677},$$(28)

where the subindices and upper indices e and i refer to electrons and ions. The pressure and mass of each species are p_e, p_i, m_e and m_i. The collisional frequency of ion-ion collisions is v_ii and of electron-ion collissions is v_ei. The quantities x_i and x_e represent the ratio between cyclotron frequency (Ω) and collision frequency (v) for electrons and ions, respectively $$x_{\text{e}}=\frac{{\Omega }_e}{v_{ei}};x_i=\frac{{\Omega }_i}{v_{ii}},$$(29)

where $${\Omega }_{\text{e}}=\frac{\left|e\right|B}{m_{\text{e}}};{\nu }_{ei}=3.7{\times }^{-6}\frac{\ln \left(\Lambda \right)n_{\text{e}}}{T_{\text{e}}^{3/2}};p_{\text{e}}=n_{\text{e}}{k}_B{T}_{\text{e}},$$(30) $${\Omega }_i=-\frac{\left|e\right|Z_{i}B}{m_i};v_{ii}=6\times {10}^{-8}\frac{1\text{n}\left(\text{Λ}\right)n_i{z}_i^4}{T_i^{3/2}};p_i=n_i{k}_B{T}_i,$$(31)

e is the electron charge, Z_i is the charge of the ion, n_e and n_i the electron and ion number density. All the definitions in SI units. The Coulomb logarithm could be approximated by $$\ln \left(\Lambda \right)=23.4-1.15{\log }_{10}\left(n_{\text{e}}\right)+3.45{\log }_{10}{T}_{\text{e}},$$(32)

for T_e < 50 eV, with n_e given in cm⁻³ and T_e in electronvolts. It follows from these expressions that the case of null magnetic field the heat flux becomes isotropic (q = −k_‖∇T) since x_e = x_i = 0 and $k_{\perp }^{\text{i,e}}={\kappa }_{\parallel }^{\text{i,e}}$.

The parallel conductivity coefficients given in Eqs. (23) and (26) have a dependence on temperature as ~T^5/2 (this follows from the dependence on temperature of the pressure and the collision frequencies, see Eqs. (30), (31)). The expression for the electron heat flux is the same as given by Spitzer, after exchanging the constant 3.203 for 3.1616 in Eq. (23). This dependence on temperature is specially problematic due to the high temperatures of the corona, when solving the thermal conductivity with an explicit integration method. We propose to combine the schemes #1 and #2 by solving the energy equation, Eq. (5), where the heat flux is given by $$\text{q=}{}_{-q_{\parallel }}\widehat{\text{b}}{}_{-{\kappa }_{\perp }}{\nabla }_{\perp }T{+}_{\kappa \times }\widehat{\text{b}}\times {\nabla }_{\perp }T,$$(33)

with the hyperbolic equation for the evolution of the parallel heat flux component q_‖, $$\frac{\partial q_{\parallel }}{\partial t}=\frac{\text{l}}{\tau }\left(f_{sat}{\kappa }_{\parallel }\left(\widehat{\text{b}}\cdot \nabla \right)T+q\parallel \right).$$(34)

In the latter expression, k_‖ is the conductivity coefficient in parallel direction to the magnetic field $${\kappa }_{\parallel }={\kappa }_{{}_{\parallel }}^e+{\kappa }_{{}_{\parallel }}^i,$$(35)

together which includes the contributions from electrons ${\kappa }_{\parallel }^e$ (Eq. (23)) and ions ${\kappa }_{\parallel }^i$ (Eq. (26)), assuming T_e = T_i = T, n_e =n_i and Z_i = 1. Accordingly, k_⊥ is the conductivity coefficient in the perpendicular direction to the magnetic field $${\kappa }_{\perp }={\kappa }_{\perp }^e+{\kappa }_{\perp }^i,$$(36)

which includes the contributions from electrons ${\kappa }_{\perp }^e$ (Eq. (24)) and ions ${\kappa }_{\perp }^i$ (Eq. (27)); κ_× is the transversal conductivity coefficient, perpendicular to both the magnetic field and the temperature gradient $${\kappa }_{\times }={\kappa }_{\times }^{\text{e}}+{\kappa }_{\times }^i,$$(37)

and includes the contributions from electrons ${\kappa }_{\times }^{\text{e}}$ (Eq. (25)) and ions ${\kappa }_{\times }^i$ (Eq. (28)).

Since q_‖ is an independent variable, its initial value should be calculated like $$q_{\parallel }\left(t=0\right)=f_{\text{sat}}{\kappa }_{\parallel }\left(\widehat{\text{b}}\cdot \nabla \right)T.$$(38)

The model above takes into account six different conductivity coefficients. It is of interest to compare the magnitude of the terms as a function of height and magnetisation of the solar stratified atmosphere. For that, the left panel of Fig. 7 presents those conductivities as functions of height, considering a stratified atmosphere with the temperature model of Vernazza et al. (1981) extended to the corona. We adopted the following approximate dependence for the magnetic field, appropriate for a quiet solar region, $$B_z=\left(100\text{G}\right)\exp \left({}_{-Z}/\left(6\text{Mm}\right)\right).$$(39)

The results will of course vary if a different strength and stratification for the magnetic field is adopted. However, the conclusions below will remain qualitatively similar, since the dependence of the conductivity coefficients on the thermodynamic parameters is much stronger.

According to the left panel of Fig. 7, the parallel conductivity $k_{\parallel }=k_{\parallel }^{\text{e}}+k_{\parallel }^i$ is the most important one at all heights. Although in the photosphere it is equal in magnitude to the perpendicular one $k_{\perp }=k_{\perp }^{\text{e}}+k_{\perp }^i$, which has its largest values in the lower atmosphere. The biggest contribution of the transverse conductivity $k_{\times }=k_{\times }^{\text{e}}+k_{\times }^i$ is before the transition region, after, it drops drastically and increases smoothly according to height. Now, to determine specifically the contributions of each conductivity due to electrons or ions, the right panel of Fig. 7 provides the ratios ${\kappa }_{\parallel }^i/{\kappa }_{\parallel }^e,{\kappa }_{\perp }^{\text{e}}/{\kappa }_{\parallel }^{\text{e}},{\kappa }_{\perp }^i/{\kappa }_{\parallel }^{\text{e}},{\kappa }_{\times }^i/{\kappa }_{\parallel }^{\text{e}}$ and ${\kappa }_{\times }^e/{\kappa }_{\parallel }^{\text{e}}$ The ratios are computed over the parallel conductivity of electrons $\left(k_{\parallel }^{\text{e}}\right)$ which is the most dominant at almost all heights. In this plot, we can see that below the transition region (z = 2.1 Mm), the perpendicular conductivity of electrons $\left({\kappa }_{\perp }^{\text{e}}\right)$ is the second more important and is larger than all the conductivities from ions. The parallel conductivity of ions $\left({\kappa }_{\parallel }^i\right)$ is proportional to the parallel conductivity of electrons $\left({\kappa }_{\parallel }^{\text{e}}\right)$ and is around one order of magnitude smaller. This might be evident since we are assuming equal temperatures for both electrons and ions to plot this figures. However, it shows that in the corona, the parallel conductivity of ions $\left({\kappa }_{\parallel }^i\right)$ is the second most dominant. Additionally, both transverse conductivities for ions and electrons $({\kappa }_{\times }^i$ and ${\kappa }_{\times }^{\text{e}})$ may contribute below the transition region. At other heights, the other components do not contribute much in this atmosphere model since they are smaller than the numerical perpendicular conductivity ratio of the code (~10⁻⁵), which was found in Sect. 4.2.

2D test with a realistic solar temperature model. We present an adaptation of the hot plate test, discussed above in Sect. 4.3, to a realistic density and temperature profile for the solar atmosphere. This is done to verify the performance of our implementation of the thermal conductivity using realistic high temperatures of the solar corona. Including conductivity in such model in the parabolic treatment decreases the time step three orders of magnitude in comparison to the advective time steps.

The initial state consists of a stratified temperature profile of Vernazza et al. (1981) smoothly joint to a constant-temperature corona at 1.2 × 10⁶ K. The pressure and densities are obtained from the hydrostatic equilibrium equation and the ideal gas law $$P\left(z\right)=p_0\exp \left(-\frac{m_{p}g}{k_{\text{B}}}{\displaystyle \underset{z0}{\overset{z}{\int }}\frac{\text{d}\tilde{z}}{T\left(\tilde{z}\right)}}\right),$$(40) $$P\left(z\right)=\frac{m_p}{k_{\text{B}}}\frac{p\left(z\right)}{T\left(z\right)},$$(41)

where m_p is the proton mass, k_B is the Boltzmann constant, p₀ = 0.1 Pa and z₀ = 10 Mm is a reference height. The magnetic field strength is taken from Eq. (39), and we set it to be inclined by 45° degrees. The numerical domain extends over [−1, 1] × [0, 7] Mm with 400 × 1400 points. We set a hot region at the top boundary between x = −0.25 Mm and x = 0.25 Mm with a temperature 10% higher than the equilibrium. The boundary conditions are periodic in the x direction and fixed to the initial values at the top and bottom.

We have evolved this setup using the scheme #1 with the time step dt = dt_TC, and with the scheme #2 we set the time steps dt = 50dt_TC and dt = 100dt_TC. The top panels of Fig. 8 display colormaps of the temperature and magnetic field lines at t = 25 s obtained by the different schemes. It can be observed that, as expected the hot temperature is advected along the magnetic field lines: No important differences are observed between the results obtained by the different schemes. The bottom left panel shows a more detailed comparison of a 1D cut of the temperature at x = −0.15 Mm, where we find that the solutions exhibit a similar tendency. Finally, the bottom right panel shows the cpu times for the different runs. The results suggest that the hyperbolic treatment can model heat conductivity under realistic coronal values, reducing significantly the computational cost. In these simulations the advection time step was 2893 times larger than the conduction one. For the simulation where the time step is dt = 50dt_tc, the efficiency ratio was of 31.38, and when dt = 100dt_tc it was of 60.1.

Fig. 6 Unstable examples of the tests using scheme #2 with time steps beyond the stability limit. Top left shows the temperature evolution in the one-dimensional test with dt = 200dtTC. Top right presents a temperature colormap for the static ring problem using dt = 33dtTC and t = 100. The bottom left panel shows a temperature colormap of the temperature for the hotplate test, using a time step dt = 500dttc. Bottom right panel shows a contour of temperature in the three-dimensional static ring problem in an early stage, t = 5 s, using a time step dt = 45dttc.

Fig. 7 Thermal conductivity coefficients for a stratified solar atmosphere with the temperature model of Vernazza et al. (1981) extended to the corona and a magnetic field for the quiet sun given by Eq. (39). Left panel: different components for ions and electrons. Right panel: ratios between different conductivity coefficients and the parallel electron conductivity coefficient.

6 Discussion and Conclusions

We have described the implementation of two numerical methods to model thermal conductivity in the code MANCHA3D. The scheme #1 corresponds to the explicit evolution of the parabolic heat flux term, which leads to very restrictive time steps. Scheme #2 is a hyperbolic treatment of the heat flux which allows the code to achieve considerably large speed-ups since it solves a separate hyperbolic equation for the parallel heat flux. Both schemes were tested in one, two and three spatial dimension by reproducing standard tests for anisotropic thermal conductivity.

We used the general heat flux expressions derived by Braginskii (1965). This formulation recovers the commonly used Spitzer expression in the limit of the strong magnetic field, but also takes into account the effects of the conductivity perpendicular and transverse to the magnetic field. The importance of the latter conductivity coefficients depends on the ratio between the collisional and cyclotron plasma frequencies, i.e. on plasma magnetization. This implementation is needed for modeling the whole solar atmosphere since it does a smooth transition between solar corona with remarkably anisotropic conductivity and the lower regions with large plasma ß and high collisionality where the heat flux becomes isotropic.

We propose to use a combination of both schemes, the parabolic treatment for the perpendicular and transverse conductivities and the hyperbolic treatment for the parallel conductivity, which is the one that imposes strict time steps due to its dependence on temperature. We show that in some regions of the atmosphere the perpendicular and transverse conductivities are not negligible. Namely, below the transition region, the perpendicular and transverse conductivities for ions and electrons might increase the heat flux. We present a two-dimensional test for the heat conduction in a stratified atmosphere, prove the robustness of the implemented methods to model thermal conduction, and confirm the significant speed up of the hyperbolic treatment over the computational cost of the simulations.

Our adaptation is similar to the one previously considered in the MURAM code (Rempel 2017) and the PENCIL code Warnecke & Bingert (2020). However, only the field aligned heat flux is considered in these two codes, which is, strictly speaking, valid only for the corona. In those works separate simulations from different zones of the atmosphere are coupled. The heat flux implemented in our code allows smooth transition from the corona to the convection zone, bringing us one step closer to realistic modeling of the solar atmosphere including the corona.

A drawback of the adaptation we have used is that the hyperbolic scheme slows down the temperature evolution in regions of high temperature and low density like the corona, which would normally slow down the code significantly, however this temperature variation should not affect much the physical results, as has been confirmed by Rempel (2017) and Warnecke & Bingert (2020). We also would like to point out that further analysis is still needed, since we have not considered complex magnetic configurations in an stratified atmosphere and we have not considered other non-ideal terms so far. In addition, an important issue that needs to be addressed in such complex scenarios is is how to increase correctly the magnitude of the timestep without reaching the unstable limit for the hyperbolic scheme. This matter will be addressed in our next project, in which we will study solar convection in a simulation box extended to the corona, focusing on the interplay between the thermal conduction, the radiative losses and other non-ideal effects. We emphasize that the scheme will allow running simulations of the whole atmosphere with reasonable time steps, limited by advection processes below the chromosphere and not by the thermal conduction in the hottest part of the computational domain. Furthermore, we are planning to extend this work with the self-consistent two-fluid model developed by Hunana et al. (2022) which overcomes the shortcomings of the Braginskii formulation in the weakly collisional regime.

Fig. 8 Top panels: colormaps of the temperature and magnetic field lines at t = 25 s obtained with the different schemes and time steps. Bottom left panel: 1D cut of the temperature. Bottom right panel: computational times.

Acknowledgments

The authors thank the support by the European Research Council through the Consolidator Grant ERC-2017-CoG-771310-PI2FA and by the Spanish Ministry of Economy and the Industry and Competitiveness through the grant PGC2018-095832-B-I00.

References

All Tables

All Figures

Fig. 1 Left panel shows the evolution of temperature in the one-dimensional thermal conductivity test computed with the scheme #1 (explicit evolution of the parabolic term) and with the scheme #2 (hyperbolic equation for heat flux) with three different time steps. Right panel shows a comparison between simulations using scheme #2, one including the saturation factor of Eq. (11) and another one without it.

In the text

	Fig. 2 Two-dimensional static ring problem. Top panels show colormaps of the temperature and magnetic field lines at the initial and final states. The bottom left panel shows a comparison of the temperature at y = 0 between the different simulations at the final stage. Bottom right panel gives a graphical representation of the computational times.
In the text

Fig. 3 Graphical representation of the method used to calculate the numerical perpendicular conductivity. The dots represent the final maximum temperature minus initial maximum temperature (∆T) in the simulations with different values of k⊥ shown at the horizontal axis. The dashed line is the difference of temperature (∆Tnum) in the simulation with κ⊥ = 0. We fit the data to a hyperbolic tangent function (solid line). The numerical perpendicular conduction coefficient к⊥,num is estimated by locating the position where the first derivative of the fit is 100 times smaller than its maximum value.

In the text

	Fig. 4 Left top to right middle plot: colormaps of the temperature, together with magnetic field lines, at t = 1.0, obtained in the different simulations, using scheme #1 and scheme #2 with the time steps dt = 10dttc, dt = 100dttc, dt = 200dttc and dt = 400dttc. The bottom left plot shows a 1D cut of the temperature at z = −0.27, and the right bottom panel shows the simulation time for each run.
In the text

Fig. 5 Results for the three-dimensional static ring problem. The top panels correspond to the contour of temperature at the initial and final stages, together with the magnetic field lines, obtained with the simulation with scheme #1. Bottom left panel shows a comparison of one snapshot of the temperature at x = −0.5 and y = 0.21 at t = 10 for different simulations using scheme #2 and different time steps. In the bottom right panel, a graphical representation of the computational times of the different runs is given.

In the text

Fig. 6 Unstable examples of the tests using scheme #2 with time steps beyond the stability limit. Top left shows the temperature evolution in the one-dimensional test with dt = 200dtTC. Top right presents a temperature colormap for the static ring problem using dt = 33dtTC and t = 100. The bottom left panel shows a temperature colormap of the temperature for the hotplate test, using a time step dt = 500dttc. Bottom right panel shows a contour of temperature in the three-dimensional static ring problem in an early stage, t = 5 s, using a time step dt = 45dttc.

In the text

	Fig. 7 Thermal conductivity coefficients for a stratified solar atmosphere with the temperature model of Vernazza et al. (1981) extended to the corona and a magnetic field for the quiet sun given by Eq. (39). Left panel: different components for ions and electrons. Right panel: ratios between different conductivity coefficients and the parallel electron conductivity coefficient.
In the text

	Fig. 8 Top panels: colormaps of the temperature and magnetic field lines at t = 25 s obtained with the different schemes and time steps. Bottom left panel: 1D cut of the temperature. Bottom right panel: computational times.
In the text