# How numerical treatments of the transition region modify energy flux into
the solar corona

T. A. Howson,1,2 C. Breu,2
1School of Design and Informatics, Abertay University, Bell Street, Dundee,
DD1 1HG, U.K.
2School of Mathematics and Statistics, University of St Andrews, St Andrews,
Fife, KY16 9SS, U.K.
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The large temperature gradients in the solar transition region present a
significant challenge to large scale numerical modelling of the Sun’s
atmosphere. In response, a variety of techniques have been developed which
modify the thermodynamics of the system. This sacrifices accuracy in the
transition region in favour of accurately tracking the coronal response to
heating events. Invariably, the modification leads to an artificial broadening
of the transition region. Meanwhile, many contemporary models of the solar
atmosphere rely on tracking energy flux from the lower atmosphere, through the
transition region and into the corona. In this article, we quantify how the
thermodynamic modifications affect the rate of energy injection into the
corona. We consider a series of one-dimensional models of atmospheric loops
with different numerical resolutions and treatments of the thermodynamics.
Then, using Alfvén waves as a proxy, we consider how energy injection rates
are modified in each case. We find that the thermodynamic treatment and the
numerical resolution significantly modify Alfvén travel times, the
eigenfrequencies and eigenmodes of the system, and the rate at which energy is
injected into the corona. Alarmingly, we find that the modification of the
energy flux is frequency dependent, meaning that it may be difficult to
compare the effects of different velocity drivers on coronal heating if they
are imposed below an under-resolved transition region, even if the
sophisticated thermodynamic adaptations are implemented.

###### keywords:

Sun: oscillations - Sun: corona – Sun: transition region

††pubyear: 2023††pagerange: How numerical treatments of the transition region
modify energy flux into the solar corona–How numerical treatments of the
transition region modify energy flux into the solar corona

##  1 Introduction

The solar corona is the largest and hottest layer of the Sun’s atmosphere and
is the subject of one of the great outstanding questions in astrophysics; the
coronal heating problem. This concerns how the surprisingly high temperatures
of the corona are maintained against energy losses, including thermal
conduction to the cooler layers below and radiative losses into Space. The
complexity of the solar atmosphere has ensured that this problem has resisted
decades of sustained effort including sophisticated observational, analytical
and numerical studies. It is widely accepted that the required energy is
ultimately injected into the atmosphere by complex convective motions at the
solar surface and a wide variety of contemporary models show how a hot corona
can be sustained as a result (e.g. Gudiksen & Nordlund, 2005; Bingert & Peter,
2011; Reale et al., 2016; Kanella & Gudiksen, 2018; Breu et al., 2022;
Kuniyoshi et al., 2023; Reid et al., 2023). However, the specific nature of
the energy dissipation mechanisms remains hotly contested. Thorough reviews of
contemporary research in this area are presented by (e.g. Reale, 2014;
Klimchuk, 2015; Van Doorsselaere et al., 2020; Viall et al., 2021).

In recent decades, significant growth in computational power has enabled the
solar atmosphere to be modelled with large scale, high-resolution MHD
simulations. Increasingly, these simulations include the full,
gravitationally-stratified atmosphere with the different layers of the
atmosphere considered within a single numerical model (e.g. Hansteen et al.,
2015; Cheung et al., 2019; Howson & De Moortel, 2022; Robinson et al., 2022;
Chen et al., 2023; Guo et al., 2023). As each layer is associated with
distinct physical processes occurring on disparate spatial and temporal
scales, incorporating the complete atmosphere within one simulation represents
a significant numerical challenge. Despite this, contemporary models attempt
to track the flux of energy and mass from convective layers at the base of the
simulation volume, through a complex and dynamic chromosphere and into the
corona, which is heated to realistic temperatures (e.g. simulations produced
with the MuRAM and Bifrost codes Vögler et al., 2005; Gudiksen et al., 2011).
The success of these models can then be assessed by generating synthetic
emission from the simulation results for comparison with real observations. As
emission from the corona is sensitive to the plasma density and temperature,
the exchange of mass between the atmospheric layers and energy dissipation
rates are important components of coronal heating models.

One particularly challenging aspect of solar atmospheric modelling concerns
the transition region. This is a narrow layer of the solar atmosphere that
sits at the interface between the relatively cool and dense chromosphere and
the hot and tenuous corona. Over the transition region, the plasma temperature
increases by more than two orders of magnitude, over a short distance. As a
result, there are very large temperature gradients which present a
considerable problem for the finite difference schemes implemented within MHD
solar atmospheric codes. Under-resolving these gradients can significantly
impair the accuracy of simulations, including vastly under-estimating the
upflow of plasma into the corona following heating events Bradshaw & Cargill
(2013) and artificially suppressing thermal non-equilibrium cycles Johnston et
al. (2019). A potentially naïve solution may be to simply increase the number
of grid points used by the numerical schemes such that the transition region
temperature gradients remain well-resolved. However, whilst this can be a
prudent strategy in 1-D codes, the computational cost in three dimensions is
prohibitive. In response, several numerical techniques have been developed
(e.g. Lionello et al., 2009; Mikić et al., 2013; Johnston et al., 2017;
Johnston & Bradshaw, 2019; Johnston et al., 2021). These are described in more
detail in Section 2, however they generally work by adapting the effects of
the transition region in order to accurately model the coronal response to
heating events (e.g. density enhancement due to the evaporation of
chromospheric and transition region plasma Hirayama, 1974; Fisher et al.,
1985; Tian & Chen, 2018) even with relatively coarse numerical resolution.
These techniques are designed to correctly track the evolution of plasma in
the corona, such that synthetic emission can be generated from heating models
for direct comparison with observations (e.g. Antolin et al., 2016; Pontin et
al., 2017; Kanella & Gudiksen, 2019; Warren et al., 2020). However, they also
broaden the transition region and thus modify the flux of other forms of
energy (e.g. Poynting flux) from the lower atmosphere into the corona. By
modifying this energy flux, the artificial TR broadening may have unintended
consequences for large scale coronal heating simulations.

As these methods are being increasingly implemented in multi-dimensional
models of the solar atmosphere (e.g. Van Damme et al., 2020; Zhou et al.,
2021; Howson & De Moortel, 2022; Keppens et al., 2023; Li et al., 2023b;
Pelouze et al., 2023), it is now essential to quantify how coronal energy
injection rates are affected by these thermodynamic treatments. To this end,
in this paper, we consider how simple propagating Alfvén waves are modified as
they propagate through different simulated transition regions. By using these
waves as a proxy for more complex atmospheric dynamics, we will discuss the
frequency and resolution-dependent effects on energy transmission. This will
allow us to estimate how mechanical energy flux is affected by transition
region modifications in larger and more complex models. The remainder of the
article is presented as follows: In section 2, we describe our simple models,
in section 3, we describe our results and in section 4, we discuss the
implications for contemporary modelling of the fully coupled solar atmosphere.

##  2 Numerical methods

For the majority of the simulations conducted within this article, we used the
Lare2D code (Arber et al., 2001; Arber, 2018). The code advances the full,
non-ideal and non-linear MHD equations given in normalised form by:

$\frac{\text{D}\rho}{\text{D}t}=-\rho\nabla\cdot{\bf{v}},$ (1)
$\rho\frac{{\text{D}{\bf v}}}{{\text{D}t}}={\bf j}\times{\bf B}-\nabla
P-\rho{\bf g}+{\bf F}_{\text{visc.}},$ (2)
$\rho\frac{{\text{D}\epsilon}}{{\text{D}t}}=-P(\nabla\cdot{\bf
v})-\nabla\cdot{\bf q}-\rho^{2}\Lambda(T)+\eta\lvert{\bf
j}\rvert^{2}+Q_{\text{visc.}}+Q_{\text{bg.}},$ (3) $\frac{\text{D}{\bf
B}}{\text{D}t}=\left({\bf B}\cdot\nabla\right){\bf v}-\left(\nabla\cdot{\bf
v}\right){\bf B}-\nabla\times\left(\eta\nabla\times{\bf B}\right),$ (4)
$P=2k_{B}nT.$ (5)

Here $\rho$ is the plasma density, ${\bf v}$ is the velocity, j is the current
density, B is the magnetic field, $P$ is the gas pressure, g is the
gravitational acceleration, $\epsilon$ is the specific internal energy
density, $\eta$ is the resistivity, $k_{B}$ is the Boltzmann constant, $n$ is
the number density and $T$ is the temperature. For numerical stability, small
shock viscosity terms are included, which contribute a frictional force to the
equation of motion (2) and an associated small heating term to the energy
equation (3). These are described in detail in (Arber, 2018; Reid et al.,
2020). By testing different dissipation coefficients, we confirmed that the
viscous effects are small (for these transport coefficients), such that any
wave damping is effectively negligible.

The energy equation includes contributions from thermal conduction,
$\nabla\cdot{\bf q}$, optically thin radiation, $\rho^{2}\Lambda(T)$, and a
background heating, $Q_{\text{bg.}}$. The radiative loss curve is described in
detail by Klimchuk et al. (2008) and the background heating term is
implemented to maintain an initial equilibirum. The magnitude of this heating
term is discussed in Sect. 2.1. The vector, ${\bf q}$, represents the heat
flux and is defined according to the Braginskii (1965) model for thermal
conduction in a magnetised plasma. In particular, q is given by

${\bf q}=\left({\bf k}\cdot\nabla T\right){\bf
n}+\nabla\cdot\left(\frac{b^{2}_{\text{min}}}{B^{2}+b^{2}_{\text{min}}}\right)\kappa\nabla
T,$ (6)

where ${\bf n}=\bf B/(B^{2}+b^{2}_{\text{min}})$ is parallel to the magnetic
field and $k=\kappa\,{\bf n}$ with $\kappa=\kappa_{0}T^{5/2}$ and
$\kappa_{0}=10^{-11}\text{ Jm}^{-1}\text{ K}^{-7/2}\text{ s}^{-1}$. The
constant $b_{\min}$ is defined as a small constant to avoid numerical issues
at magnetic null points $({\bf B}={\bf 0})$ and in the limit $B^{2}\gg
b^{2}_{\text{min}}$, equation 6 recovers the Spitzer-Härm parallel
conductivity, with efficient heat conduction along field lines but negligible
energy transfer across them. Although, no magnetic null points occur within
our simulations, for completeness, we note that if ${\bf B}\to{\bf 0}$,
isotropic conduction would be recovered.

The second spatial derivatives of the plasma temperature required for
advancing the energy equation 3 are particularly problematic in the transition
region, where the temperature changes very rapidly with height. As discussed
in Sect. 1, spatially under-resolving this region in coronal heating models
has significant consequences for plasma evolution. In this article, we
consider two numerical techniques that allow more accurate modelling of
coronal plasma even at lower resolutions. These are the L09 method (Lionello
et al., 2009; Mikić et al., 2013) and the TRAC approach (Johnston & Bradshaw,
2019; Johnston et al., 2021). The L09 method works by modifying the
thermodynamics below a fixed cut-off temperature, $T_{c}$. Specifically, at
temperatures below $T_{c}$, the parallel thermal conductivity is increased:

$\kappa^{*}_{\parallel}(T)=\begin{cases}\kappa_{0}T^{5/2},&\,\,\,T\geq
T_{c},\\\ \kappa_{0}T_{c}^{5/2},&\,\,\,T<T_{c},\\\ \end{cases}$ (7)

and the optically thin radiative losses are decreased:

$\Lambda^{*}(T)=\begin{cases}\Lambda(T),&\,\,\,T\geq T_{c},\\\
\Lambda(T)\left(\frac{T}{T_{c}}\right)^{5/2},&\,\,\,T<T_{c}.\\\ \end{cases}$
(8)

Here, the asterisks represent the modified expressions implemented for our
simulations. This will result in some broadening of the transition region, and
this effect is greater for higher cut-off temperatures.

Figure 1: Initial conditions for the 256x (dashed) and 4096x (solid)
resolution simulations with the SH (black), TRAC (blue), L09a (red) and L09b
(green) thermodynamic treatments. The upper panels show the temperature
profile over the whole loop (left) and zoomed in over the transition region
(right). The crosses in the upper right hand panel reveal the simulation
resolution (for the 256x case) by illustrating the location of the grid
points. The lower left panel shows the logarithm of the density close to the
lower foot point of the loop and the lower right panel shows the local Alfvén
speed profile over the entire loop, for each model.

The TRAC method, on the other hand, uses a variable cut-off temperature
$T_{c}$, which is continuously updated in response to simulation conditions.
In particular, it sets $T_{c}$ as low as possible such that sufficient
numerical resolution is maintained in the transition region. As such, the TRAC
method results in the minimal possible transition region broadening for a
given numerical resolution. Full details of the implementation are described
in Johnston & Bradshaw (2019) and Johnston et al. (2021). We note that, in
general, the TRAC approach is typically preferable as it eliminates
unnecessary transition region broadening and does not require a suitable cut-
off temperature to be identified a priori. In this article, we consider four
different treatments for the transition region; the unmodified Spitzer-Härm
conduction (no temperature cut-off), the TRAC treatment and two L09 cases with
fixed cut-off temperatures at $T_{c}=2.5\times 10^{5}$ K and $T_{c}=5\times
10^{5}$ K. These fixed temperature cut-offs are representative of values used
in the existing literature (e.g. Howson & De Moortel, 2022; Van Damme et al.,
2020). For brevity, we refer to these cases as SH, TRAC, L09a and L09b,
respectively.

In addition to these numerical treatments, thermal conduction can be limited
by the maximum conductive flux that the plasma is able to support (the free-
streaming limit). This occurs when the energy-transporting particles are all
travelling in the same direction at the electron thermal speed. In LareXd,
this saturated flux, $F_{s}$ is implemented as

$F_{s}=\alpha nk_{B}Tv_{\text{th, e}}$ (9)

where $\alpha$ is a user-imposed coefficient, $n$ is the number density,
$k_{B}$ is the Boltzmann constant, $T$ is the temperature and $v_{\text{th,
e}}=\sqrt{k_{B}T/m_{e}}$ is the electron thermal speed ($m_{e}$ is the
electron mass). For the classically derived expression, where the limit is
imposed by the thermal speed, a value of $\alpha=1.5$ is required (e.g.
Bradshaw & Cargill, 2006). However, it has been suggested that using
$v_{\text{th, e}}/6$ is more appropriate (e.g. $\alpha=1.5/6=0.25$). In either
case, for the relatively cool loops considered here, the maximum conductive
flux is well below the saturation threshold and thus the choice of $\alpha$
(between these two values) is moot. The choice would become significant,
however, for much hotter plasma (e.g. during solar flares). For the current
study, we selected the classical value, $\alpha=1.5$.

###  2.1 Initial conditions and numerical relaxation

We followed standard techniques for modelling closed, curved coronal loops as
straight magnetic field lines embedded in the solar photosphere at both
magnetic foot points. We consider a field line of length $l=120$ Mm (-60 Mm
$\leq y\leq y$ 60 Mm) which is assumed to be aligned with the central axis of
a coronal loop. The magnetic field is uniform, parallel with the $y$ axis and
has a field strength of 25 G. For simplicity, we neglect any expansion of the
loop cross-section with height. We have implemented a range of numerical
resolutions: $\\{256,512,1024,2048,4096\\}$ grid points along the field line.
We note that although we are only considering one-dimensional problems, for
consistency with existing multi-dimensional models of the Sun’s atmosphere, we
elected to use a two-dimensional code. However, in this paper we assume
perfect invariance across the magnetic field. The cross-field direction ($x$
axis) is simulated with a nominal number of grid points (4) and periodic
boundaries to maintain this invariance. The $z$ direction is perfectly
invariant but the magnetic and velocity fields are permitted to have non-zero
components in this direction.

In order to mimic the curved geometry of a coronal loop, we define gravity to
act in the $y$ direction with a sinusoidal profile:

$g_{\parallel}=g_{0}\sin\left(\frac{\pi y}{2y_{\text{max}}}\right),$ (10)

where $g_{0}\approx 274\text{ m s}^{-2}$ is the gravitational acceleration at
the solar surface ($y=\pm 60$ Mm). The acceleration points vertically
downwards for $y<0$ and upwards (towards y=$y_{\text{max}}$) for $y>0$ Mm.
This field-aligned component of the gravity assumes that the loop has a semi-
circular profile. We neglect any cross-field contribution due to gravity.

Coronal loops are embedded in the relatively cool and dense lower layers of
the atmosphere (photosphere and chromosphere) and extend into the hot and
tenuous corona. To include this in our initial conditions, we assume a
temperature profile, $T=T(y)$ defined by:

$T(y)=T_{\text{ch}}+\frac{T_{\text{co}}-T_{\text{ch}}}{2}\left\\{\tanh\left(\frac{y+a}{b}\right)-\tanh\left(\frac{y-a}{b}\right)\right\\},$
(11)

where $T_{\text{ch}}=2\times 10^{4}$ K is the chromospheric temperature,
$T_{\text{co}}=10^{6}$ K is the initial coronal temperature, $a=55$ Mm
controls the location of the transition region and $b=0.5$ Mm controls the
initial width of the temperature transition. In Bradshaw & Cargill (2013), the
authors detail the numerical resolution requirements for accurately
reproducing coronal evolution in gravitationally stratified one-dimensional
solar atmospheric loops. As discussed in Sect. 1, spatially under-resolving
the transition region will lead to lower evaporative upflows of dense plasma
in response to heating, and ultimately, lower coronal densities. According to
Table 1 in Bradshaw & Cargill (2013), for coronal loops with relatively cool
apex temperatures (3 MK), a grid resolution of approximately 25 km in the
transition region is sufficient. Given that the loops in our current study are
cooler (1 MK), the resolution requirements are less stringent. As our highest
resolution case has a grid size of $120\text{ Mm}/4096\approx 29$ km, it is
likely able to reproduce coronal behaviour accurately. For the lower
resolutions cases, however, accurate coronal apex densities and temperatures
can only be reproduced with the modifications in thermodynamics discussed
above.

To find the loop’s initial density profile, we assumed hydrostatic force
balance and solved the resulting differential equation using a Runga-Kutte
scheme. Whilst this ensures the system is in pressure balance, the conduction,
radiation and heating terms in the energy equation 3 mean it is not in
thermodynamic balance. In order to achieve this for our initial equilibrium we
then perform a numerical relaxation using the Lare2d code. In this article we
set the background heating to be $Q_{\text{bg.}}=5\times 10^{-6}\text{ J
m}^{-3}\text{ s}^{-1}$ in all cases, which generates apex temperatures of
approximately 1.3 MK (see upper left panel of Fig. 1). We impose a large
viscosity term to damp field-aligned flows until the system settles into a
state of numerical equilibrium and we take this as the initial conditions for
our simulations. This relaxation viscosity is then set to zero for the
subsequent wave driving simulations (Sects. 3.3 & 3.4).

The resulting density and temperature profiles are shown in Fig. 1 for the
256x (dashed lines) and 4096x (solid lines) resolution cases with each of the
thermodynamic treatments. The upper left panel shows the temperature profiles
along the full loop length. We note that the majority of the curves are very
similar, closely approximating the profile of the high resolution SH model
(solid black line). The exceptions are the 256x SH case (dashed black line)
and the two L09 simulations (green lines). The 256x SH case significantly
underestimates the coronal density (see lower left panel and Bradshaw &
Cargill, 2013) and, as a result, coronal radiative losses will be much lower
in this case (the losses scale with the square of the density, see equation
3). Even though energy losses in the corona are dominated by conduction, the
reduction in the radiative cooling rate is sufficient to allow noticeably
higher temperatures in this setup. Conversely, the thermodynamic treatment in
the L09b simulations leads to significant over-broadening of the transition
region (see green curves in right hand panel). This ultimately leads to the
lower coronal temperatures, demonstrating that cut-off temperatures which are
too high can be unhelpful for coronal loop models. However, such high cut-offs
may remain necessary for very high temperature loops where the resolution
requirements are even more demanding.

The effects of each treatment on the transition region are clear in the
zoomed-in temperature profile displayed in the upper right panel of Fig. 1. At
both resolutions (256x and 4096x) the SH treatments produce the steepest
transition regions with broadening evident for all other cases. The broadening
is almost independent of resolution for the fixed temperature cut-off cases
(L09a in red and L09b in green), however, we see much more broadening in the
low resolution TRAC case (dashed blue) than in the high resolution equivalent
(solid blue). This is because the temperature cut-off adjusts in the TRAC
treatment according to the numerical resolution and the properties of the loop
(in order to ensure coronal evolution is modelled accurately). As the SH and
TRAC transition regions are very similar for the 4096x case (compare solid
blue and black curves), we will assume that this represents sufficient
resolution for this loop and thus use the 4096x SH case (solid black) as a
benchmark simulation representing the correct evolution.

The lower left panel of Fig. 1 shows the logarithm of the plasma density in
the bottom 20 Mm of the simulation domain. Again, this shows the effects of
the thermodynamic treatment on the transition region broadening. Both SH cases
(solid and dashed black lines) and the high resolution TRAC case (solid blue
line) produce a steep transition region, where the density decreases by
approximately two orders of magnitude over a very narrow layer. In all other
simulations, this decrease occurs over a wider region, allowing for the
gradients to be better resolved in the low resolution cases. We note that
despite this broadening, all simulations (with the exception of the 256x SH
case) exhibit similar density profiles in the corona (e.g. -40 Mm $\leq y\leq$
40 Mm). These density profiles are associated with variations in the local
Alfvén speed, $v_{A}=B/\sqrt{\mu_{0}\rho}$, displayed in the lower right hand
panel of Fig. 1. We note that the significantly lower coronal density attained
in the 256x SH setup, produces much higher coronal Alfvén speeds than in all
the other cases. In comparison, all other models generate similar Alfvén speed
profiles, especially in the upper atmosphere (e.g. -40 Mm $\leq y\leq$ 40 Mm).
However, as we discuss in detail below, even the relatively small differences
that persist are sufficient to cause significantly different wave dynamics
within each model.

###  2.2 Boundary conditions

We generated Alfvén waves within our model, by transversely perturbing the
magnetic field at the lower $y$ boundary in the invariant $z$ direction. We
consider wave drivers of the form

$v_{z}(t)=v_{0}\sin\omega t,$ (12)

where $v_{0}$ is a small amplitude to generate linear waves and $\omega$ is
the wave frequency. For the simulations described below, we use
$v_{0}=100\text{ m s}^{-1}$ which is less than 1% of the minimum Alfvén speed
in each simulation and much smaller than the coronal Alfvén speeds (see lower
right panel of Fig. 1). For the frequency, $\omega$, we consider the natural
frequencies of the system and this is discussed in more detail below. At the
other foot point of the loop, we impose ${\bf v}={\bf 0}$ for the velocity and
a zero-gradient condition for all other quantities. This fixes a wave node at
this upper boundary. As we are using a two dimensional MHD code, we also
define the $x$ boundaries to be periodic. This maintains the invariance in
this direction.

Figure 2: Alfvén travel time, $\tau(y)$ (see equation 13, for the 256x (dashed
lines) and 4096x (solid lines) resolutions with the SH (black), TRAC (blue),
L09a (red) and L09b (green) thermodynamic treatments. The left-hand panel
shows the travel time over the whole loop and, for clarity, the right hand
panel restricts the domain to the region close to the upper $y$ boundary. We
note that the 4096x SH curves (solid black) is not visible as they almost
exactly coincide with the 4096x TRAC curves (solid blue).

##  3 Results

In order to illustrate why the different transition region treatments can have
a large influence on Alfvén wave dynamics (and hence, energy flux) in coronal
loops, we first begin by considering the Alfvén travel time for each of our
initial conditions. This represents how long it takes an Alfvén wave front to
propagate from the driven foot point to a particular height. In other words,
we calculate the travel time $\tau(y)$ as

$\tau(y)=\int_{0}^{y}\frac{d\mathrm{s}}{v_{A}},$ (13)

, where $v_{A}$ is the local Alfvén speed and $s$ is an infinitesimal field
line element. We show the function, $\tau(y)$, for a variety of our initial
conditions in Fig. 2. The left hand panel shows the profile over the entire
length of the loop and, for clarity, the right hand panel shows $\tau(y)$
close to the non-driven foot point. We see that the total travel time varies
across the simulations from approximately 88s (for the low resolution SH case)
to around 116 s (for the high resolution L09b case). We note that here we have
only displayed the lowest (256x) and highest (4096x) numerical resolutions and
all other cases lie within the two extremes for any given thermodynamic
treatment. We see that the right hand panel clearly demonstrates the
significant differences between the Alfvén travel times for the different
resolutions and transition region treatments. There are two pertinent points
which drive these differences. Firstly, the different Alfvén speed profiles in
each case (as discussed for the lower right panel of Fig. 1), naturally lead
to different Alfvén travel times. For example, the low coronal densities in
the 256x SH simulation (dashed black line) result in higher Alfvén speeds and
shorter travel times. We also note that even in cases where the coronal Alfvén
speeds are very similar, differences in the wave speed in the transition
region can result in discrepancies in the total travel time. For example, the
SH and L09a treatments exhibit around a 10% difference in the travel time
despite the relatively short length of the transition region (it is much less
than 10% of the volume even in the broadened cases). Differences in Alfvén
speeds in the lower atmosphere can have a significant impact due to the low
speeds and, hence relatively long time that the wave takes to propagate
through this region.

Secondly, the numerical resolution can modify the travel time simply because
the local Alfvén speed is only known at a relatively low number of grid
points. In practice, the integral in equation 13 is calculated as a discrete
summation over the simulation grid points. Therefore, in the transition
region, where the wave speed changes rapidly, the low resolution simulation
does not track the travel time well. As a result, the low resolution curves
(dashed lines) show different travel times to their higher resolution
counterparts. This effect is clearest for the fixed temperature cut-off models
(as the cut-off does not change as a function of resolution), L09a (red) and
L09b (green curves) which show very similar density and Alfvén wave profiles
(see lower row of Fig. 1) but still produce different travel times. Using sub-
grid interpolation methods for equation 2 would reduce this difference in the
calculated Alfvén travel times but this would not reflect the wave propagation
as simulated by the Lare2d code.

###  3.1 Eigenfrequencies and eigenmodes

We can use the Alfvén travel times calculated above, together with the WKB
approximation, to provide estimates for the natural frequencies of the coronal
loops. However, these will likely be inaccurate (especially for the
fundamental mode and low number harmonics (e.g. Wright & Smith, 1990) due to
the significant non-uniformity of the wave speed along the modelled loops.
Instead, we calculate the eigenfrequencies and corresponding eigenmodes of the
field lines by considering the wave equation for non-constant propagation
speed:

$\frac{\partial^{2}v_{z}}{\partial
t^{2}}=v_{A}^{2}(y)\frac{\partial^{2}v_{z}}{\partial y^{2}},$ (14)

where $v_{A}$ is the local Alfvén speed displayed in the lower right panel of
Fig. 1. Then, by considering non-trivial, oscillatory, separable solutions,
$v_{z}=Y(y)T(t)$, of this partial differential equation, we can express the
spatial variation, $Y(y)$ as

$\frac{\mathrm{d}^{2}Y}{\mathrm{d}y^{2}}+\frac{\omega^{2}}{v^{2}_{A}(y)}Y=0,$
(15)

where $\omega$ is a real constant. The Alfvén eigenmodes and corresponding
eigenfrequencies are given by functions, $v_{z}(y)$, and constants, $\omega$,
such that there are wave nodes ($v_{z}=0$) at the two field line foot points.
In order to find these, we implement a shooting method to find numerical
solutions of equation 15.

In Fig. 3, we display the eigenmodes of the first seven harmonics for the
4096x SH (benchmark) simulation. For clarity, we have normalised the maximum
of each curve to $1-n/10$, where $n$ is the harmonic number beginning with
$n=0$ for the fundamental mode, $n=1$ for the first overtone, and so on. As
such, the amplitude of each eigenmode is arbitrary. We note that due to the
relatively low Alfvén speeds in the chromosphere and lower transition region,
for higher overtones, the majority of the wave nodes are located close to the
two loop foot points. As such, low resolution simulations (e.g. 256x) will not
be able to resolve the short wavelengths in the chromosphere for higher
frequency modes. This is discussed in more detail in Sect. 3.4. We also note
that, due to the high density at these low altitudes, the magnitude of the
eigenmodes here are much smaller than in the coronal volume.

Figure 3: Eigenmodes of the fundamental (black) mode and first (purple),
second (dark blue), third (light blue), fourth (green), fifth (orange) and
sixth (red) overtones.These are eigenmodes of equation 15 calculated for the
4096x SH model. Figure 4: _Left column:_ Eigenfrequencies of the fundamental
(zeroth harmonic) and higher overtones for loops modelled with different
numerical resolutions and thermodynamic treatments. In the upper three panels,
we show the effects of resolution on loops modelled with the TRAC, L09b and SH
treatments, respectively. For these panels, the dashed black line shows the
4096x SH frequencies which we use as a benchmark. The colour used for each
resolution is consistent across these panels. The bottom left panel shows the
effects of the thermodynamic treatment on the eigenfrequencies in the 256x
(dashed lines) and 4096x (solid line) models. _Right column:_ The relative
error between the curves in the adjacent panels in the left-hand column and
the benchmark case (4096x SH).

We note that the exact forms of these eigenmodes will be sensitive to a wide
variety of factors, including but not limited to the relative sizes of the
chromosphere and corona, loop curvature, loop asymmetry and loop expansion. As
such, they may not be representative of eigenmodes in real atmospheric loops.
However, here we simply wish to compare the effects of the numerical
resolution and thermodynamic treatment on the natural frequencies of the
modelled loops. To this end, in the left hand column of Fig. 4, we display the
eigenfrequencies for the first 7 harmonics (including the fundamental mode as
the zeroth harmonic) for different numerical resolutions and thermodynamic
treatments. In the upper three panels, we show how these vary for different
numerical resolutions in the TRAC, L09b and SH cases, respectively. In each
panel, we also include the eigenfrequencies for the 4096x SH simulation as a
benchmark (dashed line). For clarity, the right hand column of Fig. 4 shows
the relative error between each calculated eigenfrequency and the benchmark
solution. It is clear that as the resolution is increased, the TRAC (first
row) and SH (third row) treatments converge to the benchmark eigenfrequencies.
However, for the resolutions typically attainable in large scale, 3D MHD codes
(e.g. corresponding to 256 or 512 grid points in the loop-aligned direction),
the error in the frequency calculation can be approximately 40%. At low
resolutions, all thermodynamic treatments produce poor estimates of the
eigenfrequency calculated for the fully resolved loop (particularly for higher
harmonics). The relative error for each of these cases will depend on how
important the transition region is for determining the eigenfrequency of a
given loop. In particular, the relative accuracy may depend on loop lengths
and temperatures and the depth of the chromosphere. We also note that for the
L09b treatment (second row), whilst the eigenfrequencies do converge at high
resolution, they do not converge to the benchmark case. This is easy to see at
higher harmonics. This behaviour arises because the fixed temperature cut-off
consistently over-broadens the transition region even at high resolutions.
This phenomenon is significantly reduced (albeit still present, particularly
for higher harmonics) in L09a case due to the reduced broadening of the
transition region.

In the bottom row of Fig. 4, we display the effects of the thermodynamic
treatment on the calculated eigenfrequencies (left panel) for the two extreme
resolutions (256x, dashed lines; 4096x, solid lines) and the relative error
(compared to the benchmark solution). At low resolution, we see that the
choice of thermodynamic treatment has significant implications for the natural
frequencies of the system (particularly at higher harmonics) and, as discussed
before, that no choice accurately reproduces the benchmark solution. At high
resolution, however, all methods (except L09b) reproduce the eigenfrequencies
reasonably accurately (with small errors at high harmonics for the L09a case,
due to the fixed transition region broadening). This is not a particularly
surprising result; modifying the Alfvén speed profile in the lower atmosphere
will certainly impact on the natural frequencies of the modelled loop.
However, this can have important consequences for energy flux in solar
atmospheric models as different resonances can be excited in each of these
setups. In Sect. 3.4, we will consider whether this can lead to systematic
errors for energy injection rates in general coronal heating models. However,
as a brief aside, we will first consider implications of these results for
seismological inversions.

Figure 5: Transverse velocity, $v_{z}(y)$, excited by a continuous, high
frequency, sinusoidal wave driver imposed at the lower $y$ bounndary. We show
results from the 4096x resolution simulations for the SH (black), TRAC (blue),
L09a (red) and L09b (green) cases. The left panel shows an early time,
$t\approx 60$ s, and the right panel shows a later time, $t\approx 550$ s,
after several wave reflections have occurred. The solid black line (SH) is the
benchmark result.

###  3.2 Seismological implications

At the broadest level, coronal seismology uses the properties of observed
waves to deduce properties of the solar corona that cannot be measured
directly (e.g. see reviews by Nakariakov & Verwichte, 2005; Andries et al.,
2009; Ruderman & Erdélyi, 2009; Arregui et al., 2012; De Moortel & Nakariakov,
2012; Nakariakov & Kolotkov, 2020). It is a well-developed field and has been
used to provide estimates of coronal magnetic field strengths (e.g. Nakariakov
& Ofman, 2001; Soler et al., 2010), the transverse density structuring of
loops (e.g. Goddard et al., 2017; Pascoe et al., 2017; Goddard et al., 2018;
Pascoe et al., 2018) and, through the use of Bayesian inference, to provide
evidence for and against different solar atmospheric models (e.g. Arregui &
Asensio Ramos, 2011; Montes-Solís & Arregui, 2019; Arregui, 2022). These
studies demonstrate how expected wave behaviour (e.g. propagation speed,
damping rates) derived from mathematical models can be used to deduce unknown
parameters (e.g. field strength) from solar observations. However, these
methods can have large uncertainties associated with them, not least because
it is often difficult to definitively identify what wave mode is being
observed.

As discussed in the previous section, the lower atmosphere can play a very
important role in establishing the natural frequency of these magnetic
structures. They are not simply coronal loops, but are embedded in the
transition region and chromosphere too. If wave nodes are established at the
upper transition region, then assuming standing oscillations are purely
coronal is likely a good approximation. This may be a reasonable view of
oscillations excited impulsively in the corona (e.g. from a nearby solar flare
Nakariakov et al., 1999; Li et al., 2023a). In this paradigm, the large
density gradients in the transition region may act as effective reflectors of
wave energy, effectively forming wave nodes in the upper transition region.
However, if these standing waves are driven by oscillatory flows in the
chromosphere/photosphere (as we are assuming in this article and are often
invoked for driving decayless oscillations, e.g. Nisticò et al., 2013;
Anfinogentov et al., 2015; Zhong et al., 2022; Petrova et al., 2023), then the
contribution of the lower atmosphere to the natural frequencies must be
accounted for.

As a very simple example of this idea in practice, let’s suppose we observe a
fundamental standing Alfvén wave (in reality we may be more likely to observe
a kink mode but the same argument applies). For our setup, we may expect to
observe a frequency of approximately $0.0563\text{ s}^{-1}$ (the fundamental
frequency for the 4096x SH model). However, these waves will typically be
interpreted as a coronal-only oscillations. Thus, if we instead consider the
eigenfrequencies of a coronal-only system (using equation 15), then we find a
fundamental frequency of $0.0947\text{ s}^{-1}$. Indeed this is a similar
value to the frequency of the second overtone of the full system which appears
to have wave nodes at the base of the corona (see dark blue line in Fig. 3).
As the observed frequency can be used to estimate the magnetic field strength,
we will obtain an estimate that is 0.0947/0.0563 $\approx$ 1.7 times too big.
This simple calculation does not consider inaccuracies in numerical wave
modelling due to the thermodynamic treatments discussed throughout this
article. However, despite this, the simple argument highlights an important
consideration for seismological inversions; the location of the nodes for
observed standing modes must be clearly identified in order to understand the
true natural frequencies of the oscillating structure.

###  3.3 Modelling propagating waves

Returning to the effects of the thermodynamic model on simulated wave
dynamics, we now consider the case of a wave that is continuously driven from
the lower $y$ boundary. Whereas in Sect. 3.1 we considered a purely analytical
description, here we now model the propagation of Alfvén waves using the
Lare2d code. We impose a driver as described in equation 12 with frequency,
$\omega=0.069\text{ s}^{-1}$. We note that this is a comparable to the
magnitude of the fundamental frequency in each case but is non-resonant.

In Fig. 5, we display the wave velocity ($v_{z}$) as a function of position
along the loop for the four high resolution (4096x) simulations. The left hand
panel shows the wave at an early time ($t\approx 60$ s), before the
propagating wave fronts have reached the opposite foot point. The right hand
panel, on the other hand, is much later, ($t\approx 550$ s), after several
reflections have occurred. In both panels, we have normalised all curves by
the maximum velocity obtained in the L09b simulation and the solid black line
corresponds to the benchmark (4097x SH) solution. We see that at early times
(left hand panel), the TRAC treatment (blue) provides very good agreement with
the benchmark case. This is unsurprising given that TRAC produces minimal
transition region broadening at this high resolution (compare solid blue and
black lines in the upper right panel of Fig, 1). The L09a (red) and L09b
(green) on the other hand both show lower mean propagation speeds, with the
effect greatest for the more significant broadening in the L09b case. This is
in agreement with the longer travel times for these setups shown in Fig. 2. We
also note that they both exhibit larger amplitudes, suggesting that more
energy is transmitted into the corona with these treatments. We will assess
this point in more detail in Sect. 3.4. The small scale oscillations that
follow the leading wave front in all cases are associated with wave
reflections excited as the front passes through inhomogeneities in the local
Alfvén speed (e.g. Asgari-Targhi et al., 2021; Pascoe et al., 2022).

In the right hand panel, we see that at later times that there is little
agreement between with L09a, L09b and benchmark solutions. After a few
reflections, the relatively small differences visible in the left panel
compound to produce very different waves. Once again, the L09b case has the
largest amplitude, suggesting that energy is injected into the corona more
efficiently in this case. The TRAC treatment, on the other hand, still
reproduces the benchmark solution with reasonable accuracy. However, there are
two important caveats to note here. Firstly, as time progresses and more wave
reflections take place, the differences between the TRAC and benchmark
solutions will become increasingly pronounced. Furthermore, this favourable
result for the TRAC method is less applicable at lower resolutions, where the
transition region broadening will be significant in the TRAC case too. In
general, the TRAC treatment is beneficial because it reduces broadening when
possible. However, if the resolution or loop parameters are such that
significant broadening is required (e.g. for high temperature loops), then the
TRAC case will perform as poorly as the fixed temperature cut-offs.

Figure 6: Pearson correlation between the wave velocity in the benchmark
solution (4096x SH) and in the models with other thermodynamic treatments. The
panels show how the correlation changes as a function of time at high
resolution (4096x; upper panel) and at low resolution (512x; lower panel). The
black curve (SH case) is not included in the upper panel because this would be
measuring the correlation between one result and itself.

In Fig, 6, we provide a measure of the accuracy of the models produced by each
thermodynamic treatment in the 4096x (upper) and 512x (lower) resolution
cases. In particular, at every time output from the simulation, we measure the
correlation (with the Pearson correlation coefficient) between the wave
velocity ($v_{z}$) in each model and the benchmark solution (4096x SH). A
correlation of 1 indicates a perfect match between the solutions, a
correlation of 0 indicates there is no relationship between solutions and a
score of -1 indicates the solutions are perfectly out-of-phase. In practice, a
correlation of close to 1 indicates a good match. For the 4096x case,
unsuprisingly we see that the TRAC solution with minimal transition region
broadening produces a good match with the benchmark. As explained above this
match steadily worsens as time progresses. In general we see that reducing the
modification of the thermodynamics is beneficial at high resolution and thus
the L09a case performs better than the L09b simulation (with both worse than
the TRAC case). Meanwhile, in the lower panel of Fig. 6, we see that for the
lower resolution cases (although this still represents a high resolution in
terms of large scale 3D simulations), all thermodynamic treatments produce a
poor match with the benchmark solution. The reduced broadening cases (TRAC,
blue and L09a, red) produce marginally better results than the SH (black) and
L09b (green) cases but in general the benchmark solution is not reproduced.

These results are concerning given that this resolution is representative of
that used in state-of-the-art models of the solar atmosphere. That said, we
know that many of our models do not accurately reproduce dynamics in the lower
atmosphere anyway (e.g. due to neglecting partial ionisation, radiative
transfer etc.). For some studies, coronal modellers may take the view that we
know there is wave energy in the corona, so, as long as we model waves with
the correct amplitudes/wavelengths etc., understanding precisely how waves are
injected into the corona is unimportant in comparison to the dynamics (e.g.
phase mixing, resonant absorption, Kelvin-Helmholtz instability) that occur as
they propagate within the corona. However, increasingly coronal models are
being driven by photospheric motions (e.g. through self-consistent convection
or through drivers imposed below the chromosphere). Additionally, precisely
tracking the energy flux through the lower atmosphere is important for
understanding the energetics of the atmosphere as a whole. Indeed if the
thermodynamic treatments permit an artificially high (or low) energy flux into
the corona, then we will obtain incorrect heating rates, for example. There is
also a more subtle issue. If the different transition region treatments are
associated with errors of different magnitudes for different driver types,
then comparing the effects of different drivers becomes problematic. In
particular, if a low frequency driver is relatively unaffected by these errors
in comparison to a high frequency driver, then making a fair comparison is
challenging. In Howson et al. (2020); Howson & De Moortel (2022), the authors
presented comparisons of heating associated with long (DC) and short (AC) time
scale driving and found DC driving to be more effective for heating the
corona. The latter study used the L09a treatment and thus it is important to
establish whether the energy injection rates are artificially enhanced for the
different thermodynamic treatments and whether the DC or AC driving is
affected more.

###  3.4 Broadband driving

In order to explore this point in more detail, we will now turn our attention
to a broadband wave driver to establish whether there are systematic errors in
the energy injection rates. In the previous sections, we have discussed how
the numerical treatment of the TR modifies the system in response to
continuous transverse driving at a single frequency. However, it is unlikely
that solar photospheric motions oscillate with a fixed frequency for long
periods of time and thus it is important to quantify the impact of the TR
treatment on energy flux for waves driven with a broadband driver. To this
end, we consider drivers, $v_{z}(t)$ defined by

$v_{z}(t)=\sum_{i=1}^{N}u_{0}\sin\left(\omega_{i}t+\psi_{i}\right).$ (16)

These broadband wave drivers are defined as the sum of $N$ sinusoidal
components with different frequencies, $\omega_{i}$ and phase shifts
$\psi_{i}$. These phase shifts are randomly selected from a uniform
distribution on $[0,2\pi]$. For this article we take $N=50$ and we restrict
our consideration of the frequency space to a range between two cut-off
frequencies, $\omega_{\text{min}}$ and
$\omega_{\text{max}}=3\omega_{\text{min}}$. In particular, the $i^{\text{th}}$
frequency is defined as

$\omega_{i}=\frac{i}{N}\left(\omega_{\text{max}}-\omega_{\text{min}}\right)+\omega_{\text{min}}.$
(17)

We consider two different frequency ranges, to represent low and high
frequency broadband drivers defined by
$\omega_{\text{min}}=0.8\omega_{f}\approx 0.041\text{ s}^{-1}$ and
$\omega_{\text{min}}=3.2\omega_{f}\approx 0.164\text{ s}^{-1}$. Here,
$\omega_{f}\approx 0.0513\text{ s}^{-1}$ is the fundamental frequency in the
4096x SH loop (see Fig. 4). As the amplitude of each component, $u_{0}$, is a
constant, the power in the broadband driver is independent of frequency over
the range $[\omega_{\text{min}},\omega_{\text{max}}]$. As such, within our
frequency range, the wave driver has a white noise profile. The temporal
profile of the driver is shown in Fig. 7.

Figure 7: Temporal variation of the wave velocity imposed at the lower $y$
footpoint by the low frequency broadband driver. Figure 8: Energy content in
the coronal portion of the domain as a function of time for low (left panels)
and high (right panels) frequency drivers. We show results from the SH
(black), TRAC (blue), L09a (red) and L09b (green) simulations at high (4096x,
upper panels) and low (512x, lower panels) resolution. We have normalised all
curves to the maximum energy content in the high frequency, high resolution
L09b simulation.

From the discussion outlined above on Figs. 5 & 6, we know that the simulated
wave dynamics will be different for each thermodynamic treatment and numerical
resolution. However, the important question here is whether more or less
energy is injected into the corona for some treatments/resolutions in
comparison to the benchmark case. In Fig. 8, we show the change in the energy
density, $E(t)$, within the range -40 Mm $y$ 40 Mm as a function of time for
high (4096x, upper panels) and low (512x, lower panels) resolution
simulations. For this analysis we have calculated,

$E(t)=\int_{y=-40\text{ Mm}}^{y=40\text{
Mm}}\left(\frac{B^{2}}{2\mu_{0}}+\frac{\rho
v^{2}}{2}+\frac{P}{\gamma-1}+\rho\Phi\right)\,\,\mathrm{d}y,$ (18)

where the terms in the integrand are the magnetic, kinetic, internal and
gravitational potential energies, respectively. For each case, we subtracted
the initial value to show the change in energy density. The $y$ range for the
integral ensures that we restrict our attention to the coronal subvolume of
the domain. The left and right columns correspond to the low and high
frequency broadband drivers, respectively. We show the results from SH
(black), TRAC (blue), L09a (red) and L09b (green) simulations and all curves
normalised by the maximum energy content attained in the high frequency, high
resolution L09b case.

In the high resolution, low frequency case (upper left panel), we see that the
TRAC simulation (blue) reproduces the coronal energy content of the benchmark
case (black curve) well. The L09a and L09b provide successively worse
estimates, with the coronal energy content being overestimated by 5x at some
points during the L09b case. This shows that artificially broadening the
transition region permits a greater average flux of wave energy into the
corona. This is because the reduced Alfvén speed gradients in the broadened
cases reduce the efficiency of wave reflections back into the lower
atmosphere. For the high frequency, high resolution cases (upper right panel),
we see that the TRAC method still reproduces the benchmark solution well,
however, the L09a cases now provides an equally poor estimate as the L09b
simulation. The higher frequency waves have shorter wavelengths which are more
liable to reflect in the transition region and hence even the weak broadening
produces significant over estimations in the energy flux. This shows that the
overestimation of energy injection rates can also be frequency dependent.

In the lower resolution simulations (lower panels), we see that the mean
energy injection rate is lower than in their high resolution counterparts.
This effect is particularly profound for the high frequency cases (lower right
panel). This is largely a consequence of relatively high numerical dissipation
rates in the lower atmosphere due to the relatively low resolution and short
wave lengths (particularly for the high frequency waves). This point
notwithstanding, we now see different energy injection rates for all
simulations and the TRAC cases no longer coincide with the SH results (or
indeed with the high resolution benchmark). For both frequency drivers, we see
that the wave energy in the coronal volume at any given point is very
sensitive to the particular thermodynamic treatment. As with the high
resolution cases, there is a tendency for broadened transition regions to
permit enhanced energy flux into the corona and again these results are
frequency dependent. These results show that accurately modelling MHD waves as
they propagate through the lower atmosphere and into the corona is extremely
challenging, particularly for high frequency modes.

##  4 Discussion and conclusions

In this article, we have considered how numerical resolution and a variety of
thermodynamic treatments (Lionello et al., 2009; Mikić et al., 2013; Johnston
& Bradshaw, 2019; Johnston et al., 2021) can modify the flux of energy from
the lower solar atmosphere into the corona. As they are well understood
mathematically, we have used Alfvén waves with a range of frequencies as a
proxy for energy injection mechanisms. We have shown that the Alfvén travel
times, eigenmodes and eigenfrequencies, and energy injection rates are all
highly sensitive to the resolution and thermal conduction models used within
simulations. Additionally, we have highlighted the importance of the lower
atmosphere on seismological inversions when using wave modes excited beneath
the transition region.

Increasingly, numerical models of the solar atmosphere are including each of
the distinct layers (photosphere, chromosphere, transition region, corona)
within the simulation volume. These contemporary studies treat the physics of
the lower atmosphere with varying degrees of completenness. For example,
simply treating the chromosphere as a mass reservoir (e.g. Van Damme et al.,
2020; Reid et al., 2021; Cozzo et al., 2023; Skirvin et al., 2023), or
including more accurate chromospheric physics such as radiative transfer (e.g.
Battaglia et al., 2021; Nóbrega-Siverio & Moreno-Insertis, 2022; Hansteen et
al., 2023; Martínez-Sykora et al., 2023). Inevitably, these different
approaches will lead to different energy densities at the top of the
chromosphere. However, the nature of the transition region can then have
further consequences for the amount of energy reaching the corona. We have
shown that the choice of numerical treatment for thermodynamics in the
transition region will modify the mechanical energy at higher altitudes.
Furthermore, and perhaps more alarmingly, we have also shown that these
treatments are associated with energy injection errors which depend on the
frequency of the driver. In light of this, significant care is required when
undertaking direct comparisons of coronal heating driven by different
photospheric convection profiles (e.g. long vs short time scale driving in
Howson et al., 2020; Howson & De Moortel, 2022).

In the present study, we have considered loops that are relatively easy to
model numerically. In particular, they are not especially hot (e.g. Wang et
al., 2003) and they are not dynamic. As such, our negative findings may be
even worse in many other situations. For example, hot loops with apex
temperatures several times larger than in our models, are commonplace in
active regions, and are likely important for identifying heating events (e.g.
Bradshaw & Klimchuk, 2011). These loops require higher numerical resolution or
enhanced transition region broadening (either automatically with TRAC or with
higher fixed temperature cut-offs for the L09 approach), and as such will
likely exacerbate our negative results. It is also important to reiterate that
short wave lengths in the chromosphere (due to the relatively low Alfvén
speeds) can lead to significant wave damping and reduced energy transmission,
particularly for high frequencies. This may be unimportant for studies
interested in wave propagation in the corona, however, as discussed above,
this does have implications for comparing the effects of different
photospheric drivers on coronal heating.

It may seem that waves reflected at the transition region may have further
opportunities to be transmitted into the corona following a second (and
subsequent) reflections at the driven boundary. However, whilst this is an
important and subtle point, this effect does not necessarily permit enhanced
coronal energy injection over long time periods. In our simulations, upon
returning to the driven boundary, reflected waves interact with the imposed
velocity field. This modifies the Poynting flux injected into the domain and
can have important consequences for the energetics of the system. If the
reflected waves are in-phase with the wave driver, then resonances will be
excited and the Poynting flux will increase. However, for non-broadband
drivers this is unlikely, and the waves will typically be non-resonant. If
non-resonant reflected waves have the same amplitude as the driver, then, on
average, the Poynting flux will be equally positive and negative over the
course of a wave period. In such a case, there will be no net energy injection
after the first reflected waves reach the boundary. However, typically the
reflection coefficient at the transition region will be less than unity and
the reflected waves will have lower amplitudes. In this case, the driver will
still inject energy into the system but at a reduced rate. As such, it may be
misleading to think that the reflected waves have multiple attempts to
propagate into the corona as they can also reduce the rate at which energy is
injected into the system.

On the basis of our results, we recommend the use of the TRAC method as the
most suitable treatment for resolutions currently attainable in large scale,
multi-dimensional coronal heating models. As the TRAC approach is associated
with the minimum possible broadening (for any given loop and numerical
resolution), it will generally produce the smallest errors in the simulated
wave dynamics. That said, at attainable numerical resolutions, it will still
provide a poor match to a fully resolved benchmark case, particularly for high
frequency waves. At this stage, we offer no concrete solution to these issues
and believe an appropriate first step will be accurately quantifying the
effects of the thermodynamic treatments for a variety of loop lengths and
geometries, magnetic field strengths, wave modes (e.g. slow waves, kink waves)
and also for long timescale driving (e.g. DC heating mechanisms). Only then
will we be able to determine how significant these issues are for contemporary
three-dimensional solar atmospheric modelling.

## Acknowledgements

The research leading to these results has received funding from the UK Science
and Technology Facilities Council (consolidated grant ST/S000402/1). The
authors would like to thank Dr J Reid for his help and comments during the
preparation of this manuscript. Finally, the authors would also like to thank
the anonymous referee for considering our work and providing helpful
suggestions to improve our article.

## Data Availability

The data from the numerical simulations and analysis presented in this paper
are available from the corresponding author upon reasonable request.

## References

  * Andries et al. (2009) Andries J., van Doorsselaere T., Roberts B., Verth G., Verwichte E., Erdélyi R., 2009, Space Science Reviews, 149, 3
  * Anfinogentov et al. (2015) Anfinogentov S. A., Nakariakov V. M., Nisticò G., 2015, Astronomy and Astrophysics, 583, A136
  * Antolin et al. (2016) Antolin P., De Moortel I., Van Doorsselaere T., Yokoyama T., 2016, The Astrophysical Journal, 830, L22
  * Arber (2018) Arber T. D., 2018, LareXd User Guide. https://github.com/Warwick-Plasma/Lare2d/tree/main/Manual, version 3.3 edn
  * Arber et al. (2001) Arber T. D., Longbottom A. W., Gerrard C. L., Milne A. M., 2001, Journal of Computational Physics, 171, 151
  * Arregui (2022) Arregui I., 2022, Frontiers in Astronomy and Space Sciences, 9, 826947
  * Arregui & Asensio Ramos (2011) Arregui I., Asensio Ramos A., 2011, The Astrophysical Journal, 740, 44
  * Arregui et al. (2012) Arregui I., Oliver R., Ballester J. L., 2012, Living Reviews in Solar Physics, 9, 2
  * Asgari-Targhi et al. (2021) Asgari-Targhi M., Asgari-Targhi A., Hahn M., Savin D. W., 2021, The Astrophysical Journal, 911, 63
  * Battaglia et al. (2021) Battaglia A. F., Canivete Cuissa J. R., Calvo F., Bossart A. A., Steiner O., 2021, Astronomy and Astrophysics, 649, A121
  * Bingert & Peter (2011) Bingert S., Peter H., 2011, Astronomy and Astrophysics, 530, A112
  * Bradshaw & Cargill (2006) Bradshaw S. J., Cargill P. J., 2006, Astronomy and Astrophysics, 458, 987
  * Bradshaw & Cargill (2013) Bradshaw S. J., Cargill P. J., 2013, The Astrophysical Journal, 770, 12
  * Bradshaw & Klimchuk (2011) Bradshaw S. J., Klimchuk J. A., 2011, The Astrophysical Journal Supplement, 194, 26
  * Braginskii (1965) Braginskii S. I., 1965, Reviews of Plasma Physics, 1, 205
  * Breu et al. (2022) Breu C., Peter H., Cameron R., Solanki S. K., Przybylski D., Rempel M., Chitta L. P., 2022, Astronomy and Astrophysics, 658, A45
  * Chen et al. (2023) Chen F., Cheung M. C. M., Rempel M., Chintzoglou G., 2023, The Astrophysical Journal, 949, 118
  * Cheung et al. (2019) Cheung M. C. M., et al., 2019, Nature Astronomy, 3, 160
  * Cozzo et al. (2023) Cozzo G., Reid J., Pagano P., Reale F., Hood A. W., 2023, arXiv e-prints, p. arXiv:2306.06047
  * De Moortel & Nakariakov (2012) De Moortel I., Nakariakov V. M., 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 3193
  * Fisher et al. (1985) Fisher G. H., Canfield R. C., McClymont A. N., 1985, The Astrophysical Journal, 289, 414
  * Goddard et al. (2017) Goddard C. R., Pascoe D. J., Anfinogentov S., Nakariakov V. M., 2017, Astronomy and Astrophysics, 605, A65
  * Goddard et al. (2018) Goddard C. R., Antolin P., Pascoe D. J., 2018, The Astrophysical Journal, 863, 167
  * Gudiksen & Nordlund (2005) Gudiksen B. V., Nordlund Å., 2005, The Astrophysical Journal, 618, 1020
  * Gudiksen et al. (2011) Gudiksen B. V., Carlsson M., Hansteen V. H., Hayek W., Leenaarts J., Martínez-Sykora J., 2011, Astronomy and Astophysics, 531, A154
  * Guo et al. (2023) Guo M., Duckenfield T., Van Doorsselaere T., Karampelas K., Pelouze G., Gao Y., 2023, The Astophysical Journal Letters, 949, L1
  * Hansteen et al. (2015) Hansteen V., Guerreiro N., De Pontieu B., Carlsson M., 2015, The Astrophysical Journal, 811, 106
  * Hansteen et al. (2023) Hansteen V. H., Martinez-Sykora J., Carlsson M., De Pontieu B., Gošić M., Bose S., 2023, The Astrophysical Journal, 944, 131
  * Hirayama (1974) Hirayama T., 1974, Solar Physics, 34, 323
  * Howson & De Moortel (2022) Howson T. A., De Moortel I., 2022, Astronomy and Astophysics, 661, A144
  * Howson et al. (2020) Howson T. A., De Moortel I., Fyfe L. E., 2020, Astronomy and Astrophysics, 643, A85
  * Johnston & Bradshaw (2019) Johnston C. D., Bradshaw S. J., 2019, The Astrophysical Journal Letters, 873, L22
  * Johnston et al. (2017) Johnston C. D., Hood A. W., Cargill P. J., De Moortel I., 2017, Astronomy and Astrophysics, 597, A81
  * Johnston et al. (2019) Johnston C. D., Cargill P. J., Antolin P., Hood A. W., De Moortel I., Bradshaw S. J., 2019, Astronomy and Astrophysics, 625, A149
  * Johnston et al. (2021) Johnston C. D., Hood A. W., De Moortel I., Pagano P., Howson T. A., 2021, Astronomy and Astrophysics, 654, A2
  * Kanella & Gudiksen (2018) Kanella C., Gudiksen B. V., 2018, Astronomy and Astrophysics, 617, A50
  * Kanella & Gudiksen (2019) Kanella C., Gudiksen B. V., 2019, Astronomy and Astrophysics, 621, A95
  * Keppens et al. (2023) Keppens R., Popescu Braileanu B., Zhou Y., Ruan W., Xia C., Guo Y., Claes N., Bacchini F., 2023, Astronomy and Astrophysics, 673, A66
  * Klimchuk (2015) Klimchuk J. A., 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140256
  * Klimchuk et al. (2008) Klimchuk J. A., Patsourakos S., Cargill P. J., 2008, The Astrophysical Journal, 682, 1351
  * Kuniyoshi et al. (2023) Kuniyoshi H., Shoda M., Iijima H., Yokoyama T., 2023, The Astrophysical Journal, 949, 8
  * Li et al. (2023a) Li D., Bai X., Tian H., Su J., Hou Z., Deng Y., Ji K., Ning Z., 2023a, arXiv e-prints, p. arXiv:2306.04973
  * Li et al. (2023b) Li X., Keppens R., Zhou Y., 2023b, The Astrophysical Journal Letters, 947, L17
  * Lionello et al. (2009) Lionello R., Linker J. A., Mikić Z., 2009, The Astrophysical Journal, 690, 902
  * Martínez-Sykora et al. (2023) Martínez-Sykora J., De Pontieu B., Hansteen V. H., Testa P., Wargnier Q. M., Szydlarski M., 2023, The Astrophysical Journal, 949, 112
  * Mikić et al. (2013) Mikić Z., Lionello R., Mok Y., Linker J. A., Winebarger A. R., 2013, The Astophysical Journal, 773, 94
  * Montes-Solís & Arregui (2019) Montes-Solís M., Arregui I., 2019, Astronomy and Astrophysics, 622, A88
  * Nakariakov & Kolotkov (2020) Nakariakov V. M., Kolotkov D. Y., 2020, ARA&A, 58, 441
  * Nakariakov & Ofman (2001) Nakariakov V. M., Ofman L., 2001, Astronomy and Astrophysics, 372, L53
  * Nakariakov & Verwichte (2005) Nakariakov V. M., Verwichte E., 2005, Living Reviews in Solar Physics, 2, 3
  * Nakariakov et al. (1999) Nakariakov V. M., Ofman L., Deluca E. E., Roberts B., Davila J. M., 1999, Science, 285, 862
  * Nisticò et al. (2013) Nisticò G., Nakariakov V. M., Verwichte E., 2013, Astronomy and Astrophysics, 552, A57
  * Nóbrega-Siverio & Moreno-Insertis (2022) Nóbrega-Siverio D., Moreno-Insertis F., 2022, The Astrophysical Journal Letters, 935, L21
  * Pascoe et al. (2017) Pascoe D. J., Goddard C. R., Anfinogentov S., Nakariakov V. M., 2017, Astronomy and Astrophysics, 600, L7
  * Pascoe et al. (2018) Pascoe D. J., Anfinogentov S. A., Goddard C. R., Nakariakov V. M., 2018, The Astrophysical Journal, 860, 31
  * Pascoe et al. (2022) Pascoe D. J., De Moortel I., Pagano P., Howson T. A., 2022, Monthly Notices of the Royal Astronomical Society, 516, 2181
  * Pelouze et al. (2023) Pelouze G., Van Doorsselaere T., Karampelas K., Riedl J. M., Duckenfield T., 2023, Astronomy and Astrophysics, 672, A105
  * Petrova et al. (2023) Petrova E., Magyar N., Van Doorsselaere T., Berghmans D., 2023, The Astrophysical Journal, 946, 36
  * Pontin et al. (2017) Pontin D. I., Janvier M., Tiwari S. K., Galsgaard K., Winebarger A. R., Cirtain J. W., 2017, The Astrophysical Journal, 837, 108
  * Reale (2014) Reale F., 2014, Living Reviews in Solar Physics, 11, 4
  * Reale et al. (2016) Reale F., Orlando S., Guarrasi M., Mignone A., Peres G., Hood A. W., Priest E. R., 2016, The Astrophysical Journal, 830, 21
  * Reid et al. (2020) Reid J., Cargill P. J., Hood A. W., Parnell C. E., Arber T. D., 2020, Astronomy and Astrophysics, 633, A158
  * Reid et al. (2021) Reid J., Cargill P. J., Johnston C. D., Hood A. W., 2021, Monthly Notices of the Royal Astronomy Society, 505, 4141
  * Reid et al. (2023) Reid J., Threlfall J., Hood A. W., 2023, Monthly Notices of the Royal Astronomy Society, 518, 1584
  * Robinson et al. (2022) Robinson R. A., Carlsson M., Aulanier G., 2022, Astronomy and Astrophysics, 668, A177
  * Ruderman & Erdélyi (2009) Ruderman M. S., Erdélyi R., 2009, Space Science Reviews, 149, 199
  * Skirvin et al. (2023) Skirvin S. J., Gao Y., Van Doorsselaere T., 2023, The Astrophysical Journal, 949, 38
  * Soler et al. (2010) Soler R., Arregui I., Oliver R., Ballester J. L., 2010, The Astrophysical Journal, 722, 1778
  * Tian & Chen (2018) Tian H., Chen N. H., 2018, The Astrophysical Journal, 856, 34
  * Van Damme et al. (2020) Van Damme H. J., De Moortel I., Pagano P., Johnston C. D., 2020, Astronomy and Astrophysics, 635, A174
  * Van Doorsselaere et al. (2020) Van Doorsselaere T., et al., 2020, Space Science Reviews, 216, 140
  * Viall et al. (2021) Viall N. M., De Moortel I., Downs C., Klimchuk J. A., Parenti S., Reale F., 2021, in Raouafi N. E., Vourlidas A., eds, Geophysical Monograph Series Vol. 1, Solar Physics and Solar Wind. p. 35, doi:10.1002/9781119815600.ch2
  * Vögler et al. (2005) Vögler A., Shelyag S., Schüssler M., Cattaneo F., Emonet T., Linde T., 2005, Astronomy and Astophysics, 429, 335
  * Wang et al. (2003) Wang T. J., Solanki S. K., Curdt W., Innes D. E., Dammasch I. E., Kliem B., 2003, Astronomy and Astrophysics, 406, 1105
  * Warren et al. (2020) Warren H. P., et al., 2020, The Astrophysical Journal, 896, 51
  * Wright & Smith (1990) Wright A. N., Smith P. R., 1990, Journal of Geophysical Research, 95, 3745
  * Zhong et al. (2022) Zhong S., Nakariakov V. M., Kolotkov D. Y., Anfinogentov S. A., 2022, Monthly Notices of the Royal Astronomical Society, 513, 1834
  * Zhou et al. (2021) Zhou Y.-H., Ruan W.-Z., Xia C., Keppens R., 2021, Astronomy and Astrophysics, 648, A29