11institutetext: School of Mathematics and Statistics, University of St.
Andrews, St. Andrews, Fife, KY16 9SS, U.K. 22institutetext: Rosseland Centre
for Solar Physics, University of Oslo, PO Box 1029 Blindern, NO-0315 Oslo,
Norway

# The effects of driving time scales on coronal heating in a stratified
atmosphere

T. A. Howson 11 I. De Moortel 1122

###### Abstract

Aims. We investigate the atmospheric response to coronal heating driven by
random velocity fields with different characteristic time scales and
amplitudes.

Methods. We conducted a series of three-dimensional magnetohydrodynamic
simulations of random driving imposed on a gravitationally stratified model of
the solar atmosphere. In order to understand differences between alternating
current (AC) and direct current (DC) heating, we considered the effects of
changing the characteristic time scales of the imposed velocities. We also
investigated the effects of the magnitude of the velocity driving.

Results. In all cases, complex foot point motions lead to a proliferation of
current sheets and energy dissipation throughout the coronal volume. For a
given driving amplitude, DC driving typically leads to a greater rate of
energy injection when compared to AC driving. This ultimately leads to the
formation of larger currents, increased heating rates and higher coronal
temperatures in DC simulations. There is no difference in the spatial
distribution of energy dissipation across simulations, however, energy release
events in AC cases tend to be more frequent and last for less time than in DC
cases. This results in more asymmetric temperature profiles for field lines
heated by AC drivers. Higher velocity driving is associated with larger
currents, higher temperatures and the corona occupying a larger fraction of
the simulation volume. In all cases, the majority of heating is associated
with small energy release events, which occur much more frequently than larger
events.

Conclusions. When combined with observational results that highlight a greater
abundance of oscillatory power in lower frequency modes, these findings
suggest that energy release in the corona is more likely to be driven by
longer time scale motions. In the corona, AC and DC driving will occur
concurrently and their effects remain difficult to isolate. The distribution
of field line temperatures and the asymmetry of temperature profiles may
reveal the frequency and longevity of energy release events and therefore the
relative importance of AC and DC heating.

##  1 Introduction

Despite extensive research over much of the past century, the specific
processes which maintain the high temperatures of the solar corona remain
poorly constrained. It is widely accepted that energy is injected into the
Sun’s atmosphere by the convective flows that are observed at the photosphere.
However, the nature of the subsequent transmission and dissipation of this
energy is still unclear. This is known as the coronal heating problem and
research continues into the energy dissipation mechanism(s) and the associated
mass cycle between different layers of the atmosphere (see reviews by Narain &
Ulmschneider, 1990, 1996; Klimchuk, 2006; Parnell & De Moortel, 2012; Reale,
2014; Sakurai, 2017).

The majority of the mechanisms proposed in response to the coronal heating
problem fall into one of two broad categories, namely; AC heating and DC
heating. This dichotomy arises according to the time scales of the velocity
flows that drive energy into the corona. In particular, if the driving time
scales are short (in relation to the Alfvén crossing time along a given
coronal structure), then we will see AC heating (e.g. see reviews by Arregui,
2015; Van Doorsselaere et al., 2020) and if the driving time scales are long,
then we will see DC heating (e.g., see review by Wilmot-Smith, 2015). With
regards to nomenclature, we note that AC and DC heating mechanisms are
sometimes referred to as _wave-based_ and _reconnection-based_ , respectively.
However, we note that MHD wave modes are able to drive magnetic reconnection
(e.g. Sakai et al., 1984; McLaughlin et al., 2009; Howson et al., 2021) and
magnetic reconnection can drive oscillations, shocks and wave heating (e.g.
Bárta et al., 2007; Hsieh et al., 2009; Kigure et al., 2010). For either AC or
DC heating, due to the high magnetic and viscous Reynolds numbers $(\gg 1)$
that are expected for typical coronal conditions, energy must be transported
to small scales in the magnetic and/or velocity fields before it can be
efficiently dissipated as heat.

Several processes have been proposed to promote a transfer of energy to small
scales. In the context of DC heating, these include braiding (e.g. Parker,
1972; Van Ballegooijen, 1988; Peter et al., 2004; Wilmot-Smith et al., 2011),
flux tube twisting (Bareford et al., 2013; Gordovskyy et al., 2016) and/or MHD
avalanches (e.g. Hood et al., 2016; Reid et al., 2018). Meanwhile, for AC
heating, these processes include resonant absorption (e.g. Ionson, 1978;
Davila, 1987; Poedts et al., 1989; Ofman et al., 1998) and phase mixing (e.g.
Hasegawa & Chen, 1974; Heyvaerts & Priest, 1983; Parker, 1991). Moreover, many
studies have highlighted the propensity of MHD waves to trigger the
development of turbulent-like flows either through the interaction of counter-
propagating modes (e.g. Iroshnikov, 1964; Matthaeus et al., 1999; van
Ballegooijen et al., 2011) or through the onset of dynamic instabilities (e.g.
Browning & Priest, 1984; Terradas et al., 2008; Antolin et al., 2014; Hillier
et al., 2019). For these cases, energy will cascade to smaller scales until it
inevitably reaches the dissipation length scale whereupon it will be converted
to heat.

The dissipation mechanism notwithstanding, it remains unclear whether coronal
oscillations carry sufficient energy to balance expected losses (for example,
see different estimates by. Tomczyk et al., 2007; De Pontieu et al., 2007;
McIntosh et al., 2011; Thurgood et al., 2014; Srivastava et al., 2017). In
particular, in closed coronal loops, the interaction between an imposed driver
with a reflected wave will result in the time-averaged energy injection rate
for simple, sinusoidal wave drivers being too low unless the driver frequency
is resonant, the amplitude is very large or transport coefficients are
significantly enhanced above classical values (e.g. Howson et al., 2019, 2021;
Prokopyszyn et al., 2019). Despite this, periodic motions with a wide range of
frequencies have been detected (e.g. Morton et al., 2015, 2016), and
significantly for this study, there is much more power in low frequency modes
(similar to DC models).

An important difference between many models of AC and DC heating arises
according to whether the system will store energy in the background coronal
magnetic field (rather than simply in the perturbed component). As MHD wave
energy is typically equipartitioned between the kinetic and magnetic
components (over the duration of the a wave period), the amount of energy in a
typical AC heating model is constrained by the wave amplitudes, which can
often be observed directly. However, DC heating models will typically induce
an increase in the background magnetic energy which will ultimately allow a
sustained Poynting flux of energy into the system. This is difficult to
measure directly and thus leads to large uncertainties in the energy content
of the corona. There is of course no reason why short time scale driving
cannot lead to the storage of magnetic energy in the background field, and
this study eliminates this difference to allow a fair comparison between AC
and DC models. This approach was also followed in Howson et al. (2020) where
similar driving was imposed on a potential coronal arcade. In this previous
study, we found that DC driving produced greater energy injection rates,
current formation and ultimately, higher temperatures.

Understanding the mass-cycle that results from the coupling of the corona to
lower layers of the atmosphere remains a significant but important challenge.
In particular, the small length scales that exist in the transition region
present difficulties for accurately reproducing the mass flux into the corona
during heating events (Bradshaw & Cargill, 2013). Overcoming this problem with
innovative numerical techniques (e.g. Lionello et al., 2009; Johnston &
Bradshaw, 2019; Johnston et al., 2021) is critical as, despite the low
plasma-$\beta$ conditions in the corona, cooling times and any synthetic
emission produced from simulations are highly sensitive to the evolution of
the coronal density (e.g. Antiochos & Sturrock, 1978; Bradshaw & Cargill,
2010; Bradshaw & Klimchuk, 2011; Winebarger et al., 2018). Furthermore, the
efficiency of wave heating models critically depends on this density profile
and it remains unclear whether the heating is able to self-consistently
generate and/or maintain the density structuring (Cargill et al., 2016; Van
Damme et al., 2020).

In this paper we investigate the effects of driving time scales on energy
release within a stratified atmosphere that includes the coupling between the
corona and the chromosphere. In Sect. 2, we outline our model, including the
initial conditions and the driving profiles. Then, in Sect. 3, we present our
results and, finally, in Sect. 4, we discuss the implications of our findings
in the context of coronal heating research.

##  2 Numerical method

The numerical simulations presented within this article were conducted with
the Lagrangian-Remap code, Lare3d (Arber et al., 2001). The scheme advances
the full, resistive, three dimensional, MHD equations in normalised form. They
are given by

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

Here, $\rho$ is the plasma density, $\@vec{v}$ is the velocity, $\@vec{j}$ is
the current density, $\@vec{B}$ is the magnetic field, $\@vec{g}$ is the
acceleration due to gravity, $P$ is the gas pressure, $\epsilon$ is the
specific internal energy density, $k_{B}$ is the Boltzmann constant and $n$ is
the number density. In the energy equation (3), the second and third terms on
the right hand side correspond to thermal conduction and optically thin
radiation, respsectively. In the conduction term, $\@vec{F}_{C}$ is the
Spitzer-Härm heat flux. In the radiation term, $\Lambda(T)$, is a piecewise
continuous function that approximates radiative losses in an optically thin
plasma. This function is described in Klimchuk et al. (2008). We impose a
temperature floor of $2\times 10^{4}$ K which prevents the plasma cooling
below this temperature. In practice, this places a lower bound on the
chromospheric temperature which is acceptable as we are not fully simulating
the chromospheric physics (instead it acts as a mass reservoir).

In order to alleviate the difficulties of numerically resolving the conductive
flux in the spatially under-resolved transition region, we employ a technique
outlined in Linker et al. (2001); Lionello et al. (2009); Mikić et al. (2013).
In particular, adapted thermal conduction and radiative loss terms are imposed
below a cutoff temperature, $T_{c}$, which we set to be $2.5\times 10^{4}$ K.
This has the effect of artificially broadening the transition region and has
been shown to accurately simulate coronal temperature and density profiles
(when compared to high resolution 1D models) in both static (Lionello et al.,
2009) and dynamic loops (Mikić et al., 2013).

In these equations, we have included the resistivity, $\eta$ and viscosity,
$\nu$ as non-ideal terms which dissipate energy from the magnetic and velocity
fields, respectively. The viscosity is a sum of contributions from two small
shock viscosity terms which are included within all simulations to ensure
numerical stability. Together, these contribute a force,
<EMAIL_ADDRESS>on the right-hand side of the equation of motion 2
and a heating term, $Q_{\text{visc}}$ to the energy equation (3). These terms
are discussed in detail in Arber (2018). Using the notation detailed in the
referenced manual, we have set $\text{visc1}=0.05$ and $\text{visc2}=0.25$.

As we are not considering the correct physical treatment of the chromosphere,
we seek to only allow significant energy dissipation within the coronal
portion of the domain (see below). As such, we set $\eta=0$ in the volume
close to the upper and lower boundaries. This has the added benefit of
minimising the slippage of magnetic foot points through the imposed driving.
More specifically, the resistivity is defined to be

$\eta=\eta_{0}\left\\{\tanh(z+\alpha)-\tanh(z-\alpha)\right\\},$ (6)

where $\alpha=z_{\text{max}}-9$ Mm and $\eta_{0}$ is a constant that yields a
characteristic magnetic Reynolds number of $5\times 10^{3}$ in the coronal
volume (for the current numerical resolution and characterstic length scales).
The scheme employed here does not enforce energy conservation and, thus, any
temperature increase is not associated with numerical dissipation.

Figure 1: Initial density (blue) and temperature (red) profiles as a function
of $z$. The simulations are initially invariant in the $x$ and $y$ directions.

###  2.1 Initial conditions

For this study, we have simplified the geometry of the corona by modelling
curved loops as straight magnetic structures which are embedded in a dense
chromosphere at both foot points. The initial conditions consisted of a
uniform magnetic field with a strength of 20 G. The field is initially
parallel to the $z$-axis. We considered a numerical domain with dimensions of
-10 Mm $\leq x\leq$ 10 Mm, -10 Mm $\leq y\leq$ 10 Mm and -30 Mm $\leq z\leq$
30 Mm. The simulations used a numerical domain of $256\times 256\times 512$
grid points, in the $x,y$ and $z$ directions, respectively. In these
simulations, we assume that the graviational acceleration, $\@vec{g}$ (see
equation 2) is field aligned, and is modified with height as would be expected
in a semi-circular loop. In particular, the parallel component of gravity is
given by

$g_{\parallel}=g_{0}\sin\left(\frac{\pi z}{2z_{\text{max}}}\right),$ (7)

where $g_{0}\approx 274\text{ m s}^{-2}$ is the acceleration due to gravity at
the solar surface. This acceleration is vertically downwards for $z<0$ and
upwards (towards the loop foot point at z=$z_{\text{max}}$) for $z>0$.

Our simulations begin with an imposed temperature profile, $T=T(z)$, given by

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

where $T_{ch}=2\times 10^{4}$ K is the chromospheric temperature,
$T_{\text{co}}=10^{6}$ K is the initial coronal temperature, $a=22.5$ Mm
controls the location of the transition region and $b=0.3$ Mm controls the
width of the transition region. The temperature profile is displayed as the
red line in Fig. 1.

We define initial conditions that are in hydrostatic equilibrium using

$\frac{\mathrm{d}P}{\mathrm{d}z}=-\rho g_{\parallel},$ (9)

to define the density, $\rho=\rho(z)$. Here we use equations (5) and (8) to
eliminate the pressure and temperature and then integrate numerically using
the 4th order Runge-Kutta scheme. The logarithm of the resulting density
profile is shown by the blue line in Fig. 1.

As there is no background heating term included in these conditions, the
system is not in an initial thermodynamic equilibrium. In particular, the
coronal volume will lose energy to the lower layers of the atmosphere through
thermal conduction and all plasma above the temperature floor ($2\times
10^{4}$ K) will cool due to optically thin radiative losses. Instead, we aim
to see whether the driving can support a hot corona and how the long-term
state differs for different imposed velocity profiles.

Figure 2: Imposed velocity drivers. _Upper left:_ Schematic of the drivers
imposed on each of the foot points of magnetic field lines. _Upper right:_
Imposed velocity field on the lower boundary of the medium amplitude AC
simulation. _Centre:_ Temporal evolution of the magnitude of the imposed
driver on the lower boundary of the medium amplitude AC simulation on the line
$y=0$ Mm. _Lower:_ Temporal evolution of the magnitude of the imposed driver
on the lower boundary of the medium amplitude DC simulation on the line $y=0$
Mm.

###  2.2 Boundary conditions

We model the convective flows that exist at the photosphere by imposing a
transverse, space- and time-dependent velocity profile at both the upper and
lower $z$ boundaries of the domain. We note that our simulations do not
include the full-complexity of chromospheric physics and as such the exact
manner in which energy is transmitted through the lower atmosphere is not
considered here. As described in Howson et al. (2020), we define a driver
using a sum of many individual two-dimensional Gaussians, each of which has a
particular amplitude, direction, length scale and time scale. In particular,
on both photospheric boundaries ($z=\pm 30$ Mm), we impose
$\@vec{v}=(v_{x},v_{y},0)$, where

$\displaystyle
v_{x}=\displaystyle\sum_{i=1}^{N}v_{i}\cos{\theta_{i}}\exp\left\\{\frac{-(r-r_{i})^{2}}{l_{i}^{2}}\right\\}\exp\left\\{\frac{-(t-t_{i})^{2}}{\tau_{i}^{2}}\right\\},$
(10) $\displaystyle
v_{y}=\displaystyle\sum_{i=1}^{N}v_{i}\sin{\theta_{i}}\exp\left\\{\frac{-(r-r_{i})^{2}}{l_{i}^{2}}\right\\}\exp\left\\{\frac{-(t-t_{i})^{2}}{\tau_{i}^{2}}\right\\}.$
(11)

Here, for each $i$ in the summation, $v_{i}$ is the amplitude of the velocity
component, $\theta_{i}$ defines the direction of the component, $r_{i}$ is the
centre of a two-dimensional Gaussian, $l_{i}$ is a parameter which defines the
spatial scales of the velocity driver (the width of the 2D Gaussian), $t_{i}$
is the time of peak amplitude for each component and $\tau_{i}$ defines the
time scale of each component.

For each velocity component, all parameters are randomly selected from some
statistical distribution. In particular, for all $i$, the $v_{i}$ are normally
distributed with mean $v_{\mu}$ and variance $v^{2}_{\mu}/25$, the
$\theta_{i}$ are uniformly distributed on the interval $\left[0,2\pi\right]$,
the $r_{i}$ are uniformly distributed over the driven boundaries, the $l_{i}$
are normally distributed with mean $L/4=2.5$ Mm and variance
$L^{2}/400=0.25\text{ Mm}^{2}$, the $t_{i}$ are uniformly distributed over the
duration of the simulation, and the $\tau_{i}$ are normally distributed with
mean $\tau_{\mu}$ and variance $\tau^{2}_{\mu}/16$. We note that smaller
values of $\tau_{\mu}$ create shorter time scales for the velocity driver (AC
driving) and larger values are associated with longer time scales (DC
driving). In equations 10 & 11, the number of terms included in the summation,
$N$, is selected to be a function of the typical time scale $\tau_{\mu}$ and
is chosen such that a similar number of components are active at all times.
This ensures that the complexity of the velocity driver is consistent between
different simulations in the parameter space.

In this article, we conduct a parameter study on the mean amplitude, $v_{\mu}$
and the velocity time scale $\tau_{\mu}$. We consider three different
amplitudes; low, medium and high, and two different characteristic time
scales; AC and DC driving. The velocity drivers imposed in the low, medium and
high amplitude simulations have a mean magnitude of approximately, 2, 4 and 6
$\text{ km s}^{-1}$, respectively. The velocity Gaussian components have an
average duration of approximately 70 s in the AC simulations and approximately
1400 s in the DC simulations. For reference, the initial Alfvén crossing time
between the two $z$ boundaries is approximately 450 s.

In Fig. 2, we show the key features of the imposed driver. In the upper left
hand panel, we show a schematic of the model with blue arrows representing the
magnetic field and red arrows representing the action of the velocity drivers
on the $z$ boundaries. In the upper right hand panel, we show a snapshot of
the velocity field on the lower boundary of the medium amplitude AC
simulation. In the lower two panels, we show the temporal evolution of the
imposed velocity on the medium amplitude AC (centre) and DC (lower)
simulations in the line $x=0$ Mm, $z=-30$ Mm. We note that the velocity
components have much longer lifetimes in the lower panel (DC) than in the
central panel. Movies of AC and DC driving are included in the files that
accompany this article.

In all simulations, the $x$ and $y$ boundaries are defined to be periodic. For
the $z$ boundaries, with the exception of the imposed velocity drivers
described above, zero-gradient conditions are enforced for all variables.

##  3 Results

We begin by describing the general evolution of the medium amplitude AC
simulation. We initially focus on the general properties of the system that
persist across all simulations. The specific differences that result from
different velocity drivers will be detailed in subsequent sections.

In all cases, the imposed driving at field line foot points stresses the
magnetic field and induces the formation of currents throughout the domain. In
the resistive subvolume (see equation 6), this then permits energy release
through Ohmic heating and magnetic reconnection, leading to an increase in the
plasma temperature. Conductive fronts then transfer this heat along field
lines to the dense transition region and chromosphere. This generates an
increase in the gas pressure and propels the evaporation of plasma into the
coronal volume. The higher temperatures/densities enhance the
conductive/radiative losses in the corona and, as a result, the temperature
does not continuously increase. Although the driving, heating and cooling
processes are inherently non-uniform in space and time, after a long time (in
terms of the Alfvén travel time and the driving time scales), a statistically
steady-state is obtained, where the coronal temperature, mass content and
energy release rate remain relatively constant in time.

Figure 3: Horizontal cuts of the modulus of the current density, the
temperature and the density at $t=10900$ s in the medium amplitude AC
simulation. We have normalised the current to the maximum value in this plane.
The temporal evolution of this figure is shown in a movie that accompanies
this article. Figure 4: The Ohmic heating ($\eta j^{2}$; left), temperature
(centre) and density (right) averaged in $y$ and shown as a function of $x$
and $z$. For the density, the logarithm is plotted and the colour bar is
saturated at $10^{-10}\text{ kg m}^{-3}$ in order to show the details of
coronal structures. The panels shown here correspond to the medium amplitude
AC simulation.

In Fig. 3, we show horizontal cuts of the modulus of the current density
(upper), the temperature (centre) and the density (lower) at $z=0$ Mm
(simulation midplane) at $t=10900$ s (approx. 3 hours). The temporal evolution
of these quantities is shown in the accompanying movie. We see that the
currents form narrow sheets that are distributed throughout the entire cross-
section. However, despite their spatial spread, the largest values are very
localised, which leads to distinct heating events. In the movie, we see that
these are also short-lived in time, producing a burtsy, impulsive heating
profile even in this AC driving case.

In this resistive regime, the currents lead to energy release and an increase
in the plasma temperature. As such, it is not surprising that the temperature
profile (central panel of Fig. 3) shows signatures of features that can be
identified in the profile of $|\@vec{j}|$. For example, we see a high current
and temperature region in the area around $x=-5$ Mm, $y=-2$ Mm and similar
ribbon structures at $x=-5$ Mm, $y=6$ Mm. However, it is clear that the
temperature features tend to be more diffuse. There are several reasons for
this. Firstly, as plasma does not cool instantaneously, the temperature will
show regions of previously heated plasma, which are no longer associated with
large currents in the upper panel. Secondly, thermal conduction plays an
important role. Whilst only negligible amounts of energy are transferred
across field lines, it will efficiently move heat along field lines. As such,
the temperature profile can show signatures of heating (and currents) at
different altitudes (that are not seen in the upper panel of Fig. 3).
Furthermore, at this time many field lines have a significant component
parallel to this plane (due to the sustained driving), and conduction will
therefore transfer heat across the plane. Finally, as the plasma-$\beta$ is
non-zero, the heating will lead to cross-field pressure forces which will
slightly modify the distribution of hot plasma.

The density profile (bottom panel of Fig. 3) exhibits different structures to
those observed in the other two panels. Indeed, there is some out-of-phase
behaviour between the temperature and density panels. In particular, the high
density region at $x=-3$ Mm, $y=4$ Mm, is at a low temperature. This is
characteristic of the formation of a coronal condensation. As the plasma is at
the field line apex, there is no gravitational force to pull the dense plasma
towards either of the magnetic foot points. Therefore, the plasma can remain
relatively stable at the midplane for a long time before eventually falling
towards one of the foot points. The radiative losses here will be enhanced
(due to the high density and cool temperature) and this increased cooling rate
will further reduce the temperature, lowering the gas pressure and drawing
more plasma into the condensation.

In Fig. 4, we display vertical profiles of the Ohmic heating rate ($\eta
j^{2}$; left), the temperature (centre) and the density (right) for this
simulation at the same time ($t=10900$ s). For all three panels we have
averaged along the $y$ axis to provide an insight into the nature of the
variables throughout the whole domain. The mean Ohmic heating is normalised to
the maximum value obtained during the simulation. For clarity, we plot the
logarithm of the mean density and note that the chromospheric density is
saturated to allow the coronal structure to be seen. The temporal evolution of
this figure is shown in the accompanying movie.

For the Ohmic heating (left panel), the $\eta=0$ regions are apparent for
$|z|\gtrsim 21$ Mm. We note that large currents do develop within these
regions but they do not contribute directly to plasma heating. We see that
currents tend to form in long narrow structures that are aligned with magnetic
field lines and that the highest Ohmic heating rates occur in small, intense
bursts. These energy release events are also short-lived in time (see
accompanying movie) and are reminiscent of nanoflares.

In the central panel we see that a corona with a temperature of approximately
0.6 MK is maintained for this simulation (more energetic simulations are
discussed in more detail below). The temperature of the chromosphere remains
at $2\times 10^{4}$ K, partly due to the lack of heating here but mainly due
to the huge radiative losses that any heated plasma will experience here. As
with the Ohmic heating profile, the temperature also shows narrow structuring
that is aligned with the magnetic field. This is partly due to the nature of
the heating but also due to thermal conduction spreading thermal energy along
(and not across) magnetic field lines. We note that it is difficult to
identify the signatures of individual heating events once the integration in
the $y$ direction has been performed.

In the right hand panel, we see that the heating in this simulation is
sufficient to sustain a coronal density of approximately $10^{-12}\text{ kg
m}^{-3}$ with some dense structures forming in response to the driving. The
atmosphere remains well stratified and, once again, as plasma evaporation does
not occur instantaneously with heating events, there is little connection
between the density and the temperature/heating profiles. We note that the
density enhancement of the visible structures in relation to the ambient
corona often exceeds an order of magnitude. This is interesting in the context
of wave heating models which typically require a significant density gradient
across coronal flux tubes. However, in the accompanying movie, we note that
the evolution of the density (as well as the other two variables) is highly
dynamic. This is a consequence of the short time scaling driving considered in
this simulation and it is unclear whether the density structuring is
sufficiently long-lived to allow for classical wave heating mechanisms (e.g.
phase mixing and resonant absorption) to release energy efficiently. Given
these structures will have longer lifetimes in the DC simulations, a
combination of AC and DC driving may produce classical wave heating processes
(such as phase mixing and resonant absorption) in these structures.

###  3.1 Poynting flux

The imposed drivers inject energy into the system at both the upper and lower
$z$ boundaries of the domain. Since no perpendicular flows are permitted
through the boundary, the only energy flux term is given by the Poynting flux.
Meanwhile, there is an energy loss term that is associated with the radiative
losses. Therefore, in terms of the volume integrated energy, $\Xi$, we have

$\frac{\mathrm{d}\Xi}{\mathrm{d}t}+\frac{1}{\mu_{0}}\int\@vec{E}\times\@vec{B}\cdot\mathrm{d}\@vec{S}=-\int
R_{L}\,\,\mathrm{d}V.$ (12)

Here, the first integral represents the Poynting flux where, $\mu_{0}$ is the
magnetic permeability in a vacuum, $\@vec{E}$ is the electric field and the
integral is calculated over both the upper and lower $z$ boundaries. The
second term on the right-hand side represents the volume integrated radiative
losses. Here we recall that the energy lost by plasma at $2\times 10^{4}$ K is
effectively zero (due to the inclusion of the temperature floor). For the
current geometry, the energy injected through the driven boundaries can be
expressed as (e.g. Parnell & De Moortel, 2012; Howson et al., 2020)

$\int_{A_{1}}\left(v_{x}B_{x}+v_{y}B_{y}\right)B_{z}\,\,\mathrm{d}A\,-\int_{A_{2}}\left(v_{x}B_{x}+v_{y}B_{y}\right)B_{z}\,\,\mathrm{d}A,$
(13)

where the two area integrals are calculated over the upper and lower
boundaries respectively.

Figure 5: Schematic of the energetic consequences for different magnetic field
and driving profiles. With the current geometry in mind, we show driving on
the lower (left panels) and upper (right panels) $z$ boundaries.

In Fig. 5, we show a schematic of the implications of equation (13) on
different driving profiles. For simplicity, we consider a case with
$B_{y}=v_{y}=0$ and show transverse driving ($v_{x}$) at the lower (left hand
panels) and upper (right hand panels) boundaries. For the current setup, at
both boundaries we have $B_{z}>0$ and in Fig. 5, we show the different
possibilities for the case $B_{x}<0$. Ultimately, at either boundary, when the
driving is increasing the magnitude of the transverse component of the field
($\left|B_{x}\right|$), it is injecting energy into the system. Otherwise, if
the magnitude of the transverse component of the field is decreased by the
driving, the imposed velocity will remove energy from the system. This
argument extends naturally to the case with $v_{y}\neq 0$.

As each of the velocity components in the driver definition (eq. 10) are
independent from each other, when any given term is switched on, it is as
likely to remove energy from the system as it is to inject new energy. It is
only the continual action of each driver component that allows a net influx of
energy. For example, a driver component acting in the positive $x$ direction
on the lower $z$ boundary will induce a negative $B_{x}$ term here (see upper
left hand panel of Fig. 5) and thus begin to inject energy into the system.
Therefore, in general, the longer each driving component acts for, the greater
the Poynting flux it will drive into the simulation domain. As such, we can
expect the long time scale cases (DC driving) to have a larger time-averaged
Poynting flux than the short time scale cases (AC driving). This phenomenon is
fundamentally responsible for the results that follow.

Figure 6: Poynting flux injected by medium amplitude AC (solid) and DC
(dashed) driving. In the upper panel we show the instantaneous flux and in the
lower panel we show the cumulative energy injection. In both cases, we have
normalised the variables by the maximum of the dashed lines.

In Fig. 6, we show the evolution of the Poynting flux for medium amplitude AC
(solid) and DC (dashed) driving. In the upper panel, we show the instantaneous
Poynting flux (equation 13) and in the lower panel we show the cumulative
Poynting flux. This second plot is the total energy injected since the start
of each simulation and is simply the time integral of the lines in the upper
panel. For both simulations, we see that after an initial increase, the
instantaneous Poynting flux reaches a relatively steady state (when viewed
over long enough time periods). As the mean value is positive in both cases,
this ensures that there is a constant flux of energy injected into each of the
simulations. However, over the experiment run-time, this does not necessarily
lead to increasingly energetic systems as the radiative loss term is effective
at removing energy, particularly in the dense chromosphere.

For the reasons discussed above, we see considerably higher energy injection
rates for the DC driving (dashed lines) than for the AC driving (solid lines).
We also see that the Poynting flux fluctuates rapidly in the AC simulations
and much more slowly in the DC cases. This is a direct effect of the temporal
variability associated with the different driving profiles. Despite this, the
cumulative Poynting flux (lower panel) increases approximately linearly for
both simulations.

Figure 7: Vertical energy flux within the simulation domain. For the medium
amplitude AC simulation, we show the temporally- and spatially- (in $x$ and
$y$) averaged total flux (black), Poynting flux (blue) and conductive flux
(red).

Once energy is injected into the simulation domains, it is transmitted through
the volume via a series of flux terms; namely, Poynting flux, kinetic energy
flux, enthalpy flux, gravitational potential flux and conductive flux. In Fig,
7, we show the combined effects of these terms on the transfer of energy in
the $z$ direction. To create this plot, we have averaged the $z$ component of
each energy flux term in the $x$ and $y$ directions and over the duration of
the simulation. We have then summed each term to give the mean combined flux
parallel to the $z$ axis (black line). Here, a positive value means energy is
flowing in the positive $z$ direction. We also show the contribution from the
Poynting flux (blue line) and the conductive flux (red line). For all three
curves, we have normalised by the maximum of the absolute value of the total
energy flux (black line).

We see that throughout much of the domain, the energy transfer is dominated by
the Poynting flux (compare black and blue lines). However, in the corona and
transition region, the conductive flux also acts as an important conduit for
transporting energy within the simulation volume. This relates to a coronal
cooling term which transfers energy from hot plasma at the loop apex to the
cooler layers of the atmosphere below. As the chromosphere is isothermal,
there is no conductive flux within this layer.

It is important to note that the rapid decrease in the absolute value of the
energy flux for $|z|\gtrsim$ 26 Mm means that much of the energy is not
transmitted to the corona in the central portion of the domain. Indeed only
about 20 % of the injected energy passes into the upper chromosphere. Some of
the remaining energy is stored in the perturbed magnetic field, however, much
is lost through radiative effects. For example, this can include slow wave
chromospheric damping. In reality, whilst a large proportion of mechanical
energy will not be passed into the corona, we cannot infer an accurate
estimate for the transmission rate from these simulations due to the lack of
chromospheric physics included in our model.

Figure 8: Volume integrated quantities for all six simulations. In the upper
two panels, we show the change in the volume integrated total energy (left)
and magnetic energy (right). In both cases, we show the percentage relative to
the initial total energy. We also show the volume integrated Ohmic heating
(lower left), normalised to the maximum rate in the high amplitude, DC
simulation. Finally, we show the maximum temperature at $z=0$ in the lower
right panel. In all figures, solid lines refer to AC cases and dashed lines
refer to DC cases. Each colour refers to a different driver amplitude.

###  3.2 Energetics and thermal evolution

The increased Poynting flux in DC driving cases provides a larger source of
energy that can be dissipated in order to heat the plasma. In this section, we
compare the effects of driving time scales and amplitudes on the components of
the volume-integrated energy and on the temperature of the coronal plasma.

In the upper left panel of Fig, 8, we show the change in the volume integrated
energy for all six of the simulations. The change is shown as a percentage of
the initial value, which is the same in all cases. For all panels in this
figure, solid lines correspond to AC simulations and dashed lines correspond
to DC cases. The different colours show the different driving amplitudes. The
red, green and blue lines show high, medium and low amplitude cases,
respectively. We note that the high amplitude DC simulation was run for a
shorter time (approximately 7000 s) than the other simulations, as the
computational time step became prohibitively short beyond this point.

As expected, we typically see that higher amplitude driving is associated with
a larger increase in the total energy. However, the random nature of the
imposed velocity means that this is not necessarily true at all times (e.g.
compare dashed blue and green lines). In most cases, the total energy
increases during the course of the simulation. Despite this, the low amplitude
AC driver (solid blue line) does not inject enough energy to balance the
thermal losses (radiation) from the simulation. All simulations show an
initial decrease in the total energy, although this occurs over a short time
period ($<500$ s) in the most energetic cases (red lines). This is because
thermal losses are important from the simulation onset, whereas it takes some
time for the drivers to begin to inject significant Poynting flux. In
agreement with the discussion in the previous section, we see that, for all
amplitudes, the DC driving consistently injects more energy than the AC
driving (compare solid and dashed lines). Further, the characteristic time
scales of each driver can be identified in this plot. In particular, the AC
curves (solid lines) exhibit higher frequency variation than the DC cases.

In the upper right panel of Fig, 8, we restrict our attention to the magnetic
component of the energy. Again, the lines are normalised by the maximum of the
initial (total) energy. By comparing the curves with the previous panel, we
see that the majority of the energy increase is stored in the magnetic field.
The initial, uniform field represents the minimum energy state for the amount
of magnetic flux that passes through the simulation domain. As no flux is
removed (or indeed added) to the simulation, the magnetic energy can only
increase from the initial state. As such, the change in magnetic energy
(relative to the initial value) is never negative in any of the simulations.
The injected magnetic energy is associated with a non-potential component of
the field which can be dissipated by resistive effects.

The lower left hand panel of Fig. 8 shows the Ohmic heating produced by the
dissipation of currents in the field. Here, we have normalised all curves by
the maximum of the dashed red curve. As a direct result of the increased
energy injection, we see more heating in higher amplitude cases and in the DC
simulations. Despite integrating over the entire simulation volume, the Ohmic
heating rate shows high temporal variation, particularly in the AC cases. This
suggests that the energy release is inherently intermittent and bursty.
However, we note that in all cases, there is a background threshold that the
heating never falls below. As such, there is a steady component to this
heating, at least when computed over a sufficiently large volume. It is likely
that this background threshold is a function of the transport coefficients and
may be much lower in higher Reynolds number plasmas.

In the lower right hand panel of Fig. 8, we show the mean temperatures in the
$z=0$ Mm plane (effectively the top of the corona) during all simulations.
Unsurprisingly, we see that the simulations with higher Ohmic heating rates
produce hotter coronae. In particular, DC driving cases support hotter
atmospheres than their AC driving counterparts. It is also of note that the
short time variations observed in the AC Ohmic heating profiles (solid lines
in the lower left panel) are reduced in the thermal evolution. As with the
comparison between the current and temperature in Fig. 3, this is largely
because the plasma does not cool instantaneously following a heating event.
Consequently, the temporal gradients are reduced.

Figure 9: The volume integrated Ohmic heating rate in $-1\leq x,y,z\leq 1$ Mm
for the medium amplitude AC (solid blue) and DC (dashed black) simulations.
Here we have normalised both lines by the maximum of the DC case (black).

As much of the local variation is lost by averaging over the entire simulation
volume, in Fig. 9, we show the local Ohmic heating rate as a function of time
in a small sub-volume given by $-1\leq x,y,z\leq 1$ Mm. For clarity, we only
show the medium amplitude AC (solid blue line) and DC (dashed black line
cases). In both cases, we have normalised the curves by the maximum of the
dashed black line. As expected, the mean of the DC heating rate is larger,
however, it only exceeds the AC heating rate at a few discrete, large scale
events. In both cases, the energy release is inherently bursty, although
again, there is a minimum threshold which effectively acts as a steady heating
term. The DC heating events tend to be larger and more long-lived but also
less frequent than the AC events. The importance of the frequency of heating
is reviewed in Reale (2010) and means we should expect that more of the
cooling phase will be observed for DC heating.

###  3.3 Field line analysis

As the thermodynamic evolution of coronal plasma is typically dominated by
field-aligned conduction, the response of the solar atmosphere to heating
events is frequently studied using 1D modelling. The complexity of these
simulations prohibits any single field line being tracked throughout the
experiments. However, we are able to provide insight into the characteristic
behaviour of different field lines by sampling many field lines at different
times. In the following, we have tracked $10^{4}$ field lines with foot points
uniformly distributed on the lower boundary at all output times during the
simulation. The starting points do not move in time and thus, as the field is
advected by the imposed drivers, we have not tracked the same field lines. The
simulations produce 400 data outputs and thus we have a sample of $4\times
10^{6}$ field lines from each experiment.

Figure 10: Mean temperature profile along field lines from all simulations.
Solid lines show AC cases and dashed lines show DC cases. Each colour refers
to different driver amplitudes. For each simulation, the averaging was
calculated for $10^{4}$ different field lines sampled every 35 s.

In Fig. 10, we show the average temperature profile along these field lines
for the different simulations. As the field lines are different lengths, we
have only considered the central 60 Mm of each field line (the minimum
possible length). The variable, $s$, measures the distance along each field
line and $s=0$, is defined to be at $z=0$ Mm. All six simulations are shown
with dashed and solid lines referring to AC and DC driving, respectively and
red, green and blue corresponding to high, medium and low amplitude driving.
The simulations with the higher apex temperatures (lower right hand panel of
Fig. 8) are hotter throughout the coronal volume but all cases retain the cool
chromospheres at the temperature floor of $2\times 10^{4}$ K. This is because
the high density here means that radiative losses are always in excess of any
heat that conducts from the corona. However, it is apparent that the hotter
simulations support an extended coronal volume as the transition regions move
closer to the simulation boundaries.

Figure 11: Temperature profile along two field lines from all the medium
amplitude AC (solid lines) and DC (dashed lines) simulations.

We note that the large temperature gradients that might be expected in the
transition regions are not visible in Fig. 10. This is partially due to the
numerical treatment of conduction employed here (see Sect. 2), but is also
caused by averaging over many field lines. Each field line has transition
regions at different $z$ values, according to its own thermal evolution. As
such, the wide spread of transition region locations produces the relatively
shallow average temperature profiles. In Fig. 11, we show the temperature
profile of field lines in the medium amplitude AC (solid lines) and DC (dashed
line) simulations. Here, the temperature jumps in the transition regions are
somewhat steeper (compare to green lines in Fig. 10) although they are still
reduced by the adjustment to thermal conduction. Although the DC simulations
are in general hotter, as we see here, this is not necessarily the case for
all field lines. As mentioned earlier, the DC heating events tend to have
lower frequencies and thus allow more time for field lines to cool than in the
AC simulations. We also note that the AC temperature profiles (solid lines)
are more asymmetric than for the DC cases. This is a consistent trend across
field lines in the simulations and is discussed in more detail below.

Figure 12: Properties of field lines sampled from all six simulations. _Upper
left:_ The mean Ohmic heating rate along field lines. _Upper right:_ The
distribution of field line lengths. _Lower left:_ The contribution of field
lines with different Ohmic heating rates to the total energy release. _Lower
right:_ The distribution of field line temperatures. In the upper left panel,
curves are normalised to the maximum of the dashed red line. The lines in all
other panels are normalised such the area under each curve is unity.

In Fig. 12, we show some general properties of the field lines traced in the
simulation results. In the upper left hand panel, we show the mean Ohmic
heating profile along all of the sampled field lines in each simulation. To
allow for different length field lines, we have normalised all lengths to
unity ($x$ axis) before computing the average. We have also normalised the
value of the mean Ohmic heating rates to the maximum of the red dashed curve.
In all cases the heating rate falls to zero close to the upper and lower
boundaries because the resistivity is set to zero here (eq. 6).
Unsurprisingly, we see that the more energetic simulations produce higher
Ohmic heating along the entire coronal section of the field lines. In all
cases, the average heating rate is approximately constant across this coronal
volume. This is largely a result of the initial uniform magnetic field and is
likely to change if the field strength changes in height (for example in an
arcade structure as used in Howson et al., 2020). Additionally, despite this
average uniformity, we note that heating on individual field lines can be
highly localised and thus non-uniform.

In the upper right hand panel of Fig. 12, we show the distribution of field
line lengths from each simulation. As the height of the box is 60 Mm and field
lines connect to opposite boundaries, the minimum length of a field line in
all cases is 60 Mm. Indeed, all field lines initially have this length.
However, as the driving progresses, and the magnetic field becomes stressed,
the length of field lines increases as different field lines take increasingly
convoluted paths between the upper and lower boundaries. Both higher amplitude
and longer time scale driving increase the mean length of field lines more
than lower amplitude and/or shorter time scale driving. This is because
individual components of the imposed velocities (see eqs. 10 & 11) will
displace magnetic foot points by a greater distance if the amplitude is larger
or they are active for longer (increased time scale). We note that this also
affects the width of the field line length distributions, with more energetic
simulations producing a greater spread. The greater field line lengths in DC
and high amplitude simulations will promote apex temperature increases due to
a decrease in thermal conductive losses to the lower atmosphere. However, as
the mean length increase is relatively small, this effect will be negligible
in comparison to the increased energy injection for these simulations.

Figure 13: Maximum temperature against the integrated Ohmic heating rate along
the length of each field line. The heating rate is normalised to the maximum
observed in the high amplitude DC simulation. Each panel shows the results of
a different simulation. Individual points each show different field lines and
the colour shows the simulation time when the field line was traced. Solid
lines show the expected scaling law assuming thermodynamic equilibrium. Dashed
lines show the scaling measured in each simulation.

In the lower left hand panel of Fig. 12, we show that the total heating is
mainly composed of heating from many field lines with weak heating rather than
the few field lines with strong heating. To generate this plot, we grouped
field lines according to the magnitude of their heating ($x$ axis normalised
by the maximum heating rate in each simulation). Then, we calculated the total
heating produced by field lines in each group ($y$ axis). Whilst the exact
shape of the distribution differs between simulations, it is clear that in all
cases, energy release is dominated by relatively weak, but numerous, events.

In the lower right hanx panel of Fig. 12, we show the distribution of the
maximum temperature of the sampled field lines. Once again, we see that the
more energetic simulations produce more field lines with higher temperatures.
Additionally, we note that the temperature distributions are positively skewed
with a large number of cool field lines and a relatively small number of
hotter field lines. This is associated with the Ohmic heating distributions
(lower left panel), which show that relatively few field lines have high
heating rates and thus not many field lines reach the very high temperatures.

Although the majority of energy release is associated with low energy events,
the highest temperatures in each simulation are only attained where the Ohmic
heating rate is very high. In Fig. 13, we show how the maximum temperature on
each of the traced field lines compares to the total Ohmic heating on that
particular field line (each point represents an individual field line). For
clarity, we only show a subsample of the field lines by selecting every tenth
snapshot from each simulation. We also omit the first 2000 seconds of
evolution as many of the field lines have very similar characteristics at
early times due to the uniformity of the initial conditions. In each panel, we
show the result of a different simulation and the colours of the points refer
to the simulation time at which the field line was traced. The solid line in
each panel shows the scaling predicted by the RTV scaling law (Rosner et al.,
1978), under the assumption that all field lines have the same length. The
dashed line shows the actual scaling found in each simulation. Whilst the
assumption on field line length is violated, the fractional variation is small
and this is not a large factor in the disagreement between the solid and
dashed lines. Instead, the discrepancy is a result of field lines not being in
thermodynamic equilibrium when the measurement is taken. In general, the DC
simulations (upper row) show better agreement with the predicted scaling law,
because their evolution is much less dynamic and thus remains closer to the
equilibrium state. In all cases, the observed power law is less steep than
expected (compare gradients of dashed and solid lines). This is largely
because previously heated field lines (where the Ohmic heating is now reduced)
can still show enhanced temperatures, as they do not cool instantaneously.

In agreement with previous plots, we see that the DC and higher amplitude
simulations exhibit hotter temperatures than other cases. We also see that the
DC simulations (except the high amplitude case with the reduced run time) show
a wider distribution in terms of the field line apex temperatures (also see
lower right hand panel of Fig. 12) than their AC counterparts. Further, as the
colour of each point refers to the simulation time, the temporal evolution can
be tracked by moving from blue to green to red points. The upper left hand
panel (high amplitude DC simulation) only contains blue points as the
simulation run time is much shorter than for other cases. In all other panels,
the difference between early and late times can be seen by comparing blue and
red points. We see that blue points (early times) generally lie at lower
temperatures and lower Ohmic heating rates as it takes time for the drivers to
produce the small scales required for efficient energy release. As such, the
simulations are cooling at early times (see also lower right hand panel of
Fig. 8). In all cases, the energy release requires complexity to be injected
by the imposed velocity drivers. Therefore, it typically requires many
different Gaussian components to have been activated (eqs. 10 & 11) before the
highest heating rates develop. This takes longer for long time scale drivers
and, as such, we see a more distinct population of blue points in the DC
panels (medium and low amplitudes) than in the AC cases. The only exception to
this behaviour is in the low amplitude AC case (lower right panel). This
simulation never reaches a heating phase as the Poynting flux is insufficient
to overcome the thermal losses.

Figure 14: Temporal evolution of the temperature in the line $x=z=0$ Mm for
the medium amplitude AC (upper) and DC (lower) simulations.

###  3.4 Distinguishing features of AC and DC heating

The different characteristic time scales of the velocity driving translate
into increased dynamism for the AC simulations within the coronal fraction of
the atmosphere. This can be seen in short time scale variation in the velocity
field and the density, for example. As such, synthetic observables generated
by these simulations can be used to easily distinguish between the two driving
time scales (Fyfe et al., 2021). This enhanced dynamism can also be observed
in the evolution of the apex temperature which is affected by the different
heating characteristics and the different advection profiles for the AC and DC
cases. In Fig. 14, we show the evolution of the temperature in the line
$x=z=0$ Mm for the medium amplitude AC (upper) and DC (lower) simulations. In
the upper panel, we see much higher frequency variation, lower mean
temperatures and no large scale, high temperature formation as can be seen in
the lower panel. This highlights the fundamental differences between the AC
and DC evolution. The infrequent but large heating events that characterise
the energy release in the DC simulations are easily identifiable in the lower
panel (e.g. for $y>0$ Mm at $t\approx 2500$ s). Whilst a wide range of
frequencies coexist in the Sun’s atmosphere, the differences between the
temperature evolution, suggests that these heating mechanisms may be
distinguishable.

Figure 15: Evolution of the mean asymmetry, $\Lambda$, of temperature profiles
along the traced field lines. Results are shown for high (red), medium (green)
and low (blue) amplitude AC (solid) and DC (dashed) simulations.

An additional characteristic which distinguishes between the heating in
different simulations is identified by considering the asymmetry of
temperature profiles along individual field lines. As mentioned previously,
Fig, 11 shows two field lines from AC simulations (solid lines) which have
more asymmetric temperature profiles than the two field lines from DC
simulations (dashed lines). In order to show that this is typical for field
lines in the different simulations, we define a measure of the temperature
asymmetry on a given field line by

$\Lambda=\frac{\left|\int T(s)-\int T(-s)\right|\,\,\mathrm{d}s}{\int
T(s)\,\,\mathrm{d}s}.$ (14)

Here, $s$ paramaterises the field line. In this equation, the numerator
represents the total difference between the temperature profile on either side
of the loop apex. Hence, it is zero if the temperature profile is symmetric.
The denominator then normalises the measure such that the asymmetry of field
lines with different mean temperatures can be compared.

In Fig. 15, we show the evolution of the mean (across all traced field lines)
of $\Lambda$ for the high (red), medium (green) and low (blue) amplitude AC
(solid) and DC (dashed) lines. We have normalised all curves by the maximum of
the solid red line. We see that the AC simulations produce more asymmetric
temperature profiles than the equivalent DC cases. After the initial
adjustment phase of the simulations (the initial conditions are symmetric),
the asymmetry remains roughly constant and the difference between the AC and
DC cases is consistent. Additionally, we also note that this measure of the
asymmetry in the field line temperature profiles is not affected by the
amplitude of the velocity drivers in each simulation. Therefore, measuring the
temperature asymmetry on coronal field lines may provide important information
about the driving time scales that are most relevant for energy release.

##  4 Discussion and Conclusions

In this article, we have presented the results of a series of large scale MHD
simulations of coronal heating in a gravitationally stratified solar
atmosphere. Energy is injected by perpendicular velocity drivers imposed at
magnetic foot points at both the upper and lower boundaries of the
computational domains. The fundamental characteristics (e.g. amplitude, length
scales, time scales) of these drivers are selected from statistical
distributions to produce random motions. These drivers are designed to mimic
the photospheric flows which inject energy into the Sun’s atmosphere.

From a parameter study on the driving time scales and on the velocity
amplitudes, we find that longer driving time scales and higher velocities
inject a greater Poynting flux into the system, produce larger currents and
ultimately dissipate more energy as heat. Typically, the simulations evolve
towards a steady state where the rate of energy dissipation in the corona is
balanced by the loss mechanisms (e.g. radiation and conduction to lower layers
of the atmosphere). Heating in DC simulations is characterised by larger scale
and longer lasting energy release events that occur more infrequently than in
AC experiments.

In contrast to many existing models of AC heating, the short time scale
driving in this setup stores large quantities of magnetic energy in the
coronal field. This typically enhances the Poynting flux and permits greater
energy injection than is possible in simulations where the driver only
perturbs the field about an initial equilibrium. As such, even in the AC
simulations (medium or high amplitude), the rate of energy injection is
sufficient to balance coronal losses. Despite this, as with all large scale 3D
MHD simulations, the viscous and magnetic Reynolds numbers that can be
obtained here are orders of magnitude smaller than might be expected in the
solar corona. As a result, it is uncertain whether the energy injected by
these drivers would be dissipated on sufficiently short time scales to be
relevant for coronal heating. Further, it remains unclear whether the relative
efficiency of AC and DC heating is constant across a large range of different
Reynolds numbers.

Another limitation of the current model is that it does not include a full and
accurate treatment of chromospheric physics. Indeed, in these simulations the
lower layers of the atmosphere simply act as a mass reservoir from which
plasma can evaporate in response to heating events. Despite this, and although
the transmission of energy through the chromosphere remains poorly understood,
coronal observations have measured oscillatory power at a range of different
frequencies (e.g. Morton et al., 2015, 2019). As such, we know that motions
with similar characteristic time scales to the driving presented here are
transmitted into the corona. In fact, these observational results detect
significantly more power at lower frequencies (similar to DC driving) than at
higher frequencies (AC driving). Whilst this may be partly due to
observational biases; e.g. power at higher frequencies being difficult to
observe due to smaller length scales and higher dissipation rates, it seems
that there is significantly more power in low frequency modes. As our
comparison, which assumes the same power for AC and DC modes, shows greater
energy dissipation in the DC simulations, these observations suggest that this
result may be enhanced further. As such, we expect DC heating mechanisms to
dominate energy release in the closed corona.

The evolution induced by the random driving is inherently complex and thus, it
is neither possible, nor particularly helpful, to analyse the nature of
individual heating events. Instead, we have focussed on volume integrated and
averaged quantities which provide information about the typical properties of
plasma parameters (e.g. temperature). One such variable which shows a clear
difference between AC and DC simulations is the asymmetry of the temperature
profiles along magnetic field lines. In particular, AC driving typically
produces temperature profiles which are, on average, more asymmetric than for
DC driving. Furthermore, this difference is independent of the magnitude of
the energy release. Therefore, whilst this parameter may be difficult to
measure in observational studies, it may provide a useful insight into the
fundamental nature of coronal heating.

_Acknowledgements._ The research leading to these results has received funding
from the UK Science and Technology Facilities Council (consolidated grant
ST/N000609/1), the European Union Horizon 2020 research and innovation
programme (grant agreement No. 647214). IDM received funding from the Research
Council of Norway through its Centres of Excellence scheme, project number
262622. This work used the DiRAC Data Analytic system at the University of
Cambridge, operated by the University of Cambridge High Performance Computing
Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This
equipment was funded by BIS National E-infrastructure capital grant
(ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC
DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National
e-Infrastructure.

## References

  * Antiochos & Sturrock (1978) Antiochos, S. K. & Sturrock, P. A. 1978, ApJ, 220, 1137
  * Antolin et al. (2014) Antolin, P., Yokoyama, T., & Van Doorsselaere, T. 2014, Astrophys. J. Letts., 787, L22
  * Arber (2018) Arber, T. D. 2018, LareXd User Guide, 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 (2015) Arregui, I. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140261
  * Bareford et al. (2013) Bareford, M. R., Hood, A. W., & Browning, P. K. 2013, A&A, 550, A40
  * Bárta et al. (2007) Bárta, M., Karlickỳ, M., Vrsnak, B., & Goossens, M. 2007, Central European Astrophysical Bulletin, 31
  * Bradshaw & Cargill (2010) Bradshaw, S. J. & Cargill, P. J. 2010, ApJ, 717, 163
  * Bradshaw & Cargill (2013) Bradshaw, S. J. & Cargill, P. J. 2013, ApJ, 770, 12
  * Bradshaw & Klimchuk (2011) Bradshaw, S. J. & Klimchuk, J. A. 2011, ApJS, 194, 26
  * Browning & Priest (1984) Browning, P. K. & Priest, E. R. 1984, A&A, 131, 283
  * Cargill et al. (2016) Cargill, P. J., De Moortel, I., & Kiddie, G. 2016, ApJ, 823, 31
  * Davila (1987) Davila, J. M. 1987, The Astrophysical Journal, 317, 514
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2007, Science, 318, 1574
  * Fyfe et al. (2021) Fyfe, L., Howson, T., & De Moortel, I. 2021, arXiv e-prints, arXiv:2110.14257
  * Gordovskyy et al. (2016) Gordovskyy, M., Kontar, E. P., & Browning, P. K. 2016, A&A, 589, A104
  * Hasegawa & Chen (1974) Hasegawa, A. & Chen, L. 1974, Phys. Rev. Lett., 32, 454
  * Heyvaerts & Priest (1983) Heyvaerts, J. & Priest, E. R. 1983, Astron. Astrophys., 117, 220
  * Hillier et al. (2019) Hillier, A., Barker, A., Arregui, I., & Latter, H. 2019, MNRAS, 482, 1143
  * Hood et al. (2016) Hood, A. W., Cargill, P. J., Browning, P. K., & Tam, K. V. 2016, ApJ, 817, 5
  * Howson et al. (2021) Howson, T., De Moortel, I., & Pontin, D. 2021, arXiv e-prints, arXiv:2109.15019
  * Howson et al. (2020) Howson, T. A., De Moortel, I., & Fyfe, L. E. 2020, A&A, 643, A85
  * Howson et al. (2019) Howson, T. A., De Moortel, I., Reid, J., & Hood, A. W. 2019, A&A, 629, A60
  * Hsieh et al. (2009) Hsieh, M. H., Tsai, C. L., Ma, Z. W., & Lee, L. C. 2009, Physics of Plasmas, 16, 092901
  * Ionson (1978) Ionson, J. A. 1978, Astrophys. J., 226, 650
  * Iroshnikov (1964) Iroshnikov, P. S. 1964, Sov. Ast., 7, 566
  * Johnston & Bradshaw (2019) Johnston, C. D. & Bradshaw, S. J. 2019, ApJ, 873, L22
  * Johnston et al. (2021) Johnston, C. D., Hood, A. W., De Moortel, I., Pagano, P., & Howson, T. A. 2021, arXiv e-prints, arXiv:2106.03989
  * Kigure et al. (2010) Kigure, H., Takahashi, K., Shibata, K., Yokoyama, T., & Nozawa, S. 2010, PASJ, 62, 993
  * Klimchuk (2006) Klimchuk, J. A. 2006, Sol. Phys., 234, 41
  * Klimchuk et al. (2008) Klimchuk, J. A., Patsourakos, S., & Cargill, P. J. 2008, ApJ, 682, 1351
  * Linker et al. (2001) Linker, J. A., Lionello, R., Mikić, Z., & Amari, T. 2001, J. Geophys. Res., 106, 25165
  * Lionello et al. (2009) Lionello, R., Linker, J. A., & Mikić, Z. 2009, ApJ, 690, 902
  * Matthaeus et al. (1999) Matthaeus, W. H., Zank, G. P., Oughton, S., Mullan, D. J., & Dmitruk, P. 1999, ApJ, 523, L93
  * McIntosh et al. (2011) McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477
  * McLaughlin et al. (2009) McLaughlin, J. A., De Moortel, I., Hood, A. W., & Brady, C. S. 2009, A&A, 493, 227
  * Mikić et al. (2013) Mikić, Z., Lionello, R., Mok, Y., Linker, J. A., & Winebarger, A. R. 2013, ApJ, 773, 94
  * Morton et al. (2015) Morton, R. J., Tomczyk, S., & Pinto, R. 2015, Nature Communications, 6, 7813
  * Morton et al. (2016) Morton, R. J., Tomczyk, S., & Pinto, R. F. 2016, ApJ, 828, 89
  * Morton et al. (2019) Morton, R. J., Weberg, M. J., & McLaughlin, J. A. 2019, Nature Astronomy, 196
  * Narain & Ulmschneider (1990) Narain, U. & Ulmschneider, P. 1990, Space Sci. Rev., 54, 377
  * Narain & Ulmschneider (1996) Narain, U. & Ulmschneider, P. 1996, Space Sci. Rev., 75, 453
  * Ofman et al. (1998) Ofman, L., Klimchuk, J. A., & Davila, J. M. 1998, ApJ, 493, 474
  * Parker (1972) Parker, E. 1972, The Astrophysical Journal, 174, 499
  * Parker (1991) Parker, E. N. 1991, ApJ, 376, 355
  * Parnell & De Moortel (2012) Parnell, C. E. & De Moortel, I. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 3217
  * Peter et al. (2004) Peter, H., Gudiksen, B. V., & Nordlund, Å. 2004, ApJ, 617, L85
  * Poedts et al. (1989) Poedts, S., Goossens, M., & Kerner, W. 1989, Solar physics, 123, 83
  * Prokopyszyn et al. (2019) Prokopyszyn, A. P. K., Hood, A. W., & De Moortel, I. 2019, A&A, 624, A90
  * Reale (2010) Reale, F. 2010, Living Reviews in Solar Physics, 7, 5
  * Reale (2014) Reale, F. 2014, Living Reviews in Solar Physics, 11, 4
  * Reid et al. (2018) Reid, J., Hood, A. W., Parnell, C. E., Browning, P. K., & Cargill, P. J. 2018, A&A, 615, A84
  * Rosner et al. (1978) Rosner, R., Tucker, W. H., & Vaiana, G. S. 1978, ApJ, 220, 643
  * Sakai et al. (1984) Sakai, J., Tajima, T., & Brunel, F. 1984, Solar physics, 91, 103
  * Sakurai (2017) Sakurai, T. 2017, Proceedings of the Japan Academy, Series B, 93, 87
  * Srivastava et al. (2017) Srivastava, A. K., Shetye, J., Murawski, K., et al. 2017, Scientific Reports, 7, 43147
  * Terradas et al. (2008) Terradas, J., Andries, J., Goossens, M., et al. 2008, ApJ, 687, L115
  * Thurgood et al. (2014) Thurgood, J. O., Morton, R. J., & McLaughlin, J. A. 2014, ApJ, 790, L2
  * Tomczyk et al. (2007) Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192
  * Van Ballegooijen (1988) Van Ballegooijen, A. 1988, Geophysical & Astrophysical Fluid Dynamics, 41, 181
  * van Ballegooijen et al. (2011) van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R., & DeLuca, E. E. 2011, ApJ, 736, 3
  * Van Damme et al. (2020) Van Damme, H. J., De Moortel, I., Pagano, P., & Johnston, C. D. 2020, A&A, 635, A174
  * Van Doorsselaere et al. (2020) Van Doorsselaere, T., Srivastava, A. K., Antolin, P., et al. 2020, Space Sci. Rev., 216, 140
  * Wilmot-Smith et al. (2011) Wilmot-Smith, A., Pontin, D., Yeates, A., & Hornig, G. 2011, Astronomy & Astrophysics, 536, A67
  * Wilmot-Smith (2015) Wilmot-Smith, A. L. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140265
  * Winebarger et al. (2018) Winebarger, A. R., Lionello, R., Downs, C., Mikić, Z., & Linker, J. 2018, ApJ, 865, 111