# Influence of the Lower Atmosphere on Wave Heating and Evaporation in Solar
Coronal Loops

Mingzhe Guo  Centre for mathematical Plasma Astrophysics, Department of
Mathematics, KU Leuven, Celestijnenlaan 200B, B-3001 Leuven, Belgium  Timothy
Duckenfield  Department of Mathematics, Physics and Electrical Engineering,
Northumbria University, Newcastle upon Tyne, NE1 8ST, UK  Tom Van Doorsselaere
Centre for mathematical Plasma Astrophysics, Department of Mathematics, KU
Leuven, Celestijnenlaan 200B, B-3001 Leuven, Belgium  Konstantinos Karampelas
Centre for mathematical Plasma Astrophysics, Department of Mathematics, KU
Leuven, Celestijnenlaan 200B, B-3001 Leuven, Belgium  Gabriel Pelouze
Université Paris-Saclay, CNRS, Institut d’astrophysique spatiale, 91405,
Orsay, France  Yuhang Gao  School of Earth and Space Sciences, Peking
University, Beijing, 100871, China  Centre for mathematical Plasma
Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B,
B-3001 Leuven, Belgium

###### Abstract

We model a coronal loop as a three-dimensional magnetic cylinder in a
realistic solar atmosphere that extends from the chromosphere to the corona.
Kink oscillations, believed ubiquitous in the solar corona, are launched in
the loop. Heating is expected due to the dissipation of wave energy at small
structures that develop from the Kelvin-Helmholtz instability induced by kink
oscillations. Increases in temperature and internal energy can be observed in
the coronal counterpart of the driven loop. With the presence of thermal
conduction, chromospheric evaporation can also be seen. Although the volume
averaged temperature and density changes seem slight ($\sim 4\%$ relative to a
non-driven loop), the enthalpy flow from the lower atmosphere redistributes
the density and temperature in the vertical direction, thus enhancing the
dissipation of wave energy in the corona. The efficient heating in the coronal
counterpart of the loop can complement the thermal conductive losses shown in
the current model and thus maintain the internal energy in the corona.

##  1 INTRODUCTION

More and more observations have confirmed the omnipresence of kink waves in
coronal loops (e.g., Nakariakov & Kolotkov, 2020; Nakariakov et al., 2021, for
recent reviews) and thus bring kink waves into the forefront of coronal
heating consideration (see De Moortel & Nakariakov, 2012; Arregui, 2015; Van
Doorsselaere et al., 2020, for recent reviews). In the heating process,
effective dissipation of wave energy can be achieved when the Lundquist number
and Reynolds number, which are typically very large ($\sim 10^{12}$) in the
solar corona, become very small (Soler & Terradas, 2015). Thus efficient
heating relies on the occurrence of small spatial scales. Generally, kink
waves are believed to suffer resonant absorption (see Goossens et al., 2011,
for a review) that transfers wave energy from collective kink modes to local
Alfvén waves with the occurrence of transverse inhomogeneity around the loop
boundary (e.g., Guo et al., 2020). The converted localized azimuthal Alfvén
waves undergo phase mixing (e.g., Heyvaerts & Priest, 1983). The subsequent
phase mixed Alfvén waves can enhance the Kelvin-Helmholtz instability (KHI,
Browning & Priest, 1984), which is induced by velocity shear between loops and
the background corona, thereby facilitating the wave dissipation by further
generating smaller scales. This scenario is first demonstrated by Ofman et al.
(1994) in coronal loop models. Recently, numerical progress confirmed that
transverse waves in coronal loops can induce the KHI, forming the TWIKH rolls
(Transverse Wave Induced Kelvin-Helmholtz rolls, see e.g., Terradas et al.,
2008; Antolin et al., 2015; Karampelas et al., 2017; Howson et al., 2017; Guo
et al., 2019a, b; Shi et al., 2021). Non-linear damping of kink waves
associated with the KHI has been investigated numerically by e.g., Magyar &
Van Doorsselaere (2016) and analytically by e.g., Van Doorsselaere et al.
(2021). During such a process, the energy of collective kink waves can
dissipate at such turbulent small structures, and quantitative evaluations of
heating effects have been made by e.g., Karampelas et al. (2017, 2019); Guo et
al. (2019a, b). Recent numerical progress shows that the heating effects
induced by kink waves can balance the radiative loss in the solar corona (Shi
et al., 2021; De Moortel & Howson, 2022), which shed new light on
understanding coronal heating from the perspective of waves.

The influence of chromospheric evaporation on the heating of coronal loops was
first considered by Ofman et al. (1998) based on resonant kink waves. However,
the aforementioned investigations associated with TWIKH rolls mainly focus on
the coronal counterpart of loops. In reality, the lower atmosphere can act as
a mass and energy reservoir for physical processes in the corona. In the
process associated with impulsive heating, nanoflares (Parker, 1988) for
instance, enthalpy flux plays an important role in the thermal evolution of
coronal loops (e.g., Bradshaw & Cargill, 2010; Cargill et al., 2012).
Regarding the AC models, Ohmic heating has been confirmed by both observations
(Van Doorsselaere et al., 2007) and numerical simulations (Karampelas et al.,
2019; Guo et al., 2019b) in coronal loops. In such a scenario, wave heating is
expected to happen near loop footpoints, where the current density that is
essential for resistive dissipation has maximum values at loop ends (Van
Doorsselaere et al., 2007). In this way, the cooler chromosphere can be heated
straightforwardly. Given this, a natural question arises: Can the heating near
loop footpoints trigger chromosphere evaporation? The scenario of evaporation
is usually discussed when chromospheric materials are heated during such
energy release processes as solar flares (e.g., Fisher et al., 1985; Milligan
et al., 2006; Tian & Chen, 2018). Although the energy dissipation of waves is
not so fierce as flares, gentle chromosphere evaporations are still expected
(Fisher et al., 1985). The influence of such evaporation flows on the
efficiency of wave heating in the corona remains to be evaluated in recent
models. Recent progress by Van Damme et al. (2020) investigated the influence
of chromosphere evaporation on the phase mixing of Alfvén waves in a two-
dimensional coronal loop model. The density shell around the loop boundary
that is essential for energy dissipation of Aflvén waves seems not
significantly changed in the modest heating process. However, wave energy
dissipation happens over almost the whole deformed cross-section of loops when
it comes to transverse oscillations (e.g., Karampelas & Van Doorsselaere,
2018; Guo et al., 2019a, b). It remains to be seen how chromosphere
evaporation influences wave energy dissipation in coronal loops deformed by
kink oscillations.

In this paper, we aim to investigate the influence of the lower atmosphere on
coronal heating effects. The paper is organized as follows. Section 2
describes the model, including the equilibrium and numerical setup. In Section
3, we present the numerical results. Section 4 summarizes our findings, ending
with some further discussion.

##  2 NUMERICAL MODEL

We use a similar magnetic flux tube employed in Pelouze et al. (2023, P23),
with the key difference that we simulate half of a closed loop with both
footpoints anchored in the chromosphere. The model is initiated from a 2D
hydrostatic equilibrium in cylindrical coordinates ($r,z$). A loop-oriented
gravity $g(z)=-g_{\odot}\sin\left[\pi(L-z)/(2L)\right]$ is considered, with
$L=100$Mm being the half length of the loop. The initial magnetic field is
$\bm{B}=B_{0}\hat{z}$, with $B_{0}$ set to 42G. The temperature profile,
derived from Aschwanden & Schrijver (2002), is given by

$\displaystyle T(r,z)=\left\\{\begin{array}[]{ll}T_{\rm ch},&z\leq z_{\rm
ch},\\\ T_{\rm ch}+\left[T_{\rm co}(r)-T_{\rm
ch}\right]\left[1-\left(\displaystyle\frac{L-z}{L-z_{\rm
ch}}\right)^{2}\right]^{0.3},&z>z_{\rm ch},\end{array}\right.$ (3)

where $T_{\rm ch}=20000$K is the temperature of the chromosphere, and $z_{\rm
ch}=4$Mm represents its thickness. The transverse temperature profile is
defined as

$\displaystyle T_{\rm co}(r)=T_{\rm e}+\displaystyle\frac{1}{2}\left(T_{\rm
i}-T_{\rm
e}\right)\left\\{1-\tanh\left[\left(\displaystyle\frac{r}{R}-1\right)b\right]\right\\},$
(4)

where $T_{\rm i}=1.2$MK ($T_{\rm e}=3.6$MK) is the temperature inside
(outside) the tube. $R=1$Mm gives the radius of the loop, and $b=10$ defines
the thickness of the inhomogeneous layer, corresponding to $l/R=0.6$ at the
loop bottom. The transverse distribution of the density in Figure 1a clearly
shows the variation of this layer thickness along the height. The initial
state is not in magnetohydrostatic balance, thus a relaxation is needed. In
this stage, we redefine the velocity rewrite layer, which was proposed in P23
to absorb the upward flows in the loop. With the presence of the velocity
rewrite layer, the large flow induced by the initial state, which is not in
magnetohydrostatic balance, is weakened. Thus the process of reaching an
equilibrium becomes more gradual and gentle, ensuring the stability of the
simulation. The velocity is modified as $v^{\prime}(r,z)=\alpha(t)v_{i}(r,z)$.
The parameter $\alpha(t)$ is defined as

$\displaystyle\alpha(t)=\left\\{\begin{array}[]{ll}0.9995+0.0005\displaystyle\frac{t}{t_{c}},&0<t\leq
t_{\rm c},\\\ 1,&t>t_{\rm c},\end{array}\right.$ (7)

where $t_{\rm c}=37.3{\rm ks}$ represents the critical time. The vertical
velocity in the loop is then suppressed to be less than $1.3{\rm kms^{-1}}$
and significant quantity changes (e.g., a drop in magnetic field strength in
P23 111In the current model, the magnetic field undergoes a slight change from
the initial value of 42G to around 40G after relaxation.) are avoided after
the relaxation. The relaxed state is then straightforwardly converted to 3D by
rotating the 2D axisymmetric results. To avoid any unphysical oscillations due
to boundary changes from the 2D cylindrical coordinate to the 3D Cartesian
domain, we allow the 3D system to relax for another 2.8ks. The initial
temperature and density distributions of the loop axis of the 2D and 3D
simulations are shown in Figure 2.

Note that this simulation is gravitationally stratified, leading to the
transverse structuring of the loop changing with height. As shown in Figure
1b, the inhomogeneous layer thickness is smaller at the loop apex ($z=100$Mm)
than at lower heights. One may expect this change in layer thickness to have
an impact on the formation of KHI (e.g., Magyar & Van Doorsselaere, 2016).

Figure 1: Transverse distribution of the density along the $r$-direction in
the 2D simulation (left) and $x$-direction in the 3D simulation (right). The
left panel represents the initial density in the 2D simulation. The right
panel shows the initial state of the 3D simulation after relaxation. Different
colors represent different heights. Figure 2: Spatial distribution of
temperature and density along the $z$-direction. The left panel displays the
initial quantities before relaxation in the 2D equilibrium. The right panel
shows the initial state of the 3D simulation after relaxation. Solid and
dashed lines represent the distribution at the loop axis ($x=0$, $y=0$) and
the external loop region ($x=8$Mm, $y=0$), respectively.

The chromosphere and transition region are included in the current model. The
transition region has been artificially broadened to ensure its resolution,
following the scheme developed by Linker et al. (2001); Lionello et al.
(2009); Mikić et al. (2013). As described in P23, a critical temperature
$T_{\rm c}=2.5\times 10^{5}{\rm K}$ is defined, below (above) which the
parallel thermal conductivity is $\kappa_{\parallel}=\kappa_{0}T_{\rm
c}^{5/2}$ ($\kappa_{\parallel}=\kappa_{0}T^{5/2}$), with $\kappa_{0}$ being
$5.6\times 10^{-7}{\rm erg~{}cm^{-1}~{}s^{-1}~{}K^{-7/2}}$. This modified
thermal conductivity changes the temperature length scale in the vertical
direction, thus leading to a significant broadening of the transition region
after relaxation, see also Van Damme et al. (2020) and P23. Therefore, a
coarser resolution in the vertical direction is allowed in the current
simulation.

To solve the 3D ideal MHD equations, we employ the PLUTO code (Mignone et al.,
2007). A piecewise linear scheme is used for spatial reconstruction. Numerical
fluxes are computed by the Roe Riemann solver, and a second-order
characteristic tracing method is used for time marching. The hyperbolic
divergence cleaning method is adopted to ensure the divergence-free nature of
the magnetic field. Anisotropic thermal conduction is included in our
simulations. In the 2D run, the computational domain is $[0,6]{\rm
Mm}\times[0,100]{\rm Mm}$. We consider 128 uniformly spaced cells in the
$r$-direction and a uniform grid of 1024 cell points in the $z$-direction. In
the 3D case, the simulation domain is $[-6,6]{\rm Mm}\times[-3,3]{\rm
Mm}\times[0,L/2]$. We consider uniform 1024 cell points in the $z$-direction
and 256 uniformly spaced cells in the $x$-direction and 128 uniform grid cells
in the $y$-direction. The resolution in the $x,y$ plane is about 46.9km.
Although increasing the transverse resolution may reveal more pronounced
heating effects, the main findings of the current work are not expected to be
influenced. On one hand, the changing width of the boundary layer allows the
current resolution to sufficiently resolve fine structures away from the loop
apex 222Even at the loop apex, we can still observe relatively larger KHI
structures, as shown in Figure 3b and related animations.. On the other hand,
the turbulent structures generated in the model will further extend the
thickness of the boundary layer, making it more feasible to reveal small
scales within it. The vertical resolution seems lower than that in the
transverse direction. Nevertheless, it is enough to resolve the broadened
transition region that is more than 2Mm wide after relaxation, considering the
broadening scheme described above.

In the 2D simulation, an axisymmetric (outflow) boundary is employed at $r=0$
($r=100$Mm). A symmetric boundary condition is used at $z=100$Mm, considering
the symmetric property of the fundamental kink mode considered in this study.
At the bottom of the loop, the density and pressure are extrapolated from the
hydrostatic equilibrium. The magnetic field is extrapolated following the zero
normal gradient condition, as described by Karampelas et al. (2019). The
vertical velocity $v_{z}$ is set to be zero, while the transverse velocity
$v_{r}$ is to be outflow. In the 3D case, all the lateral boundaries are set
to be outflow. Boundary conditions in the $z$-direction are kept the same as
in the 2D except for the transverse velocity $v_{x}$ and $v_{y}$ at $z=0$,
which are described by a continuous and monoperiodic driver (Pascoe et al.,
2010; Karampelas et al., 2017; Guo et al., 2019a). In the current study, the
driver period is 288s, which matches the eigenfrequency of the loop at the
initial state. Note that the initial eigenfrequency of the loop slightly
varies over time during the driving phase. Due to the modest changes in the
density structure in the current model, the variation in eigenfrequency is not
significant. Nonetheless, a broad-band driver (e.g., Afanasyev et al., 2019)
needs further study and discussion. The amplitude of the driver is chosen to
be 0 (non-driven model) and $4{\rm kms^{-1}}$ (driven model). These two models
are mentioned hereafter by using subscripts “0” and “4”, respectively.

Figure 3: Isosurfaces of the logarithm of density and temperature at the loop
apex ($z=100$Mm) and $y=0$ for (a) non-driven and (b) driven model at
$t=8060$s. An animation of the density and temperature isosurfaces for the
driven model is available in the online Journal. The animation proceeds from
$t=0\text{ to }9340$s.

##  3 RESULTS

Figure 4: Surface-averaged temperature, density, enthalpy flux, and mass flux
changes relative to the non-driven case are shown versus height, with zoomed-
in lower panels for the regions outlined by dashed boxes.

We first examine the dynamic process of the loop. Figure 3 shows the
isocontours of density and temperature for both models at $t=8060$s when the
loop has fully evolved. The TWIKH rolls are clearly seen from the density and
temperature distribution at the apex of the driven loop in Figure 3b. The core
of the loop cross-section is sufficiently heated at the loop apex, even though
the Ohmic dissipation is smaller than that near the loop footpoints. To
further examine the evolution in the vertical direction, surface averaged
temperature and density are examined along the $z$-direction, which are given
by

$\displaystyle\bar{T}(z,t)=\frac{1}{A}\int_{A}T{\rm
d}A,\bar{\rho}(z,t)=\frac{1}{A}\int_{A}\rho{\rm d}A.$ (8)

where $A$ represents the surface area of the region $-6{\rm Mm}\leq x\leq
6{\rm Mm},-3{\rm Mm}\leq y\leq 3{\rm Mm}$ at given heights. ${\rm d}A$ is the
surface element of the plane. Figure 4a,b show the changes in
$\left(\bar{T}-\bar{T}_{0}\right)/\bar{T}_{0}$ and
$\left(\bar{\rho}-\bar{\rho}_{0}\right)/\bar{\rho}_{0}$, which clearly reveal
the variations caused by waves in comparison to the non-driven case. To reveal
the evaporation from the lower atmosphere, the enthalpy flux and mass flux are
also considered as

$\displaystyle\bar{E}(z,t)=\frac{1}{A}\int^{t}_{0}\int_{A}\left(\frac{\gamma}{\gamma-1}p\right)v_{z}{\rm
d}A{\rm d}t,$ (9)
$\displaystyle\bar{M}(z,t)=\frac{1}{A}\int^{t}_{0}\int_{A}\rho v_{z}{\rm
d}A{\rm d}t.$ (10)

Figure 4c shows the changes of the $z$-component of the enthalpy flux
($\Delta\bar{E}=\bar{E}-\bar{E}_{0}$). From Figure 4a, we observe the expected
temperature increases near the loop footpoint. This heating induces upward
enthalpy flows, as seen in Figure 4c. The upward enthalpy flows are present
throughout the simulation, although periodic downward flows occur after about
50min. Figure 4b confirms the mass changes induced by the upward flow, as
evidenced by the density increases. A periodicity of $\sim 20$min 333This
period is close to the slow transit time in the coronal counterpart. of the
temperature profile is observed. The density and enthalpy flux exhibit a
similar periodicity as shown in Figure 4. This periodicity is associated with
the ponderomotive force and is attributed to slow waves, which has also been
observed in previous numerical studies (e.g., Magyar & Van Doorsselaere, 2016;
Karampelas et al., 2017). This has also been confirmed by Van Damme et al.
(2020) by considering an Alfvén wave driver. The excited slow oscillation can
influence the vertical distribution of the temperature profile, density
distribution, and enthalpy flux. The variations in the temperature and density
profiles relative to the non-driven case seem slight (less than $4\%$). This
is partly because the evaporation induced by waves, as demonstrated in Van
Damme et al. (2020), can only cause modest variations in the dramatically
changing atmosphere background in the vertical direction. In addition, the
surface averaged values here are naturally smaller than the local quantity
values in the loop region, which is no longer straightforward to quantify
since our current loop is fully deformed and transversely diffused than
previous coronal models. Nevertheless, the surface average procedure can still
reveal the average variations of quantities of loop cross-section at a given
height. The evolutions of quantities in the domain of $2{\rm Mm}\leq z\leq
8{\rm Mm}$ are zoomed in and presented in the lower panels of Figure 4. We see
upward enthalpy flows in this domain, indicating the occurrence of
chromospheric evaporation. The density drops and temperature increases around
$z=4$Mm are the consequence of this evaporation. Additionally, the lower
region presents a decrease (increase) in temperature (density) variations
around $z=2$Mm. Seen from Figure 2, the boundary between the chromosphere and
the transition region shifts downward to about $z=1$Mm, while the lower
boundary of the transition region outside the loop is around $z=2$Mm. This
means that denser and cooler masses from the external loop region are involved
when computing the surface averaged quantities around $z=2$Mm.

Heating effects induced by kink waves can be quantified by energy analysis
(e.g., Karampelas et al., 2017; Guo et al., 2019a, b). The continuous driving
at the loop footpoint changes the equilibrium compared with the non-driven
case, leading to the transition region moving with time. We thus consider an
upper region starting from $9$Mm as the coronal counterpart in the following
analysis. The volume-averaged Poynting flux and internal energy are given by

$\displaystyle
S(t)=-\displaystyle\frac{1}{V}\int^{t}_{0}\int_{A^{\prime}_{1}}\displaystyle\frac{1}{\mu_{0}}\left[(\bm{v}\times\bm{B})\times\bm{B}\right]\cdot{\rm
d}\bm{A^{\prime}}_{1}{\rm d}t,$ (11) $\displaystyle
I(t)=\displaystyle\frac{1}{V}\int_{V}\displaystyle\frac{p}{\gamma-1}{\rm
d}V-\displaystyle\frac{1}{V}\int^{t}_{0}\int_{A^{\prime}_{2}}\left(\displaystyle\frac{\rho
v^{2}}{2}+\rho\Phi+\displaystyle\frac{\gamma
p}{\gamma-1}\right)\bm{v}\cdot{\rm d}\bm{A^{\prime}}_{2}{\rm d}t,$ (12)

where $V$ represents the volume of the concerned domain ($-6{\rm Mm}\leq x\leq
6{\rm Mm},-3{\rm Mm}\leq y\leq 3{\rm Mm},9{\rm Mm}\leq z\leq 100{\rm Mm}$).
${\rm d}\bm{A}^{\prime}_{1}$ (${\rm d}\bm{A}^{\prime}_{2}$) is the normal
vector of the bottom (lateral) surface. We also incorporate the energy fluxes
of lateral boundaries into the internal energy variations, following a similar
procedure in Karampelas et al. (2019). The input energy is defined as the
Poynting flux of the bottom interface. As aforementioned, an enthalpy flux can
be seen through the bottom boundary ($z=9$Mm) of the concerned domain.
However, this enthalpy flux is much smaller than the Poynting flux here, thus
has been neglected when considering the input energy from the bottom. Figure 5
shows the volume averaged energy density changes relative to the initial
state, namely $S(t)-S(0)$ and $I(t)-I(0)$. In the non-driven model, the
internal energy ($I_{0}$) drops with time due to the thermal conductive loss
in the corona, while the internal energy variation in our driven model
($I_{4}$) stays positive before $120$min and approaches zero after about
$120$min. This means that the internal energy first increases and then
asymptotically approaches the initial equilibrium state by the end of the
simulation. Therefore, the heating induced by waves can overcome or balance
the thermal conductive loss in the corona, for the duration of the driving in
the current model.

Figure 5: Volume averaged input energy (black) and internal energy density
(red) changes relative to the initial state in the coronal counterpart
($-6{\rm Mm}\leq x\leq 6{\rm Mm},-3{\rm Mm}\leq y\leq 3{\rm Mm},9{\rm Mm}\leq
z\leq 100{\rm Mm}$). Subscript “0” (“4”) denotes the non-driven (driven)
model. Dot-dashed line shows the difference between the two models in the
concerned domain.

Heating efficiency can be evaluated by considering the ratio between the
relative internal energy increase and the input energy. Figure 5 also shows
the difference between the driven and non-driven models in the same domain.
From the Poynting flux and the difference in internal energy changes, we can
compute the heating efficiency of the driven model. A rapid increase in
Poynting flux before 40min can be observed, leading to an average input energy
flux of about $78.1{\rm Jm^{-2}s^{-1}}$ over the considered coronal
counterpart of the loop. The average internal energy increasing rate is about
$32.8{\rm Jm^{-2}s^{-1}}$, leading to a heating rate of $42\%$. In previous
studies (e.g., Guo et al., 2019a; Karampelas et al., 2019), we report the
expansion of the loop cross-section when small scales and turbulence
developed, inducing a decrease in the local magnetic field. This decrease can
also be traced down from the corona to lower heights with decreasing strength.
The variations of the magnetic field in the lower interface between the corona
and the transition region cause a slower increase of the input energy from
$t\sim 50$min, leading to an average input energy flux of about $30.7{\rm
Jm^{-2}s^{-1}}$ from $50{\rm min}\leq t\leq 100{\rm min}$. The corresponding
internal energy growth rate is $16.8{\rm Jm^{-2}s^{-1}}$ and the heating rate
increases to $55\%$. In the last 20 minutes of the simulation, the loop is
fully evolved and the expansion is saturated, we see the input energy flux
increase to $45.1{\rm Jm^{-2}s^{-1}}$. The corresponding internal energy
increasing rate of this phase is $34.5{\rm Jm^{-2}s^{-1}}$, leading to a
heating rate of $76.5\%$. This shows a significant heating efficiency of the
current model. Comparing with our previous models (e.g., Karampelas et al.,
2017, 2019; Guo et al., 2019a, b; Shi et al., 2021), the presence of thermal
conduction allows enthalpy flows from the lower atmosphere to the corona. Such
flows change the vertical density and temperature distribution in the coronal
counterpart of the loop (Figure 4). The energy dissipation at the loop cross-
section becomes more efficient, as indicated by the temperature increase in
Figure 3. Therefore, the upward evaporation flows lead to a more efficient
heating to balance the thermal conductive losses in the current model.

##  4  DISCUSSION AND SUMMARY

We model a three-dimensional magnetic cylinder in a realistic solar atmosphere
from the chromosphere to the corona with thermal conduction included.
Frequently observed kink oscillations are excited in this magnetic flux tube.
Based on previous studies, heating effects are expected due to the presence of
the Kelvin-Helmholtz instability induced by kink oscillations. Temperature and
internal energy increases are indeed observed in the coronal counterpart of
our driven model, with respective to a non-driven case. With the inclusion of
thermal conduction, chromospheric evaporation can be observed, leading to an
increase in temperature and a decrease in density. Although the spatially
averaged density and temperature changes are gentle ($\sim 4\%$ relative to
the non-driven case), the corresponding enthalpy flow from the lower
atmosphere can still influence the energy dissipation in the corona by
redistributing the vertical density structuring. The efficient heating in the
coronal counterpart of the loop can balance the thermal conductive loss in the
current model, and thus maintain the internal energy in the corona.

In the current simulation, we adopt an approximate model derived from
Aschwanden & Schrijver (2002) to describe the vertical distribution of the
solar atmosphere. This temperature distribution seems a coarser approximation
than the VAL model (Vernazza et al., 1981) or its upgraded versions (Avrett &
Loeser, 2008). Especially, it lacks sufficient details in the lower
atmosphere. Nevertheless, it describes the main properties of the parameters
(e.g., temperature and density) of the chromosphere and the transition region.
In addition, note that the initial analytical profiles of temperature and
density change after relaxation and a new equilibrium is achieved, as shown in
Figure 2b. The initial expressions of the VAL model of the atmosphere would
also be changed after the relaxation.

Explicit resistivity and viscosity are not included in the current ideal MHD
computation. In the solar corona, the magnetic Reynolds number is of the order
of $10^{12}$, leading to an extremely small resistivity if considering typical
spatial scales and Alfvén speed in the corona. In practical coronal models,
however, even the numerical resistivity is several orders larger than the
realistic value. Given the numerical resolution limitation, we should treat
such models, especially large 3D models, as an ideal approximation to the real
corona. Nevertheless, to examine the influence of explicit
resistivity/viscosity on wave energy dissipation is still necessary. Howson et
al. (2017) found that enhanced resistivity and viscosity that are larger than
numerical values can suppress the onset of TWIKH rolls, thus reducing wave
energy dissipation. The influence of explicit dissipation coefficients on wave
heating has been quantitatively examined by Karampelas et al. (2019).
Comparing with ideal simulations with only numerical dissipation, the heating
effects are enhanced but the heating locations are the same in the
resistive/viscous models in Karampelas et al. (2019). This confirms that the
numerical dissipation can also play an effective role in such modeling
associated with waves heating effects.

In previous studies (e.g., Ofman et al., 1998), resonant absorption and phase
mixing are believed to happen in resonant layers near a loop boundary. Cargill
et al. (2016) believed that the resonant layer is not sustained due to the
temporal evolution of the density gradient induced by Alfvén waves. In our
current results, however, TWIKH rolls spread over the whole cross section of
the loop, associated with chromospheric evaporation not only happening near
the boundary layer, but rather in the entire loop. Nevertheless, due to the
slight density and temperature changes, the evaporation has no significant
effect on the resonance conditions of the current model.

In this simulation we are continuously driving waves at the eigenfrequency of
the standing mode, and so the resultant evaporations may be viewed as an upper
limit to the evaporations which would result from realistic, small amplitude
(decayless) kink waves which are not so carefully driven. The limitation of
current-day 3D MHD simulations having a magnetic Reynolds number several
orders of magnitude off that for the corona further reinforces our results as
an upper limit.

The realistic thermal conductive loss in the quiet solar corona is about
$200{\rm J~{}m^{-2}~{}s^{-1}}$ (Withbroe & Noyes, 1977). The input energy in
our current model does not seem to be able to complement such a high
conductive loss in the real corona. However, it should be noted that the
magnetic field strength of the current model is still vertically uniform. In
reality, if the loop footpoint is anchored in the photosphere, the magnetic
field should be much (or even several orders of magnitude) larger than that in
the corona. This probably indicates that the real energy input from the
photosphere is sufficient. In fact, observations have confirmed that the
energy flux is indeed sufficient to balance the energy losses even in the
active region (e.g., Fujimura & Tsuneta, 2009).

We have simulated a more realistic coronal loop by including a dense and cold
chromosphere, and despite the additional dynamics and energetics, it is
encouraging to see that the results about wave heating inferred from coronal-
only simulations still seems to hold. Examples include the formation of
extended TWIKH rolls along the whole loop; and the increased efficiency of
wave heating once the kinetic energy saturates at the later stages of the
simulation and a turbulent density profile has developed.

We thank the referee for helpful comments that improved the manuscript. The
authors acknowledge the funding from the European Research Council (ERC) under
the European Unions Horizon 2020 research and innovation program (grant
agreement No. 724326). TVD was also supported by the C1 grant TRACEspace of
Internal Funds KU Leuven, and a Senior Research Project (G088021N) of the FWO
Vlaanderen. K.K. acknowledges support by an FWO (Fonds voor Wetenschappelijk
Onderzoek-Vlaanderen) postdoctoral fellowship (1273221N). Y.G. acknowledges
the support from the China Scholarship Council (CSC) under file No.
202206010018.

## References

  * Afanasyev et al. (2019) Afanasyev, A., Karampelas, K., & Van Doorsselaere, T. 2019, ApJ, 876, 100. doi:10.3847/1538-4357/ab1848
  * Antolin et al. (2015) Antolin, P., Okamoto, T. J., De Pontieu, B., et al. 2015, ApJ, 809, 72. doi:10.1088/0004-637X/809/1/72
  * Arregui (2015) Arregui, I. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140261. doi:10.1098/rsta.2014.026110.48550/arXiv.1501.06708
  * Aschwanden & Schrijver (2002) Aschwanden, M. J. & Schrijver, C. J. 2002, ApJS, 142, 269. doi:10.1086/341945
  * Avrett & Loeser (2008) Avrett, E. H. & Loeser, R. 2008, ApJS, 175, 229. doi:10.1086/523671
  * Bradshaw & Cargill (2010) Bradshaw, S. J. & Cargill, P. J. 2010, ApJ, 710, L39. doi:10.1088/2041-8205/710/1/L39
  * Browning & Priest (1984) Browning, P. K. & Priest, E. R. 1984, A&A, 131, 283
  * Cargill et al. (2012) Cargill, P. J., Bradshaw, S. J., & Klimchuk, J. A. 2012, ApJ, 752, 161. doi:10.1088/0004-637X/752/2/161
  * Cargill et al. (2016) Cargill, P. J., De Moortel, I., & Kiddie, G. 2016, ApJ, 823, 31. doi:10.3847/0004-637X/823/1/31
  * De Moortel & Nakariakov (2012) De Moortel, I. & Nakariakov, V. M. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 3193. doi:10.1098/rsta.2011.064010.48550/arXiv.1202.1944
  * De Moortel & Howson (2022) De Moortel, I. & Howson, T. A. 2022, ApJ, 941, 85.doi:10.3847/1538-4357/aca072
  * Fisher et al. (1985) Fisher, G. H., Canfield, R. C., & McClymont, A. N. 1985, ApJ, 289, 414. doi:10.1086/162901
  * Fisher et al. (1985) Fisher, G. H., Canfield, R. C., & McClymont, A. N. 1985, ApJ, 289, 425. doi:10.1086/162902
  * Fujimura & Tsuneta (2009) Fujimura, D. & Tsuneta, S. 2009, ApJ, 702, 1443. doi:10.1088/0004-637X/702/2/1443
  * Goossens et al. (2011) Goossens, M., Erdélyi, R., & Ruderman, M. S. 2011, Space Sci. Rev., 158, 289. doi:10.1007/s11214-010-9702-7
  * Guo et al. (2019a) Guo, M., Van Doorsselaere, T., Karampelas, K., et al. 2019, ApJ, 870, 55. doi:10.3847/1538-4357/aaf1d0
  * Guo et al. (2019b) Guo, M., Van Doorsselaere, T., Karampelas, K., et al. 2019, ApJ, 883, 20. doi:10.3847/1538-4357/ab338e
  * Guo et al. (2020) Guo, M., Li, B., & Van Doorsselaere, T. 2020, ApJ, 904, 116. doi:10.3847/1538-4357/abc1df
  * Heyvaerts & Priest (1983) Heyvaerts, J. & Priest, E. R. 1983, A&A, 117, 220
  * Howson et al. (2017) Howson, T. A., De Moortel, I., & Antolin, P. 2017, A&A, 602, A74. doi:10.1051/0004-6361/201630259
  * Karampelas et al. (2017) Karampelas, K., Van Doorsselaere, T., & Antolin, P. 2017, A&A, 604, A130. doi:10.1051/0004-6361/201730598
  * Karampelas & Van Doorsselaere (2018) Karampelas, K. & Van Doorsselaere, T. 2018, A&A, 610, L9. doi:10.1051/0004-6361/201731646
  * Karampelas et al. (2019) Karampelas, K., Van Doorsselaere, T., & Guo, M. 2019, A&A, 623, A53. doi:10.1051/0004-6361/201834309
  * Karampelas et al. (2019) Karampelas, K., Van Doorsselaere, T., Pascoe, D. J., et al. 2019, Frontiers in Astronomy and Space Sciences, 6, 38. doi:10.3389/fspas.2019.00038
  * Linker et al. (2001) Linker, J. A., Lionello, R., Mikić, Z., et al. 2001, J. Geophys. Res., 106, 25165. doi:10.1029/2000JA004020
  * Lionello et al. (2009) Lionello, R., Linker, J. A., & Mikić, Z. 2009, ApJ, 690, 902. doi:10.1088/0004-637X/690/1/902
  * Magyar & Van Doorsselaere (2016) Magyar, N. & Van Doorsselaere, T. 2016, A&A, 595, A81.doi:10.1051/0004-6361/20162901010.48550/arXiv.1609.06883
  * Mikić et al. (2013) Mikić, Z., Lionello, R., Mok, Y., et al. 2013, ApJ, 773, 94. doi:10.1088/0004-637X/773/2/94
  * Milligan et al. (2006) Milligan, R. O., Gallagher, P. T., Mathioudakis, M., et al. 2006, ApJ, 638, L117. doi:10.1086/50055510.48550/arXiv.astro-ph/0509664
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228. doi:10.1086/513316
  * Nakariakov & Kolotkov (2020) Nakariakov, V. M. & Kolotkov, D. Y. 2020, ARA&A, 58, 441. doi:10.1146/annurev-astro-032320-042940
  * Nakariakov et al. (2021) Nakariakov, V. M., Anfinogentov, S. A., Antolin, P., et al. 2021, Space Sci. Rev., 217, 73. doi:10.1007/s11214-021-00847-2
  * Ofman et al. (1994) Ofman, L., Davila, J. M., & Steinolfson, R. S. 1994, Geophys. Res. Lett., 21, 2259. doi:10.1029/94GL01416
  * Ofman et al. (1998) Ofman, L., Klimchuk, J. A., & Davila, J. M. 1998, ApJ, 493, 474. doi:10.1086/305109
  * Parker (1988) Parker, E. N. 1988, ApJ, 330, 474. doi:10.1086/166485
  * Pascoe et al. (2010) Pascoe, D. J., Wright, A. N., & De Moortel, I. 2010, ApJ, 711, 990. doi:10.1088/0004-637X/711/2/990
  * Pelouze et al. (2023) Pelouze, G., Van Doorsselaere, T., Karampelas, K., et al. 2023, arXiv:2301.03100
  * Shi et al. (2021) Shi, M., Van Doorsselaere, T., Guo, M., et al. 2021, ApJ, 908, 233. doi:10.3847/1538-4357/abda54
  * Soler & Terradas (2015) Soler, R. & Terradas, J. 2015, ApJ, 803, 43. doi:10.1088/0004-637X/803/1/43
  * Terradas et al. (2008) Terradas, J., Andries, J., Goossens, M., et al. 2008, ApJ, 687, L115. doi:10.1086/593203
  * Tian & Chen (2018) Tian, H. & Chen, N.-H. 2018, ApJ, 856, 34. doi:10.3847/1538-4357/aab15a
  * Van Damme et al. (2020) Van Damme, H. J., De Moortel, I., Pagano, P., et al. 2020, A&A, 635, A174. doi:10.1051/0004-6361/201937266
  * Van Doorsselaere et al. (2007) Van Doorsselaere, T., Andries, J., & Poedts, S. 2007, A&A, 471, 311. doi:10.1051/0004-6361:20066658
  * Van Doorsselaere et al. (2020) Van Doorsselaere, T., Srivastava, A. K., Antolin, P., et al. 2020, Space Sci. Rev., 216, 140. doi:10.1007/s11214-020-00770-y
  * Van Doorsselaere et al. (2021) Van Doorsselaere, T., Goossens, M., Magyar, N., et al. 2021, ApJ, 910, 58. doi:10.3847/1538-4357/abe630
  * Vernazza et al. (1981) Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635. doi:10.1086/190731
  * Withbroe & Noyes (1977) Withbroe, G. L. & Noyes, R. W. 1977, ARA&A, 15, 363. doi:10.1146/annurev.aa.15.090177.002051