Time step constraints

The time evolution of the variables is calculated through a third-order Runge-Kutta scheme.

The Runge-Kutta method is used to solve general differential equations of the form,

dyi

dt =Fi(yk,t). (4.9)

Knowing theyi at a specific timet₀we can calculate theyiatt₀+δtby evaluatingFi. The Kutta method follows this procedure in multiple sub-steps. The third-order Runge-Kutta used in our code solves the time evolution problem by employing three sub-steps for one time step, hence the name third order. The solution will have the form,

t y

t0 y0 k1 =δtF(y0,t0)

t1= t0+ ₁₅⁸δt y1 =y0+ ₁₅⁸k1 k2 =δtF(y₁,t1) t₂= t₁+ ²₃δt y₂ =y₀+ ¹₄k₁+ ₁₂⁵k₂ k₃ =δtF(y₂,t₂)

t=t0+δt y=y0+ ¹₄k1+0k₂+ ³₄k3+O(δt⁴)

This method, in general, requires less memory than other similar schemes, thus, making our code much more time-efficient.

The time step will generally be defined by the Courant time step given by, δt= min wherec_δt,c_δt,v,c_δt,sare the Courant-Friedrichs-Lewy (CFL) coefficients. The CFL coeffi -cients are smaller than unity, and they are generally used for convergence while solving differential equations numerically.δxmindefines the minimum grid space in the numerical box andU_maxdenotes the maximum velocity in our box as,

U_max=max

u+ q c²_s+u²_A

, (4.11)

wherecsis the sound speed anduAis the Alfvén speed. The maximum velocity in our box will be either the advective velocity, the sound speed, or the Alfvén velocity. Furthermore, D_maxdenotes the maximum diffusivity,

D_max =max

ν, γχ, η

(4.12) whereνis the kinematic viscosity,χis the thermal diffusivity, andηis the magnetic resis-tivity. Finally, theHmax represents the right-hand side of the energy equation. However, since our simulations are dominated by thermal diffusivity, this term will not be signifi-cant. In conclusion, the time step will be constrained by the thermal diffusivity, limiting the computational speed significantly. To boost the speed, we have to employ various techniques discussed next.

4.2.3.1 Non-Fourier treatment of heat flux

The efficiency of the Spitzer heat conduction puts constraints on the time step. To increase the time step and speed up our simulations, we use a non-Fourier description of the heat flux (Warnecke and Bingert 2020),

∂q

∂t = − 1

τ_Spitzer(q+K∇T), (4.13)

where τ_Spitzer is the heat flux relaxation time, and K is the Spitzer tensor with K = K0T^5/2bˆbˆ operating along the magnetic fieldlines. Here ˆbis the unit vector along the mag-netic fieldlines. From Eq. (4.13) we solve for qand then this is substituted in Eq. (4.5).

With this approach we replace the Spitzer time stepδt_Spitzer = δx²/χ_Spitzer with two new time steps, Here δx is the grid spacing and χSpitzer = |K|/ρcV is the diffusivity of the Spitzer heat conduction. The Eq. (4.13) can be expressed as a wave equation, therefore theδt₁ is also related to the wave propagation speedc_Spitzer. Notice here thatδt₁has a linear dependence on the grid sizeδx, which leads to a substantial increase of the computational time as the spatial resolution is increased.

Even though theτ_Spitzercan be adjusted manually, it is preferable to use the automatic adjustment implemented in the code by Warnecke and Bingert (2020). We first set the time step of the heat conduction to be the dominant, i.e. the smallest. The next dominant time scale in the box will be the Alfvén crossing timeδtA = δx/uA withuA ∝ B/√

ρthe Alfvén speed. The idea is that theτSpitzeris limited by the Alfvén time scale and we set,

c_Spitzer = √

2uA, → τ_Spitzer = χ_Spitzer

2u²_A . (4.15)

One problem that arises here is that if τ_Spitzer evolves freely, then it might get too large.

In that case, the Spitzer heat conduction will not be dominant anymore and might lead to non-physical solutions. For that reason, we set an upper limit to the time stepτSpitzer. In addition, the Spitzer coefficientχ_Spitzerreduces to low values in parts of low temperature and high density in the chromosphere or transition region. Therefore, this will reduce τSpitzer in those regions. Even though the Spitzer heat conduction is not dominating in those lower parts, the reducedτ_Spitzerwill decrease the overall time step tremendously. To avoid that problem, we set a lower limit of the time step to be in between of the Alfvén, and the advective time step,

wherecsis the sound speed,uis the advective velocity and we setτ_Spitzer,0 =100 s. Using the non-Fourier treatment of the heat flux Warnecke and Peter (2019b) report a boost in the time step by a factor of two.

4.2.3.2 Semi-relativistic Boris correction

The Alfvén speed approaches the speed of light in regions of high magnetic field strength and low density, i.e. above active regions (Rempel 2017). As a result, two main problems arise. First of all, the MHD approximation of neglecting the displacement current for non-relativistic velocities discussed in Chap. 2 is no longer valid. In addition, a significant increase of the Alfvén speed will lead to a very small time step. For those reasons, it is implemented in the code a semi-relativistic correction to the Lorentz force (Boris 1970;

Gombosi et al. 2002). The Lorentz force in this framework is replaced by,

0 10 20 30 40 50

Figure 4.1: Initial temperature and density profile implemented in our numerical exper-iments. Solid line show the initial temperature profile. Dashed line shows the initial density stratification based on hydrostatic equilibrium.

j×B special relativity. For the non-relativistic limit, the factorγA ' 1, and we retrieve the usual Lorentz force. For relativistic velocities, uA ' cthe above correction reduces the Lorentz force. The definition of the unit vector along the magnetic fiedlines ˆb = B/|B|

explains the factor BB/B² in Eq. (4.17). This method can also be used to increase the time step and prevent a crash of our simulations. The details of this method are discussed in Warnecke and Bingert (2020) and the implementation in the Pencil Code are discussed in Chatterjee (2018). The idea is to reduce the speed of light in the regions where the Alfvén speed is large. Following the work of Rempel (2017) the Alfvén time step will be,

dtA → dtA

wherecA is the limited speed of light. In our simulations we choose cA = 10000 km/s which leads to a time step ofdtA ' 20 ms. With this method, we ensure that the Alfvén time step will have a lower limit and it will not affect the overall time step in our simu-lations. As discussed in Warnecke and Peter (2019b), a combination of the non-Fourier heat flux scheme and the Boris correction leads to an increase of the time step by a factor of 10 or even a factor of 30. In conclusion, the Boris correction and the non-Fourier treat-ment of the heat flux combined will significantly boost the computational speed of our simulations. This allows performing a large number of simulations for a parameter study possible.