# Scalable explicit implementation of anisotropic diffusion with Runge-Kutta-
Legendre super-time-stepping

Bhargav Vaidya1,2, Deovrat Prasad3, Andrea Mignone2, Prateek Sharma3, Luca
Rickler2
1Centre of Astronomy, Indian Institute of Technology Indore, Simrol, Khandwa
Road, Indore 453552, India
2Dipartimento di Fisica Generale, Universita degli Studi di Torino , Via
Pietro Giuria 1, 10125 Torino, Italy
3Joint Astronomy Programme and Department of Physics, Indian Institute of
Science, Bangalore 560012, India
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

An important ingredient in numerical modelling of high temperature magnetised
astrophysical plasmas is the anisotropic transport of heat along magnetic
field lines from higher to lower temperatures.Magnetohydrodynamics (MHD)
typically involves solving the hyperbolic set of conservation equations along
with the induction equation. Incorporating anisotropic thermal conduction
requires to also treat parabolic terms arising from the diffusion operator. An
explicit treatment of parabolic terms will considerably reduce the simulation
time step due to its dependence on the square of the grid resolution ($\Delta
x$) for stability. Although an implicit scheme relaxes the constraint on
stability, it is difficult to distribute efficiently on a parallel
architecture. Treating parabolic terms with accelerated super-time stepping
(STS) methods has been discussed in literature but these methods suffer from
poor accuracy (first order in time) and also have difficult-to-choose tuneable
stability parameters. In this work we highlight a second order (in time) Runge
Kutta Legendre (RKL) scheme (first described by Meyer et al. 2012) that is
robust, fast and accurate in treating parabolic terms alongside the hyperbolic
conversation laws. We demonstrate its superiority over the first order super
time stepping schemes with standard tests and astrophysical applications. We
also show that explicit conduction is particularly robust in handling
saturated thermal conduction. Parallel scaling of explicit conduction using
RKL scheme is demonstrated up to more than $10^{4}$ processors.

###### keywords:

methods: numerical – (magnetohydrodynamics) MHD – conduction – instabilities –
galaxies: clusters: intra-cluster medium

##  1 Introduction

Since most baryons in the universe are in a magnetised plasma state, magnetic
fields play a crucial role in the dynamics and thermodynamics of astrophysical
objects — ranging from stars and interstellar medium to the intra-cluster and
intergalactic medium. Magnetohydrodynamic (MHD) simulations have matured
(Evans & Hawley 1988; Tóth 2000; Balsara 2001) and have contributed to several
breakthroughs in our understanding, from accretion to the interstellar medium
(Hawley et al. 1995; Korpi et al. 1999). Magnetic fields not only produce
forces and stresses in a plasma, they also affect transport properties by
predominantly allowing diffusion of heat and momentum along field lines and
suppressing transport across them. Anisotropic transport affects fundamental
properties such as convection/buoyancy in stratified plasmas (Balbus 2000;
Quataert 2008), and thermal instability and condensation of cold gas out of
the hot phase (Field 1965; Sharma et al. 2010b).

While the numerical solution of the ideal MHD equations can be carried out
very accurately in highly nonlinear regimes across thousands of processors
(Beresnyak 2011; Federrath et al. 2011), simulations with anisotropic thermal
conduction (and, likewise, similar diffusive processes) have not yet reached
the same level of fidelity and scalability.

Taking into account diffusion processes changes the mathematical structure of
the underlying system of conservation laws from purely hyperbolic to mixed
hyperbolic-parabolic type. The numerical discretisation of such systems can
then pose a more restrictive limitation on the choice of the time step which,
for a standard explicit scheme, must scale with the square of the grid size
rather than with $\Delta x$ alone. In addition, when saturation effects are
considered, thermal conduction itself becomes a mixed hyperbolic/parabolic
operator. Two approaches are commonly employed to circumvent the stability
time step constraint.

In the first one, the parabolic part of the equations is solved using fully-
or semi-implicit methods (Press et al. 1986), allowing for a time step much
longer than the explicit one (see e.g., Sharma & Hammett 2011). However,
scalable implementation of implicit schemes on massively parallel clusters is
technically challenging (Kannan et al. 2015) as it requires solving large
sparse matrices.

In the second approach, a particular stability polynomial can be employed to
construct an explicit multi-stage time stepping scheme that has extended
stability properties. If $s$ is the number of stages, the length of the super-
step that one can take (asymptotically) can be shown to scale as $\sim
s^{2}\Delta t_{\rm p}$, where $\Delta t_{\rm p}$ is the standard explicit
parabolic time step. This technique in general is classified as super-time-
stepping (STS). This technique offers therefore an overall speedup $\propto s$
and, being explicit, it can be easily parallelized. van Der Houwen & Sommeijer
1980; Verwer et al. 1990 introduced Runge-Kutta Chebyshev (RKC-STS) methods. A
variant of this method was popularized by Alexiades et al. 1996 (AAG-STS).
More recently, super-time-stepping based on Legendre polynomials, known as
Runge-Kutta Legendre STS (RKL-STS), was proposed by Meyer et al. (2012, 2014).

In this paper we review and compare the various STS schemes in the context of
anisotropic thermal conduction in a MHD plasma. We show that RKL-STS method is
an attractive option for simulating anisotropic diffusion on massively-
parallel supercomputers and offers enhanced stability with respect to its
predecessors, in particular, AAG-STS. We note that RKL-STS can also be used to
model other isotropic and anisotropic transport processes such as plasma
viscosity (e.g., Dong & Stone 2009), non-ideal terms in Ohm’s law
(e.g.,O’Sullivan & Downes 2007) , and cosmic ray diffusion (e.g., Sharma et
al. 2009; Pakmor et al. 2016). There are also applications in other diverse
areas such as image processing (e.g., Weickert 1998).

The paper is structured as follows. In section 2 we review the fundamental
equations of MHD with thermal conduction. In Section 3 we review the numerical
approaches to anisotropic thermal conduction and also discuss the monotonicity
properties in multiple dimensions. Selected numerical benchmarks are presented
in 4. The parallel scalability of explicit conduction is demonstrated in 5 and
conclusions are drawn in section 6.

##  2 Governing Equations

Here we consider the magnetohydrodynamic (MHD) equations in presence of
conduction:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bm{v})$
$\displaystyle=0,$ (1a) $\displaystyle\frac{\partial}{\partial
t}\left(\rho\bm{v}\right)+\nabla\cdot\left(\rho\bm{v}\bm{v}-\bm{BB}+p_{t}\mathsf{I}\right)$
$\displaystyle=0,$ (1b) $\displaystyle\frac{\partial E}{\partial
t}+\nabla\cdot\left[\left(E+p_{t}\right)\bm{v}-(\bm{B}\cdot\bm{v})\bm{B}\right]$
$\displaystyle=-\nabla\cdot\bm{F}_{c},$ (1c)
$\displaystyle\frac{\partial\bm{B}}{\partial
t}-\nabla\times(\bm{v}\times\bm{B})$ $\displaystyle=0,$ (1d)

where $\rho$ is mass density, $\bm{v}$ is fluid velocity, $\bm{B}$ is magnetic
field, $p_{t}=p+B^{2}/2$ is the total (gas + magnetic) pressure,
$E=\rho\epsilon+\rho v^{2}/2+B^{2}/2$ is the total energy density while
$\mathsf{I}$ is the identity tensor.

Thermal conduction effects, important for dilute gases and plasmas, are
accounted for by the additional term on the right hand side of the total
energy equation (1c). The thermal conduction flux $\bm{F}_{c}$ smoothly varies
between classical and saturated regimes,

$\bm{F}_{\rm c}=\frac{q}{q+F_{\rm class}}\bm{F}_{\rm class},$ (2)

where $F_{\rm class}$ and $q$ are, respectively, the magnitudes of the
classical thermal conduction flux (Cowie & McKee, 1977; Balsara et al., 2008)

$\bm{F}_{\rm class}=-\kappa_{\parallel}\bm{\hat{b}}(\bm{\hat{b}}\cdot\nabla
T)$ (3)

and the saturated flux

$q=5\phi\rho c_{\rm iso}^{3}\,.$ (4)

In Eq. (3), $\kappa_{\parallel}$ is the thermal conductivity parallel to the
direction of the magnetic field (a function of temperature; for expressions
see pp. 38 of Huba 2004) and $\bm{\hat{b}}\equiv\bm{B}/B$ is the magnetic
field unit vector. We do not consider thermal conduction perpendicular to
magnetic field lines because it is typically much smaller, but including it is
straightforward. In absence of magnetic fields, the classical heat flux is
simply given by (with $\kappa=\kappa_{\parallel}$)

$\bm{F}_{\rm class}=-\kappa\nabla T\,.$ (5)

Saturation of heat flux (Eq. 4) comes into play when the mean free path of the
electrons is large compared to the temperature scale-height $L_{T}=T/|\nabla
T|$ (Cowie & McKee, 1977). The saturated heat flux is given by
$q\bm{\hat{F}}_{\rm class}$ with $q$ specified in Eq. (4), while
$\bm{\hat{F}}_{\rm class}=\bm{F}_{\rm class}/F_{\rm
class}=-\bm{\hat{b}}~{}{\rm sgn}(\bm{\hat{b}}\cdot\nabla T)$ is a unit vector
along the local magnetic field line but down the temperature gradient. In Eq.
(4) the parameter $\phi$ is an uncertainty factor of the order of unity (taken
to be 0.3) while $c_{\rm iso}\equiv(p/\rho)^{1/2}$ is the isothermal sound
speed. In the large temperature gradient limit, thermal conduction is
described by a hyperbolic operator. On the contrary, for $F_{\rm class}/q\to
0$, thermal conduction is described by a parabolic (diffusion) operator. The
flux given in Eq. (2) is a harmonic mean of the fluxes in the two regimes
(weighted toward the smaller of the two) and reflects the mixed
parabolic/hyperbolic mathematical nature of the underlying differential
operators. Eqs. 2.10 & 2.10 in Balsara et al. (2008) gives some other ways of
combining saturated and classical heat fluxes. For more details on our
implementation of saturated conduction see Appendix A of Mignone et al.
(2012).

Note that thermal conduction described by Eq. (2) is very similar to the
cosmic ray streaming equation (see Eqs. 1.1, 1.2 & 6.1 in Sharma et al. 2010a)
in the sense that at temperature extrema, where gradient is zero and heat flux
is not saturated, the behaviour is diffusive but at other points (in the long
mean free path regime) it is hyperbolic. The harmonic mean used in Eq. (2) is
analogous to the regularisations of the cosmic ray streaming equation at
extrema proposed in Sharma et al. (2010a).

##  3 Numerical Approaches for Anisotropic Thermal Conduction

In this section we review the various methods for implementing anisotropic
thermal conduction.

We shall use the PLUTO code (Mignone et al., 2007, 2012) to solve Eqs. (1a-1d)
with the thermal conduction flux given by Eq. (2). In PLUTO, STS schemes have
been implemented in an operator split fashion so that we first advance Eqns.
(1a-1d) without the conduction flux using standard Godunov-type methods,
followed by the solution of

$\frac{\partial(\rho\epsilon)}{\partial t}=-\nabla\cdot\bm{F}_{c},$ (6)

by means of the STS/explicit approach.

A conservative approach is used during both steps so that the same flux is
used for adjacent grid cells sharing a face. During the parabolic step, Eq.
(6) is discretised using standard second-order finite differences:

$-\nabla\cdot\bm{F}_{c}\approx-\sum_{d}\frac{(\bm{\hat{e}}_{d}\cdot\bm{F}_{c})_{i+\frac{1}{2}}-(\bm{\hat{e}}_{d}\cdot\bm{F}_{c})_{i-\frac{1}{2}}}{\Delta
h_{d}}$ (7)

where $d=x,y,z$ spans across dimensions, $\Delta h_{d}=(\Delta x,\,\Delta
y,\,\Delta z)$ and $\bm{\hat{e}}_{d}$ are unit vectors in the three
directions. The interface flux is computed using central differences for the
temperature gradient and an upwind scheme for the saturated heat term, as
shown in appendix A of Mignone et al. (2012).

###  3.1 Explicit methods

In the standard explicit approach (forward in time centred in space) the time
step required to advanced the full system (1a)-(1d) is given by the shortest
timescale:

$\Delta t=\min\Big{(}\Delta t_{\rm h},\,\Delta t_{\rm p}\Big{)}\,,$ (8)

where

$\Delta t_{\rm
h}=\frac{C_{h}}{\displaystyle\max\left(\frac{\lambda_{x}}{\Delta
x}+\frac{\lambda_{y}}{\Delta y}+\frac{\lambda_{z}}{\Delta z}\right)}$ (9)

is the MHD/hydro (hyperbolic) time step whereas

$\Delta t_{\rm p}=\frac{C_{p}}{\displaystyle\max\left(\frac{\chi_{x}}{\Delta
x^{2}}+\frac{\chi_{y}}{\Delta y^{2}}+\frac{\chi_{z}}{\Delta z^{2}}\right)}$
(10)

is the diffusion time scale. In Eqs. 9 and 10, the maximum is taken over the
whole computational domain, $C_{h}\leq 1$ is the hyperbolic Courant number,
$C_{p}\leq 1/2$ is the parabolic Courant number,
$\lambda_{d}=|v_{d}|+c_{d,f}$, $c_{d,f}$ is the fast magneto-sonic speed in
the $d$ direction, $\Delta x,\Delta y$ and $\Delta z$ are the mesh spacings in
the three directions while $\chi_{x},\chi_{y}$ and $\chi_{z}$ are the thermal
diffusivities with the dimensions of $[{\rm length}]^{2}/[{\rm time}]$.
Adopting an ideal EoS ($p=[\gamma-1]\rho\epsilon$) one obtains, from Eq. (6)
together with Eq. (5) or Eq. (3):

$\chi_{d}=\left\\{\begin{array}[]{ll}\displaystyle\kappa\frac{\gamma-1}{\rho}&\qquad\mathrm{(HD)}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\kappa_{\parallel}b_{d}\cos\theta\frac{\gamma-1}{\rho}&\qquad\mathrm{(MHD)}\end{array}\right.$
(11)

where $b_{d}=B_{d}/B$ is the component of the magnetic field unit vector in
the direction $x_{d}$ and $\theta$ is the angle between $\bm{B}$ and $\nabla
T$ while $d=x,y,z$ (note that, owing to operator splitting, density can be
considered constant during the diffusion step).

Various explicit methods are discussed in detail in section 2 of Sharma &
Hammett (2007). A notable symmetric explicit scheme with a very small
numerical diffusion perpendicular to the local field lines is described in
Günter et al. (2005). This method has been used widely in modelling fusion
plasmas, in which the temperature and field gradients are relatively gentle.
However, for highly nonlinear problems with large temperature gradients, this
method gives unphysical temperature oscillations (see top right panel in Fig.
6 of Sharma & Hammett 2007).

It must emphasised that, for large grid resolutions and/or thermal
conductivities, the diffusion timescale can become much shorter than the
MHD/hydro time step (e.g., in galaxy clusters; Wagh et al. 2014; Yang &
Reynolds 2015), that is $\Delta t_{\rm p}\ll\Delta t_{\rm h}$. A fully
explicit method can thus become very inefficient. In these cases the ability
to take larger steps than the stability limit (Eq. 8) is highly desirable.

A slight improvement over explicit time-stepping is obtained via sub-cycling,
in which the MHD module is evolved with $\Delta t_{\rm h}$ but the conduction
module is applied for a number of sub-steps ($N_{\rm sub}\approx\Delta t_{\rm
h}/\Delta t_{\rm p}$; assuming $\Delta t_{\rm p}<\Delta t_{\rm h}$), each
using the stable time step $\Delta t_{\rm p}$.While this approach is faster
than the fully explicit approach, as we show in section 3.3, we can do much
better than this using super-time-stepping.

###  3.2 Implicit methods

Implicit methods are very attractive because their time step is not limited by
the stability limit (Eq. 8). A larger time step only leads to a slight loss in
accuracy. Also, for very high resolution simulations (i.e., with very small
$\Delta x$), the CFL time step (Eq. 9) becomes much longer than the stability
time step with diffusion (Eq. 10). For this reason implicit and semi-implicit
schemes for diffusion are very popular (e.g., Sharma & Hammett 2011; Dubois &
Commerçon 2016; Kannan et al. 2015; Pakmor et al. 2016).

The key problem with implicit methods is that they involve solving a sparse
matrix equation, which is difficult to solve in parallel in a scalable way.
Even a tridiagonal matrix equation which results from the semi-implicit method
of Sharma & Hammett (2011) is non-trivial to solve in parallel. The MHD codes
using implicit diffusion typically use sparse matrix libraries such as
Hypre111http://acts.nersc.gov/hypre/ or PETSc222http://www.mcs.anl.gov/petsc/
to solve the matrix equation iteratively. This global solution approach is
fundamentally different from the explicit update in which information only
propagates across nearest grid cells in a single time step. Most importantly,
even the scalable implicit methods for diffusion with MHD do not scale on more
than $\sim$1000 processors on massive distributed memory supercomputers (e.g.,
see Fig. 11 in Caplan et al. 2017). In contrast, the STS methods that we
advocate (in particular, RKL-STS) in this paper, being explicit, are trivial
to parallelise and show good scaling up to tens of thousands of processors
(see section 5).

###  3.3 Super-time-stepping: RKC, RKL, AAG

Super-time-stepping (STS) is a way of choosing explicit time steps that are
on-average much longer than the explicit stability time step (Eq. 10). The
anisotropic diffusion equation has non-positive eigenvalues; i.e., none of the
Fourier modes grow with time. Following Meyer et al. 2014 (see their section
2), we can write the diffusion equation in a semi-discrete ODE form

$\frac{du}{dt}=\mathsf{M}u,$ (12)

where, for a linear problem, $\mathsf{M}$ is a symmetric constant coefficient
matrix which results from the discretization of the parabolic operator. The
eigenvalues $\lambda$ of such an operator are non-positive and real.

The STS method can indeed be thought of as a multi-stage Runge-Kutta method in
which the intermediate stages are chosen for stability, rather than for a
higher order accuracy. The update of an $s-$stage scheme is quantified in
terms of an amplification factor, $R_{s}$, defined as

$u(t+\tau)=R_{s}(\tau\mathsf{M})u(t).$ (13)

We now define a stability polynomial for an $s$-stage STS scheme,

$R_{s}(\tau\lambda)\equiv\prod_{i=1}^{s}(1+\lambda\Delta t_{i})$

with $\sum_{j=1}^{s}\Delta t_{j}=\tau$, and we impose
$|R_{s}(\tau\lambda)|\leq 1$ for all values of $\lambda$ between the largest
negative eigenvalue of $\mathsf{M}$ and 0 in order to ensure stability.
Temporal accuracy is achieved by matching terms in the stability polynomial
with the expansion of the analytic solution of Eq. (12),

$u(t+\tau)=\left(1+\tau\mathsf{M}+\frac{1}{2}(\tau\mathsf{M})^{2}+...\right)u(t).$
(14)

For first (second) order STS the stability polynomial only matches the above
analytic expression up to the linear (quadratic) term.

The sub-step sequence in STS was originally chosen by writing the stability
polynomial, $R_{s}$ (a polynomial of degree $s$ in $\lambda$), in form of a
shifted Chebyshev polynomial (CP; e.g., Verwer et al. 1990); we refer to this
scheme as RKC (Runge-Kutta Chebyshev). The sub-steps in a super-step are
chosen to exploit the recursion properties of CPs. Since the absolute value of
CPs is always $\leq 1$ if their argument lies in $[-1,1]$, the magnitude of
amplification factor after a super-time step is $\leq 1$. Since CPs attain
unity within $[-1,1]$, RKC-STS requires damping for the method to work for
practical diffusion problems.

Recently, Meyer et al. (2012, 2014) have proposed Legendre polynomials (LPs)
as the basis for the construction of robust STS schemes. For a general
$s$-stage RKL scheme, the stability polynomial is chosen to be of the form

$R_{s}(\tau\lambda)=a_{s}+b_{s}\mathcal{P}_{s}\left(w_{0}+w_{1}(\tau\lambda)\right)$
(15)

where $w_{0}=1$ is chosen for all RKL schemes. The LPs,
$\mathcal{P}_{j}\left(x\right)$ obey a three point recursion relationship
given by

$\mathcal{P}_{j}\left(x\right)=\left(\frac{2j-1}{j}\right)x\mathcal{P}_{j-1}\left(x\right)-\left(\frac{j-1}{j}\right)\mathcal{P}_{j-2}\left(x\right).$
(16)

RKL-2 method enforces the stability polynomial at sub-stage $j$ to be
$a_{j}+b_{j}\mathcal{P}_{j}\left(\tau\mathsf{M}\right)$ and uses the property
of LPs that their absolute value is bounded by unity if their argument lies in
(-1,1). Following Meyer et al. (2014), in Appendix A we illustrate the second
order temporal accuracy of the RKL-2 scheme along with explicit formulae for
an $s=3$ stage scheme.

The key advantage of RKL schemes over RKC is that they are more robust because
LPs, unlike CPs, are always smaller than unity in magnitude if their argument
lies in (-1,1). Thus, there is no need of an explicit damping parameter ($\nu$
in Eqs. 9 & 10 of Alexiades et al. 1996). Moreover, Meyer et al. (2014) show
that RKL has superior linear stability and monotonicity properties compared to
RKC.

For stability, the argument of the LP in Eq. (15) should be $\geq-1$ (note
that for a parabolic operator the eigenvalues are non-positive); i.e.,
$1-w_{1}\tau|\lambda|_{\rm max}\geq-1$ or $\tau\leq 2/(|\lambda|_{\rm
max}w_{1})$, where $|\lambda|_{\rm max}$ is the eigenvalue of $\mathsf{M}$
with largest absolute value. Note that $w_{1}=2/(s^{2}+s)$ and $4/(s^{2}+s-2)$
for RKL-1 and RKL-2, respectively. Given a super-time-step that can be
conveniently chosen as the hyperbolic CFL time step $\Delta t_{\rm h}$ (Eq.
9), the number of sub-stages in RKL-STS is given by (see Eq. 20, 21 in Meyer
et al. 2014)

$\begin{array}[]{lcll}\Delta t_{\rm h}&=&\displaystyle\Delta t_{\rm
p}\frac{(s^{2}+s)}{2}&\text{for RKL-1,}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\Delta t_{\rm h}&=&\displaystyle\Delta t_{\rm
p}\frac{(s^{2}+s-2)}{4}&\text{for RKL-2}.\end{array}$ (17)

where, abbreviations RKL-1 and RKL-2 stand for the first and second order
Runge-Kutta-Legendre schemes.We assume here that $s>1$; $s=1$ corresponds to a
standard forward Euler update. In this paper, RKL without a qualifier refers
to RKL-2. We explicitly mention if we are using the first order method. It is
possible to make the scheme even faster by taking multiple $\Delta t_{\rm h}$
as the size of a super-step but of course one needs to worry about accuracy.

Likewise, for AAG-STS, the largest number of first order Euler sub-steps $s$
satisfies (see Eq. 10 in Alexiades et al. 1996)

$\Delta t_{\rm h}=\Delta t_{\rm
p}\frac{s}{2\sqrt{\nu}}\left[\frac{(1+\sqrt{\nu})^{2s}-(1-\sqrt{\nu})^{2s}}{(1+\sqrt{\nu})^{2s}+(1-\sqrt{\nu})^{2s}}\right].$
(18)

The maximum $s$ for STS schemes gives a speed-up of the conduction module
scaling as $(\Delta t_{\rm h}/\Delta t_{\rm p})^{1/2}$. The speed-up is
marginally better for RKC (compare Eqs. 18, 19 in Meyer et al. 2012) although
RKL has better stability (as shown in section 4), and hence is much more
attractive.

It is instructive to compute the effective parabolic CFL number
$C_{p}=\chi\Delta t/\Delta x^{2}=\Delta t_{\rm h}/2\Delta t_{\rm p}$ as a
function of the number of sub-steps $s$. In the limit of large $s$ ($s\gg 1$),
the RKL schemes (17) yield

$C_{p}\to\frac{s^{2}}{4}\quad\mathrm{(RKL-1)}\,,\qquad
C_{p}\to\frac{s^{2}}{8}\quad\mathrm{(RKL-2)}.$ (19)

However, in the same limit, AAG-STS yields an effective gain which depends on
$\nu$:

$C_{p}\to\left\\{\begin{array}[]{lcl}\displaystyle\frac{s^{2}}{2}&\quad\mathrm{for}&s\lesssim
1/(2\sqrt{\nu}).\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\frac{s}{4\sqrt{\nu}}&\quad\mathrm{for}&s\gtrsim
1/(2\sqrt{\nu}).\end{array}\right.$ (20)

This equation sets an upper limit on the gain that can be reached using the
traditional AAG-STS scheme: no substantial gain is obtained when $s\gtrsim
1/(2\sqrt{\nu})$.

Table 1: Expected Flops relative to MHD evolution

MHD | explicit | sub-cycling | STS |  implicit†
---|---|---|---|---
1 | $(1+r^{-1})t^{-1}$ | $1+t^{-1}r^{-1}$ | $1+t^{-1/2}r^{-1}$ | $1+r^{-1}$

Notes: $r$ is the ratio of number of Flops needed for the MHD module and the
conduction module; $t=\Delta t_{\rm p}/\Delta t_{\rm h}$ (taken to be $<1$;
Eqs. 10 & 9).
† here it is assumed that a single implicit update for conduction takes the
same number of Flops as a single explicit update; the implicit update is
likely to be more expensive and the number of expected Flops is larger.

As mentioned earlier, the key advantage of STS methods is their simplicity,
parallel scaling on massively distributed clusters, and their large speed-up
compared to the explicit method. Table 1 shows the expected scaling of Flops
(equivalently, run-time) using various conduction schemes relative to pure MHD
evolution. Here, $t=\Delta t_{\rm p}/\Delta t_{\rm h}$ (assumed to be $<1$)
and $r$ is the ratio of Flops (floating point operations) for MHD and
conduction modules. The typical value of $r$ is a few (say 5; this is
relatively independent of the physical set-up and dimensionality) but $t$ can
be smaller than 0.01 (this depends on the physical set-up). For these typical
values the ratios of Flops required in Table 1 (for MHD, explicit, sub-
cycling, STS, implicit) become: 1, 120, 21, 3, 1.2. Therefore, STS promises a
good speed-up compared to the explicit method, and is competitive relative to
implicit methods (especially given the better parallel scaling of the former).

###  3.4 Monotonicity in Multi-Dimensions

A serious problem for anisotropic diffusion, in presence of large gradients,
is that the temperature can behave non-monotonically. The use of limiters for
interpolating transverse temperature gradients has proved useful in preventing
negative temperatures in such cases (Sharma & Hammett 2007; see Appendix B for
details).

The use of limiters to calculate the transverse temperature gradients (Eq. 51)
has shown to improve the robustness of all schemes (e.g., see Sharma & Hammett
2007, 2011). However, we show in section 4.1.3 that limiters reduce the
accuracy of STS schemes for a large number of sub-steps. Given this, we keep
limiters as an option to be used only when temperature gradients are large and
the wrong sign of the heat flux gives a negative temperature in some grid
cells.

##  4 Numerical Tests: effectiveness of RKL-STS

In the following we present selected numerical benchmarks. We adopt an ideal
equation of state so that

$\rho\epsilon=\frac{p}{\gamma-1}$ (21)

where $\gamma$ is the ratio of specific heats which equals $5/3$ unless
otherwise stated.

###  4.1 Scalar Diffusion Equation

In this section we consider numerical tests based on the solution of the
scalar diffusion equation

$\frac{\partial(\rho\epsilon)}{\partial t}=-\nabla\cdot\bm{F}_{c},$ (22)

where $\bm{F}_{c}$ is given by Eq (2). Other fluid variables are not evolved
in time.

In all computations, we assign the super-step $\Delta t_{\rm h}$ and compute
the number of sub-steps by solving either Eq. 17 (for RKL) or Eq. 18 (for AAG-
STS) for $s$ and then rounding the solution to the next larger integer, i.e.,
$s\,\to\,1+{\rm floor}(s)$.

####  4.1.1 Gaussian diffusion in 1-D

A standard one dimensional diffusion test is useful to compare the stability
and accuracy of different methods described in section 3.3. For this purpose,
we solve the standard heat equation (without saturation) for temperature:

$\frac{\partial}{\partial t}T(x,t)=\frac{\partial}{\partial
x}\left(\kappa\frac{\partial T(x,t)}{\partial x}\right),$ (23)

using a constant diffusion coefficient $\kappa=1$. In the PLUTO code, this is
equivalent to solving Eq. 6 by setting $\rho=1$ and $\gamma=2$, so that
$\rho\epsilon=p=T$ in code units.

The initial condition consists of a Gaussian temperature profile
$T(x,0)=\exp(-x^{2}/2\sigma^{2})$ inside the domain $x\in[-2,2]$. Then Eq. 23
admits the well-known analytical solution

$T^{\rm ref}(x,t)=\frac{1}{\sqrt{1+\frac{2\kappa
t}{\sigma^{2}}}}\exp\left[-\frac{x^{2}}{2\sigma^{2}\left(1+\frac{2\kappa
t}{\sigma^{2}}\right)}\right]\,.$ (24)

Boundary conditions are specified using the exact solution.

Figure 1 shows the computed temperature profile (symbols) versus the
analytical solution (Eq. 24; solid line) at different times using RKL-2 on
$N_{x}=128$ grid points and a parabolic CFL number $C_{p}~{}(\equiv\Delta
t_{\rm h}/2\Delta t_{\rm p})=10$. As expected, the initial peak of the
Gaussian temperature profile spreads symmetrically about $x=0$ with time.

Figure 1: Temperature profiles (blue square symbols) for the 1D Gaussian
diffusion problem at $t=0,\,0.08,\,0.22,\,0.45$ obtained with RKL-2 using
$N_{x}=128$ grid zones and a parabolic CFL number $C_{p}=10$. The exact
analytical solution (Eq. 24) is given by the black solid line.

Figure 2: Error analysis for the 1D Gaussian test problem using different
explicit schemes. Left panel: errors in L1 norm as functions of the grid
resolution $N_{x}$ for AAG-STS (red), RKL-1 (green) and RKL-2 (blue) methods.
The two dotted lines give the ideal scalings for first and second order
schemes. Middle panel: total number of steps (integration steps $\times$
number of sub-steps) for the explicit 1st order scheme (yellow) and AAG-STS,
RKL-1 and RKL-2 (red, green and blue respectively). The computational workload
follows the expected scaling of $N_{x}^{2}$ (explicit) and $N_{x}^{3/2}$ (STS)
shown by dotted lines. Right panel: error in L1 norm as a function of the
parabolic CFL number $C_{p}=\kappa\Delta t_{\rm h}/\Delta x^{2}$ at the fixed
resolution $N_{x}=4096$.

In order to verify the order of accuracy of the various methods considered
here, we first perform a resolution study by doubling the grid resolution from
$N_{x}=64$ up to $N_{x}=16384$. We quantify the accuracy for different
integration methods, AAG-STS (we use $\nu=10^{-3}$ unless otherwise stated),
RKL-1, RKL-2, by estimating the error in $L1$ norm as

$\epsilon_{L1}=\frac{1}{N_{x}}\sum_{i}\left|T^{n}_{i}-T^{\rm
ref}(x_{i},t^{n})\right|.$ (25)

Since we use standard centred differencing for the second derivative in Eq. 23
with local error $O(\Delta x^{2})$, we must ensure that the error in our
simulation is not dominated by spatial discretisation. To this end, we adjust
the time step such that the ratio of the time step and the grid spacing is
fixed; i.e.,

$\frac{\Delta t_{\rm h}}{\Delta x}=\frac{C_{p}\Delta
x}{\kappa}=\frac{C_{p0}}{N_{x0}\kappa}={\rm const},$ (26)

where $N_{x0}=64$ and $C_{p0}=0.45$ are, respectively, the number of zones and
the parabolic CFL number at the lowest resolution. This yields a local
truncation error (for $T$) of $O(\Delta t^{m})+O(\Delta t^{3})$ where $m=2,3$
for a first or second order scheme, respectively. Notice also that the
parabolic CFL number increases linearly with the spatial resolution.

The left panel in Fig. 2 shows the errors obtained for AAG-STS, RKL-1 and
RKL-2 methods together with the expected first and second-order convergence
rates. From the figure one can verify that the order of accuracy for AAG-STS
and RKL-1 is that of a $1^{\rm st}$ order method whereas RKL-2 converges as a
$2^{\rm nd}$ order method, as expected. We remark that integration with AAG-
STS shows the occurrence of negative temperatures during sub-steps for
$N_{x}\geq 8192$ (corresponding to $s\gtrsim 30$), although the solution
remains positive at the end of the super-step. This can pose serious
difficulties in more complex applications in which the diffusion coefficient
is a nonlinear function of the temperature (e.g., Spitzer conductivity
$\kappa\propto T^{5/2}$) requiring $T$ to be non-negative at all times. On the
contrary, RKL-1 and RKL-2 never exhibit such a behaviour (i.e., the solution
remains positive for all sub-steps; this property is enforced by construction
in the RKL scheme as described at the end of section 2.2 of Meyer et al.
2014).

In the middle panel of Figure 2 we plot the computational time (measured as
the total number of steps and sub-steps) as a function of the number of grid
points $N_{x}$. For an explicit scheme, let $N_{t}\propto N_{x}^{2}$ be the
total number of steps required to reach some final time step. Then, from Eqs.
17 & 20, we expect the number of sub-steps using STS schemes to scale roughly
as $s\propto\sqrt{\Delta t_{\rm h}/\Delta t_{\rm p}}\propto\sqrt{N_{x}}$. The
total number of steps (integration steps $\times$ number of sub-steps) is
therefore expected to be $N_{t}\propto N_{x}s\propto N_{x}\sqrt{N_{x}}$. This
behaviour is verified in the middle panel of Figure 2 from which we conclude
that STS techniques provide an effective asymptotic gain over standard
explicit time-stepping proportional to the square root of the number of grid
points.

Next, in the right panel of Figure 2, we plot the L1-norm errors by changing
the parabolic CFL number $C_{p}=\kappa\Delta t_{\rm h}/\Delta x^{2}$ for a
fixed grid resolution $N_{x}=4096$. Although comparable for $C_{p}\lesssim 1$,
the errors grow linearly with the CFL number for the first order schemes like
AAG-STS and RKL-1. However, for $C_{p}\gtrsim 200$, AAG-STS becomes unstable
and integration is no longer possible, while RKL-1 remains stable without
showing any significant undershoot even during each cycle sub-step. The second
order scheme maintains approximately the same accuracy for $C_{p}\lesssim
10^{3}$ and the error starts to increase more rapidly for larger CFL although
the solution remains well-behaved and positive at all times.

Figure 3: Two-dimensional contour plot of the logarithm of the L1 error (Eq.
25) for the 1D Gaussian test problem as a function of the parabolic CFL number
$C_{p}$ and the $\nu$ parameter for the AAG-STS scheme with fixed resolution
$N_{x}=1024$. Blue colours indicate smaller errors while red, orange and
yellow denote increasingly larger errors and unstable behaviour. The number of
sub-steps $s$ is over-plotted using white contour lines. The value of $s$
depends on the CFL number which varies from $4$ to $4096$. The thick solid
line represents $\nu=1/(8C_{p})$, the efficiency limit above which the
parabolic CFL number grows linearly with the number of sub-steps (Eq. 20).
Notice that all computations fail when $s\gtrsim 50$.

Finally we investigate the stability properties of the AAG-STS method alone,
by varying both the parabolic CFL number and the $\nu$ parameter at the fixed
grid resolution of $N_{x}=1024$ zones. Fig. 3 shows the logarithm of the error
as a function of $C_{p}$ and $\nu$. The corresponding number of sub-steps $s$
is also over-plotted using white contour lines. For sufficiently small values
of the damping parameter ($\nu\lesssim 10^{-4}$), we see that the error
rapidly increases when $C_{p}$ exceeds $\sim 200$ ($s\sim 20$) and
computations eventually become unstable for $C_{p}\gtrsim 10^{3}$ ($s\gtrsim
45$), irrespective of the value of $\nu$. For $s\lesssim 1/(2\sqrt{\nu})$, Eq.
20 shows that $C_{p}\approx s^{2}/2$ is independent of $\nu$. For larger
values of the damping parameter ($\nu\gtrsim 10^{-4}$), however, the number of
steps required to complete the calculation at a given CFL number increases and
therefore we observe a loss of efficiency that reduces the limiting CFL number
from $C_{p}\approx 10^{3}$ (at $\nu\approx 10^{-4}$) to $C_{p}\approx 200$
(for $\nu\approx 10^{-2}$). Of course, increasing the number of sub-steps at a
given CFL number leads to larger stability at the cost of extra computational
work. Finally, we point out that our temperature profiles are marginally
affected by the numerical resolution, the effect of which is that of
triggering (in case of numerical instability for $s\gtrsim$ 40) the growth of
Nyquist ($k=\pi/\Delta x$) mode.

####  4.1.2 Saw-tooth profile with saturated heat flux

In the next example we compare the performance of the selected integration
schemes by also including saturation of the heat flux. More specifically, we
consider the initial sawtooth temperature profile

$T(x,0)=10+20\left[\frac{x}{10}-\mathrm{floor}\left(\frac{x}{10}+\frac{1}{2}\right)\right]$
(27)

and solve Eq. 6 by setting $\rho=1$ and $\gamma=2$ so that $\rho\epsilon=p=T$
in code units, as in the previous example. The thermal conduction flux is
given by Eq. 2 with $\bm{F}_{\rm class}$ and $q$ defined in Eq. 5 and Eq. 4,
respectively, with $\kappa=40$, $\phi=0.3$. Computations are carried out on
the one-dimensional domain $x\in[-L_{x}/2,L_{x}/2]$ using $N_{x}=400$ zones
and periodic boundary conditions. The total number of steps to reach some
final time $t_{\rm stop}$ is therefore $N_{\rm step}=t_{\rm stop}/\Delta
t_{\rm h}$, with $\Delta t_{\rm h}$ computed from the parabolic CFL number

$N_{\rm step}=\frac{t_{\rm
stop}\kappa}{C_{p}}\left(\frac{N_{x}}{L_{x}}\right)^{2}=\frac{1600}{C_{p}}\,,$
(28)

where we have used $t_{\rm stop}=0.1$.

Figure 4: Numerical solution of the sawtooth temperature profile with heat
saturation at $t=0.1$. Top panel: computations obtained with parabolic CFL
number $C_{p}=2$ for AAG-STS ($\nu=10^{-3},\,10^{-2}$ orange and red,
respectively), RKL-1 (green) and reference solution (black dashed line).
Bottom panel: same as before but using a parabolic CFL number $C_{p}=50$.
Spurious oscillations appear using AAG-STS with $\nu=10^{-3}$.

We perform two sets of computations corresponding to $C_{p}=5$ and $C_{p}=50$
and compare the results with a reference solution obtained on a much finer
grid. Results are shown in Fig. 4. In the case of AAG-STS we employ
$\nu=10^{-3}$ (orange) and $\nu=10^{-2}$ (red).

Away from extrema, the evolution is initially dominated by the contribution of
the saturated flux only (since $\partial^{2}T/\partial x^{2}=0$ everywhere
except where the first derivative is discontinuous). The diffusive part of the
flux acts in those region where a change in slope is present. For $C_{p}=5$,
all methods yield well-behaved solutions with comparable errors whereas for
$C_{p}=50$, high frequency spurious oscillations appear when using AAG-STS
with a lower values of the damping parameter ($\nu=10^{-3}$) as shown by the
orange line in Figure 4. Oscillations originate in proximity of the maxima
where the second derivative of $T$ does not vanish and propagate downstream as
the system evolves. Increasing the grid resolution or the parabolic CFL number
tends to amplify this unstable behaviour. We note that oscillations disappear
if saturation is not included.

On the other hand, computations remain stable when using RKL-1 (green) or
RKL-2 (not shown) even for larger parabolic CFL numbers ($C_{p}\lesssim 200$).
For $C_{p}\gtrsim 200$ we observe larger numerical errors in the solution even
with RKL. The production of temperature oscillations at maxima is analogous to
numerical oscillations produced at extrema when numerically solving the cosmic
ray streaming equation (see Fig. 2.1 in Sharma et al. 2010a). Unlike in cosmic
ray streaming, oscillations are not produced at temperature minima because of
a smaller streaming speed ($\propto\sqrt{T}$). We have verified this
dependence on streaming speed by running with a higher initial temperature (T
= 30, instead of T = 10) and noticing oscillatory behavior both at temperature
maxima and minima in the case of AAG-STS with $\nu=10^{-3}$, while other cases
show stable solutions.

####  4.1.3 Ring diffusion

The 2-D Cartesian ring diffusion test problem, presented in Parrish & Stone
(2005); Sharma & Hammett (2007), is useful to study the monotonicity
properties of various numerical schemes for anisotropic diffusion in presence
of temperature discontinuities. Temperature discontinuities are fairly common
in astrophysical fluids and plasmas, and an ideal numerical scheme should not
lead to negative temperatures in presence of large temperature gradients. For
an explicit update with $\Delta t\leq\Delta t_{\rm p}$, Sharma & Hammett
(2007) showed that temperature monotonicity is maintained if limiters are used
to interpolate traverse temperature gradients (see Appendix B for details);
using arithmetic mean for interpolation leads to non-monotonicity in general.
The use of limiters with an implicit/semi-implicit update does not strictly
maintain monotonicity but improves monotonicity substantially compared with
arithmetic averaging (see Figs. 4 & 5 in Sharma & Hammett 2011).

We numerically solve the anisotropic thermal diffusion equation, Eq. 6 with
$\kappa_{\parallel}=1$; saturation of heat flux is ignored for this test
problem. We set $\rho=1$, $\gamma=2$, and $\rho\epsilon=p=T$ in code units.
The magnetic field lines are circular with $B_{x}=y/(x^{2}+y^{2})^{1/2}$ and
$B_{y}=-x/(x^{2}+y^{2})^{1/2}$, and the fluid is static. All variables except
$T$ ($\rho$, $\bm{B}$, $\bm{v}$) are held fixed; $T$ evolves only because of
anisotropic diffusion. The computational domain is $[-1,1]\times[-1,1]$, equal
number of grid points ($N$) are used in the two directions, and periodic
boundary condition is imposed on $T$. The initial condition on $T$ is

$T=\left\\{\begin{array}[]{ll}\displaystyle
12&\mathrm{for}\quad\displaystyle\frac{11\pi}{12}\leq\theta\leq\frac{13\pi}{12},\quad
0.5\leq r\leq 0.7,\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle
10&\mathrm{otherwise},\end{array}\right.$ (29)

where $r=(x^{2}+y^{2})^{1/2}$ and $\theta$ ($\tan\theta=y/x$) are polar
coordinates. The magnetic field lines are not aligned with the Cartesian grid
and the transverse heat flux (Eq. 51) is non-zero. The anisotropic diffusion
equation is evolved until $t=1$ using different schemes.

Figure 5: Temperature (in code units) contour plots for the ring diffusion
test problem at $t=1$ evolved using RKL-2 STS without (left panel) and with
(right panel) limiters for calculating the transverse temperature gradients in
Eq. 49. The grid resolution is $(50,100,200)$ and the number of RKL stages is
$5,20,50$. Without limiters, the minimum temperature is always smaller than
10, the initial minimum. With limiters, the evolution is more diffusive and
the temperature is below 10 only for $s=50$ and $N=50$.

Figure 5 shows the temperature at $t=1$ using RKL-2 scheme with different
resolutions ($N=50,~{}100,~{}200$) and different number of STS stages
($s=20,~{}50,~{}100$); left (right) panel shows results without (with)
limiters for interpolating transverse temperature gradients (see Appendix B).
For all cases without limiters the minimum temperature at $t=1$ is less than
10, the initial minimum temperature.333For the present test problem with a
maximum and minimum temperature of 12 and 10, the minimum temperature at later
times is not negative but it will become negative if the initial minimum and
maximum temperatures are say 0.1 and 10. A temperature ratio of 100 is
commonplace in multiphase astrophysical flows, such as the interstellar
medium. Once temperature becomes negative the MHD solver breaks down because
of an imaginary sound speed. With limiters, only $s=50$ and $N=50$ run shows a
temperature below 10. Therefore, monotonicity is better maintained with the
use of limiters. A comparison of left and right panels shows that the use of
limiters introduces large numerical diffusion perpendicular to field lines,
especially for smaller resolution and larger number of stages (i.e., for
smaller $N$ and larger $s$). A similar behaviour is observed for other schemes
such as AAG-STS, RKC-STS and RKL-1. For very large $\Delta t_{\rm h}/\Delta
t_{\rm p}$, corresponding to $s\sim\sqrt{\Delta t_{\rm h}/\Delta t_{\rm
p}}\gtrsim 20$, it is better not to use the limiters because of excessive
transverse diffusion. (A very large $s$, although numerically stable, also
leads to a loss of accuracy.) However, limiters are necessary for preventing
negative temperatures in presence of large gradients. For practical
implementation we recommend the use of limiters only when temperature becomes
smaller than a reasonable floor value in the anisotropic conduction step.

###  4.2 MHD Tests

Unlike section 4.1, the test problems described in this section solve the full
set of MHD equations in presence of thermal conduction (Eqs. 1a-1d). These
tests demonstrate the coupling of the MHD hyperbolic conservation laws with
the parabolic update of thermal conduction.

####  4.2.1 Supernova blast-wave

We consider an MHD blast wave in cylindrical co-ordinates with initial
parameters similar to the $L1$ model of Meyer et al. (2012) but without
radiative cooling. For this 2D axisymmetric test, we solve the standard set of
ideal MHD equations, taking into account anisotropy in thermal conduction flux
along with saturated conduction with $\phi=0.3$ (see Eqs. 2, 3 and 4).

The initial magnetic field of $0.3~{}\mu$G is oriented along the $z-$
axis.Energy of $10^{51}\rm erg$, equivalent to a supernova explosion, is
injected in a region of spherical radius ($[r^{2}+z^{2}]^{1/2}$) $7.8$ pc
around the origin of a $300$ pc domain in cylindrical geometry. A high
pressure is set inside this spherical region (resolved with 10 zones) such
that one-third of the supernova energy goes as thermal energy. The radial
velocity of this hot material is set by equating the kinetic energy to the
remaining two-third of energy. The temperature in the ambient medium is set to
be $8000$ K and the initial density is set to be $0.7m_{u}$ cm-3 ($m_{u}$ is
atomic mass unit) throughout the domain. Axisymmetric boundary conditions are
imposed at $r=0.0$ and initial conditions are imposed at the vertical boundary
$z=0.0$. The conditions at the outer boundaries in both $r-$ and $z-$
directions are imposed to be outflow.

Figure 6 shows the comparison of logarithmic value of temperature in Kelvin
for runs with and without Thermal Conduction at time $t=0.1$ Myr. The RKL runs
with thermal conduction use Spitzer conduction along field lines and no
conduction across the field lines. The black lines shown in each panel depict
the magnetic field lines. Since thermal conduction flux $\bm{F_{\rm c}}$ is
only along the field lines, this anisotropy leads to an asymmetric expansion
of the inner hot bubble along $z-$ axis. The evolution of temperature and
magnetic field lines in Figure 6 evidently shows that the diffusion in the
transverse direction is suppressed for the run with anisotropic thermal
conduction. While in absence of thermal conduction the outer shock maintains
its spherical structure.

Figure 6: Comparison of logarithmic value temperature (in K) for the 2D
magnetised blast wave test problem without (left panel) and with (right panel)
anisotropic thermal conduction. The RKL-2 method is used for the run with
conduction. The snapshot shown is at time $t=0.1$ Myr. The black lines
indicate the magnetic field lines.

Further, we compare the execution times for this test using the standard
explicit time-stepping and the RKL-2 method. This comparison for varying grid
resolutions is listed in table 2. The RKL-2 method is faster than the explicit
method by one order of magnitude for our highest resolution run with a grid
resolution of $1536^{2}$.

Table 2: A comparison of the execution time for the blast wave problem $N$ |  Wall-clock time $[{\rm s}]$
---|---
explicit | RKL-2
$192$ | $292$ | $112$
$384$ | $4448$ | $1276$
$768$ | $62700$ | $11828$
$1536$ | $1.11\times 10^{6}$ | $1.45\times 10^{5}$

This observation is consistent with the fact that STS techniques like RKL-2
method provide an asymptotic gain over standard explicit time-stepping methods
as demonstrated in the middle panel of Fig. 2 for the 1D Gaussian diffusion
test (Sec. 4.1.1).

####  4.2.2 Local thermal instability

In this section we describe the results from 2-D MHD simulations of local
thermal instability (TI) with anisotropic conduction (Sharma et al. 2010b).
Our computational domain is a 40 kpc $\times$ 40 kpc periodic Cartesian box
with a mean initial electron density of $n_{e}=0.1$ cm-3 (mean mass per
particle/electron, $\mu/\mu_{e}=0.62/1.17m_{p}$) and a uniform temperature of
0.7 keV (typical cool cluster core parameters). Only classical thermal
conduction with the conductivity given by the Spitzer value (Eq. 11 in Sharma
et al. 2010b) is included; i.e., $\bm{F}_{\rm c}=\bm{F}_{\rm class}$ instead
of Eq. 2. We use a tabulated cooling function corresponding to the plasma of a
third solar metallicity (Sutherland & Dopita 1993); the cooling function is
set to zero for $T<10^{6}$ K. Following Sharma et al. (2010b), we include a
spatially uniform heating rate density which is equal to the average cooling
rate density over the computational domain. With this, the computational box
is in global thermal balance, and this mimics the observed rough thermal
balance inferred in cool core clusters. Moreover, this setup shows the
exponential linear growth of the local thermal instability. With these
parameters, the initial cooling time (which is approximately equal to the
growth timescale of local TI; e.g., see Eq. 19c in McCourt et al. 2012) is
0.095 Gyr. A magnetic field of 5 $\mu$G is initialised at 45 degrees to the
Cartesian box. Homogeneous, isotropic, isobaric random density perturbations
are initialised to seed the local TI. The density perturbations
($\delta\equiv[\rho-\rho_{0}]/\rho_{0}$) are given by

$\delta(x,y)=\sum_{|k|=2,|l|=2}^{10}a_{k,l}\cos\left(\frac{2\pi(kx+ly)}{L}+\phi_{k,l}\right),$
(30)

where $k,l$ are mode labels, $a_{k,l}=1.5\times 10^{-3}r(k^{2}+l^{2})^{-1/2}$
and $\phi_{k,l,m}=2\pi r$ ($r$ is a random number uniformly distributed
between -0.5 to 0.5, which is different for the amplitude and phase and for
different modes), and $L=40$ kpc is the box size. These choices give a maximum
over-density amplitude ${\rm max}(\delta)\approx 0.003$. The setup is very
similar to but not identical as Sharma et al. (2010b).

We run the local TI test problem using different methods for anisotropic
thermal conduction: (a) fully explicit evolution in which both the MHD and
conduction modules are evolved using a time step $\Delta t_{\rm p}$ (Eq. 10;
this is typically shorter than the MHD CFL step $\Delta t_{\rm h}$; Eq. 9);
(b) sub-cycling of conduction module in which MHD module uses $\Delta t_{\rm
h}$ but the conduction module is sub-cycled and applied $\Delta t_{\rm
h}/\Delta t_{\rm p}$ times using a time step of $\Delta t_{\rm p}$; (c) MHD
module is evolved using $\Delta t_{\rm h}$ and conduction module is evolved
using AAG-STS with $\nu=0.01$ and $s\sim(\Delta t_{\rm h}/\Delta t_{\rm
p})^{1/2}$ stages (see Eq. 18); and (d) MHD module is evolved using $\Delta
t_{\rm h}$ and conduction module is evolved using RKL-2 with $s\sim(\Delta
t_{\rm h}/\Delta t_{\rm p})^{1/2}$ stages (see Eq. 17). The wall-time taken
for different methods to run the $512\times 512$ TI test problem until 0.87
Gyr (9.16 cooling times) is listed in Table 3. All the methods use the
monotonized-centered (MC) limiter to calculate the transverse terms in the
anisotropic heat flux ($F_{x}^{T}$ in Eq. 51 and analogous expression for
$F_{y}^{T}$). If we do not use limiters for interpolating transverse
temperature gradients (and instead use simple averaging), all the different
runs (a-d) blow up at some point in nonlinear evolution due to negative
temperature somewhere in the computational domain. This test highlights the
importance of using limiters for robustness in presence of large temperature
gradients.

Table 3: Wall-time for the local TI test run till 0.87 Gyr

explicit | sub-cycling |  AAG-STS† | RKL-2
---|---|---|---
6 h 1 m 29 s | 43 m 45 s | 9 m | 16 m 6 s

The grid resolution is $512\times 512$. † AAG-STS run crashes at 0.87 Gyr, but
others do not.

We note that the AAG-STS run crashes at 0.87 Gyr even with limiters because of
a large $C_{p}=\Delta t_{\rm h}/2\Delta t_{\rm p}$ (see Fig. 3). Other runs
could go for much longer without numerical problems. Both RKL-2 and AAG-STS
methods clearly show a significant speed-up relative to sub-cycling and
explicit methods. On comparing the STS methods, we see that AAG-STS is faster
than RKL-2 due to a smaller number of computations per cycle. However, the
first-order accurate AAG-STS scheme exhibits unstable behavior while RKL-2
maintains stability during the entire integration as described below.

Figure 7: Contour plots of log10 temperature (keV) for RKL-2 (top panels) and
AAG-STS (bottom panels) runs at linear (0.30 Gyr) and nonlinear (0.45 Gyr and
0.86 Gyr) stages of the instability. The arrows shows the local magnetic field
unit vector. The non-linear stage starts at around 0.45 Gyr and the magnetic
field lines no longer remain aligned at 450 to the x-axis. The evolution of TI
differs substantially with the two methods at late times in the non-linear
phase. The AAG-STS run crashes at 0.87 Gyr.

Figure 7 shows temperature snapshots at different stages of TI evolution using
RKL-2 (top panels) and AAG-STS (bottom panel) schemes. In the linear stage all
schemes give a similar evolution. In the saturated nonlinear state AAG-STS and
RKL-1 start to deviate quantitatively from each other. Unlike AAG-STS, the
RKL-2 temperature snapshots at 0.86 Gyr are very similar to the ones obtained
using explicit and sub-cycling methods (not shown in Fig. 7); this time is
very close to the time when the AAG-STS run blows up (at 0.87 Gyr) due to
numerical instability.

Figure 8: The evolution of various quantities like kinetic energy density
[ergs cm-3], magnetic energy density [ergs cm-3], the minimum and maximum
temperature in the computational domain ($T_{\rm min}$ and $T_{\rm max}$) with
time using different methods for anisotropic conduction to study the local TI
. Explicit (red solid line), sub-cycling (green solid line), and RKL-2 (black
solid line) methods show similar time evolution but AAG-STS (blue dashed line)
starts to deviate in the non-linear stage. Notice the numerous spikes in
$T_{\rm max}$ and $T_{\rm min}$ for AAG-STS, which are symptoms of numerical
instability which blows up the code at 0.87 Gyr. We impose a numerical
temperature floor when the temperature becomes negative.

Figure 8 shows the time evolution of various quantities as a function of time
using different methods for anisotropic thermal conduction. Top left panel of
Figure 8 shows that the evolution of kinetic energy in the box is similar for
all the runs. The top right panel shows the average magnetic energy evolution.
Here, AAG-STS deviates from the other runs in the non-linear stage. Similar
deviations are seen in AAG-STS for $T_{\rm min}$ and $T_{\rm max}$ (minimum
and maximum temperature in the computational domain), but not for RKL-2. The
numerical fragility of AAG-STS due to an imperfect stability parameter $\nu$
is evident from spikes in $T_{\rm min}$ and $T_{\rm max}$. The AAG-STS run
blows up at 0.87 Gyr due to the numerical instability of AAG-STS for a large
number of sub-stages.

##  5 Parallel scaling of RKL conduction

In this section, we demonstrate the parallel scaling of isotropic and
anisotropic thermal conduction in PLUTO code. For this purpose, we use the
blast wave test problem in three dimensions on a Cartesian grid with size
$L=30$ pc. To set up a blast wave we initialise an ambient static medium with
density $\rho_{0}$ and pressure $p_{0}$. A blast region with high density and
pressure is set within a spherical radius $r_{0}$ around the origin. The
density and pressure in this region is a factor 10 and 1000 times larger than
the ambient medium respectively. We smooth the transition between the blast
region and the ambient medium using a smoothing function as follows

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\rho_{0}(1+9\Lambda[r])$
$\displaystyle p$ $\displaystyle=$ $\displaystyle
p_{0}(1+999\Lambda[r])/\gamma$ (31)

where, $\Lambda[r]=(1+\exp([r-r_{0}]/[0.1~{}{\rm pc}]))^{-1}$ is a smoothing
function, $r=(x^{2}+y^{2}+z^{2})^{1/2}$ is the spherical radius,
$p_{0}=2.295\times 10^{-9}$ dyn cm-2, $\rho_{0}=2.1\times 10^{-23}$ g cm-3,
$r_{0}=1.0$ pc, and $\gamma=5/3$. Additionally, for MHD runs the initial
magnetic field is along $y$ axis, $\bm{B}=B_{0}\hat{\bm{y}}$ where
$B_{0}=\sqrt{2p/\beta}$ with $\beta=1$. Thus, the field strength inside the
high pressure blast region is about 30 times larger than the ambient medium.
We use Spitzer conduction along field lines and no conduction across the field
lines for the MHD run. The HD run uses isotropic Spitzer conduction. For
comparison, we also include a HD run without conduction.

We carried out strong scaling studies using a $512^{3}$ Cartesian grid on CPU-
only nodes of IISc Cray XC40 cluster
SahasraT444http://www.serc.iisc.in/facilities/cray-xc40-named-as-sahasrat/.
Unlike the 2-D blast wave test problem in section 4.2.1, we choose a 3-D blast
wave problem for testing parallel strong scaling on more than ten thousand
processors. In strong scaling, the problem size remains fixed but the number
of processors is increased progressively. A large enough problem size is
needed (that is why 3-D instead of 2-D test problem is chosen) for
communication not to dominate over computation even on largest number of
processors. We run the blast wave test problem for $t=0.5$ kyr on processors
ranging from 256 to 16384. On average, the number of sub-steps $s$ ranges
between 5-12 for the HD run and between 3-8 for MHD.

Figure 9: Strong scaling test for the 3D blast wave test with 5123 cells. The
left panel shows the wall clock time for the simulations run till $t=0.5$ kyr
using different number of processors. The right panel shows the comparison of
numerical speed-up with the ideal scaling (black dashed line). Scaling study
using pure HD (red circles), HD+RKL (green triangles) and MHD+RKL (magneta
stars) shows a close to perfect scaling up to 104 processors.

The left panel of Figure 9 shows the wall clock time as a function of
processor count for our strong scaling studies. The wall clock time behaves as
expected, with the MHD-RKL run taking longer than the HD-RKL run. Both runs
with conduction show close to the inverse linear scaling with the processor
count. The right panel of Figure 9 shows that the speed-up improves with
increasing core count, but drops to 80% for more than $10^{4}$ processors. The
same trend is also seen for pure HD run, indicating that the addition of
thermal conduction has not resulted in any performance degradation. Our strong
scaling tests show that RKL methods yield scalable algorithm on PetaFlop
facilities.

This feature is usually not shared by implicit methods which, as already
mentioned in Sec. 3.2, require inverting large sparse matrices, an operation
difficult to achieve efficiently on multi-core systems (Botchev & van der
Vorst, 2001). Nevertheless, the debate on which approach can be more efficient
in massive parallel computations has yet to be settled. In the work by Pakmor
et al. (2016), for instance, a semi-implicit scheme is employed to solve the
cosmic-ray transport equation (coupled to the MHD equations) on unstructured
moving mesh in the context of galaxy dynamics. Their semi-implicit solver
requires solving a single linear system of equations per timestep and weak
scaling tests claim good parallel efficiency, up to 480 cores. However, more
complex scenarios including nonlinear systems of equations (e.g Ohmic
diffusion, ambipolar diffusion, etc.) may introduce additional complexities
making the applicability of implicit scheme less efficient. Although a
comparison with implicit methods is outside the scope of this paper, we point
out to the recent work by Caplan et al. (2017), where explicit STS schemes are
compared with implicit Krylov solvers in the context of magnetised solar
corona. The left panels of Fig. 11 in their paper shows the superiority of
explicit schemes over implicit ones for their choice of problem and model. The
comparison also shows that the strong parallel scaling of implicit methods
saturates and deviates sharply from ideal behavior at large number of cores
($>1000$).

The explicit RKL schemes presented in this work confirm this prediction (the
scaling performance with explicit RKL method shown in Fig. 9 is able to
achieve a high efficiency of 80% even up to 1000 cores) and can be naturally
extended to more complex systems of equations with with minimal modifications.

##  6 Conclusions

In this paper we have discussed various numerical methods implementing
anisotropic thermal conduction coupled with the standard set of MHD equations.
In particular, we have described the second order (in time) accurate Runge-
Kutta Legendre (RKL-2) super-time stepping (STS) method implemented in PLUTO
code. We have then compared these numerical methods on simple test problems,
and also on astrophysical test problems like blast wave and thermal
instability in which MHD evolution is coupled with conduction. The major
conclusions of our paper are:

  1. 1.

Using the 1-D Gaussian diffusion test, we show that the RKL-2 scheme is second
order accurate in time, in comparison to the standard AAG-STS (Alexiades et
al. 1996) scheme which is first order.

  2. 2.

Super-time stepping schemes based on Chebyshev polynomials such as AAG-STS
require an ad-hoc damping parameter ($\nu$), which has to be big enough for
numerical stability. Figure 3 shows that AAG-STS becomes unstable for number
of substages $s\gtrsim$ few 10s. Moreover, for $2s\sqrt{\nu}\gtrsim 1$ there
is no real speed-up using AAG-STS (see Eq. 20). In absence of such a parameter
and with sufficient inherent damping, RKL-STS schemes are more robust and can
better exploit the super-time stepping strategy (see the right panel of Fig.
2). The robustness of RKL-STS is also useful for implementing saturated
conduction (e.g., see Fig. 4).

  3. 3.

The Cartesian ring diffusion test problem (Fig. 5) shows that all STS schemes
break down if we use a large number of stages with low resolution. Moreover,
the use of limiters to interpolate the transverse temperature gradient leads
to somewhat larger perpendicular diffusion. Limiters do help prevent non-
monotonicity of temperature in presence of temperature discontinuities. For
practical implementation of anisotropic diffusion, limiters need not be used
as a default option, but may be recommended at locations and times at which
the temperature falls below a floor value. The local thermal instability test
problem in section 4.2.2 demonstrates the utility of limiters in presence of
temperature discontinuities.

  4. 4.

As demonstrated in Tables 1, 2 & 3, and also indicated by the middle panel of
Figure 2, the run-time with super-time stepping is substantially shorter than
with explicit update or even sub-cycling.

  5. 5.

Last, but perhaps most importantly as shown in Figure 9, the explicit schemes
such as RKL-STS shows an excellent scaling (efficiency of $\sim 80\%$ up to
104 processors) on modern distributed PetaScale supercomputers.

## Acknowledgments

The authors would sincerely like to thank the referee for the valuable
comments which have played a significant role in improving the paper. BV would
like to thank the support provided by University of Torino and also would like
to acknowledge the hospitality of IISc during the research visit in 2016. This
work is partly supported by the DST-India grant no. Sr/S2/HEP-048/2012 and an
India-Israel joint research grant (6-10/2014[IC]). We thank the SERC-IISc
staff for facilitating our use of SahasraT cluster for the parallel scaling
runs; only with their intervention we could finish our runs in a reasonable
time. PS acknowledges the hospitality of KITP where this paper was completed.
This research was supported in part by the National Science Foundation under
Grant No. NSF PHY-1125915.

##  Appendix A Illustration of RKL-2 scheme

This Appendix is closely based on Meyer et al. (2012, 2014) and is included
here for completeness. Second order accuracy for the RKL-2 scheme can be
achieved by matching the first three terms in Eqs. (14) & (15); i.e., by
imposing $R_{s}(0)$ = 1, $R_{s}^{{}^{\prime}}(0)$ = 1 and
$R_{s}^{{}^{\prime\prime}}(0)$ = 1. These three requirements can be used to
estimate the values of three coefficients in Eq.(15) as,

$\displaystyle a_{s}$ $\displaystyle=$ $\displaystyle 1-b_{s}$ (32)
$\displaystyle b_{s}$ $\displaystyle=$
$\displaystyle\frac{\mathcal{P}^{{}^{\prime\prime}}_{s}(1)}{(\mathcal{P}^{{}^{\prime}}_{s}(1))^{2}}=\frac{s^{2}+s-2}{2s(s+1)}\,\,\,\,\,\,\,s\geq
2$ (33) $\displaystyle w_{1}$ $\displaystyle=$
$\displaystyle\frac{\mathcal{P}^{{}^{\prime}}_{s}(1)}{\mathcal{P}^{{}^{\prime\prime}}_{s}(1)}=\frac{4}{s^{2}+s-2}\,\,\,\,\,\,\,s\geq
2$ (34)

For $s<2$, we can choose $b_{0}=b_{1}=b_{2}=1/3$ (Meyer et al., 2014). Here we
have used various properties of the Legendre polynomials and their
derivatives, i.e.,

$\displaystyle\mathcal{P}_{s}\left(1\right)$ $\displaystyle=$ $\displaystyle
1$ (35) $\displaystyle\mathcal{P}^{{}^{\prime}}_{s}(1)$ $\displaystyle=$
$\displaystyle\frac{s(s+1)}{2}$ (36)
$\displaystyle\mathcal{P}^{{}^{\prime\prime}}_{s}(1)$ $\displaystyle=$
$\displaystyle\left(\frac{s^{2}+s-2}{4}\right)\mathcal{P}^{{}^{\prime}}_{s}(1)$
(37)

Like Eq. (15) for the stability polynomial at the end of $s$ stages, the
stability polynomial till $j$ stages is chosen to be
$a_{j}+b_{j}\mathcal{P}_{j}\left(1+w_{1}(\tau\lambda)\right)$.

The RKL-2 scheme with $s$ stages can be expressed as,

$\displaystyle Y_{0}$ $\displaystyle=$ $\displaystyle u(t_{0})$ $\displaystyle
Y_{1}$ $\displaystyle=$ $\displaystyle
Y_{0}+\tilde{\mu}_{1}\tau\mathsf{M}Y_{0}$ $\displaystyle Y_{j}$
$\displaystyle=$
$\displaystyle\mu_{j}Y_{j-1}+\nu_{j}Y_{j-2}+(1-\mu_{j}-\nu_{j})Y_{0}$
$\displaystyle+\tilde{\mu}_{j}\tau\mathsf{M}Y_{j-1}+\tilde{\gamma}_{j}\tau\mathsf{M}Y_{0}\hskip
45.52458pt2\leq j\leq s$ $\displaystyle u(t_{0}+\tau)$ $\displaystyle=$
$\displaystyle Y_{s}$ (38)

where the coefficients, $\mu_{j}$, $\tilde{\mu}_{1}$, $\tilde{\mu}_{j}$,
$\nu_{j}$ and $\tilde{\gamma}_{j}$ can be obtained using the recursion
relation for the Legendre polynomials (see Eq.(16)) and rearrangement of
terms. The expressions for these coefficients are given by,

$\displaystyle\mu_{j}$ $\displaystyle=$
$\displaystyle\frac{2j-1}{j}\frac{b_{j}}{b_{j-1}}$ (39)
$\displaystyle\tilde{\mu}_{1}$ $\displaystyle=$ $\displaystyle b_{1}w_{1}$
(40) $\displaystyle\tilde{\mu}_{j}$ $\displaystyle=$
$\displaystyle\mu_{j}w_{1}$ (41) $\displaystyle\nu_{j}$ $\displaystyle=$
$\displaystyle-\frac{j-1}{j}\frac{b_{j}}{b_{j-2}}$ (42)
$\displaystyle\tilde{\gamma}_{j}$ $\displaystyle=$ $\displaystyle-
a_{j-1}\tilde{\mu}_{j}$ (43)

To get a better sense of the scheme, we explicitly list the value of the above
coefficients for a representative small value of $s=3$. For this three stage
scheme we have,

$\displaystyle\mu_{2}$ $\displaystyle=$
$\displaystyle\frac{3}{2};\,\,\,\mu_{3}=\frac{25}{12}$
$\displaystyle\tilde{\mu}_{1}$ $\displaystyle=$
$\displaystyle\frac{2}{15};\,\,\,\tilde{\mu}_{2}=\frac{3}{2}w_{1}=\frac{3}{5};\,\,\,\tilde{\mu}_{3}=\frac{25}{12}w_{1}=\frac{5}{6}$
$\displaystyle\nu_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2};\,\,\,\nu_{3}=-\frac{2}{3}\frac{b_{3}}{b_{1}}=-\frac{5}{6}$
$\displaystyle\tilde{\gamma}_{2}$ $\displaystyle=$
$\displaystyle(b_{1}-1)\tilde{\mu}_{2}=-\frac{2}{5};\,\,\,\tilde{\gamma}_{3}=(b_{2}-1)\tilde{\mu}_{3}=-\frac{5}{9}$

On substituting the above values of various coefficients for $s=3$ in Eq.(A)
we get;

$\displaystyle Y_{0}$ $\displaystyle=$ $\displaystyle u(t_{0})$ (44)
$\displaystyle Y_{1}$ $\displaystyle=$ $\displaystyle
Y_{0}+\frac{2}{15}\tau\mathsf{M}Y_{0}$ $\displaystyle=$ $\displaystyle
Y_{0}+\frac{2}{15}\tau\lambda Y_{0}$ $\displaystyle=$
$\displaystyle\left(1+\frac{2}{15}(\tau\lambda)\right)Y_{0}$ $\displaystyle=$
$\displaystyle\left(\frac{2}{3}+\frac{1}{3}\mathcal{P}_{1}\left(1+\frac{2}{5}(\tau\lambda)\right)\right)Y_{0}$
$\displaystyle=$
$\displaystyle\left(a_{1}+b_{1}\mathcal{P}_{1}\left(1+w_{1}(\tau\lambda)\right)\right)Y_{0}=R_{1}(\tau\lambda)Y_{0}$
$\displaystyle Y_{2}$ $\displaystyle=$
$\displaystyle\frac{3}{2}Y_{1}-\frac{1}{2}Y_{0}+\frac{3}{5}\tau\mathsf{M}Y_{1}-\frac{2}{5}\tau\mathsf{M}Y_{0}$
$\displaystyle=$
$\displaystyle\left(\frac{3}{2}+\frac{3}{5}(\tau\lambda)\right)Y_{1}-\left(\frac{1}{2}+\frac{2}{5}(\tau\lambda)\right)Y_{0}$
$\displaystyle=$
$\displaystyle\left(1+\frac{2}{5}(\tau\lambda)+\frac{2}{25}(\tau\lambda)^{2}\right)Y_{0}$
$\displaystyle=$
$\displaystyle\left(\frac{2}{3}+\frac{1}{3}\mathcal{P}_{2}\left(1+\frac{2}{5}(\tau\lambda)\right)\right)Y_{0}=R_{2}(\tau\lambda)Y_{0}$
$\displaystyle Y_{3}$ $\displaystyle=$
$\displaystyle\frac{25}{12}Y_{2}-\frac{5}{6}Y_{1}-\frac{1}{4}Y_{0}+\frac{5}{6}\tau\mathsf{M}Y_{2}-\frac{5}{9}\tau\mathsf{M}Y_{0}$
$\displaystyle=$
$\displaystyle\left(1+(\tau\lambda)+\frac{1}{2}(\tau\lambda)^{2}+\frac{1}{15}(\tau\lambda)^{3}\right)Y_{0}$
$\displaystyle=$
$\displaystyle\left(\frac{7}{12}+\frac{5}{12}\mathcal{P}_{3}\left(1+\frac{2}{5}(\tau\lambda)\right)\right)Y_{0}=R_{3}(\tau\lambda)Y_{0}$
$\displaystyle u(t_{0}+\tau)$ $\displaystyle=$ $\displaystyle Y_{3}$ (48)

where $\lambda$ is the eigenvalue of parabolic operator matrix $\mathsf{M}$.
Note that the first three terms of final solution of $Y_{3}$ (Eq. A) match the
Taylor expansion of the exponential function (Eq.(14). Thus this scheme is
second order-accurate with a leading order in error given by the fourth term,
i.e., $(\tau\lambda)^{3}/15$.

Finally to estimate the value of the super-step $\tau$, we should have
$|R_{3}(\tau\lambda)|\leq 1$. The stability condition will be satisfied only
if the argument of LP is $\geq-1$ (note that $\lambda$ is non-positive), i.e.,
$-5\leq\tau\lambda\leq 0$. While the similar stability condition for one step
Euler-scheme requires $-2\leq\tau\lambda\leq 0$. This clearly shows the $s=3$
stage RKL-2 method allows to choose a large time-step as compared to the Euler
scheme. Higher values of sub-steps $s$, will result in a wider range of stable
time steps giving an obvious advantage over standard explicit schemes.

##  Appendix B Limiting transverse temperature gradients for monotonicity

Simple finite differencing of Eq. 6 accentuates temperature extrema in
presence of large temperature gradients expected naturally in astrophysical
plasmas (e.g., Sharma & Hammett 2007). The heat flux can be decomposed into
normal and transverse components (Eqs. 50, 51). Expressing the heat flux ,
$\bm{F}_{\rm c}$ with components normal and transverse to the field lines we
obtain the following explicit formulation of Eq. 6 in 2D Cartesian geometry
(generalisation to 3-D and non-Cartesian coordinates is straightforward),

$\frac{3}{2}nk_{B}\frac{\partial T}{\partial t}=-\frac{\partial}{\partial
x}(F_{\rm c,x}^{N}+F_{\rm c,x}^{T})-\frac{\partial}{\partial y}(F_{\rm
c,y}^{N}+F_{\rm c,y}^{T}).$ (49)

For simplicity we consider the classical limit (i.e., ignore saturated
conduction). On comparing Eq. 49 with Eq. 3, we get

$\displaystyle F_{\rm c,x}^{N}$ $\displaystyle=$
$\displaystyle-\kappa_{\parallel}b_{x}^{2}\frac{\partial T}{\partial x},$ (50)
$\displaystyle F_{\rm c,x}^{T}$ $\displaystyle=$
$\displaystyle-\kappa_{\parallel}b_{x}b_{y}\frac{\partial T}{\partial y},$
(51)

as the normal and transverse components of the heat flux (analogous
expressions can be constructed for $F_{\rm c,y}^{N}$ and $F_{\rm c,y}^{T}$).

The normal component of the heat flux is naturally face-centred and always
carries heat from higher to lower temperatures. However, the transverse
temperature gradients are not naturally located at the faces, and need to be
interpolated there. As the transverse flux can have any sign (Eq. 51), without
special treatment, this component can lead to negative temperatures in regions
with large temperature gradients (this also applies to saturated conduction in
Eqs. 2 & 4; therefore, the transverse temperature gradient required to
evaluate ${\rm sgn}[\bm{\hat{b}}\cdot\nabla T$] for saturated flux should also
use limiters for robustness). Sharma & Hammett (2007) introduced limiters
(similar to those used in the reconstruction step in finite volume methods;
LeVeque 2002) to interpolate the transverse temperature gradients at the cell
faces and showed that the resulting explicit scheme preserves temperature
extrema.

## References

  * Alexiades et al. (1996) Alexiades V., Amiez G., Gremaud P.-A., 1996, Communications in numerical methods in engineering, 12, 31
  * Balbus (2000) Balbus S. A., 2000, ApJ, 534, 420
  * Balsara (2001) Balsara D. S., 2001, Journal of Computational Physics, 174, 614
  * Balsara et al. (2008) Balsara D. S., Tilley D. A., Howk J. C., 2008, MNRAS, 386, 627
  * Beresnyak (2011) Beresnyak A., 2011, Physical Review Letters, 106, 075001
  * Botchev & van der Vorst (2001) Botchev M. A., van der Vorst H. A., 2001, Journal of Computational and Applied Mathematics, 137, 229
  * Caplan et al. (2017) Caplan R. M., Mikić Z., Linker J. A., Lionello R., 2017, in Journal of Physics Conference Series. p. 012016 (arXiv:1610.01265), doi:10.1088/1742-6596/837/1/012016
  * Cowie & McKee (1977) Cowie L. L., McKee C. F., 1977, ApJ, 211, 135
  * Dong & Stone (2009) Dong R., Stone J. M., 2009, ApJ, 704, 1309
  * Dubois & Commerçon (2016) Dubois Y., Commerçon B., 2016, A&A, 585, A138
  * Evans & Hawley (1988) Evans C. R., Hawley J. F., 1988, ApJ, 332, 659
  * Federrath et al. (2011) Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, ApJ, 731, 62
  * Field (1965) Field G. B., 1965, ApJ, 142, 531
  * Günter et al. (2005) Günter S., Yu Q., Krüger J., Lackner K., 2005, Journal of Computational Physics, 209, 354
  * Hawley et al. (1995) Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742
  * Huba (2004) Huba J. D., 2004, NRL: Plasma Formulary, Naval Research Laboratory, Washington, DC 20375-5320
  * Kannan et al. (2015) Kannan R., Springel V., Pakmor R., Marinacci F., Vogelsberger M., 2015, arXiv preprint arXiv:1512.03053
  * Korpi et al. (1999) Korpi M. J., Brandenburg A., Shukurov A., Tuominen I., Nordlund Å., 1999, ApJ, 514, L99
  * LeVeque (2002) LeVeque R. J., 2002, Finite volume methods for hyperbolic problems. Vol. 31, Cambridge university press
  * McCourt et al. (2012) McCourt M., Sharma P., Quataert E., Parrish I. J., 2012, MNRAS, 419, 3319
  * Meyer et al. (2012) Meyer C. D., Balsara D. S., Aslam T. D., 2012, Monthly Notices of the Royal Astronomical Society, 422, 2102
  * Meyer et al. (2014) Meyer C. D., Balsara D. S., Aslam T. D., 2014, Journal of Computational Physics, 257, 594
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Mignone et al. (2012) Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS, 198, 7
  * O’Sullivan & Downes (2007) O’Sullivan S., Downes T. P., 2007, ] 10.1111/j.1365-2966.2007.11429.x, 376, 1648
  * Pakmor et al. (2016) Pakmor R., Pfrommer C., Simpson C. M., Kannan R., Springel V., 2016, MNRAS, 462, 2603
  * Parrish & Stone (2005) Parrish I. J., Stone J. M., 2005, ApJ, 633, 334
  * Press et al. (1986) Press W. H., Flannery B. P., Teukolsky S. A., 1986, Numerical recipes. The art of scientific computing
  * Quataert (2008) Quataert E., 2008, ApJ, 673, 758
  * Sharma & Hammett (2007) Sharma P., Hammett G. W., 2007, Journal of Computational Physics, 227, 123
  * Sharma & Hammett (2011) Sharma P., Hammett G. W., 2011, Journal of Computational Physics, 230, 4899
  * Sharma et al. (2009) Sharma P., Chandran B. D. G., Quataert E., Parrish I. J., 2009, ApJ, 699, 348
  * Sharma et al. (2010a) Sharma P., Colella P., Martin D. F., 2010a, SIAM Journal on Scientific Computing, 32, 3564
  * Sharma et al. (2010b) Sharma P., Parrish I. J., Quataert E., 2010b, ApJ, 720, 652
  * Sutherland & Dopita (1993) Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
  * Tóth (2000) Tóth G., 2000, Journal of Computational Physics, 161, 605
  * Verwer et al. (1990) Verwer J., Hundsdorfer W., Sommeijer B., 1990, Numerische Mathematik, 57, 157
  * Wagh et al. (2014) Wagh B., Sharma P., McCourt M., 2014, MNRAS, 439, 2822
  * Weickert (1998) Weickert J., 1998, Anisotropic diffusion in image processing. Vol. 1, Teubner Stuttgart
  * Yang & Reynolds (2015) Yang H.-Y., Reynolds C. S., 2015, arXiv preprint arXiv:1512.05796
  * van Der Houwen & Sommeijer (1980) van Der Houwen P. J., Sommeijer B. P., 1980, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 60, 479