© The Authors 2024

1 Introduction

Our closest star continues to offer a wealth of enigmatic and mystifying phenomena that forge an exciting laboratory for physics relating to plasma, particles, nuclear, fluid dynamics, and thermodynamics, to name a few. Since the early 1950s, researchers have persevered to connect theory to observations, and vice versa. The development of faster and larger supercomputers, as well as better observational equipment over the past decade, has further accelerated this progress, and yet some of the oldest questions, such as the coronal heating problem, still remain unanswered. If such questions are to be answered, the need to capture the complexity of varying physical processes that occur across the full extent of the solar atmosphere must be met.

To that end, Bifrost is a 3D massively parallel magnetohydrodynamics (MHD) code, aiming to simulate stellar atmospheres from the upper convection zone through to the lower corona in detail. First and foremost, Bifrost’s aim is realism, and as such the precise implementation of additional physics such as radiative transfer, the generalised Ohm’s law, and thermal conduction is a key attribute in the development of the code. We focus on thermal conduction, which causes additional transport of heat in the upper solar atmosphere (chromosphere and corona). The contribution manifests itself in the internal energy equation, $$\frac{\partial e}{\partial t}=-\nabla \cdot eu-p\nabla \cdot u+Q,$$(1)

where Q represents the divergence of the heat flux from thermal conduction.

Spitzer & Seeger (1963) showed that the heat flux behaves as a non-linear diffusion equation in the direction of the magnetic field, such that $$Q=-{\nabla }_{\Vert }\cdot (q),$$(2) $$q=-{\kappa }_0{T}^{5/2}{\nabla }_{\Vert }T,$$(3)

where κ₀ ~ 10⁻⁶ erg cm⁻¹ s⁻¹ K^−7/2 is the Spitzer conduction coefficient. A full decomposition of the conductive flux, including perpendicular and transverse components, can be found in Braginskii (1965). In the corona, where the temperatures are largest, and hence the conductive flux has the most effect on the internal energy, the perpendicular and transverse terms are found to be much smaller than the parallel terms. In the chromosphere, however, the perpendicular contribution becomes comparable to the parallel conductive flux (Navarro et al. 2022). Nevertheless, perpendicular and transverse values are on the same scale as numerical diffusion within Bifrost and so these components are not currently implemented (Gudiksen et al. 2011).

When discretised explicitly using a numerical scheme, for example as $$\frac{\partial e[i,n]}{\partial t}=\frac{e[i,n]-e[i,n-1]}{\mathrm{\Delta }t},$$(4) $$\nabla q[i,n]=\frac{q[i,n]-q[i-1,n]}{\mathrm{\Delta }x},$$(5)

in the 1D case, the second-order diffusive operator limits the Courant condition to scale with Δx². This greatly limits the suitable time-stepping (Δt) for finer grids, which, even using parallel codes with high-performance computing, requires large amounts of memory and CPU hours to compute. Simulating the long-term evolution of large-scale events, with timescales on the order of solar hours, or even days, becomes time-consuming and error-ridden. Instability within explicit methods is also difficult to manage and requires large artificial diffusive parameters to maintain it.

Bifrost employs two other methods of solving thermal conduction. The first uses the Crank-Nicholson implicit method described in Gudiksen et al. (2011), $$\frac{e[n+1]-e^{*}}{\mathrm{\Delta }t}=\theta {\nabla }_{\Vert }\cdot q[n+1]+(1-\theta ){\nabla }_{\Vert }\cdot q[n].$$(6)

The quantities at time-step n are calculated before the MHD time-step, starred quantities are produced after the MHD timestep, and [n + 1] quantities are deduced implicitly, using a multi-grid method (Malagoli et al. 1996). Finally, the energy term is overwritten altogether. This method is known to be unconditionally stable, but still subject to other numerical artefacts and error build-up. Nonetheless, it is a more robust method than the explicit version, offering larger time-steps with stability, but at the cost of more computing time due to the complex machinery of the method.

The second method implements the hyperbolic PDE approximation proposed by Rempel (2017), $$\frac{\partial q}{\partial t}=\frac{1}{\tau }\left(-q-{\kappa }_0{T}^{5/2}{\nabla }_{\Vert }T\right),$$(7) $$\frac{\partial e}{\partial t}=\dots -\nabla \cdot q,$$(8)

which converges on the usual parabolic Spitzer Equation (3) over time. For τ > 0, this is the equivalent of solving a wave equation for T , and as such the solution is a wave, with a propagation speed of $c=\sqrt{\kappa /\tau }$, when κ and τ are constant. The wave method, as it will now be referred to, is robust, reducing any instabilities that occur, and faster to compute than the above implicit method, as it is still calculated explicitly. It must be noted that this method only converges on the correct flux values in time, relative to the convergence timescale τ. Care must be taken with the choice of τ. Since the Courant condition now also depends on 1/τ, the convergence rate must not become too small or the system will become unstable.

It can be seen that each method has its positives and negatives. The following sections discuss the utilisation of these methods in the context of MHD codes. Section 2 discusses the structure of the numerical schemes and specifics of the machinery in the Spitzer module within Bifrost. Sections 3 and 4 rigorously test the capabilities of each method, through two idealised tests and comparing them to analytical solutions. We conclude with a discussion on how to use these methods optimally when applied to more complex MHD simulations and remark on the development of this module for the future.

2 Numerical schemes

Since thermal conduction is not a driver of the evolution of any MHD quantities except internal energy, its contribution may be calculated independently of the main MHD calculations in a simulation. For the explicit and wave methods, the calculated conduction component is simply added to the current $\frac{\partial e}{\partial t}$ before the next time-step. For the implicit method, the final internal energy, e, is overwritten to include the effect of the conductive flux. Thus, whilst Bifrost’s main calculation works on an explicit Hyman or Runge-Kutta method, the thermal conduction calculations are free to be calculated in whichever method suits the simulation. It is important to understand how the methods’ different behaviours affect the outcome of the Spitzer module, and ultimately the heat transfer in any simulation.

2.1 Explicit method

The explicit method is a first-order Euler method. The greatest cause of instability for this method is the Courant condition, $C~\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}$. We can split this into three parts for each derivative calculated; two spatial (∇T, ∇ ⋅ q) and one temporal $\left(\frac{\partial e}{\partial t}\right)$. Firstly, in order to ensure stable handling of high temperatures within the code (important for the corona, where T can reach above 10 MK), we calculate the temperature gradient of log(T), such that the flux term becomes $$q=\kappa (T){\nabla }_{\Vert }\log T={\kappa }_0{T}^{7/2}\widehat{b}(\widehat{b}\cdot \nabla )\log T,$$(9)

where b̂ = B/|B| is the unit vector for the magnetic field. This also helps to lessen the instability within the first spatial derivative due to steep gradients. The second derivative essentially concerns the rate of change of the temperature gradient, ∇² T. This can cause issues around, for example, shock regions, as is seen in Sect. 4. Finally, the time derivative must heed that the total change in energy per time-step does not exceed the amount of energy in the system, or else the numerical model will break and the solution will become unphysical. Bifrost monitors this relation through the condition $$\kappa (T)={\kappa }_0{T}^{7/2}<\beta \left(\frac{\mathrm{\Delta }{x}^2}{\mathrm{\Delta }t}\right)e,$$(10)

for β є [0, 1). This is observed before the final spatial derivative is taken such that problems due to large gradients in T can also be monitored this way. In Bifrost, β = 0.7, which allows for energy loss from other processes.

2.2 Implicit method

This method is an adaptation and extension to 3D of the nonlinear multi-grid scheme presented in Malagoli et al. (1996). Firstly, the hydrodynamic energy equation is calculated, alongside an explicit approximation, very similar to above, for κ(T)∇T at the current time-step, which is treated as a source term for the implicit calculation. The weight of this first step to the total conductive flux contribution is decided by θ є [0, 1]. For θ = 0, we recover the explicit method. As such, it is important that θ is suitably large that instabilities are completely smoothed for larger time-steps. If θ ~ 1 however, the implicit method will diffuse the solution too much unnecessarily. Here, we used θ = 0.9, although the results for the following tests do not change significantly for θ > 0.5, where the instability of the explicit part begins to dominate for θ < 0.5.

The Crank–Nicholson implicit calculation uses a non-linear multi-grid solver for the inversion of the operators, which accelerates the convergence of the relaxation iterations. It performs corrections to the calculated solution on a selection of fine to coarser grid patterns, in order to reduce errors with wavelengths comparative to each grid size (Millar & Bancroft 2003). Once the coarsest resolution has been reached, the resolution is once again increased back to the original resolution, in a pattern known as a V-cycle, due to its implied shape. For each grid resolution, h, the interaction between it and the next coarsest level, H = 2h, uses the full approximation scheme (see Brandt & Livne (2011) for a detailed explanation):

1. At level h, the solution to the implicit non-linear problem Lu = f, where u is unknown, is denoted as L^hu^h = f^h.

2. After a few Gauss-Seidel relaxation sweeps, the smoothed solution is approximated as ũ^h.

3. A correction term, v^h, is needed, such that v^h + ũ^h = u^h.

4. The residual, r^h = L^h(ũ^h + v^h) − L^hũ^h, is injected into the coarser grid level, H .

5. Using an interpolation scheme, $I_h^H$ , the following values can be interpolated to the new resolution:

   • $u\to u^H=I_h^H{\tilde{u}}^h+v^h$

   • $f\to f^H=L^H\left(I_h^H{\tilde{u}}^H\right)+{\widehat{I}}_h^H{\tilde{u}}^h$, where I and Î pertain to separate interpolation processes,

   such that the linear problem becomes L^Hu^H = f^H.

6. An approximation, ũ^H, is calculated for the above using Gauss-Seidel sweeps. The new correction term can be approximated by $v^H={\tilde{u}}^H-I_h^H{\tilde{u}}^h$.

7. The coarse correction term is injected back into the finer grid scheme through interpolation, and used to update the finer solution approximation, ${\tilde{u}}_{\text{NEW}}^h={\tilde{u}}^h+{\widehat{I}}_H^h{v}^H$.

Between steps 6 and 7, it is possible to include more levels of coarser grids, by injecting the new residual to lower grids and repeating steps 5 and 6, before updating each solution approximation in the finer grids. By the coarsest grid, the calculations must be processed on only one core in order to keep communication costs low in relation to the computations. Each level is iterated for a select number of repetitions and the full cycle is also repeated as desired. The more repetitions, the faster the solution will converge, and so it is beneficial to iterate each level multiple times. The coarser grids may be repeated many times without a large computational cost, and as such tend to be repeated more than the finer grids. Bifrost has four levels, with each level half the resolution of the level above it. For the purpose of this paper, the iterations of each level are 4, 4, 4 and 16, respectively, with the whole cycle being repeated twice.

2.3 Wave method

By rearranging Equation (7), we recover the hyperbolic equation for calculating q, given by q = qb̂, proposed by Rempel (2017), $$\frac{\partial q}{\partial t}=\frac{1}{\tau }\left(-\frac{q^{*}}{f_{\text{sat}}}-q\right),$$(11)

and where q* = κ₀ T^5/2(b̂ ·∇) T recovers the original Spitzer equation. Rempel discusses appropriate values for the saturation rate , f_sat, and convergence rate, τ, as $$f_{\text{sat}}=1+\frac{|q*|}{1.5\rho C_S^3},$$(12) $$\tau ={\left(f_{\text{CFL}}\frac{\mathrm{\Delta }{x}_{\text{min}}}{\mathrm{\Delta }t}-\left|v\right|\right)}^{-2}\frac{f_{\text{sat}}{\kappa }_0{T}^{7/2}}{e},$$(13)

where $C_S=\sqrt{\gamma p/\rho }$ is the ideal gas sound speed. It is possible to omit the saturation factor, f_sat, in Equation (11) for faster convergence, but this also decreases the stability of the method (Navarro et al. 2022), as is discussed in Sect. 4. The choice of the equation for τ is flexible and somewhat depends on the code. For example, Navarro et al. (2022) instead opt for a constant τ = 4Δt throughout, whereas Warnecke & Bingert (2020) relate τ to the Alfvén velocity and the conductivity. In the latter case, lower and upper limits are both imposed on τ for similar reasons to the ones discussed below. f_CFL was taken to be 0.7 for our tests. However, the choice, f_CFL є [0,1], has a negligible effect in these tests, since we have ignored dynamics (e.g spatially varying velocities) and because the dominating term is proportional (Δt)⁻².

This method creates an additional calculation that must be taken at each time step, whereby the new flux term, q, is calculated from the time derivative in Equation (11). Therefore, for time step [n], there are two operations:

1. Add the flux calculated at the previous time to the energy equation, $$e[n+1]=e[n]+\mathrm{\Delta }t(\dots -\nabla \cdot (q[n])).$$(14)

2. Calculate the flux to be used at the new time using Equation (11), $$q[n+1]=q[n]+\mathrm{\Delta }t\left[\frac{1}{\tau }\left(-\frac{q^{*}[n]}{f_{\text{sat}}}-q[n]\right)\right].$$(15)

Therefore, the stability of a fourth derivative calculation, in addition to those already discussed in the explicit method, must be considered. From Equation (13), it is clear that τ scales with (Δt)², and as such, the Courant condition for Equation (15) scales with 1/Δt. In other words, the system becomes unstable for smaller time-steps. Therefore, it is important that the wave implementation contains a ‘throttle’ in the form of a lower limit on τ, $${\tau }_{\min }=\alpha \mathrm{\Delta }t,$$(16)

where the choice of α is discussed more in Sect. 4.

2.4 Adaptive time-stepping

By default, Bifrost varies the time-step of the current running simulation, by monitoring the overall Courant condition. When additional physics modules, such as thermal conduction, are involved, it is important that the local Courant condition, and any other conditions for the chosen method, are also considered. The limit in Equation (10) for the explicit method is an example of this. We can compare the tests described below for both constant and variable time-stepping, in order to understand where such limiting factors may come in to play within a simulation.

3 Methodology

We now present two tests with which to validate and analyse the methods. The first is a 1D convergence test, described before by Rempel (2017) and also used by Navarro et al. (2022) in their testing of the Mancha3D code. Thus, this test is suitable for validating the implementation and efficacy of the wave and explicit methods. We note that the implicit multi-grid solver is tuned for at least a 2D set-up, so was not included for the 1D test. The second test calculates the diffusion of a Gaussian temperature profile in a 2D set-up. The numerical results are then compared to an analytical first-order approximation by Furuseth et al. (2024). Both tests require a purely thermal setup, such that all MHD time derivatives, except energy, remain zero. Examples of these tests can be seen in Figs. 1 and 2, respectively.

Fig. 1 Diffusion of the wave method with and without saturation for Δt = 10−6. The explicit solution for Δt = 10−7 is indicated by dashed lines in the final two panels. The solid lines indicate evolution at times t = 0.0, 0.125, 0.25, 0.5, 1.0, and 2.0. The dotted line references the analytical solution for convergence.

Fig. 2 Evolution of the self-similar diffusion test for Gaussian and selfsimilar initial conditions. At t = 1.40 s, the profile lies in the χ t >> 1 regime and both initial profiles show a self-similar diffused profile, albeit for different ϕ0 . The solutions are run with the explicit method for Δt = 10−5 s. We note that different scales are used for the colour maps between t = 0 s and t = 1.4 s.

3.1 1D convergence test

We considered the 1D non-linear heat diffusion equation, $$\frac{\partial T}{\partial t}={\kappa }_0\frac{\partial }{\partial x}\left(T^{5/2}\frac{\partial T}{\partial x}\right),$$(17)

with an initial temperature profile of $$T_0(x)=0.1+0.9x^5,$$(18)

and fixed boundaries such that T(0) = 0.1 and T(1) = 1 for all t. For t >> 0, the temperature solution converges ($\frac{\partial T}{\partial t}=0$) to the asymptotic stationary solution $$T={\left[{0.1}^{7/2}+\left(1-{0.1}^{7/2}\right)x\right]}^{2/7},$$(19)

where κ₀ controls the speed of convergence.

Our domain contained 256 grid points, extending across x = [0, 1]. We normalised the set-up to unity, providing a constant 1D magnetic field, B_x = 1.0, and density, p = 1.0. Internal energy and pressure were calculated in the ideal gas regime, such that they became directly proportional to T . We note that the above solution pertains to a time evolution in T, $$\frac{\partial T}{\partial t}=-{\nabla }_{\Vert }\cdot q(T),$$(20)

whereas our implementation was directly solved in the energy equation, and so pertained to a time evolution in e (see Equation 1), $$\frac{\partial e}{\partial t}=\frac{\rho }{\gamma -1}\frac{\partial T}{\partial t}=-{\nabla }_{\Vert }\cdot q(T),$$(21)

in the ideal gas and normalised regime. Therefore, given a constant density, the speed, κ_e, at which the test within Bifrost converged to the solution (19) is related to the Spitzer coefficient κ₀ by $${\kappa }_e=\frac{\rho }{\gamma -1}{\kappa }_0,$$(22)

where we have normalised the Boltzmann constant, k_B , ion mass, m_H, and mean molecular weight, µ, to 1.0. Furthermore, we set γ = 5/3 and κ₀ = 1.0 to agree with the tests in Rempel (2017) and Navarro et al. (2022).

3.2 Self-similar diffusion test

The self-similar solution to a non-linear diffusion equation with concentration-dependent coefficients in up to three dimensions, has been studied analytically by Pattle (1959). Analytical descriptions of the evolution and convergence from an instantaneous point source with a zero background concentration, and finite boundary conditions (T(r) = 0 for finite r) are provided. In the case of thermal conduction, this concentration is taken to be temperature, T . Taking the dimension as one, since diffusion is only along the magnetic field, this leads to the following behaviour:

• Evolution of a point source with initial central value T₀ and width R₀ at the foot-points $$T(r,t)=\{\begin{array}{ll}{T}_0{(1+\chi t)}^{-2/9}{\left(1-\frac{r^2}{R{(t)}^2}\right)}^{2/5},\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}r<R(t)\hfill \\ 0,\hfill & \text{\hspace{0.17em}otherwise}.\hfill \end{array}$$(23) $$R(t)=R_0{(}^1,$$(24) $$\chi =\frac{9}{5}\frac{2{\kappa }_0{T}_0^{5/2}}{R_0^2},$$(25)

  where T₀ and R₀ are determined by the total surface integral ${\varphi }_0=\frac{T_0{R}_0}{B}$ for the Beta function $B=B\left(\frac{1}{2},\frac{7}{2}\right)$.

• Convergence when χ^t >> 1 of any distribution with total surface integral ϕ₀ $$T(r<R(t),t)\to T_0{\left(\frac{5}{9}\frac{{\left(B{\varphi }_0\right)}^2}{2{\kappa }_{0}t}\right)}^{2/9}{\left(1-\frac{r^2}{R{(t)}^2}\right)}^{2/5}$$(26) $$R(t)\to {\left(\frac{9}{5}{\kappa }_0{\left(B{\varphi }_0\right)}^{5/2}t\right)}^{2/9}.$$(27)

Equations (26) and (27) provide an important result of the solution, that the converged self-similar solution does not depend on the initial central value of the distribution, and as such, any distribution that has the same ϕ₀ as the given instantaneous point source, and becomes zero in a finite domain, will in fact converge on the same shape.

This solution has been used already in Bifrost to test and validate the implementation of ambipolar diffusion (Moreno-Insertis et al. 2022). However, since the physical quantity we considered is temperature, the system became unphysical at 0 Kelvin. We therefore considered a small background temperature, T_∞, such that T_tot(r, 0) = T_∞ + T(r, 0), and compared to a first-order modification of the self-similar solution in the regime T_∞ << T₀ (Furuseth et al. 2024): $$\max \left[T_{\text{Tot}}^{(1)}(t)\right]=\{\begin{array}{ll}{T}_{\infty }+T_0{\left[1+\chi t{\left(1+\frac{T_{\infty }}{T_0}\right)}^{5/2}\right]}^{-2/9}\hfill & ,\text{\hspace{0.17em}if\hspace{0.17em}}\chi t\ll 1\hfill \\ T_{\infty }+T_0{\left[1+\chi t{\left(\frac{{\varphi }_{\text{Tot}}}{{\varphi }_0}\right)}^{5/2}\right]}^{-2/9}\hfill & ,\text{\hspace{0.17em}if\hspace{0.17em}}\chi t\gtrsim 1\hfill \\ \to T_{\infty }+{\left(\frac{5}{9}\frac{{\left(B{\varphi }_{\text{Tot}}\right)}^2}{2{\kappa }_{0}t}\right)}^{2/9}\hfill & ,\text{\hspace{0.17em}if\hspace{0.17em}}\chi t\gg t.\hfill \end{array}$$(28) $$R^{(1)}(t)=\{\begin{array}{ll}{R}_0{\left[1+\chi t{\left(1+\frac{T_{\text{∞}}}{T_0}\right)}^{5/2}\right]}^{-2/9}\hfill & ,\text{\hspace{0.17em}if\hspace{0.17em}}\chi t\ll 1\hfill \\ R_0{\left[1+\chi t{\left(\frac{{\varphi }_{\text{Tot}}}{{\varphi }_0}\right)}^{5/2}\right]}^{2/9}\hfill & ,\text{\hspace{0.17em}if\hspace{0.17em}}\chi t\gtrsim 1\hfill \\ \to {\left(\frac{9}{5}2{\kappa }_0{\left(B{\varphi }_{\text{Tot}}\right)}^{5/2}t\right)}^{2/9}\hfill & ,\text{\hspace{0.17em}if\hspace{0.17em}}\chi t\gg 1.\hfill \end{array}$$(29)

Our set-up took a Gaussian temperature profile in a 2D domain, with an initial background temperature of 1% of the peak total temperature T_Tot, $$T=9.9\times {10}^5\exp \left(-\frac{x^2+z^2}{{0.01}^2}\right)+1\times {10}^4.$$(30)

It should be noted that the analytical solutions pertain to an initial condition with a self-similar shape (see Equation (23)), but we show (see Fig. 2) that in the regime important to us, χt >> 1, the Gaussian initial condition converges to the selfsimilar solution with the same perturbed area, ϕ₀ , as with the converged zeroth-order Equations (26) and (27). In our case, the self-similar profile has R₀ = 0.1 Mm, T₀ = 9.9 × 10⁵ K, and T_∞ = 1 × 10⁴ K.

In our (1 × 1) Mm domain, centred about (0, 0), and split into (256 × 256) grid points, we imitated coronal conditions (Carlsson et al. 2019) with low density, ρ = 10⁻¹⁵ g cm⁻², following the ideal gas regime, and weak magnetic field in one direction, B_z = 0.316 G, B_x,y = 0textG. Whereas the magnitude of the magnetic field is arbitrary, the density value affects the speed of diffusion, as discussed in the previous test, but this time without unity normalisation. Therefore, when comparing to the analytical solution, we modified the analytical κ₀ → κ_α, $${\kappa }_a=\frac{(\gamma -1)\mu m_h}{k_B\rho }{\kappa }_0,$$(31)

where k_B = 1.38 × 10⁻¹⁶ erg K⁻¹ is the Boltzmann constant, m_h = 1.67 × 10⁻²⁴ g is the ion mass, and µ = 0.5 relates to full hydrogen ionisation, along with γ = 5/3. The plasma was initially at rest (u = 0), and remained at rest throughout, relating to pure thermal diffusion. The boundary conditions for the domain were periodic. Where there is a need to analyse diffusion over a longer time period (see for example 7), the domain was extended to (2 × 2) Mm with (512 × 512) grid points, such that the resolution remained the same and the boundaries never directly interacted with the profile.

4 Results

4.1 Accuracy

Figures 1 and 2 give examples of the two tests described above. Bifrost’s explicit solution to the 1D test, agrees completely with the results of Rempel (2017), for the case τ = 0, which corresponds to the normal parabolic equation. We validate that the wave method without saturation factor f_sat (right panel) converges quickly to the parabolic solution, sometime before t = 0.125, since the results are comparable to those of the explicit method. In contrast, with added saturation (middle panel), the wave method lags behind the explicit solution throughout, but still achieves convergence by t = 2. The benefits of the saturated method become apparent in Figs. 3 and 4 for variable timestepping. Where the unsaturated version (right panel) breaks down by t = 0.5 due to large instabilities, the saturated version is able to recover and maintain stability¹ by the convergence time of t = 2.0. We can see in Fig. 4 that the Δt values in both versions steadily rise above 8 × 10⁻⁴ . The larger instability within the unsaturated method causes the time-step to rise slower than for the saturated version, which reaches the maximum stable Δt (≈ 9 × 10⁻⁴) by t = 0.5. The unsaturated method, can no longer control the instability however, and the simulation crashes just after t = 0.5. This method would require a stricter control condition, lowering the maximum possible Δt for stable runs. For this reason, Bifrost defaults to using the saturated version, since the main benefit of the wave method is larger time-stepping over the explicit method, and a faster calculation time than the implicit method, which, as is described later (see Sect. 4.4), has a similar maximum time-stepping available as the saturated wave method.

From now on, we shall consider primarily the self-similar test. Figure 2 evidently shows that a Gaussian profile produces acceptable results to compare to a self-similar solution in the χ^t >> 1 regime. In the direction of the magnetic field, the tails steepen to form shock-fronts, as in the self-similar case. The other direction remains completely unchanged in its profile. Therefore, the Gaussian shape provides a better-behaving alternative to the self-similar initial condition, since the profile is naturally continuous everywhere, and does not have as large initial temperature gradients, so less ringing occurs within the solution. Figure 5 shows that the relative error, $$\text{\hspace{0.17em}Relative\hspace{0.17em}error\hspace{0.17em}}=\frac{T_{\text{numerical\hspace{0.17em}}}-T_{\text{analytical\hspace{0.17em}}}}{T_{\text{analytical\hspace{0.17em}}}},$$(32)

between the two initial conditions is sufficiently similar for longer timescales. The largest differences between the numerical method and the analytical solutions are at the peak and in the tails. The majority of the error at the tails is a consequence of neglecting higher orders in the approximation (Furuseth et al. 2024), so the best area for analysis and validation of the numerical methods comes from the temperatures near the peak of the profile. It is clear to see the larger instabilities within the selfsimilar initial condition’s evolution, and this leads to a slight decrease in the converged peak temperature, where the Gaussian converges towards a relative error of zero at the peak.

Fig. 3 Diffusion of the wave method with and without saturation for non-constant Δt. The explicit solution is indicated by dashed lines.

Fig. 4 Evolution of Δt in time for the adaptive time-stepping wave method in the 1D test, starting at Δt = 1e − 6.

4.2 Time-stepping

Figure 5 is the result of a stable explicit run, with constant time-step, Δt = 10⁻⁵ seconds. This is close to the upper limit for a stable explicit Δt, due to the condition derived in Equation (10). On the other hand, the wave method stays stable for time-steps over two orders of magnitude larger, becoming unstable in the self-similar test for constant Δt >> 3 × 10⁻³ seconds. As Δt increases, the method becomes increasingly wavy² in the initial iterations, and so the magnitude of the relative error is much larger at early times than for the explicit method. Figure 6 around time t = 0.5 s is an example of this, where the wave method with adaptive time-stepping reaches larger Δt than the constant time-step run, and as such, the relative error stays larger for longer. This is also partly due to the larger time-stepping, relating to a larger value of τ, and so a slower rate of convergence. In Fig. 7, we see that the adaptive time-step runs all have a period of time where the time-step decreases in order to manage the instability from initial waviness. These smaller Δt also allowed the adaptive run in Fig. 6 to converge slightly faster than the constant run. The erratic characteristics of the relative error later on in the same run (such as around t ≈ 1.7 s) are due to additional waviness as the Δt increases to near the stable limit again.

As with the 1D convergence test, the wave solution lags far behind the analytical solution initially (see Sect. 4.2), and a considerable amount of time passes before the method reaches errors of the same order as the explicit method; roughly 6 seconds, as is shown in Fig. 7. Due to such short time-steps needed, the explicit method would need to perform around 70 000 iterations to reach the point where the saturated wave method becomes similarly accurate, which is computationally expensive and timeconsuming. Excluding saturation in the wave method, as is seen in Fig. 6, performs almost identically regarding the relative error to the explicit method for constant time-stepping. The same cannot be said for the adaptive time-stepping version, which performs the worst of all the methods, due to the stricter conditions for stability needed, as was discussed before. Therefore, although the unsaturated method offers a time-step of an order larger than the explicit method, the requirement for constant time-stepping is still a large time-constraint, and most likely not compatible within larger simulations that need adaptive time-stepping to ease computational cost.

Another alternative would be to use the implicit method, since the implicit error in Fig. 6 starts in a similar fashion to the explicit case by quickly decreasing in magnitude from the initial overestimate (on average). However, the solution then overshoots and increases again, this time as an underestimate, before slowly decreasing with time to become an overestimate again. In the long term, the implicit error continues to grow as an overestimate, rather than converging on the analytical solution (see Fig. 8). This could highlight an issue within the multi-grid, such as over-smoothing or improper weighting between the before- and after-time-step subroutines. It could also just be the buildup in convergence errors if not enough iterations have been used within the multi-grid or relaxation regimes. However, in the period of time before the wave solution has converged, the implicit method provides a robust and stable alternative method which can run at larger Δt. There is some sacrifice in run time, due to the much larger computational time needed for implicit calculations, but since the maximum Δt is over an order greater than for the explicit method, the speed-up is still considerable.

The outcome of the variable time-stepping set-up depends significantly on the starting Δt value. As is discussed in Sect. 2, the value of τ is directly proportional to Δt². Choosing a smaller time-step to start will therefore allow the method to converge on the parabolic solution faster, but also expose the solution to more instability. The time-step must increase to more stable Δt before it can control this instability, and this takes more iterations for smaller starting Δt, meaning the method has unstable tendencies for a larger number of oscillations which can grow and crash the run. Bifrost contains a series of hyperdiffusion parameters (see Gudiksen et al. 2011) to combat such situations. With appropriate hyper-diffusion, the initial instability is sufficiently damped and lower starting time-steps, such as Δt = 10⁻⁵ s in Fig. 7 can be achieved. Then, the method converges quickly towards zero relative error due to lower τ values initially. Notably, the hyper-diffusion completely quashes the secondary error peak seen between 1.4 and 3 seconds, allowing the adaptive Δt to continue to increase to the maximum value over 4 simulation seconds faster than without hyper-diffusion (see Fig. 7, bottom panel). We can see in the top panel of Fig. 7, that the hyper-diffusive runs starting at Δt = 10⁻⁴ and 10⁻⁵ s both converge very well within 7 seconds, with the lower timestep proving the most accurate throughout, whereas the larger starting time-step of 5 × 10⁻⁴ s converges slower due to the overall smaller convergence rate, and larger instabilities.

Fig. 5 Evolution of the explicit solution for Δt = 10−5 s (solid), compared to the first- order analytical approximation (dashed). Here, the solution for every 0.4 seconds until 3.2 seconds is shown with lines from dark to light. The relative error of the Gaussian profile quickly becomes comparable to that of the self-similar profile for solutions in the χt >> 1 regime.

Fig. 6 Average relative error around the peak (T > 105 K) for various set-ups.

Fig. 7 Average relative error around the peak (T > 105 K) of variable time-stepping in the wave method, using hyper-diffusive terms (dashed).

Fig. 8 Long-term error convergence for the wave and implicit methods, with Δt = 5 × 10−4 seconds.

Fig. 9 Areas of application to the limit on τ within the wave method for various limiting values. Solutions are shown for every 0.8 seconds.

4.3 Safety limits

In Equation (13), we can see the various terms that can affect the size of the wave convergence factor, 1/τ. It should be noted that τ scales with T^5/2. This means at the peak of the profile, where T is large, the convergence will be at its slowest, whilst in the tails where T is smaller, the convergence rate can become very high (→ ∞). Here is where we see the activation of the throttle on τ, described in Equation (16). In Fig. 9, we test the activation of this limit for α = (0.7, 2.0, 4.0). As expected, the larger the limit, the more of the temperature profile has this limit applied to it. However, the relative error remains similar, independent of how much of the profile the limit is applied to. Below α = 0.7, both the 1D and self-similar simulations crash, due to instability travelling down the steep side of the profiles towards T_∞. In the α = 0.7 case, the instability exists in all iterations, but never becomes unmanageable. In the α = 2.0 case, however, the initial instability is dampened out by time t = 2.4 s, guaranteeing stability from there on. Since there seems to be no distinct loss in accuracy by using a larger α, Bifrost by default uses α = 2.0. At very small Δt, τ_min will be much greater than the calculated τ for the whole profile, and so the limit is applied everywhere. This will significantly reduce the rate of convergence and thus the short- and mid-term accuracy in those Δt regimes. However, for variable time-stepping, the proportion of time, both computationally and physically, spent at such small Δt is negligible and the accuracy is likely to be recovered quickly as Δt increases.

In contrast, what limits the time-stepping for the explicit method is the peak of the profile, where T, and hence the heat flux, is largest. As explained in Sect. 2 it is important that the heat flux does not become larger than the total amount of energy within in the system at any time. If a similar limit would be applied to κ in the explicit method, as it was with τ in the wave method, an issue arises as the Δt increases, since then κ_max, $${\kappa }_{\max }=\left(0.7\frac{\mathrm{\Delta }{x}^2}{\mathrm{\Delta }t}\right)e$$(33)

from Equation (10), becomes smaller and so the limit must be applied to more of the profile. Unlike the limit on τ, Fig. 10 shows that the application of this limit has a significant effect on the accuracy. This is since this limit directly effects the diffusion rate, whereas the limit on τ merely dampens the rate at which the method converges to the correct solution. Therefore, this must be the dominant factor that controls the time-step such that the flux is maintained appropriately, and this creates a strict limit on Δt, as is seen in Fig. 11.

4.4 Performance

Figure 11 compares the adaptive time-stepping for all methods. It is clear to see that the wave method has the advantage over the explicit method, with an average difference of over one order magnitude between the two time-steps at any point in time. The relationship between the implicit and wave methods is a little more complex. We see a dip in Δt for the wave method around t = 0.4 s, in order to control the initial waviness. The implicit time-step, however, is free to continue increasing, such that fewer iterations might be needed for the implicit method than the wave method for small t. Tables 1 and 2 provide details for simulation runs using different set-ups and methods at times t = 1.6 s and t = 3.0 s, respectively. These former time relates to the time at which the time-step of the wave method overtakes that of the implicit method, and the latter to final time in Fig. 11. We can see that the drop in Δt in the wave method does allow the implicit method to have fewer iterations for a period of time, as is seen in Table 1, but due to the wave method’s faster calculation time, the implicit method always requires more time than the wave method, even with fewer iterations. We can see from the constant time-step runs in Table 2, that the complex machinery of the implicit method makes it around three times slower than the calculations of the wave method. For t > 1.7 seconds, the adaptive time-step wave method also reaches larger timesteps than the implicit method, and so the disparity between computation time will continue to grow for longer runs.

The large jump in Δt for all methods in the first half a second is the main advantage of having adaptive time-stepping. Table 2 shows that running these methods with constant time-stepping needs a lot more iterations to complete the simulation since the strict conditions at the initial profile would apply throughout. The speed-up for each method is between 30 and 50 times faster using adaptive time-stepping, with the greatest speed-up found in the implicit method. The wave method clearly dominates in terms of time to completion, and the inclusion of higher hyperdiffusion allows for further acceleration of the time-step, leading to fewer total iterations. It is also the most perceptive method to changes in the initial time-step, due to the dependency on Δt within the convergence rate, τ. As expected, starting with a higher initial time-step leads to a faster run overall, but this is where it is important to balance the priority between computation speed and accuracy, since we have discussed before that the higher starting time-steps lead to a longer convergence with respect to the simulation time. For example, we can calculate from the wave test run with initial Δt = 5 × 10⁻⁴ s, and the hyper-diffusive wave test run with initial Δt = 1 × 10⁻⁴ s, using Fig. 7 and Table 2, that for a 14% increase in computation time, we gain an approximate 87% reduction in the average relative error at t = 3.0 simulated seconds. The initial time-step has much less influence on the implicit method, and since the method converges on the correct solution much faster, the choice of Δt at the start is less significant, and will most likely depend on external factors in more complex simulations.

Fig. 10 Application of a limit on the heat flux in the explicit method. Solutions are shown every 0.8 seconds.

Fig. 11 Δt in the adaptive time-stepping version of each method, with a starting Δt of 10−5 seconds.

5 Discussion

The analysis above has been carried out on pure thermal tests. In other words, all velocities, densities, and magnetic field remain constant throughout. However, the Spitzer equation is, in reality, an addition to the energy equation, which has knock-on effects on the entire magnetohydrodynamics system. It is important to understand how a change in heat flux, such as the wave method’s initial waviness, affects an evolving simulation. Similarly, it is important to understand how the MHD quantities interact with and influence the heat flux term. Testing this is outside the scope of this paper, but we can speculate that the inclusion of more physics will most likely create natural stabilisers for conduction, such that the time-stepping range increases and conditions for stability become dependent on factors outside of the Spitzer module. For example, with the inclusion of a non-stationary fluid (i.e. that v ≠ 0), the initial peak will diffuse faster, relaxing the conditions for stability quicker. Within the wave method, τ will be greater where the velocities are larger, which ensures that the wave method can still converge appropriately where the characteristics of the plasma may be changing quicker. However, additional physics also introduces more opportunities for shocks, and other extreme phenomena which may negatively impact the thermal conduction code’s performance. This would need to be tested in more complex simulations.

We have discussed the various methods of calculating thermal conductivity within the parallel MHD solar atmosphere code, Bifrost. Any simulation may be run using three different numerical methods for solving the heat flux term. The explicit method, as is expected, provides an accurate but computationally slow solution, since the conditions for stability require very small time-steps. This method is likely only to be used for short simulation runs where the precision of the heat flux term is an important characteristic. As such, it is unlikely to be used for long, realistic solar simulations, as the time-step limit is too taxing.

The most beneficial alternative is the saturated wave method from Rempel (2017). The maximum Δt available (allowing stability) is over an order of magnitude larger than for the explicit method, and there is a very low relative error after an initial convergence period. There are two disadvantages to the saturated wave method. Firstly, that the code can become unstable for small Δt. It is often useful to start or restart a simulation with a very small time-step to let changes in set-up, such as an increased magnetic field, settle to an equilibrium quickly. We have shown that with the help of hyper-diffusion, this instability can be maintained and no accuracy seems to be lost in the long run. The second disadvantage is the convergence time. There is a significant period of time in which the accuracy of the wave method is up to five times worse than the implicit and explicit methods. In realistic simulations, it is important to be aware of this convergence time when analysing the heat flux in the early iterations. This is usually not a problem since most simulations will be run for a set period of time before analysis starts, in order that the simulation be completely self-driven and any artificial phenomena from the initial conditions, of which the waviness of this method could be attributed to, has reduced to minimum levels. We show that a limit on the convergence speed τ must be implemented for stability, but that has no effect on the accuracy of the method, and negligible effect on the long-term convergence time. We suggest a limit of τ = 2.0, slightly higher than the minimum possible value of τ, τ_min = 0.7, in order to avoid unnecessary ringing, but found no reason to increase it further to the extent of Rempel (2017) and Navarro et al. (2022), who both use 4. Our limit is still higher than that of Warnecke & Bingert (2020), where τ is intrinsically related to the Alfvén speed. Nevertheless, our implementation of the wave method proves to be suitably dynamic and robust, and outperforms both alternative methods in the long term. It is important to consider the sensitivity of this method to the initial time-step, as considerable changes in accuracy may be possible with small alterations to this choice.

The implicit method provides an important solution for those needing more accurate medium-range heat flux values. This is the best option for analysis in the region where the wave method is not near convergence, and the explicit method would take too long with short time-steps to calculate. The larger machinery of the implicit method still has a time cost, but this is mostly counteracted by the larger time-steps than are possible for the explicit method, leaving the difference between computation times in the wave and implicit methods as manageable for low numbers of iterations. The simplicity of these tests allowed the multi-grid levels to be set with only a small number of repetitions each and still converge. For more complex simulations, it is possible that the repetition count would need to increase, which would in turn increase the computation time and cause the disparity between the wave and implicit methods to grow. With this in mind, the wave method becomes the most desirable choice for most simulations, with opportunity for the development of a new method that could be used alongside it for analysis in the convergence period.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska- Curie [grant agreement no. 945371], and is also supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622.

References

All Tables

All Figures

	Fig. 1 Diffusion of the wave method with and without saturation for Δt = 10−6. The explicit solution for Δt = 10−7 is indicated by dashed lines in the final two panels. The solid lines indicate evolution at times t = 0.0, 0.125, 0.25, 0.5, 1.0, and 2.0. The dotted line references the analytical solution for convergence.
In the text

	Fig. 2 Evolution of the self-similar diffusion test for Gaussian and selfsimilar initial conditions. At t = 1.40 s, the profile lies in the χ t >> 1 regime and both initial profiles show a self-similar diffused profile, albeit for different ϕ0 . The solutions are run with the explicit method for Δt = 10−5 s. We note that different scales are used for the colour maps between t = 0 s and t = 1.4 s.
In the text

	Fig. 3 Diffusion of the wave method with and without saturation for non-constant Δt. The explicit solution is indicated by dashed lines.
In the text

	Fig. 4 Evolution of Δt in time for the adaptive time-stepping wave method in the 1D test, starting at Δt = 1e − 6.
In the text

	Fig. 5 Evolution of the explicit solution for Δt = 10−5 s (solid), compared to the first- order analytical approximation (dashed). Here, the solution for every 0.4 seconds until 3.2 seconds is shown with lines from dark to light. The relative error of the Gaussian profile quickly becomes comparable to that of the self-similar profile for solutions in the χt >> 1 regime.
In the text

	Fig. 6 Average relative error around the peak (T > 105 K) for various set-ups.
In the text

	Fig. 7 Average relative error around the peak (T > 105 K) of variable time-stepping in the wave method, using hyper-diffusive terms (dashed).
In the text

	Fig. 8 Long-term error convergence for the wave and implicit methods, with Δt = 5 × 10−4 seconds.
In the text

	Fig. 9 Areas of application to the limit on τ within the wave method for various limiting values. Solutions are shown for every 0.8 seconds.
In the text

	Fig. 10 Application of a limit on the heat flux in the explicit method. Solutions are shown every 0.8 seconds.
In the text

	Fig. 11 Δt in the adaptive time-stepping version of each method, with a starting Δt of 10−5 seconds.
In the text