Context Thermal conductivity provides important contributions to the energy evolution of the upper solar atmosphere behaving as a nonlinear concentrationdependent diffusion equation. Recently different methods have been offered as bestfit solutions to these problems in specific situations but their effectiveness and limitations are rarely discussed. Aims. We rigorously test the different implementations of solving the conductivity flux in the massivelyparallel magnetohydrodynamics code Bifrost with the aim of specifying the best scenarios for the use of each method. Methods We compare the differences and limitations of explicit versus implicit methods and analyse the convergence of a hyperbolic approximation. Among the tests we use a newly derived 1storder selfsimilar approximation to compare the efficacy of each method analytically in a 1D purethermal test scenario. Results We find that although the hyperbolic approximation proves the most accurate and the fastest to compute in longrunning simulations there is no optimal method to calculate the midterm conductivity with both accuracy and efficiency. We also find that the solution of this approximation is sensitive to the initial conditions and can lead to faster convergence if used correctly. Hyperdiffusivity is particularly useful in aiding the methods to perform optimally. Conclusions We discuss recommendations for the use of each method within more complex simulations whilst acknowledging the areas of opportunity for new methods to be developed.
['Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Mathematical Physics']
Implementation of thermal conduction energy transfer models in the Bifrost Solar atmosphere MHD code
{'Implementation of thermal conduction energy transfer models in the Bifrost solar atmosphere MHD code': 'George Cherry ⋆ , Mikolaj Szydlarski, and Boris Gudiksen  \n- 1 Rosseland Centre for Solar physics,Universitetet i Oslo, Sem Sælands vei 13, 0371,Oslo, Norway\n- 2 Institutt for Teoretisk Astrofysikk, Universitetet i Oslo, Sem Sælands vei 13,0371,Oslo, Norway  \nNovember 4, 2024', 'ABSTRACT': 'Context. Thermal conductivity provides important contributions to the energy evolution of the upper solar atmosphere, behaving as a non-linear concentration-dependent di ff usion equation. Recently, di ff erent methods have been o ff ered as best-fit solutions to these problems in specific situations, but their e ff ectiveness and limitations are rarely discussed.  \nAims. We have rigorously tested the di ff erent implementations of solving the conductivity flux, in the massively parallel magnetohydrodynamics code, Bifrost, with the aim of specifying the best scenarios for the use of each method.  \nMethods. We compared the di ff erences and limitations of explicit versus implicit methods, and analyse the convergence of a hyperbolic approximation. Among the tests, we used a newly derived first-order self-similar approximation to compare the e ffi cacy of each method analytically in a 1D pure-thermal test scenario.  \nResults. We find that although the hyperbolic approximation proves the most accurate and the fastest to compute in long-running simulations, there is no optimal method to calculate the mid-term conductivity with both accuracy and e ffi ciency. We also find that the solution of this approximation is sensitive to the initial conditions, and can lead to faster convergence if used correctly. Hyperdi ff usivity is particularly useful in aiding the methods to perform optimally.  \nConclusions. We discuss recommendations for the use of each method within more complex simulations, whilst acknowledging the areas where there are opportunities for new methods to be developed.  \nKey words. Conduction - Magnetohydrodynamics (MHD) - Methods: numerical - Methods: analytical - Sun: corona - Sun: chromosphere', '1. Introduction': "Our closest star continues to o ff er 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.  \nTo 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  \nequation,  \n∂ e ∂ t = -∇ · e u -p ∇ · u + Q , (1)  \nwhere Q represents the divergence of the heat flux from thermal conduction.  \nSpitzer & Seeger (1963) showed that the heat flux behaves as a non-linear di ff usion equation in the direction of the magnetic field, such that  \nQ = -∇ ∥ · ( q ) , (2) q = -κ 0 T 5 / 2 ∇ ∥ T , (3)  \nwhere κ 0 ∼ 10 -6 erg cm -1 s -1 K -7 / 2 is the Spitzer conduction coe ffi cient. 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 e ff ect 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 di ff usion within Bifrost and so these components are not currently implemented (Gudiksen et al. 2011).  \nWhen discretised explicitly using a numerical scheme, for example as  \n∂ e [ i , n ] ∂ t = e [ i , n ] -e [ i , n -1] ∆ t , (4)  \n∇ q [ i , n ] = q [ i , n ] -q [ i -1 , n ] ∆ x , (5)  \nin the 1D case, the second-order di ff usive operator limits the Courant condition to scale with ∆ x 2 . 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 di ffi cult to manage and requires large artificial di ff usive parameters to maintain it.  \nBifrost employs two other methods of solving thermal conduction. The first uses the Crank-Nicholson implicit method described in Gudiksen et al. (2011),  \ne [ n + 1] -e ∗ ∆ t = θ ∇ ∥ · q [ n + 1] + (1 -θ ) ∇ ∥ · q [ n ] . (6)  \nThe quantities at time-step n are calculated before the MHD time-step, starred quantities are produced after the MHD time-step, 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, o ff ering larger time-steps with stability, but at the cost of more computing time due to the complex machinery of the method.  \nThe second method implements the hyperbolic PDE approximation proposed by Rempel (2017),  \n∂ q ∂ t = 1 τ GLYPH<16> -q -κ 0 T 5 / 2 ∇ ∥ T GLYPH<17> , (7)  \n∂ e ∂ t = ... - ∇ · q , (8)  \nwhich 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 = √ κ/τ , 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.  \nIt 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 ∂ e ∂ t before the next time-step. For the implicit method, the final internal energy, e , is overwritten to include the e ff ect 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 a ff ect 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 ∼ ∆ t ∆ x 2 . We can split this into three parts for each derivative calculated; two spatial ( ∇ T , ∇· q ) and one temporal ( ∂ e ∂ t ). 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  \nq = κ ( T ) ∇ ∥ log T = κ 0 T 7 / 2 ˆ b ( ˆ b · ∇ ) log T , (9)  \nwhere ˆ 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, ∇ 2 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  \nκ ( T ) = κ 0 T 7 / 2 < β ∆ x 2 ∆ t ! e , (10)  \nfor β ∈ [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 non-linear 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 di ff use 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.  \nThe 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 = 2 h , uses the full approximation scheme (see Brandt & Livne (2011) for a detailed explanation):  \n- 1. At level h , the solution to the implicit non-linear problem Lu = f , where u is unknown, is denoted as L h u h = f h .  \n2. After a few Gauss-Seidel relaxation sweeps, the smoothed  \nh  \n- 4. The residual, r h = L h (˜ u h + v h ) -L h ˜ u h , is injected into the coarser grid level, H .\n- solution is approximated as ˜ u . 3. A correction term, v h , is needed, such that v h + ˜ u h = u h .\n- 5. Using an interpolation scheme, I H h , the following values can be interpolated to the new resolution:\n- -u → u H = I H h ˜ u h + v h ,\n- -f → f H = L H ( I H h ˜ u H ) + ˆ I H h ˜ u h , where I and ˆ I pertain to separate interpolation processes,  \nsuch that the linear problem becomes L H u H = f H .  \n- h 7. The coarse correction term is injected back into the finer grid scheme through interpolation, and used to update the finer solution approximation, ˜ u h NEW = ˜ u h + ˆ I h H v H .\n- 6. An approximation, ˜ u H , is calculated for the above using Gauss-Seidel sweeps. The new correction term can be approximated by v H = ˜ u H -I H ˜ u h .  \nBetween 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 = q ˆ b , proposed by Rempel (2017),  \n∂ q ∂ t = 1 τ -q ∗ f sat -q ! , (11)  \nand where q ∗ = κ 0 T 5 / 2 ( ˆ b ·∇ ) T recovers the original Spitzer equation. Rempel discusses appropriate values for the saturation rate , f sat, and convergence rate, τ , as  \nf sat = 1 + | q ∗ | 1 . 5 ρ C 3 S , (12)  \nτ = f CFL ∆ x min ∆ t - | v | ! -2 f sat κ 0 T 7 / 2 e ,  \n- (13)  \nwhere CS = p γ p /ρ 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 e ff ect in these tests, since we have ignored dynamics (e.g spatially varying velocities) and because the dominating term is proportional ( ∆ t ) -2 .  \nThis 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:  \n- 1. Add the flux calculated at the previous time to the energy equation,  \ne [ n + 1] = e [ n ] + ∆ t ( ... - ∇ · ( q [ n ])) . (14)  \n- 2. Calculate the flux to be used at the new time using equation 11,  \nq [ n + 1] = q [ n ] + ∆ t " 1 τ -q ∗ [ n ] f sat -q [ n ] !# . (15)  \nTherefore, 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 ) 2 , 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 τ ,  \nτ min = α ∆ t , (16)  \nwhere 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 e ffi cacy 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  \nFig. 1. Di ff usion 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.  \n<!-- image -->  \ntest. The second test calculates the di ff usion 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. 1and 2, respectively.', '3.1. 1D convergence test': 'We considered the 1D non-linear heat di ff usion equation,  \n∂ T ∂ t = κ 0 ∂ ∂ x T 5 / 2 ∂ T ∂ x ! , (17)  \nwith an initial temperature profile of  \nT 0( x ) = 0 . 1 + 0 . 9 x 5 , (18)  \nand fixed boundaries such that T (0) = 0 . 1 and T (1) = 1 for all t . For t >> 0, the temperature solution converges ( ∂ T ∂ t = 0) to the asymptotic stationary solution  \nT = h 0 . 1 7 / 2 + (1 -0 . 1 7 / 2 ) x i 2 / 7 , (19)  \nwhere κ 0 controls the speed of convergence.  \nOur domain contained 256 grid points, extending across x = [0 , 1]. We normalised the set-up to unity, providing a constant 1D magnetic field, Bx = 1 . 0, and density, ρ = 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 ,  \n∂ T ∂ t = -∇ ∥ · q ( T ) , (20)  \nwhereas our implementation was directly solved in the energy equation, and so pertained to a time evolution in e (see equation 1),  \n∂ e ∂ t = ρ γ -1 ∂ T ∂ t = -∇ ∥ · q ( T ) , (21)  \nin the ideal gas and normalised regime. Therefore, given a constant density, the speed, κ e , at which the test within Bifrost  \nArticle number, page 4 of 11  \nconverged to the solution (19) is related to the Spitzer coe ffi cient κ 0 by  \nκ e = ρ γ -1 κ 0 , (22)  \nwhere we have normalised the Boltzmann constant, kB = , ion mass, mH , and mean molecular weight, µ , to 1.0. Furthermore, we set γ = 5 / 3 and κ 0 = 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 di ff usion equation with concentration-dependent coe ffi cients 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 di ff usion is only along the magnetic field, this leads to the following behaviour:  \n- -Evolution of a point source with initial central value T 0 and width R 0 at the foot-points  \nT ( r , t ) =            T 0 (1 + χ t ) -2 / 9 1 -r 2 R ( t ) 2 ! 2 / 5 , if r < R ( t ) 0 , otherwise . (23)  \nR ( t ) = R 0 (1 + χ t ) 2 / 9 , (24)  \nχ = 9 5 2 κ 0 T 5 / 2 0 R 2 0 , (25)  \nwhere T 0 and R 0 are determined by the total surface integral ϕ 0 = T 0 R 0 B for the Beta function B = B ( 1 2 , 7 2 ).  \n- -Convergence when χ t >> 1 of any distribution with total surface integral ϕ 0  \nFig. 2. Evolution of the self-similar di ff usion 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 di ff used profile, albeit for di ff erent ϕ 0. The solutions are run with the explicit method for ∆ t = 10 -5 s. We note that di ff erent scales are used for the colour maps between t = 0 s and t = 1 . 4 s.  \n<!-- image -->  \nT ( r < R ( t ) , t ) → T 0 5 9 ( B ϕ 0) 2 2 κ 0 t ! 2 / 9 1 -r 2 R ( t ) 2 ! 2 / 5 (26)  \nR ( t ) → 9 5 κ 0( B ϕ 0) 5 / 2 t ! 2 / 9 . (27)  \nEquations 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 ϕ 0 as the given instantaneous point source, and becomes zero in a finite domain, will in fact converge on the same shape.  \nThis solution has been used already in Bifrost to test and validate the implementation of ambipolar di ff usion (MorenoInsertis 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  \nT ∞ << T 0 (Furuseth et al. 2024):  \nmax[ T (1) Tot ( t )] =                      T ∞ + T 0 GLYPH<20> 1 + χ t GLYPH<16> 1 + T ∞ T 0 GLYPH<17> 5 / 2 GLYPH<21> -2 / 9 , if χ t ≪ 1 T ∞ + T 0 GLYPH<20> 1 + χ t GLYPH<16> ϕ Tot ϕ 0 GLYPH<17> 5 / 2 GLYPH<21> -2 / 9 , if χ t ≳ 1 → T ∞ + GLYPH<18> 5 9 ( B ϕ Tot) 2 2 κ 0 t GLYPH<19> 2 / 9 , if χ t ≫ t . (28)  \nR (1) ( t ) =                    R 0 GLYPH<20> 1 + χ t GLYPH<16> 1 + T ∞ T 0 GLYPH<17> 5 / 2 GLYPH<21> -2 / 9 , if χ t ≪ 1 R 0 GLYPH<20> 1 + χ t GLYPH<16> ϕ Tot ϕ 0 GLYPH<17> 5 / 2 GLYPH<21> 2 / 9 , if χ t ≳ 1 → GLYPH<16> 9 5 2 κ 0( B ϕ Tot) 5 / 2 t GLYPH<17> 2 / 9 , if χ t ≫ 1 . (29)  \nOur 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,  \nT = 9 . 9 × 10 5 exp -x 2 + z 2 0 . 01 2 ! + 1 × 10 4 . (30)  \nIt 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 self-similar solution with the same perturbed area, ϕ 0, as with the converged zeroth-order equations, 26 and 27. In our case, that self-similar profile uses R 0 = 0 . 1 Mm, T 0 = 9 . 9 × 10 5 K, and T ∞ = 1 × 10 4 K  \nIn 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 -15 g cm -2 , following the ideal gas regime, and weak magnetic field in one direction, Bz = 0 . 316 G, Bx ,y = 0 textG . Whereas the magnitude of the magnetic field is arbitrary, the density value a ff ects the speed of di ff usion, as discussed in the previous test, but this time without unity normalisation. Therefore, when comparing to the analytical solution, we modified the analytical κ 0 → κα ,  \nκ a = ( γ -1) µ mh kB ρ κ 0 , (31)  \nwhere kB = 1 . 38 × 10 -16 erg K -1 is the Boltzmann constant, mh = 1 . 67 × 10 -24 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 di ff usion. The boundary conditions for the domain were periodic. Where there was a need to analyse di ff usion over a longer time period, 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.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, some time before  \nFig. 3. Di ff usion of the wave method with and without saturation for non-constant ∆ t . The explicit solution is indicated by dashed lines.  \n<!-- image -->  \nt = 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 time-stepping. 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 1 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 -4 . 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 -4 ) 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 Fig. 11), has a similar maximum time-stepping available as the saturated wave method.  \nFrom 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,  \nRelative error = T numerical -T analytical , T analytical (32)  \nbetween the two initial conditions is su ffi ciently similar for longer timescales. The largest di ff erences 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 self-similar 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.  \nFigure 5 is the result of a stable explicit run, with constant time-step, ∆ t = 10 -5 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 -3 seconds. As  \n<!-- image -->  \nt  \nFig. 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': '5  \nFig. 5. Evolution of the explicit solution for ∆ t = 10 -5 s (solid), compared to the firstorder 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.  \n<!-- image -->  \n∆ t increases, the method becomes increasingly wavy 2 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.  \nFig. 6. Average relative error around the peak ( T > 10 5 K) for various set-ups.  \n<!-- image -->  \nAs with the 1D convergence test, the wave solution lags far behind the analytical solution initially (see Fig. 9), 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 time-consuming. 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 o ff ers 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.  \nFig. 7. Average relative error around the peak ( T > 10 5 K) of variable time-stepping in the wave method, using hyper-di ff usive terms (dashed).  \n<!-- image -->  \nFig. 8. Long-term error convergence for the wave and implicit methods, with ∆ t = 5 × 10 -4 seconds.  \n<!-- image -->  \nAnother 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 build-up 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.  \nThe 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 2 . 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 hyper-di ff usion parameters (see Gudiksen et al. 2011) to combat such situations. With appropriate hyper-di ff usion, the initial instability is su ffi ciently damped and lower starting time-steps, such as ∆ t = 10 -5 s in Fig. 7 can be achieved. Then, the method converges quickly towards zero relative error due to lower τ values initially. Notably, the hyper-di ff usion 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-di ff usion (see Fig. 7, bottom panel). We can see in the top panel of Fig. 7, that the hyperdi ff usive runs starting at ∆ t = 10 -4 and 10 -5 s both converge very well within 7 seconds, with the lower time-step proving the most accurate throughout, whereas the larger starting time-step of 5 × 10 -4 s converges slower due to the overall smaller convergence rate, and larger instabilities.', '4.3. Safety limits': 'In equation 13, we can see the various terms that can a ff ect 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  \nT 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.  \nIn 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,  \nκ max = 0 . 7 ∆ x 2 ∆ t ! e (33)  \nfrom 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 e ff ect on the accuracy. This is since this limit directly e ff ects the di ff usion 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 di ff erence 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 di ff erent 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  \nFig. 9. Areas of application to the limit on τ within the wave method for various limiting values. Solutions are shown for every 0.8 seconds.  \n<!-- image -->  \nFig. 10. Application of a limit on the heat flux in the explicit method. Solutions are shown every 0.8 seconds.  \n<!-- image -->  \nFig. 11. ∆ t in the adaptive time-stepping version of each method, with a starting ∆ t of 10 -5 seconds.  \n<!-- image -->  \nmore 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 time-steps than the implicit method, and so the disparity between computation time will continue to grow for longer runs.  \nThe 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-50 times faster  \nusing 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 hyperdi ff usion 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 -4 s, and the hyperdi ff usive wave test run with initial ∆ t = 1 × 10 -4 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.", '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 e ff ects on the entire magnetohydrodynamics system. It is important to understand how a change in heat flux, such as the wave method's  \nTable 1. Completion statistics at time t = 1 . 6 seconds .  \nNotes. Completion statistics for the numerical methods at t = 1 . 6 seconds with adaptive time-stepping, run on a 2D grid with [256 × 256] grid points. Starred methods have a higher hyper-di ff usivity applied. All times are taken for a simulation run on 16 cores in parallel, and exclude processing of output data.  \nTable 2. Completion statistics at time t = 3 . 0 seconds.  \nNotes. Completion statistics for the numerical methods at t = 3 . 0 seconds, run on a 2D grid with [256 × 256] grid points. Starred methods have a higher resistivity applied. All times are taken for a simulation run on 16 cores in parallel, and exclude processing of output data.  \ninitial waviness, a ff ects 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 di ff use 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.  \nWe have discussed the various methods of calculating thermal conductivity within the parallel MHD solar atmosphere code, Bifrost . Any simulation may be run using three di ff erent 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.  \nThe 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-di ff usion, 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 e ff ect on the accuracy of the method, and negligible e ff ect 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  \nτ 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.  \nThe 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 di ff erence 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.  \nAcknowledgements. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie [grant agreement Nº 945371], and is also supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622.", 'References': 'Braginskii, S. I. 1965, Rev. Plasma Phys, 1, 205  \nBrandt, A. & Livne, O. E. 2011, Multigrid Techniques (Society for Industrial and Applied Mathematics) Carlsson, M., De Pontieu, B., & Hansteen, V. H. 2019, araa, 57, 189 Furuseth, S. V., Cherry, G., & Martínez-Sykora, J. 2024, On Thermal Conduction in the Solar Atmosphere: An Analytical Solution for Nonlinear Di ff usivity without Compact Support Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154 Malagoli, A., Dubey, A., Cattaneo, F., & Levine, D. 1996, in Parallel Computational Fluid Dynamics 1995, ed. A. Ecer, J. Periaux, N. Satdfuka, & S. Taylor (Amsterdam: North-Holland), 553-560 Millar, J. & Bancroft, J. C. 2003, CREWES Research Report, 15, 54.1 Moreno-Insertis, F., Nóbrega-Siverio, D., Priest, E. R., & Hood, A. W. 2022, A&A, 662, A42 Navarro, A., Khomenko, E., Modestov, M., & Vitas, N. 2022, A&A, 663, A96 Pattle, R. E. 1959, Q. J. Mech. Appl. Math., 12, 407 Rempel, M. 2017, ApJ, 834, 10, aDS Bibcode: 2017ApJ...834...10R Spitzer, Lyman, J. & Seeger, R. J. 1963, AJP, 31, 890  \nWarnecke, J. & Bingert, S. 2020, Geophys. Astrophys. Fluid Dyn., 114, 261, aDS Bibcode: 2020GApFD.114..261W'}