Non-Fourier description of heat flux evolution in 3D MHD simulations of the solar corona

ABSTRACT

The hot loop structures in the solar corona can be well modelled by three-dimensional magnetohydrodynamic simulations, where the corona is heated by field line braiding driven at the photosphere. To be able to reproduce the emission comparable to observations, one has to use realistic values for the Spitzer heat conductivity, which puts a large constraint on the time step of these simulations and make them therefore computationally expensive. Here, we present a non-Fourier description of the heat flux evolution, which allows us to speed up the simulations significantly. Together with the semi-relativistic Boris correction, we are able to limit the time step constraint of the Alfvén speed and speed up the simulations even further. We discuss the implementation of these two methods to the Pencil Code  and present their implications on the time step, and the temperature structures, the ohmic heating rate and the emission in simulations of the solar corona. Using a non-Fourier description of the heat flux evolution together with the Boris correction, we can increase the time step of the simulation significantly without moving far away from the reference solution. However, for values of the Alfvén speed limit of 3000$\mathrm{km}/\mathrm{s}$ and below, the simulation moves away from the reference solution and produces much higher temperatures and much structures with stronger emission.

KEYWORDS:

• Magnetohydrodynamics
• solar corona
• sun
• magnetic fields
• coronal heating

1. Introduction

The solar corona can be described as a low β plasma at low densities and high temperatures. With the presence of coronal magnetic fields, this leads to plasma, where the magnetic pressure is higher than the gas pressure. Therefore, the plasma motions are dominated by the magnetic field, and the plasma can organise itself in accordance to the geometry of the magnetic field, e.g. closed loop structures. The hot plasma in the corona emits radiation in extreme UV and X-ray emission, making it observable from space-based telescopes. One of the major open questions concerning the solar corona is its heating mechanism, i.e. why is the solar corona typically more than 100 times hotter than the photosphere. One of the ideas explaining coronal heating is the field-line braiding model by Parker (1972, 1988), in which magnetic energy is released in form of nanoflares. In this model, the magnetic footpoints of the loops are irreversibly moved by the small-scale photospheric motions, get braided in chromosphere and corona, where the reconnecting field lines release magnetic energy through ohmic heating and contribute to the thermal energy budget.

Three-dimensional magnetohydrodynamic simulations modelling the solar corona were able to show that with this nanoflare heating mechanism the basic temperature structure and its dynamics can be reproduced (e.g. Gudiksen and Nordlund 2002, 2005b, Bingert and Peter 2011). These models are able to describe energy transport to the corona consistent with the nanoflare model (Bingert and Peter 2013). These types of simulations are further used to synthesise coronal emission comparable with actual observations of the corona. From these synthesised emissions, one finds that these models are able to reproduce the average Doppler shifts to some extent (Peter et al. 2004, 2006, Hansteen et al. 2010) and the formation of coronal loops, when using a data driven model with an observed photospheric magnetic field (Bourdin et al. 2013, 2014, Warnecke and Peter 2019). Furthermore, these models were used to show that the coronal magnetic field structure is close to a potential field (Gudiksen and Nordlund 2005b, Bingert and Peter 2011, Bourdin et al. 2018), and therefore nearly force free. However, the force-free approximation, broadly used to obtain coronal magnetic field with field extrapolations (for a review, we refer to Wiegelmann 2008), turns out to be not always valid (Peter et al. 2015) and fails to describe complex current structures in coronal loops above emerging active regions (Warnecke et al. 2017). Recently, Rempel (2017) showed that the solar corona can be heated by a small-scale dynamo operating in the near-surface region of the convection zone braiding the magnetic footpoints in the photosphere. Therefore, these types of models are able to reproduce the main properties of the solar corona on the resolved scale (e.g. Peter 2015). One of the most important ingredients is the vertical Poynting flux at the bottom of the corona (e.g. Galsgaard and Nordlund 1996, Bingert and Peter 2011, Bourdin et al. 2015).

Currently there are only a limited number of codes available which are used for this kind of simulations. One of the most used codes to simulate the solar corona is the Bifrost code (Gudiksen et al. 2011), which is based on earlier work of Gudiksen and Nordlund (2002, 2005b, 2005a) and the Stagger code (Galsgaard and Nordlund 1996). In these simulations, the near-surface convection is self-consistently included and produces realistic photospheric velocities. Furthermore, the Bifrost code includes a realistic treatment of the chromosphere using a non-local thermal equilibrium description. Another code is the MuRAM code (Vögler et al. 2010, Rempel 2014) that has been recently extended to the upper atmosphere (Rempel 2017). Also, there, the photospheric motions are driven by near-surface convection. Apart of these codes there are other codes used for realistic modelling of the solar corona (e.g. Mok et al. 2005, 2008, Abbett 2007, van der Holst et al. 2014).

In this paper, we present an extension to the coronal model of the Pencil Code¹ that has been used successfully to describe the solar corona using either observed magnetograms and a velocity driver mimicking the photospheric motions (Bingert and Peter 2011, 2013, Bourdin et al. 2013) or flux emergence simulations (Chen et al. 2014, 2015) as input at the lower boundary instead of simulating the near-surface convection. However, Chatterjee (2018) developed a 2D model, where the near-surface convection is included with a realistic treatment of the solar corona. Simplified two-layer simulations of the convection zone and the corona of the Sun and stars using the Pencil Code  have been successfully used to investigate the dynamo-corona interplay (Warnecke and Brandenburg 2014, Warnecke et al. 2016a), to self-consistently drive current helicity ejection into the corona (Warnecke and Brandenburg 2010, Warnecke et al. 2011, 2012a, 2013a) and the formation of sunspot-like flux concentrations (Warnecke et al. 2013b, 2016b, Losada et al. 2019).

To be able to compare the simulations of the solar corona with observations of emissivities, one needs to use a realistic value of the Spitzer heat conductivity. However, this puts a major constraint on the time step in these simulations. For simulations with a grid spacing of around 200 km the time step due to the Spitzer heat conductivity is around 1 ms. However, this can be significant lower, if one does not limit the diffusion speed by the speed of light. If one wants to study the dynamics on smaller scales and being able to reduce the fluid and magnetic diffusivities, one needs to use a higher resolution. The smaller grid spacing leads to even lower values of the time step. As the time step decreases quadratically with the grid spacing, the simulations become unfeasible for very high resolutions. To circumvent this, Chen et al. (2014), for example, used a sub-stepping scheme and Rempel (2017) used a non-Fourier scheme, where the hyperbolic equation for the heat transport is solved. Similar approaches have also been used in the dynamo community to describe the non-local evolution of the turbulent electromagnetic force (Brandenburg et al. 2004, Hubbard and Brandenburg 2009, Rheinhardt and Brandenburg 2012, Brandenburg and Chatterjee 2018). We present here a non-Fourier description of the Spitzer heat flux that has been recently implemented to the Pencil Code, see section 2.2. We compare the outcome of the simulations obtained with and without the non-Fourier scheme, see section 3. Furthermore, we also compare these simulations to those using the semi-relativistic Boris correction (Boris 1970) to the Lorentz force that has been also recently implemented to the Pencil Code (Chatterjee 2018) to limit the time step constraint due to the Alfvén speed, see section 2.3.

2. Setup

The setup of the simulations is based on the model of Bingert and Peter (2011, 2013), therefore a detailed description will not be repeated here. We model a part of the solar corona in a Cartesian box (x,y,z) of $100\times 100\times 60$  Mm${}^3$ using a uniform grid. The z = 0 layer represents the solar photosphere. We use $128\times 128\times 256$ grid points, corresponding to a resolution of $781\mathrm{km}$ in the horizontal $234\mathrm{km}$ in the vertical direction. We solve the compressible magnetohydrodynamic equations for the density ρ, the velocity $\mathit{u}$, the magnetic vector potential $\mathit{A}$ and the temperature T. (1) $\begin{array}{r}\frac{\mathrm{D}\ln \rho }{\mathrm{D}t}=-\mathbf{\nabla }\mathbf{\cdot }\mathit{u},\end{array}$(1) (2) $\begin{array}{r}\frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}=-\frac{\mathbf{\nabla }p}{\rho }+\mathit{g}+\frac{\mathit{J}\times \mathit{B}}{\rho }+\frac{1}{\rho }\mathbf{\nabla }\mathbf{\cdot }2\nu \rho \mathsf{\text{S}},\end{array}$(2) (3) $\begin{array}{r}\frac{\mathrm{D}\ln T}{\mathrm{D}t}+(\gamma -1)\mathbf{\nabla }\mathbf{\cdot }\mathit{u}=\frac{1}{c_{\mathrm{V}}\rho T}[{\mu }_0\eta {\mathit{J}}^2+2\rho \nu {\mathsf{\text{S}}}^2-\mathbf{\nabla }\mathbf{\cdot }\mathit{q}+L],\end{array}$(3) where we use a constant gravity $\mathit{g}=(0,0,-g)$ with g = 274 m/s${}^2$, a rate of strain tensor $\mathsf{\text{S}}=1/2(u_{i,j}+u_{j,i})-1/3{\delta }_{ij}\mathbf{\nabla }\mathbf{\cdot }\mathit{u}$ and a constant viscosity ν throughout the domain. Additionally we use a shock viscosity to resolve shocks formed by high Mach number flows (see Haugen et al. 2004, Gent et al. 2013 for details regarding its implementation). The pressure $p=(k_{\mathrm{B}}/\mu {\mathrm{m}}_{\mathrm{p}})\rho T$ is given by the equation of state of an ideal gas, where $k_{\mathrm{B}}$, μ and $m_{\mathrm{p}}$ are the Boltzmann constant, the molecular weight and the proton mass, respectively. The corresponding adiabatic index $\gamma =c_{\mathrm{P}}/c_{\mathrm{V}}$ is 5/3 for a fully ionised gas, with the specific heats at constant pressure $c_{\mathrm{P}}$ and constant volume $c_{\mathrm{V}}$. The heat flux $\mathit{q}$ is given by anisotropic Spitzer heat conduction (4) $\mathit{q}=-K_0{\left(\frac{T}{[K]}\right)}^{5/2}\frac{\mathit{B}\mathit{B}}{{\mathit{B}}^2}\mathbf{\nabla }T\equiv -\mathcal{K}\mathbf{\nabla }T,$(4) which only gives a contribution aligned with the magnetic field and $K_0=2\times {10}^{-11}$ W(mK)${}^{-1}$ is the value derived by Spitzer (1962) assuming a constant Coulomb logarithm. In general, the Coulomb logarithm and therefore $K_0$ depends weakly on the coronal plasma density. We limit the heat conductivity tensor such that the corresponding heat diffusion speed $\mathrm{d}x/(|\mathcal{K}|/\rho c_{\mathrm{P}})$ is 10% of the speed of light with $\mathrm{d}x$ being the grid spacing. For some of the runs, we replaced this equation by the hyperbolic equation of the non-Fourier heat flux, see section 2.2. Additionally to the anisotropic Spitzer heat conduction, we apply an isotropic numerical heat conduction, which is proportional to $|\mathbf{\nabla }\ln T|$ and a heat conduction with a constant heat diffusivity $\chi =K/c_{\mathrm{P}}\rho$. These additions are used to describe the heat flux in the lower part of the simulation, where the temperature is significantly lower and therefore the Spitzer heat conductivity is significantly smaller than in the corona. It also makes the simulation numerically more stable.

The radiative losses due to the optically thin part of the atmosphere are described by $L=-n_{\mathrm{e}}{n}_{\mathrm{H}}Q(T)$, where $n_{\mathrm{e}}$ and $n_{\mathrm{H}}$ are the electron and hydrogen particle densities. $Q(T)$ describes the radiative losses as a function of temperature following the model of Cook et al. (1989), for details see Bingert (2009).

To fulfil the exact solenoidality of the magnetic field $\mathit{B}=\mathbf{\nabla }\times \mathit{A}$ at all times, we solve for the induction equation in terms of the vector potential $\mathit{A}$. (5) $\frac{\mathrm{\partial }\mathit{A}}{\mathrm{\partial }t}=\mathit{u}\times \mathit{B}+\eta {\mathrm{\nabla }}^2\mathit{A},$(5) where we use the resistive gauge, i.e. arbitrary scalar field φ, which divergence can be added to the induction equation is chosen to be $\varphi =\eta \mathbf{\nabla }\mathbf{\cdot }\mathit{A}$. The currents are given by $\mathit{J}=\mathbf{\nabla }\times \mathit{B}$ and η is the magnetic diffusivity.

2.1. Initial and boundary conditions

At the lower boundary, we use for the vertical magnetic field the line-of-sight magnetic field from the active region AR 11102, observed on the 30th of August with the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) onboard of the Solar Dynamics Observatory (SDO), see figure 1 for an illustration. As an initial condition, we use a potential field extrapolation to fill the whole box with magnetic fields. For the temperature, we use an initial profile of a simplified representation of the solar atmosphere, similar as in Bingert and Peter (2011). The density is calculated accordingly using hydrostatic equilibrium. The velocities are initially set to zero.

Figure 1. Initial vertical magnetic field $B_z$ at the photospheric layer z = 0 (colour online).

The simulations are driven by a prescribed horizontal velocity field at the lower boundary mimicking the pattern of surface convection. As discussed in Gudiksen and Nordlund (2002, 2005b), Bingert (2009) and Bingert and Peter (2011), such a surface velocity driver is able to reproduce the observed photospheric velocity spectrum in space and time. To avoid the destruction of the magnetic field pattern caused by the photospheric velocities, we apply the following to stabilise the field: (i) we lower the magnetic diffusivity in the two lowest grid layers by a factor of 800 using cubic step function, (ii) we apply a quenching of velocities by a factor of 2, when magnetic pressure is larger than the gas pressure and (iii) we interpolate between the current vertical magnetic field and the initial one $B_z^{\mathrm{int}}$ at z = 0 layer following (6) $\frac{\mathrm{\partial }{B}_z}{\mathrm{\partial }t}=\frac{1}{{\tau }_b}(B_z^{\mathrm{int}}-B_z),$(6) where ${\tau }_b=10$ min is the relaxation time. The quenching of photospheric velocities mimics the suppression of convection in magnetised regions as observed on the solar surface (see detailed discussion in Gudiksen and Nordlund 2002, 2005b, Bingert 2009, Bingert and Peter 2011). We apply a potential field boundary condition at the bottom and top boundary of box for the magnetic field. The temperature and density are kept fix at the bottom boundary. The temperature is kept constant and the heat flux is set to zero at the top boundary allowing the temperature to vary in time. At the top boundary, we set all velocity components to zero to prevent mass leaving or entering the simulation and to suppress all flows near the top boundary. The density in the lower part is high enough to serve as a mass reservoir. All quantities are periodic in horizontal directions.

For the viscosity, we choose $\nu ={10}^{10}{\mathrm{m}}^2/\mathrm{s}$ similar to the Spitzer value for typical coronal temperatures and densities. We set $\eta =2\times {10}^{10}{\mathrm{m}}^2/\mathrm{s}$ motivated by the numerical stability of the simulations. In the solar corona, the magnetic Prandtl number $P{r}_{\mathrm{M}}=\nu /\eta$ is around ${10}^{10}-{10}^{12}$ and not 0.5 as in our simulations.

2.2. Non-Fourier heat flux scheme

To reduce the time step constraints due to the Spitzer heat conductivity, we use a non-Fourier description and solve for the heat flux $\mathit{q}$ (7) $\frac{\mathrm{\partial }\mathit{q}}{\mathrm{\partial }t}=-\frac{1}{{\tau }_{\mathrm{Spitzer}}}(\mathit{q}+\mathcal{K}\mathbf{\nabla }T),$(7) where ${\tau }_{\mathrm{Spitzer}}$ is the heat flux relaxation time, i.e. e-folding time for $\mathit{q}$ to approach $-\mathcal{K}\mathbf{\nabla }T$. $\mathcal{K}$ is the Spitzer heat conductivity tensor, which has contributions only along the magnetic field. This approach enables us to use a different time stepping constrain to solve our equations. Instead of using the time step of Spitzer heat conduction $\mathrm{d}{t}_{\mathrm{Spitzer}}=\mathrm{d}{x}^2/\gamma {\chi }_{\mathrm{Spitzer}}$ with ${\chi }_{\mathrm{Spitzer}}=|\mathcal{K}|/\rho c_{\mathrm{P}}$, we find two new time step constraints (8) $\mathrm{d}{t}_1=\mathrm{d}x\sqrt{\left(\frac{{\tau }_{\mathrm{Spitzer}}}{\gamma {\chi }_{\mathrm{Spitzer}}}\right)}\equiv \frac{\mathrm{d}x}{c_{\mathrm{Spitzer}}}\mathrm{and}\mathrm{d}{t}_2={\tau }_{\mathrm{Spitzer}},$(8) where $\mathrm{d}{t}_1$ comes from the wave propagation speed $c_{\mathrm{Spitzer}}$. To see this more clearly, we can rewrite (7(7) $\frac{\mathrm{\partial }\mathit{q}}{\mathrm{\partial }t}=-\frac{1}{{\tau }_{\mathrm{Spitzer}}}(\mathit{q}+\mathcal{K}\mathbf{\nabla }T),$(7) ) in one dimension x, with q and K being the one-dimensional counterparts of $\mathit{q}$ and $\mathcal{K}$ as (9) $\frac{\mathrm{\partial }q}{\mathrm{\partial }t}=-\frac{1}{{\tau }_{\mathrm{Spitzer}}}(q+K\frac{\mathrm{\partial }T}{\mathrm{\partial }x}).$(9) Then together with a simplified one-dimensional version of (3(3) $\begin{array}{r}\frac{\mathrm{D}\ln T}{\mathrm{D}t}+(\gamma -1)\mathbf{\nabla }\mathbf{\cdot }\mathit{u}=\frac{1}{c_{\mathrm{V}}\rho T}[{\mu }_0\eta {\mathit{J}}^2+2\rho \nu {\mathsf{\text{S}}}^2-\mathbf{\nabla }\mathbf{\cdot }\mathit{q}+L],\end{array}$(3) ), where we only consider the heat flux term (10) $\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=-\frac{1}{c_{\mathrm{V}}\rho }\frac{\mathrm{\partial }q}{\mathrm{\partial }x},$(10) we can construct a wave equation for the temperature, namely (11) $\frac{{\mathrm{\partial }}^{2}T}{\mathrm{\partial }{t}^2}=-\frac{1}{{\tau }_{\mathrm{Spitzer}}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}+\frac{\gamma {\chi }_{\mathrm{Spitzer}}}{{\tau }_{\mathrm{Spitzer}}}\frac{{\mathrm{\partial }}^{2}T}{\mathrm{\partial }{x}^2},$(11) in which $c_{\mathrm{Spitzer}}$=$\sqrt{\gamma {\chi }_{\mathrm{Spitzer}}/{\tau }_{\mathrm{Spitzer}}}$ emerges as the propagation speed. The two new time step constraints emerge from the pre-factors of the terms on the right-hand side.

By certain choices of ${\tau }_{\mathrm{Spitzer}}$, we can significantly increase the time step. Furthermore, because $\mathrm{d}{t}_1$ depends linear on the grid spacing $\mathrm{d}x$, instead of quadric as $\mathrm{d}{t}_{\mathrm{Spitzer}}$, the speed-up ratio $\mathrm{d}{t}_1/\mathrm{d}{t}_{\mathrm{Spitzer}}$ grows with higher resolutions, which leads to a computational gain. Both time step constraints are included in the CFL condition to calculate the time step of the simulation. $\mathrm{d}{t}_1$ enters the time step calculation through the advective time $\mathrm{d}{t}_{\mathrm{advec}}$ step using (12) $\mathrm{d}{t}_{\mathrm{advec}}=\frac{\mathrm{d}x}{u_{\mathrm{advec}}}\mathrm{with}{u}_{\mathrm{advec}}=max(|\mathit{u}|+\sqrt{c_{\mathrm{s}}^2+v_{\mathrm{A}}^2}+c_{\mathrm{Spitzer}}),$(12) where $u_{\mathrm{advec}}$ is the advection speed and $c_{\mathrm{s}}$ the sound speed.

The major part of the heat flux is concentrated in the transition region, where the temperature gradient is high. This can lead to strong gradients in the heat flux $\mathit{q}$ itself. We, therefore, normalise $\mathit{q}$ by the density ρ to decrease the heat flux in the lower part of the transition region compared to the upper part. The main motivation is to gain a better numerical stability and be able to resolve stronger gradients in $\mathit{q}$ better. This results in a new set of equations, where (13) $\tilde{\mathit{q}}=\mathit{q}/\rho .$(13) We basically solve now for the energy flux per unit particle instead of the energy flux density. (14) $\frac{\mathrm{\partial }\tilde{\mathit{q}}}{\mathrm{\partial }t}=\frac{1}{\rho }\frac{\mathrm{\partial }\mathit{q}}{\mathrm{\partial }t}-\tilde{\mathit{q}}\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }t}=-\frac{1}{{\tau }_{\mathrm{Spitzer}}}(\tilde{\mathit{q}}+\frac{\mathcal{K}}{\rho }\mathbf{\nabla }T)+\tilde{\mathit{q}}(\mathit{u}\mathbf{\cdot }\mathbf{\nabla }\ln \rho +\mathbf{\nabla }\mathbf{\cdot }\mathit{u}),$(14) where we use the continuity equation to derive the last term. The term in the energy equation changes correspondingly (15) $\frac{\mathrm{\partial }\ln T}{\mathrm{\partial }t}=-\frac{1}{T{c}_{\mathrm{V}}}(\mathbf{\nabla }\mathbf{\cdot }\tilde{\mathit{q}}+\tilde{\mathit{q}}\mathbf{\cdot }\mathbf{\nabla }\ln \rho )+\cdots .$(15) This formulation does not change the time step constraints shown in (8(8) $\mathrm{d}{t}_1=\mathrm{d}x\sqrt{\left(\frac{{\tau }_{\mathrm{Spitzer}}}{\gamma {\chi }_{\mathrm{Spitzer}}}\right)}\equiv \frac{\mathrm{d}x}{c_{\mathrm{Spitzer}}}\mathrm{and}\mathrm{d}{t}_2={\tau }_{\mathrm{Spitzer}},$(8) ).

Instead of choosing ${\tau }_{\mathrm{Spitzer}}$ as a constant value in time and space, we also implemented an auto-adjustment, where ${\tau }_{\mathrm{Spitzer}}$ can vary in space and time. This allows the simulation to be more flexible and to be able to optimise the time step. The main idea to choose a reasonable value for ${\tau }_{\mathrm{Spitzer}}$ is that we set the time scale of the heat diffusion to be the smallest of all relevant time scales in this problem, i.e. the heat diffusion is the fastest process. The next bigger time scale is typically the Alfvén crossing time $\mathrm{d}{t}_{\mathrm{vA}}=\mathrm{d}x/v_{\mathrm{A}}$ with the Alfvén speed $v_{\mathrm{A}}=B/\sqrt{{\mu }_0\rho }$. We want to keep the hierarchy of the time steps of each process in place while lowering the time step as much as possible. So we choose the time step of the heat diffusion to be always a bit lower than the Alfvén time step, therefore the heat diffusion is still the fastest process, but slower as before. For a fixed ratio between the $\mathrm{d}{t}_1$ and $\mathrm{d}{t}_{\mathrm{vA}}$, we “tie” ${\tau }_{\mathrm{Spitzer}}$ to $v_{\mathrm{A}}$ and we set (16) $\mathrm{d}{t}_1=\frac{\mathrm{d}{t}_{\mathrm{vA}}}{\sqrt{2}}\to c_{\mathrm{Spitzer}}=\sqrt{2}{v}_{\mathrm{A}}\to {\tau }_{\mathrm{Spitzer}}=\frac{\gamma {\chi }_{\mathrm{Spitzer}}}{2v_{\mathrm{A}}^2}.$(16)

On one hand, ${\tau }_{\mathrm{Spitzer}}$ would become very small in regions below the corona, because there ${\chi }_{\mathrm{Spitzer}}$ has very low values due to the low temperature and high density values. However, in these regions the heat transport is mainly due to the isotropic heat transport. Low values of ${\tau }_{\mathrm{Spitzer}}$ in these regions would cause a very small time step, even though the Spitzer heat flux is not important for the heat transport in these regions. Therefore, we choose the lower limit to be the advective time step, which assures that ${\tau }_{\mathrm{Spitzer}}$ will not affect the time step in these regions. On the other hand, we want to avoid ${\tau }_{\mathrm{Spitzer}}$ becoming too large and therefore the heat transport getting less efficient, i.e. $\mathit{q}$ is still sufficiently close to $-\mathcal{K}\mathbf{\nabla }T$. So, we choose ${\tau }_{\mathrm{Spitzer}}^{\max }=100$ s as a limit for ${\tau }_{\mathrm{Spitzer}}$: (17) $min(\mathrm{d}{t}_{\mathrm{vA}},\frac{\mathrm{d}x}{\sqrt{c_{\mathrm{s}}^2+{\mathit{u}}^2}})\le {\tau }_{\mathrm{Spitzer}}\le {\tau }_{\mathrm{Spitzer}}^{\max }.$(17)

To use the non-Fourier heat flux description in the Pencil Code, one has to add HEATFLUX=heatflux to src/Makefile.local and set the parameters in name list heatflux_run_pars in run.in. The relaxation time ${\tau }_{\mathrm{Spitzer}}$ can be either chosen freely and the inverse is set by using tau_inv_spitzer or one can switch on the automatically adjustment by using ltau_spitzer_va=T, then tau_inv_spitzer sets the value of 1/${\tau }_{\mathrm{Spitzer}}^{\max }$.

2.3. Semi-relativistic Boris correction

Above an active region, the magnetic field strength can be high while the density is low leading to Alfvén speeds comparable to the speed of light (e.g. Chatterjee and Fan 2013, Rempel 2017). This causes two major issues. On one hand, the MHD approximation assuming non-relativistic phase speeds is not valid anymore, i.e. we cannot neglect the displacement current. On the other hand, the high values of the Alfvén speed reduce the time step significantly. To address these two issues, we use a semi-relativistic correction of the Lorentz force following the work of Boris (1970) and Gombosi et al. (2002), where we apply a semi-relativistic correction term to the Lorentz force. This has been used and successfully tested for the MuRAM code in Rempel (2017). Here, we use the implementation discussed by Chatterjee (2018), who added this correction term to the Pencil Code. There, the Lorentz force transforms to (18) $\frac{\mathit{J}\times \mathit{B}}{\rho }\to {\gamma }_{\mathrm{A}}^2\frac{\mathit{J}\times \mathit{B}}{\rho }+(1-{\gamma }_{\mathrm{A}}^2)(\mathit{I}-{\gamma }_{\mathrm{A}}^2\frac{\mathit{B}\mathit{B}}{{\mathit{B}}^2})(\mathit{u}\mathbf{\cdot }\mathbf{\nabla }\mathit{u}+\frac{\mathbf{\nabla }p}{\rho }-\mathit{g}),$(18) where ${\gamma }_{\mathrm{A}}^2=1/(1+v_{\mathrm{A}}^2/c^2)$ is the relativistic correction factor. We note here that the correction term used here and in Chatterjee (2018) is slightly different from the one used by Rempel (2017), because Chatterjee (2018) finds a more accurate way to approximate the inversion of the enhanced inertia matrix. This leads to an additional ${\gamma }_{\mathrm{A}}^2$ in front of $\mathit{B}\mathit{B}/{\mathit{B}}^2$. If $v_{\mathrm{A}}\ll c$ and ${\gamma }_{\mathrm{A}}^2\approx 1$, we retain the normal Lorentz force expression. For $v_{\mathrm{A}}\le c$, the Lorentz force is reduced and the inertia is reduced in the direction perpendicular to the magnetic field. As the enhance inertia matrix (Rempel 2017) is originally on the right-hand side of the momentum equation, i.e. under the time derivative and it is just approximated by a correction term on the left-hand side, the semi-relativistic Boris correction does not change the stationary solution of the system and therefore does not lead to further correction terms in the energy equation. To switch on the Boris correction in Pencil Code, one sets the flag lboris_correction=T in the name list magnetic_run_pars.

The Boris correction describes the modification of the Lorentz force in the situation, where the Alfvén speed becomes comparable to the speed of light. In other words, the speed of light is a natural Alfvén speed limiter and the Boris correction describes the modification close to this limiter. We can artificially decrease the value of the limiter to a value of our choice and the Boris correction takes care of the corresponding modifications. This can significantly reduce the value of the Alfvén speed in our simulations and allow us to enhance the Alfvén time step. Unlike in Chatterjee (2018), we use the Boris correction indeed to increase the Alfvén time step, similar to what has been done by Rempel (2017). As shown by Gombosi et al. (2002), the propagation speed can be quite complicated, we choose a similar time step modification as in Rempel (2017) (19) $\mathrm{d}{t}_{\mathrm{vA}}\to \mathrm{d}{t}_{\mathrm{vA}}\sqrt{1+{\left(v_{\mathrm{A}}^2/c_{\mathrm{A}}^2\right)}^2},$(19) where $c_{\mathrm{A}}$ is the limiter. We choose for Set B $c_{\mathrm{A}}=\mathrm{10,000}\mathrm{km}/\mathrm{s}$, which corresponds to a time step of $\mathrm{d}{t}_{\mathrm{vA}}\approx 20$ ms for our simulations. The limiter $c_{\mathrm{A}}$ can be set by using va2max_boris in the name list magnetic_run_pars. The Boris correction can be used together with the automatic adjusted relaxing time ${\tau }_{\mathrm{Spitzer}}$ in the non-Fourier heat flux calculation: if one sets va2max_tau_boris in heatflux_run_pars to the same value as va2max_boris in magnetic_run_pars, then the code modifies the Alfvén speed and the Alfvén time step used in (16(16) $\mathrm{d}{t}_1=\frac{\mathrm{d}{t}_{\mathrm{vA}}}{\sqrt{2}}\to c_{\mathrm{Spitzer}}=\sqrt{2}{v}_{\mathrm{A}}\to {\tau }_{\mathrm{Spitzer}}=\frac{\gamma {\chi }_{\mathrm{Spitzer}}}{2v_{\mathrm{A}}^2}.$(16) ) and (17(17) $min(\mathrm{d}{t}_{\mathrm{vA}},\frac{\mathrm{d}x}{\sqrt{c_{\mathrm{s}}^2+{\mathit{u}}^2}})\le {\tau }_{\mathrm{Spitzer}}\le {\tau }_{\mathrm{Spitzer}}^{\max }.$(17) ) accordingly.

3. Results

We present here the results of three sets of runs, where we use different values of the heat flux relaxation time ${\tau }_{\mathrm{Spitzer}}$ in combination with and without the Boris correction. In the first set, containing only Run R, we use the normal treatment of the Spitzer heat flux without using the non-Fourier heat flux evolution and without the Boris correction. In the second set, containing 4 runs (Set H), we use the non-Fourier heat flux evolution with ${\tau }_{\mathrm{Spitzer}}$ between 10 and 1000 ms and the automatically adjustment, see section 2.2. In the third set, containing 7 runs (Set B), we use the semi-relativistic Boris correction with $c_{\mathrm{A}}=\mathrm{10,000}\mathrm{km}/\mathrm{s}$ and the non-Fourier heat flux evolution with ${\tau }_{\mathrm{Spitzer}}=$ 10–1000 ms and the automatically adjustment. We also use one run (Ba2) with even lower Alfvén speed limit of $c_{\mathrm{A}}=3000\mathrm{km}/\mathrm{s}$. An overview of the runs can be found in table 1.

Table 1. Summary of the runs. ${\tau }_{\mathrm{Spitzer}}$ is the relaxation time for non-Fourier heat flux description, see section 2.2, ${\tau }_{\mathrm{Spitzer}}=\mathrm{\infty }$ stands for the use of standard Fourier heat flux, see (4(4) $\mathit{q}=-K_0{\left(\frac{T}{[K]}\right)}^{5/2}\frac{\mathit{B}\mathit{B}}{{\mathit{B}}^2}\mathbf{\nabla }T\equiv -\mathcal{K}\mathbf{\nabla }T,$(4) ).

Runs	${\tau }_{\mathrm{Spitzer}}$ [ms]	$c_{\mathrm{A}}$ [km/s]	$\mathrm{d}t$ [ms]	$\mathrm{d}{t}_{v_{\mathrm{A}}}$ [ms]	$\mathrm{d}{t}_1$ [ms]	$\mathrm{d}{t}_2$ [ms]	$t_{\mathrm{cpu}}$ [μs]	$\mathrm{\Delta }{T}_{\mathrm{cor}}$
R	∞	∞	1.5	2.7	1.7	1.7	4.2$\times {10}^{-2}$	0
H001	10	∞	1.1	2.8	1.2	9.0	4.6$\times {10}^{-2}$	−5%
H005	50	∞	2.4	3.0	4.4	45.0	4.5$\times {10}^{-2}$	−14%
H1	1000	∞	4.5	5.2	13.0	900.0	4.6$\times {10}^{-2}$	−18%
Ha	auto	∞	2.0	2.6	2.5	2.0	4.6$\times {10}^{-2}$	−3%
B001	10	10,000	0.5	21.2	0.5	9.0	4.6$\times {10}^{-2}$	14%
B002	20	10,000	2.8	21.2	2.8	18.0	4.6$\times {10}^{-2}$	−13%
B005	50	10,000	3.4	21.2	3.7	45.0	4.5$\times {10}^{-2}$	−7%
B01	100	10,000	4.4	21.2	4.7	90.0	4.5$\times {10}^{-2}$	4%
B03	300	10,000	5.4	21.2	6.1	270.0	4.5$\times {10}^{-2}$	−0.3%
B1	1000	10,000	8.0	21.2	10.0	900.0	4.5$\times {10}^{-2}$	−4%
Ba	auto	10,000	15.4	21.2	15.0	19.5	4.7$\times {10}^{-2}$	14%
Ba2	auto	3000	47.6	70.6	49.9	64.9	4.7$\times {10}^{-2}$	18%

Notes: $c_{\mathrm{A}}$ is the Alfvén speed limit, used for the Boris correction, see section 2.3; $c_{\mathrm{A}}=\mathrm{\infty }$ stands for no Boris correction. $\mathrm{d}t$ indicates the averaged time step, $\mathrm{d}{t}_{v_{\mathrm{A}}}$ the averaged Alfvén time step and $\mathrm{d}{t}_1$ and $\mathrm{d}{t}_2$ the average time step due to the heat flux evolution, see (8(8) $\mathrm{d}{t}_1=\mathrm{d}x\sqrt{\left(\frac{{\tau }_{\mathrm{Spitzer}}}{\gamma {\chi }_{\mathrm{Spitzer}}}\right)}\equiv \frac{\mathrm{d}x}{c_{\mathrm{Spitzer}}}\mathrm{and}\mathrm{d}{t}_2={\tau }_{\mathrm{Spitzer}},$(8) ). For Run R, $\mathrm{d}{t}_1=\mathrm{d}{t}_2=\mathrm{d}{t}_{\mathrm{Spitzer}}$. All these quantities are determined as an average in the quasi-stationary state. $t_{\mathrm{cpu}}$ is wall clock time per time step per mesh point. For the timing we use the SISU Cray XC40 supercomputing cluster at CSC. $\mathrm{\Delta }{T}_{\mathrm{cor}}=(⟨T{⟩}_{\mathrm{runs}}-⟨T{⟩}_R)/⟨T{⟩}_R$ is the mean temperature deviation from the reference runs, taking as a horizontal and height ($z=$ 20–40 mm) average.

3.1. Time steps

As a first step, we look at the time steps of all the runs in table 1. In Run R, the averaged time step in the saturated stage is around 1.5 ms. This time step is constrained by the Spitzer time step $\mathrm{d}{t}_{\mathrm{Spitzer}}$, which is shown as $\mathrm{d}{t}_1=\mathrm{d}{t}_2$ in table 1. The Alfvén time step $\mathrm{d}{t}_{\mathrm{vA}}$ is around twice as large. In the Set H, the code additionally solves the non-Fourier heat flux equation that leads to an increased time step. However, the time step is actually limited by the low Alfvén time step and therefore the time step cannot be increased by a large factor. In Set H, the largest speed-up factor is around 3. For Run H001, the value of ${\tau }_{\mathrm{Spitzer}}$ is low enough to have a time step constraint of $\mathrm{d}{t}_1$ instead of $\mathrm{d}{t}_{\mathrm{vA}}$. However, the runs reach a lower time step than in Run R. For values of the relaxing time ${\tau }_{\mathrm{Spitzer}}=$ 50–1000 ms (Runs H005 and H1), the time step due to the heat flux is larger than the Alfvén time step. This means that the physical process of heat redistribution is even slower than the Alfvén speed. This leads in Run H1 to higher densities resulting in a lower Alfvén speed and a higher $\mathrm{d}{t}_{\mathrm{vA}}$, see discussion in section 3.4. Furthermore, Run H1 only runs stable, if we increase the shock viscosity to 10 times higher values than in the other runs. This will certainly lead to some additional differences independent of the direct influence of the non-Fourier heat flux description. When applying the auto-adjustment of ${\tau }_{\mathrm{Spitzer}}$ (Run Ha), the time steps $\mathrm{d}{t}_1$ and $\mathrm{d}{t}_2$ are slightly smaller than $\mathrm{d}{t}_{\mathrm{vA}}$ and limits the time step. There the speed up is less than a factor of 2, but the heat distribution is the fastest process in the system. Using the non-Fourier heat flux description leads usually to higher peak temperatures, because the temperature diffusion is less efficient. For the calculation of $\mathrm{d}{t}_1$ and $\mathrm{d}{t}_2$, the code uses the CFL pre-factors of 0.9 for both time steps, this results in $\mathrm{d}{t}_2=0.9{\tau }_{\mathrm{Spitzer}}$. As $\mathrm{d}{t}_1$ enters via (12(12) $\mathrm{d}{t}_{\mathrm{advec}}=\frac{\mathrm{d}x}{u_{\mathrm{advec}}}\mathrm{with}{u}_{\mathrm{advec}}=max(|\mathit{u}|+\sqrt{c_{\mathrm{s}}^2+v_{\mathrm{A}}^2}+c_{\mathrm{Spitzer}}),$(12) ), $\mathrm{d}t$ is often lower than $\mathrm{d}{t}_1$ and $\mathrm{d}{t}_{\mathrm{vA}}$ in our simulations.

To increase the time step further, we use the semi-relativistic Boris correction in all runs of Set B. As shown in table 1, $\mathrm{d}{t}_{\mathrm{vA}}$ significantly increases to 21.2 ms for Runs B001-Ba and to 70.6 ms for Ba2. This leads to a much larger speed-up factor of 10 for Run Ba and more than 30 for Run Ba2. For Run B001 to Run B1 with ${\tau }_{\mathrm{Spitzer}}=$ 10–1000 ms, $\mathrm{d}{t}_1$ is lower than $\mathrm{d}{t}_{\mathrm{vA}}$ and the time step can be significantly reduced, while the heat distribution is the fastest process in the system. For Run B1, we achieve a speed up of more than five, however, we need to use a comparable large value of ${\tau }_{\mathrm{Spitzer}}$, which as discussed in section 3.4 can lead to artefacts. For Runs Ba and Ba2, the auto-adjustment of ${\tau }_{\mathrm{Spitzer}}$ takes care that $\mathrm{d}{t}_1<\mathrm{d}{t}_{\mathrm{vA}}$. As discussed below, Run Ba shows a good agreement with Run R, whereas Run Ba2 tends to produce higher temperatures in the corona.

To get a better understanding of the calculation of the time step, we plot in figure 2 various contributions to the time steps for Run Ba. Without the non-Fourier heat flux description and the Boris correction, the time step is dominated by Alfvén time step $\mathrm{d}{t}_{\mathrm{vA}}$ and the Spitzer time step $\mathrm{d}{t}_{\mathrm{Spitzer}}$. The Boris correction reduces $\mathrm{d}{t}_{\mathrm{vA}}$ to $\mathrm{d}{t}_{\mathrm{vA}}^{\mathrm{Boris}}$ mostly in the regions between 5 and 30  Mm. The auto-adjustment of ${\tau }_{\mathrm{Spitzer}}$ sets $\mathrm{d}{t}_1$ to be always slightly lower than $\mathrm{d}{t}_{\mathrm{vA}}^{\mathrm{Boris}}$. Only below z = 5  Mm, $\mathrm{d}{t}_{\mathrm{vA}}^{\mathrm{Boris}}$ is small, because there the temperature diffusion is dominated by the other heat diffusion mechanism described in section 2. It is clearly visible that the $\mathrm{d}{t}_1$ is significantly higher than $\mathrm{d}{t}_{\mathrm{Spitzer}}$ (green line) and $\mathrm{d}{t}_{\mathrm{vA}}$ (red) without the Boris correction. However, we note here that because of the non-Fourier heat flux description we find higher peak temperatures in the simulation. This results in a decrease of $\mathrm{d}{t}_{\mathrm{Spitzer}}$ in comparison with runs without the non-Fourier heat flux description. In Run R, $\mathrm{d}{t}_{\mathrm{Spitzer}}$ is around 1.6 ms, where in Run Ba, it is around a factor of eight lower. Such a factor can be explained by change in temperature by a factor of 2.3.

Figure 2. Vertical distribution of time step constraints for Run Ba at time $t= 200$ min. We plot the time step due to Spitzer heat conductivity $\mathrm{d}{t}_{\mathrm{Spitzer}}$ (green), due to the heat flux $\mathrm{d}{t}_1$ (black) and $\mathrm{d}{t}_2$ (blue), due to the Alfvén speed $\mathrm{d}{t}_{\mathrm{vA}}$ (red) and reduced Alfvén speed with the Boris correction $\mathrm{d}{t}_{\mathrm{vA}}^{\mathrm{Boris}}$. The horizontal averaged values are shown with a solid line and the minimum values at each height with a dashed line (colour online).

Using the non-Fourier heat flux evolution requires to solve (7(7) $\frac{\mathrm{\partial }\mathit{q}}{\mathrm{\partial }t}=-\frac{1}{{\tau }_{\mathrm{Spitzer}}}(\mathit{q}+\mathcal{K}\mathbf{\nabla }T),$(7) ) or (14(14) $\frac{\mathrm{\partial }\tilde{\mathit{q}}}{\mathrm{\partial }t}=\frac{1}{\rho }\frac{\mathrm{\partial }\mathit{q}}{\mathrm{\partial }t}-\tilde{\mathit{q}}\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }t}=-\frac{1}{{\tau }_{\mathrm{Spitzer}}}(\tilde{\mathit{q}}+\frac{\mathcal{K}}{\rho }\mathbf{\nabla }T)+\tilde{\mathit{q}}(\mathit{u}\mathbf{\cdot }\mathbf{\nabla }\ln \rho +\mathbf{\nabla }\mathbf{\cdot }\mathit{u}),$(14) ) meaning three additional equations. However, the computational extra calculation time is around 10%, which is very small compared to the gain in time step reduction. Using the semi-relativistic Boris-Correction does not seem to increase the computation time significantly. Only if we use the auto-adjustment of ${\tau }_{\mathrm{Spitzer}}$ together with the Boris correction we find an additional 2–3% increase in the computation time, as shown in the last row of table 1.

3.2. Alfvén velocity with Boris correction

Next, we look at the influence of the semi-relativistic Boris correction on the Alfvén velocity $v_{\mathrm{A}}$. In figure 3, we plot 2D histograms of $v_{\mathrm{A}}$ for Runs R, Ba, Ba2. For Run R, the maximum speed reaches $v_{\mathrm{A}}=\mathrm{80,000}\mathrm{km}/\mathrm{s}$ at the lower part of the corona, where the density has decreased significantly with height, but the magnetic field is still strong. The median (yellow line) has its maximum at the same location with a value around $v_{\mathrm{A}}=\mathrm{18,000}\mathrm{km}/\mathrm{s}$. In Run Ba, we have applied the Boris correction with $c_{\mathrm{A}}=\mathrm{10,000}\mathrm{km}/\mathrm{s}$. Even though, this value is lower than the averaged and mean value in the region of z = 5–20  Mm, the velocity distribution does not change significantly in comparison to Run R. As a main effect of the Boris correction, the peak velocity at the top of the distribution is reduced, therefore the distribution becomes more compact. This can be also seen from the changes in the mean and median velocity. While the maximum of the mean is reduced from above $v_{\mathrm{A}}=\mathrm{20,000}\mathrm{km}/\mathrm{s}$ of Run R to nearly $v_{\mathrm{A}}=\mathrm{15,000}\mathrm{km}/\mathrm{s}$, the median changes just slightly. Also, the area between the 25 and 75 percentiles of the Alfvén velocity population moves only slightly towards lower values. This make us confident that the Boris correction with $c_{\mathrm{A}}=\mathrm{10,000}\mathrm{km}/\mathrm{s}$ does only reduce the peak velocities and not the overall velocity structure; most of the points are unaffected by the correction.

Figure 3. 2D histograms of the Alfvén speed over height z for Runs R, Ba, Ba2. We plot the mean value with red solid line and the median with a yellow solid. The dashed white-blue lines show the 25 and 75 percentiles, i.e. half of the data points are in between these lines. The black dashed line indicate the Alfvén speed limit $c_{\mathrm{A}}$ for the Boris correction (colour online).

For Run Ba2, we reduce the Alfvén speed limit to $c_{\mathrm{A}}=3000\mathrm{km}/\mathrm{s}$. This makes the velocity distribution even more compact. The maximum values are significantly reduced to $v_{\mathrm{A}}=\mathrm{35,000}\mathrm{km}/\mathrm{s}$, and the mean and median values are also lower than in Runs R, Ba. However, setting $c_{\mathrm{A}}=3000\mathrm{km}/\mathrm{s}$ does not mean that all the velocities are lower than this value, it can be understood as a significant reduction of the peak velocities and a transfer of the velocity distribution to a much more compact form.

3.3. Structure of temperature and ohmic heating

Next, we look at the horizontal averaged temperature profile over height. Even though the non-Fourier description of the heat flux can lead to higher peak temperatures, the overall temperature structure should remain roughly the same. In figure 4, we plot the horizontal averaged temperature profile over height for the reference Run R in panel (a) and a comparison with the other runs in panel (b). The horizontal averaged temperature structure in Run R shows a typical behaviour of corona above an active region with medium magnetic field strengths. The plasma above z = 10 Mm is heated self-consistently to averaged temperatures of around 1 million kelvin. This temperature profile is very similar to results of earlier work with the Pencil Code  (e.g. Bingert 2009, Bingert and Peter 2011, 2013, Bourdin et al. 2013) and other groups (e.g. Gudiksen and Nordlund 2002, 2005b, 2005a, Gudiksen et al. 2011). When comparing with the temperature profiles of the other runs, we find no large differences. For most of the runs the deviation is not more than 10%. For some runs the largest difference occurs in the transition region, where the temperature has a large gradient. Higher temperature values in this region simply mean a slightly lower transition region and lower values mean a slightly higher transition region. Nearly all runs develop a lower or similar transition region location as in Run R. Only Runs H005, H1, Ba2 develop a higher transition region. This can be explained either by sub-dominance of the heat flux time step (Runs H005, H1) or the too low limit for the Alfvén speed, see discussion below. Only in Runs B001, Ba and Ba2 the plasma is heated to 20% higher temperature in the upper corona in comparison with Run R. For Run B001, this high temperature only occur at the end of the simulation, see figure 5. In these runs, the heat diffusion might be not efficient enough to transport heat to lower layers.

Figure 4. (a) Averaged temperature $⟨T⟩$ as a function of height z for Run R. (b) Ratio of the averaged temperature profile of all runs and Run R $⟨T{⟩}_{\mathrm{runs}}/⟨T{⟩}_R$ as a function of height z. The temperatures are averaged horizontal as well as in time for the last quarter (1 h) of the simulation. The colour of the lines indicates the run names in terms of ${\tau }_{\mathrm{Spitzer}}$, the solid lines are for runs of Sets R and H, and dashed lines for Set B (colour online).

Figure 5. Time evolution of the horizontal averaged temperature $⟨T⟩$ (a) and of the horizontal averaged heating rate $⟨{\mu }_0\eta {\mathit{j}}^2⟩$ (b) at $z=$ 18–22 mm. Colour coding is the same as in figure 4 (colour online).

When we look at the temperature evolution over time, as plotted in figure 5(a), we find that each run shows a large variation in time even though we have averaged horizontally and over 18–20  Mm. This can be explained by the non-linear behaviour of the system. Because of this reason temporal variations occurring in the other runs appear not at the same time for all runs. The difference between the runs is comparable with the time variation of each run. Therefore, to be able to compare the runs, we should look at the time averaged quantities as done throughout this work.

Next, we look at the ohmic heating rate in all the runs. The ohmic heating is the main process in this type of simulations to heat the coronal plasma up to million K. Also, here, we plot the horizontal averaged profile of Run R in panel (a) of figure 6 and compare it with the other runs in figure 6(b). The profile of the ohmic heating rate shows the typical behaviour of an exponential decrease corresponding to two scale heights. Below the corona the scale height is roughly 0.5  Mm, while in the corona the scale height is around 5  Mm. Also, this is consistent with earlier finding with this kind of simulations by many groups (e.g. Gudiksen and Nordlund 2002, 2005b, 2005a, Bingert 2009, Gudiksen et al. 2011, Bingert and Peter 2011, 2013, Bourdin et al. 2013). By comparing with the other sets of runs, we find that these agree well with Run R. Only Run B001 shows a large heating rate in the lower corona, which comes here also from the last part of the simulation. Runs B01, Ba and Ba2 develop a higher heating rate in the upper corona resulting in higher temperatures at this location (see figure 4). Small changes either in the scale height of the coronal heating or in the location in the transition region can explain most of the differences we find in the comparison with Run R. This explains also the temporal changes of the heating rate at constant height, as shown in figure 5(b). The large variations in time of the heating rate can be attributed to non-linear behaviour of the system. Even in Run R, these variations are large compared to the average. Small local changes in temperature and density can also affect the heating rate. As the field is very close to a potential field, the currents are due to small perturbations from the potential field. These perturbations can easily be affected by changes in the plasma flow due to temperature and density fluctuations. Furthermore, in such dynamical non-linear systems, changes for example in the time step can affect also the realisation of the velocity solution. Even when solutions are the same on a statistical level, this can cause variations in the ohmic heating. For these types of models, large variations in time of the ohmic heating rate are a common feature (e.g. Bingert and Peter 2011, 2013) as small changes in local scale height will lead to a large change in the heating rate. Overall, the vertical horizontally averaged temperature and heating structure of all runs agree well with Run R.

Figure 6. (a) Averaged ohmic heating rate $⟨{\mu }_0\eta {\mathit{j}}^2⟩$ as a function of height z for Run R. (b) Ratio of the averaged ohmic heating rate of all runs and Run R $⟨T{⟩}_{\mathrm{runs}}/⟨T{⟩}_R$ as a function of height z. The ohmic heating rate is averaged horizontal as well as in time for the last quarter (1 h) of the simulation. The colour coding the same as in figures 4 and 5 (colour online).

3.4. Emission signatures

To further test how well the non-Fourier description of the heat flux reproduce the Fourier description, we synthesise coronal emissivities corresponding to the 171 Å channel of Atmospheric Imaging Assembly (AIA; Boerner et al. 2012) on board of SDO. We choose this AIA channel because it can be potentially compared with observations and represents well the plasma structure of around 1 million K by convolving the temperature and density structures. This can work as a good test, weather or not coronal emission structures are affected by the choice of heat flux description. For this we calculate the emission following optical thin radiation approximation, (20) $ϵ=n_e^{2}G(T),$(20) where $G(T)$ is the response function of the particular filter, we want to synthesise. Because we compare our simulations among each other and not to observation, we simplify $G(T)$ using a gaussian distribution around a mean temperature ${\log }_{10}{T}_0$, (21) $G(T)\propto \exp [-{\left(\frac{{\log }_{10}T-{\log }_{10}{T}_0}{\mathrm{\Delta }{\log }_{10}{T}_0}\right)}^2],$(21) where $\mathrm{\Delta }{\log }_{10}{T}_0$ is the temperature width used to mimic the temperature response function. We use ${\log }_{10}{T}_0=6{\log }_{10}$ K and $\mathrm{\Delta }{\log }_{10}{T}_0=0.2{\log }_{10}$ K for synthesising the emission of the AIA 171 Å channel. To calculate the emission emitted from a certain direction, we perform an integration along this direction. For the discussion below, we apply an integration along the y and z directions, respectively.

In figure 7, we plot the temperature as a side view (xz) averaged over y and in time (180–240  min) together with the synthesised emission integrated over the y and z directions also averaged in time (180–240  min) representing the AIA 171 channel for Runs R, Ha, Ba, Ba2. For these runs, we expect a good agreement with the reference run R, because the value of ${\tau }_{\mathrm{Spitzer}}$ is regulated automatically and therefore the time step is controlled by the heat flux, i.e. $\mathrm{d}{t}_1$. We find agreement between the Runs R, Ha and Ba, but we find slightly stronger emission structures in Run Ba and slightly hotter temperatures in Run Ha. To illustrate the variation in time we show in figure 8 the time evolution of the emission in a small region of the simulation box. We find a good agreement between Runs R and Ha with variation in time which are comparable with their difference. Run Ba takes a bit longer to saturate, but at around 220 min it also settles to values similar to Runs R and Ha. Run Ba2 seems to saturate to a much higher emission level than the other runs, which is already seen in figure 7.

Figure 7. Temperature and emission structure for Runs R, Ha, Ba, Ba2. We show the temperature averaged over the y direction and in time $t=$ 180–240 min (left panel) together with the synthesised emission comparable to the AIA 171 channel, representing emission at around 1 MK, integrated in the y direction (side view, middle panel) and in z direction (top view, right panel). The emission values represent the count rate of the AIA instrument and has been averaged in time $t=$ 180–240 min. The red square indicates the region which is used to calculate the temporal evolution in figure 8 (colour online).

Figure 8. Time evolution of emission averaged over a small region for Run R, Ha, Ba and Ba2. We plot the emission of the AIA 171 channel in the y direction averaged over a small region ($x=$ 23–27 mm, $z=$ 18–22 mm) as indicated by red boxes in figure 7 middle panels. The inlay shows the time evolution of the averaged emission from 150 to 300 min on a linear scale instead of logarithmic. The colour and style of the lines are the same as in figure 4 (colour online).

For Run Ba2, as pointed out in section 3.3 and shown in first column of figure 7, the corona is heated to higher temperatures, i.e. the heat transport is less efficient. We find larger temperatures mostly at the top of the corona inside the loop structures. This leads also to higher emission in the AIA 171 channel than in the Run R. This might be an artefact from the low limit of the Alfvén speed through the Boris correction in this run. Even though the Alfvén speed limiter does not effect the heat flux directly, it increases the heat flux time step and makes the heat transport less efficient.

In figure 9, we show a few other runs, which are either dominated by the Alfvén time step (Runs H005 and H1) or use a constant value of ${\tau }_{\mathrm{Spitzer}}$ (Runs B03 and B1). Run H005 shows a similar emission structure than the Runs R, Ha, Ba, however, the emission is slightly larger in the legs of the loop. Because also here the temperatures are not significantly higher, the difference is due to the slightly higher density in these regions. For Run H1, ${\tau }_{\mathrm{Spitzer}}$ is large and the time step is controlled by the Alfvén speed instead of the heat flux. This leads to larger temperatures and therefore higher emission. However, also here the density in the corona loops is larger than in Run R, leading not only to higher emission but also to a larger Alfvén time step, see table 1. Furthermore, the high shock viscosity needed to keep the run stable will also have an influence on the solution. In contrast, the time steps in Runs B03 and B1 are controlled by the time step of the heat flux ($\mathrm{d}{t}_1$). There, as expected, we find similar emission loop structures as in Run R, Ha and Ba. They are slightly larger in Run B1 than in Run B03. This means that simulations using either the automatic adjustment or a constant value of ${\tau }_{\mathrm{Spitzer}}$ reproduce the emission structure of Run R well, as long as the time step is still controlled by the heat flux time step $\mathrm{d}{t}_1$, however, the emission tends to be slightly larger. However, for too low values of the Alfvén limiter ($c_{\mathrm{A}}=3000\mathrm{km}/\mathrm{s}$) the emission and temperature become much higher than in Run R. We note here that the AIA 171 channel is relatively broad filter around the mean temperature and therefore hide some of the differences between the runs. A more narrow filters for example used on Hinode/EIS might reveal larger differences.

Figure 9. Temperature and emission structure for Runs H005, H1, B03, B1. The emission values represent the count rate of the AIA instrument. We show the temperature averaged over the y direction and in time (180–240 min, left panel) together with the synthesised emission comparable to the AIA 171 channel, representing emission at around 1 MK, integrated in the y direction (side view, middle panel) and in z direction (top view, right panel). The emission values represent the count rate of the AIA instrument and has been averaged in time (180–240 min) (colour online).

4. Discussion and conclusion

In this work, we present the new implementation of a non-Fourier description of the heat flux to the Pencil Code. We discuss the advantages and the limitations using the example of 3D MHD simulations of the solar corona. The implementation of the auto-adjustment of ${\tau }_{\mathrm{Spitzer}}$ is slightly different from the implementation used in Rempel (2017) in the sense that we ensure the heat flux time step to be always by a square root of two smaller than the Alfvén time step, whereas in Rempel (2017) there is not such a factor. Even though a detailed comparison was not conducted here, we see indications that our choice leads to a better stability of the simulations. We find that using the non-Fourier description of the heat flux alone allows for a small speed up, because in our case the time constraint of the Alfvén speed is large. For simulations with a lower magnetic field strength, we would expect a larger speed up. If we choose a constant ${\tau }_{\mathrm{Spitzer}}$, so that the heat flux time step is four times higher than the Alfvén time step, the temperatures and the emission are significantly larger than in the other runs. This seems to be an artefact of this choice of ${\tau }_{\mathrm{Spitzer}}$.

We further test the implementation of the semi-relativistic Boris correction (Boris 1970) as a limiter for the Alfvén speed. The implementation to the Pencil Code  is slightly different from the one used by Rempel (2017) and Gombosi et al. (2002), see Chatterjee (2018) for details. The Boris correction does not quench the Alfvén speed at all locations to the limit chosen, it actually reduces the peak velocities, which are not very abundant. Therefore, this correction makes the velocity distribution much more compact. The lower the limit, the more compact is the velocity distribution. Using the Boris correction allows for a significant speed up of around 10. For higher speed up, i.e. lower limit for Alfvén speed, the simulation develops higher temperatures and emission signatures than the reference run. The auto-adjustment together with Boris correction works very well to reproduce the temperatures and emission structures of the reference run with a speed up of around 10 (Run Ba). These results convince us that we can use the non-Fourier heat flux description together with the Boris correction to acquire a significant speed up of the simulation without losing a correct representation of the physical processes within the solar corona in a statistical sense. We find some differences between the solution with and without non-Fourier heat flux description and the Boris correction. However, we are not interested if the non-Fourier heat flux description is identical to Fourier heat flux description in every time step at every specific location. Instead, we are interested if the non-Fourier heat flux description reproduced the Fourier heat flux description on a statistical level. On the statistical level, we find a very good agreement.

In the future, we are planning to use these implementations to perform large-scale active region simulations similar as done by Bourdin et al. (2013, 2014), which can be then run for a much longer time and allowing the study of hot core loop formations. A first attempt is already published (Warnecke and Peter 2019). Furthermore, this implementation allows us to perform parameter studies to investigate the coronal response to different types of active regions on the Sun and also on other stars. Finally, through these improvements, we get closer to the possibility to simulate a more realistic convection-zone-corona model as started in Warnecke et al. (2012b, 2013a, 2016a).

Acknowledgments

We thanks the anonymous referees for the useful suggestions. We also thank Piyali Chatterjee for the implementation of the Boris correction to the Pencil Code  and detailed discussions about it. We furthermore thank Hardi Peter for discussion about the non-Fourier heat flux description and comments to the manuscript. The simulations have been carried out on supercomputers at GWDG, on the Max Planck supercomputer at RZG in Garching, in the facilities hosted by the CSC-IT Center for Science in Espoo, Finland, which are financed by the Finnish Ministry of Education.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Jörn Warnecke http://orcid.org/0000-0002-9292-4600

Sven Bingert http://orcid.org/0000-0001-9547-1582

Additional information

Funding

Financial support from the Max Planck Princeton Centre for Plasma Physics (JW).

Notes

1 http://github.com/pencil-code.