© ESO 2020

1. Introduction

Whether the energy carried by magnetohydrodynamic (MHD) waves in the solar corona can contribute to maintaining the million-degree plasma against its radiative losses remains a puzzle (e.g. Parnell & De Moortel 2012; Arregui 2015). Magnetohydrodynamic waves have been detected in the solar corona for more than a decade (e.g. Tomczyk et al. 2007; McIntosh et al. 2011; Mathioudakis et al. 2013) and are reported to carry a significant amount of energy, of the order of 50−200 W m⁻² (e.g. McIntosh et al. 2011; Thurgood et al. 2014; Morton et al. 2015) in comparison with the energy requirements to balance the continuous radiative losses of the coronal plasma. However, we still have no observational confirmation, nor a full modelling corroboration that this wave energy can be efficiently and effectively converted into thermal energy. These investigations are primarily carried out observing and modelling coronal loops that are arch-like dense and hot magnetic structures in the solar corona (Reale 2010). Although coronal loops do not fill up the entire coronal volume, as they are denser than the surrounding corona, they are responsible for most of the coronal emission and thus most of its radiative losses.

Moreover, recent observations of the solar corona have further constrained the properties of propagating transverse waves. In a series of work, the power spectrum of these oscillations has been derived and confirmed to be only marginally dependent on time and location, confirming that waves are an inherent and ubiquitous property of the solar corona (Morton et al. 2016, 2019). Although several mechanisms have been proposed to convert wave energy into thermal energy, realistic models of these mechanisms in coronal loops have so far not shown that they can indeed support the thermal structure of the corona. An example of such models is the mode coupling and subsequent phase mixing of Alfvén waves (Heyvaerts & Priest 1983). Pascoe et al. (2010) showed that the mode-coupling of kink and Alfvén azimuthal modes, and their subsequent phase mixing, can lead to the concentration of energy into thin boundary shells, across which the Alfvén speed varies. In this model, the boundary shell is assumed to be pre-existent to the wave propagation. Subsequent works (Pascoe et al. 2011, 2013) confirmed the robustness of these modelling results, showing that such an energy concentration naturally occurs when a boundary shell is present. Although this mode-coupling is effective in concentrating the wave energy, when the dissipation of waves is included in the model, the heating power derived from the conversion of the wave energy into thermal energy is largely insufficient to balance the radiative losses. Indeed, Pagano & De Moortel (2017) and Pagano et al. (2018) showed that there is insufficient heating of the coronal plasma to counteract the plasma radiative losses in a simplified model. The same conclusion was reached in the work of Pagano & De Moortel (2019) when propagating MHD waves were excited in the domain through the buffeting of the coronal loop footpoint with a driver derived from the multi-frequency spectrum observed by Morton et al. (2016). Other modelling efforts have reached similar conclusions (e.g. Karampelas et al. 2017, 2019), and also addressed how the mix of cold and hot plasma can show apparent heating. In contrast, other studies claim that MHD waves can efficiently heat the coronal plasma (Srivastava et al. 2017; López Ariste & Facchin 2018). These models either rely on very high-frequency oscillations, or they only focus on the energy stored in the waves, without addressing the conversion mechanism. However, analysing SDO/AIA observations, Morton et al. (2019) found that although high-frequency oscillations have larger amplitudes, they do not appear to occur regularly in the corona. Hence, their time-averaged power is small, implying they probably make little contribution to the energy supply.

At the same time, another line of research based on both modelling and observations has focused on the effect of propagating or standing transverse waves on the pre-existing coronal loop structure. Several papers have shown that standing waves cause Kelvin-Helmholtz instabilities (KHIs) at the boundary shell when a density contrast is present between the interior and the exterior region of a waveguide, leading to fragmentation of the boundary shell (e.g. Browning & Priest 1984; Terradas et al. 2008; Antolin et al. 2015; Okamoto et al. 2015). These results have been shown to be quite robust as KHIs develop even when resistivity and viscosity are included in the model (Howson et al. 2017), when standing modes are setup by the reflection of propagating waves (Karampelas et al. 2017), when gravity stratification is included (Antolin et al. 2018; Karampelas et al. 2019), and KHIs can occur in conjunction with parametric instabilities (Hillier et al. 2019). Similar structures were found in the simulations of Pagano & De Moortel (2017, 2019) with purely propagating waves, and Magyar et al. (2017) provided a first step to describing the dynamics in this setup, referred to as “uniturbulence”. Either way, the presence of MHD waves in dense loops seems to lead to the development of small-scale structures and the erosion of the density enhancement of the loop. Finally, these small-scale structures across which the density and magnetic field intensity vary rapidly could be favourable places to enhance phase mixing of waves, or where the magnetic field is sufficiently entangled (e.g. Reale et al. 2016; Pontin et al. 2017) such that nanoflares (Parker 1988) can occur.

In summary, so far, 3D numerical models seem to conclude that the direct dissipation of MHD waves is not able to supply enough energy to counteract the losses through radiation and thermal conduction in active-region coronal loops. However, waves could play a key role in shaping coronal loops and developing the small-scale structures required for their own dissipation or other heating mechanisms to become more efficient. The present paper builds on the work of Pagano & De Moortel (2019). We expanded our investigation of the heating induced by propagating transverse waves by considering a number of configurations for the loop density structuring, open or closed footpoints, and the extent of the driver. To do so, we drove an observed spectrum of Alfvén waves in the corona into a loop structure and we investigated the subsequent heating and boundary-shell evolution. We investigated how various parameters can affect (i) the heating deposited from the dissipation of waves, and (ii) how the boundary shell evolves and whether it becomes more effective. In particular, we focussed on loop structures, which are either homogeneous or change along the field-aligned direction, to model pre-existing loops that have been completely filled by plasma, or loops where the filling through evaporation is still ongoing (Reale et al. 2019), or which are short-lived structures such as spicules. We note here, though, that the process of evaporation or the motion of spicules are not included self-consistently. Conditions at the top boundary were modified to model both open and closed structures, and the horizontal extent of the boundary driver was altered to examine how the spatial extent of the footpoint motions compared to the extent of the flux tube affects the energy budget in the loop. Finally, we considered different power spectra in addition to the observed spectrum to show how model results are sensitive to the power spectrum used.

The paper is structured as follows. In Sect. 2, we illustrate the loop models and the driver setup, and in Sect. 3, we show the result for one representative simulation. In Sect. 4, we address the problem of the boundary-shell evolution, followed by Sect. 5, where we explain how different configurations impact the heating from waves. Finally, in Sect. 6, we show how the wave heating depends on the input power spectrum. Conclusions are presented in Sect. 7.

2. Model and driver

In order to study how the density structure of the loop affects the heating by transverse MHD waves in the solar corona, we ran a number of 3D numerical simulations where we solved the non-ideal MHD equations using the AMRVAC code (Porth et al. 2014). In these experiments, we varied (i) the prescribed initial condition for the loop structure, (ii) the boundary conditions that govern how transverse MHD waves are introduced or retained in the computational domain, (iii) the dissipation coefficients, and (iv) the driver we use to generate these waves. In this section, we describe in detail how these different experiments are constructed.

2.1. Loop structure

We treated the coronal loop as a magnetised cylinder consisting of a dense interior region embedded in a less dense environment. The loop interior is surrounded by a boundary shell across which the density decreases gradually until it matches the exterior density. The Alfvén speed is therefore uniform in the interior region and the exterior of the loop, but varies through the boundary shell. We neglect the effects of gravity and curvature, and focus on the coronal segment of the loop, without modelling the chromosphere and transition region.

We use a Cartesian reference frame, where z is the direction along the cylinder axis. The cylinder configuration we use is defined between z_min = −20 Mm and z_max = 180 Mm. The centre of the loop footpoint is placed at the origin of the axes. The cylinder has a radius a, and the radius of the dense interior region is given by b. Hence, the boundary shell is the interface region between the interior and the exterior, i.e. between radii b and a.

In order to study the effect of the initial coronal loop structure, we considered two different density distributions. In the uniform one, the interior region radius, b, is uniform along the length of the cylinder, and we took a = 1 Mm and b = 0.5 Mm. The density, ρ, decreases across the boundary shell and is defined as the following function of ρ_e, ρ_c, a, and b:

$$\begin{array}{c}\hfill {\displaystyle \rho ({\rho }_{\mathrm{e}},{\rho }_{\mathrm{c}},a,b)={\rho }_{\mathrm{e}}[1+\left(\frac{{\rho }_{\mathrm{c}}-1}{2}\right)[1-tanh\left(\frac{e}{a-b}[r-\frac{b+a}{2}]\right)]],}\end{array}$$(1)

where $r=\sqrt{x^2+{\mathit{y}}^2}$ is the radial distance from the centre of the cylinder, ρ_e = 1.16 × 10⁻¹⁶ g cm⁻³ is the external density, and ρ_c is the density enhancement between the exterior and interior region, which we set at ρ_c = 4. Figure 1a shows a density isosurface for this uniform configuration.

Fig. 1. Isosurfaces of density (top row) and Alfvén speed gradient (bottom row) for three different loop configurations: uniform (left-hand side), z0 = 130 Mm (centre), and z0 = 46 Mm (right-hand side).

In addition, we also considered a non-uniform cylinder where the interior region is only present near the lower boundary (for z ≤ z₀), defining its radius as

$$\begin{array}{c}\hfill b(z\le z_0)=0.5{(1-\frac{z-z_{\min }}{z_0-z_{\min }})}^4,\end{array}$$(2)

and b = 0 for z >  z₀. Using this profile of b(z ≤ z₀), the interior region radius rapidly shrinks above z = z_min, and it disappears at z = z₀ (red curves in Fig. 2). Additionally, at z = z₀ the density contrast quickly drops as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{c}}=1+\frac{3}{2}[tanh\left(\frac{z_0-z}{\mathrm{\Delta }z}\right)+1].\end{array}$$(3)

Fig. 2. Profile of density contrast ρc (blue curves) and interior region radius b (red curves), as a function of z for two density profiles of the loop with z0 = 46 Mm (dashed curves) and z0 = 130 Mm (continuous curves).

This profile for the density contrast is close to ρ_c = 4 for z <  z₀ and then decreases to ρ_c = 1 over a length Δz near z = z₀ (blue curves in Fig. 2). Figure 1b shows the density contour for this non-uniform configuration, where we use z₀ = 130 Mm and Δz = 5 Mm and Fig. 1c for z₀ = 46 Mm. The density structure narrows above z_min because the interior region gets thinner and the loop structure disappears above z = z₀, where the density contrast becomes ρ_c = 1. It should be noted that the density distribution at the cross-section z = −20 Mm is the same for both the uniform and non-uniform setup.

Such non-uniform configurations resemble the scenario where the density enhancement of the coronal loop occurs as a consequence of evaporation, which lifts plasma from the dense chromosphere into the corona. We assume that this local density enhancement is present on a time scale comparable to the wave propagation along the loop (see e.g. Antolin et al. 2018 for a similar model of a spicule).

The magnetic field, B, in the cylinder is initially uniform and aligned with the z direction. With an initial magnetic field strength of B₀ = 5.68 G, we have a plasma β = 0.02. The initial plasma temperature T is set by the equation of state:

$$\begin{array}{c}\hfill {\displaystyle p=\frac{\rho }{0.5m_{\mathrm{p}}}{k}_{\mathrm{B}}T,}\end{array}$$(4)

where m_p is the proton mass and k_B is the Boltzmann constant. The initial temperature ranges between 0.34 MK (interior) and 1.35 MK (exterior).

2.2. Driver

To reproduce the observed power distribution across different frequencies, we used a boundary driver to trigger propagating, transverse MHD waves along our magnetised cylinder, which was implemented as a 2D velocity field (along x and y directions) at the lower boundary of the domain, z = −20 Mm. In our reference setup, the velocity field is applied within a radius, a (the outer radius of the cylinder), of the centre of the lower boundary.

Our 2D velocity driver is based on the average observed spectrum of transverse oscillations in the solar corona as measured by Morton et al. (2019), which is obtained as a spatial and time average of the observed transverse oscillations in the low corona observed by the Coronal Multi-channel Polarimeter (CoMP, Tomczyk et al. 2008). This observed spectrum consists of a power law with a power index of α_p and an additional Gaussian power enhancement distributed around ν₀ = 4 mHz with maximum Δp = 0.002 in normalised units and standard deviation σ_p = 0.2 Hz. Our power spectrum is constructed such that the time integral of the velocity is zero, so the net displacement over the duration of the simulation is zero. In this study, we used the average observed spectrum as a reference for our investigation but also considered two modified spectra in order to analyse the role of the key features of the observed power spectrum (see Table 1 for the parameters used for the spectra). In the first modified spectrum, the power enhancement at ν₀ = 4 mHz was removed and uniformly redistributed across the spectrum to help us to understand the role of waves generated by the localised power enhancement on the loop structure and heating. For the second modified spectrum, we used a smaller power law index to prescribe a different distribution, with relatively more power at high frequencies than in the observed spectrum, in order to verify the efficiency of the heating if more high-frequency waves propagate into the corona from lower layers of the solar atmosphere.

Figure 3a shows the power distribution for our reference study based on the average observed spectrum (black), the spectrum without the local power enhancement (red), and the redistributed power spectrum (blue).

Fig. 3. a: spectra for transverse waves used in this study. The black curve represents the averaged observed spectrum, the red curve represents a spectrum where the power enhancement at ν0 = 4 mHz is removed, and the blue curve represents a spectrum where we used a shallower power-law index. b: velocity vx time series for the three spectra in normalised units with the same coloured legend. c: displacement associated with the velocity time series in normalised units and with the same colour legend. The black dashed curve represents the displacement derived from an additional velocity profile obtained from the same observed power spectrum.

To construct the 2D velocity driver, we derived two random time series from the given power spectra, one for the x-component of the velocity, v_x(t), and one for the y-component, v_y(t). Figure 3b shows the x-component, v_x(t), for the three cases in Table 1, and Fig. 3c shows the corresponding radial displacement of the centre of the cylinder obtained from the time integration of v_x(t) and v_y(t). We notice that the observed spectrum leads to significantly higher displacements than other spectra and we tested this result through a number of randomised time series. For example, the black dashed line represents a different displacement profile derived from the same observed power spectrum, and its displacement is still larger than the ones from the two modified spectra. In order to be able to draw conclusions on the relation between the driver spectrum and the MHD evolution, we ran several simulations with the same spectra but varied the derived time series.

In Pagano & De Moortel (2019), we modelled a velocity time series based the same observed power spectrum using a large number of pulses of different periods and amplitudes. In addition, in Pagano & De Moortel (2019), we assumed a uniform distribution of the period of the pulses, which, combined with the power spectrum, meant that low-frequency pulses had a substantially larger amplitudes than high-frequency ones, which contradicts more recent observational results by Morton et al. (2019). The approach presented in this current study is more general as it is not based on deconstructing the velocity driver into single pulses, and hence no assumptions are made about any potential correlation between period and amplitude.

2.3. MHD simulation

To study the evolution of this system, we solved the MHD equations numerically, where thermal conduction, magnetic diffusion, and joule heating are included as source terms as follows:

$$\begin{array}{cc}& {\displaystyle \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·(\rho \mathit{v})=0,}\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\displaystyle \frac{\partial \rho \mathit{v}}{\partial t}+\mathbf{\nabla }·(\rho \mathit{v}\mathit{v})+\mathrm{\nabla }p-\frac{\mathit{j}\times \mathit{B}}{c}=0,}\hfill \end{array}$$(6)

$$\begin{array}{cc}& {\displaystyle \frac{\partial \mathit{B}}{\partial t}-\mathbf{\nabla }\times (\mathit{v}\times \mathit{B})=\eta \frac{c^2}{4\pi }{\mathrm{\nabla }}^2\mathit{B},}\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\displaystyle \frac{\partial e}{\partial t}+\mathbf{\nabla }·\left[(e+\frac{{\mathit{B}}^2}{8\pi }+p)\mathit{v}\right]=-\eta j^2-\mathrm{\nabla }·{\mathit{F}}_{\mathrm{c}},}\hfill \end{array}$$(8)

where t is time, v velocity, η the magnetic resistivity, c the speed of light, $j=\frac{c}{4\pi }\mathrm{\nabla }\times \mathit{B}$ the current density, and F_c the conductive flux (Spitzer 1962). The total energy density e is given by

$$\begin{array}{c}\hfill {\displaystyle e=\frac{p}{\gamma -1}+\frac{1}{2}\rho {\mathit{v}}^2+\frac{{\mathit{B}}^2}{8\pi },}\end{array}$$(9)

where γ = 5/3 denotes the ratio of specific heats. In Pagano & De Moortel (2019), we described in detail why using values for the magnetic resistivity derived from plasma theory is currently unfeasible in state-of-the-art MHD simulations. Therefore, in this set of numerical experiments, we used an anomalous resistivity, parametrised to meet the following two competing requirements. On one hand, the dissipation of electric currents must be sufficiently efficient to convert a noticeable amount of energy into heating. On the other hand, the dissipation of currents must not damp the wave propagation too rapidly, otherwise they would not be observed. An anomalous magnetic resistivity is designed to meet these requirements, and thus

$$\begin{array}{c}\hfill \eta ={\eta }_0(|\mathit{j}|>j_0),\end{array}$$(10)

where we use j₀ = 2.5 G s⁻¹ as the threshold current, and η₀ is 10⁸ η_S, where η_S is the magnetic resistivity according to Spitzer (1962) at T = 2 MK. Additionally, this description of the dissipation mechanism allows us to run the simulations on a coarser grid. Generally, a rather small grid size is required to fully resolve the phase mixing or resonant absorption of MHD waves (e.g. Antolin et al. 2014; Pagano & De Moortel 2017; Karampelas et al. 2017; Howson et al. 2019). Such a high resolution allows for the description of the development of small-scale structures and the conversion of the wave energy into heating. However, in this work, we do not focus on the detailed description of the phase mixing, already addressed in Pascoe et al. (2010) and Pagano & De Moortel (2017), but instead are interested in the amount of wave energy that this can convert. Here, the threshold, j₀, is used to activate the dissipation mechanism only at the locations where the phase mixing happens and the high value of the magnetic resistivity allows for a significant and rapid conversion of wave energy into heating. Moreover, even if all the numerical experiments presented here are inevitably affected by numerical diffusion, with the chosen values of j₀ and η₀ the effects of the magnetic resistivity dominate on the effect of the numerical diffusion.

The computational grid has a uniform resolution of Δx = Δy = 0.15 Mm and Δz = 0.78 Mm. The simulation domain extends from z = −20 Mm to z = 180 Mm in the direction of the initial magnetic field and horizontally from x = −3 Mm to x = 3 Mm and from y = −3 Mm to y = 3 Mm, in order to contain the loop structure presented in Sect. 2.1. The boundary conditions are treated with a system of ghost cells, where we implemented periodic boundary conditions at both the x and y boundaries. For the upper z boundary, we used outflow boundary conditions, except when we wanted to focus on the dynamics occurring in closed loops, for which we used reflective boundary conditions (i.e. where the magnetic field components, density, and energy are symmetrically copied in the ghost cells and the velocity components are symmetrically copied and the sign is changed). The driver is set as a boundary condition at the lower z boundary. In addition, we added a damping layer between z = −20 Mm and z = 0 Mm, which only allowed the propagation of transverse oscillations, by multiplying the field-aligned component of the velocity, v_z, by a damping factor 0.75 after each time integration iteration in this part of the domain. We applied this improvement in order to prevent the propagation of slow waves, as this can significantly alter the energy budget. For this reason, we will only show the numerical domain for z ≥ 0, as we regard the solution that we find for z <  0 as non-physical. Finally, we scaled the amplitude of the driver to obtain transverse velocities of the order of v_x ∼ v_y ∼ 15 km s⁻¹ beyond the damping layer (i.e. above z ≥ 0), which is the order of magnitude of transverse displacements observed in the solar corona (Threlfall et al. 2013). We ran all our numerical experiments for 2224 s of physical time, which is sufficient for the driving to finish, plus an additional 80 s, which allows the last oscillations driven at the lower boundary to travel into the domain.

3. MHD evolution

We started by analysing the simulation with a uniform loop structure. The driver is applied within the radius, a, at the lower boundary and where we set open boundary conditions at the upper z-boundary. This is not necessarily the setup that best describes the conditions we find in the solar corona, but it is a useful reference simulation to compare with others to measure the effect of the density structure and the driver on the wave heating and the evolution of the boundary shell.

As the driver is applied at the lower boundary, the cylinder is displaced from its initial position and transverse waves propagate along the loop, with velocities of the order of 10 km s⁻¹. In addition, some remnant longitudinal perturbations of the order of 1 km s⁻¹ still propagate past the damping layer as well. The transverse waves propagate at different speeds across the boundary shell, and the standard phase-mixing features develop. At the same time, the displacement of the magnetised cylinder leads to the interaction between the denser interior and boundary-shell regions with the background generating compression and rarefaction just outside the boundary shell. Figure 4 shows the evolution of density and electric currents at an earlier stage (t = 617 s) and at the end of the simulation (t = 2224 s). We find that the density structure evolves significantly throughout the simulation and the initial circular cylindrical structure is replaced by a more fragmented one that shows signatures of eddies at the boundary shell. Hence, our loop structure is strongly disrupted as a result of the transverse waves and filamentary structures that are created. Standing transverse waves are known that are unstable to the KHI (Terradas et al. 2008), and propagating transverse waves show uniturbulence (Magyar et al. 2017). Since we are considering propagating waves, the occurrence of uniturbulence is more likely here, but the connection or difference between KHI and uniturbulence is currently not well understood.

Fig. 4. Density and electric current density isosurfaces at two different times (t = 617 s and t = 2224 s) for the simulation with a uniform density structure of the loop. A movie of the density evolution is available online.

Electric currents are generated across the boundary shell and dissipated by the effect of magnetic resistivity. Figure 4 shows the contour of the modulus of the electric current at the level of the threshold current at which the anomalous resistivity is triggered. Electric currents are generated by both the phase mixing of Alfvén waves and the shearing between different regions. While near the beginning of the simulation the electric currents are more uniformly distributed around the magnetised cylinder, their distribution becomes significantly more fragmented later on. Pagano & De Moortel (2019) provided a discussion on how different frequencies of the boundary-driven oscillations trigger different kinds of electric currents.

At the end of the simulation, the magnetised structure has been altered substantially and the density structure is very fragmented, having lost any resemblance of a cylinder. For the purpose of this study, we focus on the evolution of the boundary shell and the distribution of the plasma heating in time and space.

3.1. Boundary shell

As the density structure is altered by the propagating transverse waves, the boundary shell evolves in time. In order to follow this evolution, we need an operative definition of this region. For the purpose of this study, the key property of the boundary shell is the gradient of the Alfvén speed. Here, we define the boundary shell as the region where ∇V_A >  0.5 s⁻¹, a threshold chosen to match the boundary shell at t = 0 s.

Figure 5 shows the density distribution at the middle of the loop, z = 79.6, at t = 0 s, and at two times during the evolution. Contours for ∇V_A = 0.5 s⁻¹ (green contour) and ∇V_A = 1 s⁻¹ (blue contour) are over-plotted, as well as velocity vectors where |v| > 5 km s⁻¹. At t = 0 s, the boundary shell is a ring around the loop interior (region between the green contours), but it has significantly transformed at t = 741.4 s, when it is both displaced and distorted, although it still surrounds a denser (interior) region. Later, the boundary shell becomes a filamentary structure (t = 2241 s), indicating the presence of uniturbulence eddies around the loop structure. At this stage, it is no longer possible to identify a proper loop interior. Similar dynamics are present along the entire loop structure.

Fig. 5. Density cross-sections at z = 79.6 Mm at three different times (t = 0, t = 741.4 s, and t = 2241 s). Over-plotted are velocity vectors where |v| > 5 km s−1, and contours of the gradient of the Alfvén speed ∇VA = 0.5 s−1 (green contour) and ∇VA = 1 s−1 (blue contour).

In order to follow the distortion of the boundary shell, we continued to measure where ∇V_A >  0.5 s⁻¹, and, consequently, this is an evolving structure that can also appear or disappear where the gradients of the Alfvén speed vary significantly, thus a dynamic boundary shell similar to those of Antolin et al. (2015), Howson et al. (2017), Karampelas et al. (2017), and Magyar et al. (2017). It is not immediately clear whether such evolution increases the efficiency of the boundary shell in terms of providing favourable conditions for the dissipation of transverse waves. To gain more insight, we investigated three characteristics of the boundary shell: (i) its spatial extent, (ii) its efficiency (i.e. the Alfvén speed jump across the boundary shell), and (iii) its consistency (i.e. the coherence of the boundary shell in the field-aligned direction).

First, we focus on the area that is covered by the boundary shell (i.e. its spatial extent) at different z-coordinates along the loop. The upper panel of Fig. 6 shows the area expansion of the boundary shell as a function of the z-coordinate and time, normalised to the initial surface area. We find that the area covered by the boundary shell remains roughly unchanged for the first 600 s throughout the cylinder, apart from near the lower (driven) boundary, where it starts expanding earlier. After t = 600 s, the boundary shell starts expanding everywhere along the magnetised cylinder, doubling its surface near the end of the simulation. This is the time when the uniturbulence is fully developed (see Fig. 4).

Fig. 6. Map of expansion of boundary shell as a function of z and t for the simulation with a uniform loop (upper panel). Evolution of the boundary-shell efficiency in the same simulation (lower panel).

The volume of the boundary shell follows a similar evolution, remaining roughly constant for the first 600 s and then linearly expanding as soon as the uniturbulence sets in, reaching double its initial volume by the end of the simulation.

To examine the boundary-shell efficiency, we considered the total Alfvén speed variation within the boundary shell, ∫∇V_A(t)dV, normalised to the value of the same integral at t = 0 s (bottom panel Fig. 6). The larger the Alfvén speed variation, the faster the phase mixing of Alfvén waves can occur (see e.g. Heyvaerts & Priest 1983). We find that this value is up to six times larger by the end of the simulation, after a more efficient boundary shell starts to develop from around t = 600 s. Hence, after the first transitory phase of the simulation, the boundary shell expands and becomes more efficient, in principle creating more favourable conditions for wave heating to occur.

Finally, we investigated the consistency of the boundary shell (∫∇V_Adz). In order for phase mixing of the waves to develop, neighbouring propagating waves need to travel over some distance whilst a gradient in the Alfvén speed remains present. Figure 7 (top left) shows a cut of the gradient of the Alfvén speed along a vertical plane at y = 0 at t = 0. The gradient of the Alfvén speed is uniform in the z-direction and waves travelling from the lower boundary will phase mix over the full longitudinal extension of the loop. However, this configuration significantly changes over the course of the simulation, as shown in Fig. 7 (top right), at t = 1755 s. As the propagation of the waves alters the loop structure, the gradient of the Alfvén speed changes significantly in the z-direction, and transverse waves have a smaller distance over which to consistently develop phase mixing. This dynamic is specific to this footpoint driver, which continuously changes, hence breaking any initial invariance along the z-direction. In such a setup, the boundary shell constantly evolves in different ways, at different heights along the loop, leading to a reduced consistency. We now define the projected area onto the xy-plane of the region, where the boundary shell covers at least 90% of its initial z-extension as a measure of the boundary-shell consistency. The lower panel of Fig. 7 shows that the boundary-shell consistency drops to as low as 10% of its initial value during the simulation, increasing again to about 20% near the end, in the time interval when the simulation runs after the driver is switched off.

Fig. 7. Cross-section of the gradient of the Alfvén speed at y = 0, at the initial condition t = 0 s (a), and at t = 1977 s (b). Evolution of the boundary-shell consistency for the simulation with uniform loop.

In conclusion, the propagation of the transverse waves along the loop alters the initial boundary shell significantly throughout the simulation. The boundary shell expands and the gradient of the Alfvén speed increases over time, making the boundary shell larger and more efficient, in principle increasing the efficiency of phase mixing. However, at the same time, the boundary shell becomes less consistent, reducing the opportunity for waves to phase-mix. In the next section, we focus on the plasma heating in the boundary shell and compare it with the expected radiative losses.

3.2. Heating

During the simulation, most of the electric currents form in the boundary shell. As we state in Eqs. (8) and (10), when electric currents grow above a designated threshold, they are dissipated and ohmic heating occurs, where magnetic energy is transformed into internal energy of the plasma.

Figure 8 shows the distribution of this ohmic heating at two representative times in the simulation, namely at t = 617 s, just after the uniturbulence has developed, and at the end of the simulation (t = 2224 s), when the lower boundary driver has already stopped. We find that at earlier times, the heating is mostly uniformly distributed around the boundary shell, although the random footpoint motions have already created some fragmentation. In contrast, at later stages, the heating is distributed in a more heterogeneous way and the heating region is now very fragmented.

Fig. 8. Isosurface of ohmic heating at H = 108 erg cm−2 s−1 at two different times (t = 617 s and t = 2224 s) in the simulation with a uniform loop. A movie of the heating evolution is available online.

In this experiment, the magnetic resistivity is sufficient to allow for the dissipation of the waves, but it is not strong enough to damp the wave propagation fully. Figure 9 shows the heating averaged over the boundary shell cross-section as a function of the z coordinate along the loop and time. By comparing the evolution of the wave energy and thermal energy of this simulation with a corresponding ideal simulation, we confirm that the effect of the anomalous resistivity is significantly larger than the effect of the numerical dissipation.

Fig. 9. Map of average heating in boundary shell along the loop as a function of time for the simulation with a uniform loop. Although the heating is generally within the range shown here, the maximum heating is 6 × 109 erg s−1 cm2.

We find that several distinct heating events can be identified, even if waves are continuously travelling along the loop. Each event leaves a distinct signature in the time-distance map, as the heating at different z-coordinates occurs at different times. Similarly, it is clear that a large portion of the heating is deposited closer to the footpoint, as this is where stronger currents are generated. While the majority of the heating events start from the footpoint and are continued along the loop structure as the waves propagate, there are some that only form at a higher z-coordinate.

The discrete nature of the heating is also evident from the time evolution of the average heating in the boundary shell at different z coordinates, as shown in Fig. 10. We find that the time evolution of the heating can be described as a sequence of distinct events, where different amounts of energy are released for each event. While Fig. 9 clearly illustrates that these single peaks are part of more elongated heating structures, it should be noted that this pattern is not evident when focusing on fixed z-coordinates. Figure 10 illustrates again that the heating near the driven footpoint (black curve) is larger than the heating occurring further along the loop.

Fig. 10. Average heating in boundary shell as a function of time for the simulation with a uniform loop, at three different coordinates along the loop.

Finally, we compare this heating deposition with an estimate of the radiative losses, in order to understand to what extent this heating can maintain the million-degree corona. In Fig. 11, we compare the radiative losses with the heating as functions of the position along the loop and time. Both quantities are integrated over the boundary shell at each z coordinate. In particular, we compare the time integral of the average heating and average radiative losses in a cross-section (z) of the boundary shell between two simulation snapshots around the time, t. For the purpose of the time integral, the heating is approximated as a step function between the simulation snapshots. The time integral of the radiative losses is computed for a plasma whose temperature and density linearly evolve from the average density and temperature in the boundary shell of one simulation snapshot to the other. For this purpose, the radiative losses from a plasma were computed using the piecewise continuous function by Klimchuk et al. (2008). It should also be noted that the temperature variations are rather small in this simulation and the variation in the plasma radiative losses are mostly due to density variations. The radiative losses decrease in time because the boundary shell becomes less dense as it evolves. In contrast, the heating evolves in a less monotonic way, but in this simulation the heating power from the transverse waves always remains only a small percentage of the power required to balance the radiative losses (of the order of 1% of the radiative losses, with one short-lived peak of 10%).

Fig. 11. Map of ratio of average heating and average radiative losses in the boundary shell as functions of loop position and time, for the simulation with a uniform loop.

Driving the lower boundary with a different velocity time series derived from the same power spectrum but with a smaller displacement (dashed line in Fig. 3b), we find very similar results. Although the maximum contribution of the wave heating goes up to 13% of the expected radiative losses, overall it remains a small percentage of them, indicating that our conclusions do not significantly depend on the chosen time series.

When we compare the rate of change of the thermal energy contained in the computational domain in this simulation with a corresponding ideal simulation (η = 0), we find that, as expected, the increase of thermal energy in the non-ideal simulation initially exceeds the thermal energy increase in the ideal simulation because of the ongoing heating. However, after some time, the situation reverses, because the thermal pressure excess in the non-ideal simulation leads to a gradient of the thermal pressure in the z direction and consequently, mass flows across the upper boundary of the domain. These flows lead to a mass loss, and, hence, a thermal energy loss through this boundary. Such dynamics are compatible with the wave heating of open field lines in the corona and the triggering of plasma outflows. In our simulations, we find a mass loss of 6 × 10⁷ g over a time of 2200 s, corresponding to an average mass-loss rate of 2.5 × 10⁴ g s⁻¹ across the 36 Mm² surface of the outer boundary. If we consider a characteristic density of 10⁻¹⁶ g cm⁻³, we need a steady outflow of ∼0.01 km s⁻¹ to account for this mass loss.

4. Boundary shell evolution

In this section, we investigate how different initial loop structures affect the evolution of the boundary shell. In particular, we consider the non-uniform density distributions (with z₀ = 130 Mm and z₀ = 46 Mm) presented in Fig. 1, where both the density along the loop and the extent of the interior region decrease with height.

Figure 12 shows the final contour of the boundary shell for these two simulations. Compared to their initial structure in Fig. 1, we find that in both cases the boundary shell has expanded in the x and y directions and extended further upwards in the z direction. Although the expansion is already present at the lower boundary, it is most significant at higher z coordinates, where the relatively shallow boundary shell was previously not present. Additionally, the contours in Fig. 12 extend further in the z direction compared to the initial configuration, as the propagation of transverse waves along the loop structure leads to (i) steeper gradients of the Alfvén speed, and (ii) the generation of a boundary shell where this was initially not present. The first mechanisms occur quickly, as the gradient of the Alfvén speed steepens as soon as the waves propagate (∼120 s for both simulations) and the entire loop structure up until z₀ is surrounded by a cylindrical boundary shell. In contrast, the second process is slower and becomes visible only after about 1000 s, by which time, uniturbulence has already started affecting the boundary shell.

Fig. 12. Isosurfaces of Alfvén speed gradient at t = 2224 s for the simulations with z0 = 130 Mm (left-hand side) and z0 = 46 Mm (right-hand side).

The expansion of the boundary shell in these simulations (Fig. 13) follows a similar pattern to that in Fig. 6 (where the loop structure is present from the beginning), as it slowly expands to cover a surface twice as large at the initial one. However, the expansion appears much larger where no initial boundary shell was present, and this is illustrated by the red region (colour saturation) above z₀. This effect is uniform from z₀ to the end of the domain for the simulation with z₀ = 130 Mm after t = 500 s. In contrast, in the simulation with z₀ = 46 Mm we find a time evolution where the boundary shell expands above z₀ from t = 500 s at about 0.2 Mm s⁻¹. In summary, the simulations with a uniform boundary shell and with z₀ = 130 Mm both show a final volume expansion of about 2 by the end of the simulation, whereas the simulation with z₀ = 46 Mm shows a maximum expansion of 2.7. Interestingly, this process also leads to a more efficient boundary shell with respect to the phase mixing, as our measure of the efficiency of the boundary shell is about 35% higher when we consider a non-uniform loop. The time evolution pattern of the efficiency is very similar for all simulations (Fig. 13). Additionally, these boundary shells also extend in the z direction while expanding in the x − y plane, and therefore their consistency is affected too. Figure 14 compares the consistency evolution for the three simulations, relative to their respective initial consistency. We find that the consistency of the boundary shell of a uniform loop quickly drops as the simulation starts. The simulation with z₀ = 130 Mm undergoes a smaller initial decrease, followed by an increase to a value close to 1 for about 500 s, before decreasing again to smaller values, and finally ramping up again when the loops settle to a new equilibrium near the end of the simulation, after the boundary driver stops. Although fluctuating, some consistency is clearly present for a significant fraction of the time evolution. Finally, the simulation with z₀ = 46 Mm shows a nearly constant boundary-shell consistency throughout the simulation. Although it evolves, the extension in the z direction remains relatively constant. Near the end of the simulation, because the original boundary shell has extended in the z direction and because of the newly reached equilibrium, the consistency measure is twice as large as at t = 0.

Fig. 13. Map of boundary-shell expansion as a function of z and t (upper row) and evolution of the boundary-shell efficiency (lower row). for the simulations with z0 = 130 Mm (left-hand side) and z0 = 46 Mm (right-hand side).

Fig. 14. Evolution of boundary-shell consistency for the simulation with uniform loop, with z0 = 130 Mm, and z0 = 46 Mm.

In conclusion, the evolution of the boundary shell is significantly affected by the initial structure of the loop. We find that the presence of an interior region as a central uniform density enhancement is not required to form loop-like structures. A shallow initial boundary shell can be sufficient to develop a larger, steeper, and more consistent boundary shell around the loop. In fact, we find that an initially small density enhancement (here represented by the simulation with z₀ = 46 Mm) develops a relatively more favourable boundary shell compared to a fully formed uniform loop.

5. Heating the corona

In this section, we analyse how the heating resulting from the dissipation of transverse waves depends on a number of parameters. In particular, we investigated the role of the initial loop structure, the extent of the footpoint driver region, the reflection of the waves at the upper footpoint and the dissipation coefficients. In Sect. 3, we show that the dissipation of MHD waves contributes to at least 1% of the energy budget needed to maintain the thermal structure of the corona by comparing the ohmic heating in the simulation with the estimated radiative losses in the boundary shell of the loop. This finding is in agreement with our previous theoretical studies (Pagano & De Moortel 2017, 2019; Pagano et al. 2018). However, the estimates are sensitive to the specific modelling setup, and there are several physical reasons why the heating contribution could be enhanced. In this section, we examine when the heating contribution from the dissipation of MHD waves can be more effective using the simulation described in Sect. 3 as a reference.

5.1. Density structure

The total heating input is not significantly affected by the initial loop structure, however its distribution is. Figure 15 shows the logarithm of the average heating in the boundary shell as a function of the coordinate along the loop and time, normalised by the same heating quantity in our reference simulation, for the two simulations with z₀ = 130 Mm and z₀ = 46 Mm. For the most part, we find that the heating is not significantly different with respect to the reference simulation (i.e. with a uniform boundary shell); noticeable differences are only found in the portion of the loop above z₀, where the heating is substantially smaller during the time when the uniturbulence is not yet effective (t <  60 s). Near the end of the simulation, it settles at similar values. At the same time, the density is lower in the non-uniform density structures, and hence radiative losses are smaller in comparison to the reference simulation. Therefore, the heating profile can now contribute to up to 8% of the radiative losses in these non-uniform configurations.

Fig. 15. Maps of average-heating boundary-shell logarithm in the simulation with (a) z0 = 130 Mm or (b) z0 = 46 Mm normalised to the reference simulation and plotted as a function of position along the loop and time.

5.2. Role of the footpoint motions

We have so far used a driver that only affects the central part of the lower boundary (where the loop structure is initially situated) and outflow boundary conditions at the upper boundary. This scenario represents a situation where only one footpoint of a loop is subject to transverse oscillations and any wave energy not dissipated along the loop flows out at the opposite footpoint (or the structure is an open-field coronal structure). Figure 16 shows the heating profiles normalised to the reference simulation (as in Fig. 15) for two simulations where we changed these assumptions (i.e. where the outer boundary is reflective or where the driver acts on the entire lower boundary).

Fig. 16. Maps of average-heating boundary-shell logarithm in the simulation with (a) a closed loop, or (b) a wider driver normalised to the reference simulation and plotted as a function of position along the loop and time.

In the simulation with reflective upper-boundary conditions, the transverse waves cannot leave the domain from that boundary, and the propagating waves now interact with counter-propagating ones. This can make the phase mixing more efficient and further accelerate the expansion of the boundary shell. Indeed, we find that the volume of the boundary shell in this simulation expands by a factor of 3, which is larger than previously found. In particular, the boundary-shell expansion now also happens closer to the upper boundary. From an energetic point of view, we find that the thermal energy included in the domain in this simulation is always larger than the thermal energy enclosed in the respective ideal simulation, showing that the excess in thermal energy is due to the heating from the dissipation of waves. Moreover, the heating pattern differs from the one found in the simulation with an open boundary condition, as stronger heating signatures are now also present at the upper footpoint (see Fig. 16a). In this simulation, the time-integrated heating is enhanced by a factor of 2.2.

At the same time, when we consider a different extension of the boundary driver, the heating output differs significantly. Figure 16b shows the heating deposition for the simulation where the boundary driver is applied to the entire lower boundary. The heating is now always larger than in the reference simulation, by a factor of at least 5 and with peaks of up to a factor of 50. For reference, in this configuration, about eight times more energy is injected into the domain by extending the driver to the full lower boundary, in particular including the lower-density loop exterior. There are two key differences with respect to the reference simulations. Firstly, the loop interior and exterior move in phase near the lower boundary, and there is no damping of the oscillations of the flux tube due to its interaction with the background. This allows more wave energy to be converted into thermal energy via phase mixing and thus leads to higher ohmic heating. Additionally, as already shown in Pascoe et al. (2010), because of the mode coupling between the kink and Alfvén modes, wave energy flows as Poynting flux from the exterior to the boundary shell, where it can then be dissipated through phase mixing. Figure 17 compares the Poynting flux component across the external surface of the boundary shell for our reference simulation and for the simulation where the driver is applied to the entire lower boundary at the same time, t = 988 s (i.e. when the boundary shell has not yet been fragmented). We find that the Poynting flux is up to two orders of magnitude larger in the latter simulation.

Fig. 17. Isosurface of Alfvén speed gradient for the reference simulation (left-hand side) and for the wider driver simulation (right-hand side) at t = 988 s, coloured by the modulus of the component of the Poynting flux across the isosurface.

5.3. Higher velocity driver

In our reference simulations, we adjusted the amplitude of the driver in order to have transverse velocities of the order of v_x ∼ v_y ∼ 15 km s⁻¹ beyond the damping layer (i.e. past z = 0 Mm). However, because of the time variation of our driver, the speeds can drop to only a few km s⁻¹. As larger transverse velocities are regularly observed in the corona, we also ran a simulation where the speed of the driver was increased by a factor of 2.5 with respect to the simulation described in Sect. 3. Figure 18 shows the normalised heating for this simulation, which is uniformly higher than the corresponding heating in the reference simulation by a factor of at least 2, with peaks of up to 15. In this simulation, the time-integrated heating is enhanced by a factor of about 4.5, compared to an increase of a factor of 6.25 in the kinetic energy. Assuming we can extrapolate this scaling further (we require the enhancement of the heating by a factor of 100) would imply that we need to increase the amplitude of the driven waves by a factor of 11. However, non-linear effects could certainly play a role with larger amplitudes and such a dramatic increase might not be necessary in order to reach a similar heating enhancement. Either way, substantially larger amplitudes are needed to supply the entire heating requirement via transverse wave dissipation.

Fig. 18. Map of average heating boundary-shell logarithm in the simulation with a stronger driver, normalised to the reference simulation and plotted as a function of position along the loop and time.

5.4. Higher dissipation coefficient

Higher values of the resistivity coefficient or lower values for the critical current at which the resistivity is triggered might lead to more effective heating. We therefore ran two further simulations: one where η₀ = 10¹⁰ η_S (i.e. enhanced by two orders of magnitude compared to the simulations so far) and one where this higher value of the resistivity coefficient is used uniformly, without a critical current trigger.

Figure 19 shows the normalised heating for these two new simulations. The simulation where we simply enhance the resistivity coefficient (Fig. 19a) shows a larger overall heating compared to the reference simulation, particularly near the footpoint, whereas the simulation without the critical current threshold (Fig. 19b) shows stronger heating along the whole loop. For both simulations, the heating increases by more than one order of magnitude, where this is confined to the region near the footpoint when a critical current to trigger the resistivity is used.

Fig. 19. Maps of average-heating boundary-shell logarithm in the simulation with (a) η = 1010 ηS and (b) η = 1010 ηS and j0 = 0, normalised to the reference simulation and plotted as a function of position along the loop and time.

This higher dissipation and heating happens at the expense of the wave energy. However, the occurrence of heating leads to the generation of additional transverse velocities from local transverse pressure gradients, and therefore some of the plasma kinetic energy cannot be ascribed to the propagation of MHD waves. However, we can use the property of equipartition of kinetic and magnetic energies in MHD waves to identify the energy associated with transverse waves in the simulation defined as:

$$\begin{array}{cc}\hfill E_{\mathrm{W}}=& \frac{1}{2}[\frac{1}{2}\rho (v_x^2+v_y^2)+\frac{1}{8\pi }\rho (B_x^2+B_y^2)]\hfill \\ \hfill & -\frac{1}{2}\sqrt{{[\frac{1}{2}\rho (v_x^2+v_y^2)-\frac{1}{8\pi }\rho (B_x^2+B_y^2)]}^2}.\hfill \end{array}$$(11)

Figure 20 shows the average wave energy in the boundary shell along the loop as a function of time for the reference simulation and for the simulation with η₀ = 10¹⁰ η_S and j₀ = 0.

Fig. 20. Maps of average wave energy in the boundary shell as a function of time and position along the loop for (top) the reference simulation and (bottom) the simulation with η0 = 1010 ηS and j0 = 0.

We find that while the wave energy decreases only marginally along the loop in the reference simulation, it does so significantly in the other simulation, implying that the wave amplitude is damped substantially along the loop by the effect of the resistivity.

5.5. Heating summary

We can use our analysis of a number of different physical configurations to shed light on the conditions under which MHD waves can contribute significantly to the heating of the solar corona. Generally, we find that the dissipation of MHD waves via phase mixing can certainly contribute to about 1% of the energy that is radiated by the solar corona, but in certain configurations, the heating can be significantly increased.

To show the heating enhancement, in Fig. 21a we plot the average heating in the boundary shell as function of time for some of the key simulations in this study. We report in the figure how much the heating is enhanced, on average, with respect to the reference simulation. We find that the heating is, on average, enhanced by a factor of at least 4.5 when the extent of the footpoint driver is larger than the size of the loop interior or when a higher resistivity coefficient without a critical trigger current is used, or for a higher amplitude driver. One could argue that if several of these enhancements were to occur simultaneously, the wave heating may become sufficient to counteract the radiative losses.

Fig. 21. a: average heating in boundary shell as a function of time for a number of simulations listed in the colour legend, with the average heating enhancement reported. b: average heating in the boundary shell as a function of position along the loop for the same simulations normalised to the average heating at z = 0.

These different configurations also have implications for the nature of the wave generation and propagation in the solar corona. In particular, different configurations lead to different heating locations. Figure 21b shows the time-averaged heating as a function of the coordinate along the loop, normalised to its value at z = 0 Mm. When we have higher resistivity, the heating is highly localised at the footpoint. Similarly, for the simulation with a closed loop, most of the heating takes place near the reflective footpoint. Other simulations show less concentrated heating that generally decreases with z, with the exception of the simulation with a large-area driver, where the continuous inflow of wave energy in the boundary shell from the exterior region leads to a distinct heating profile, peaking at the centre of the loop rather than the footpoints. Key aspects of the different heating profiles for the different simulations are summarised in Table 2.

6. Relevance of the observed power spectrum

Finally, we investigate how using different power spectra for the boundary driver can affect the heating output. In particular, we used the different power spectra in Fig. 3a to investigate the importance of the key features associated with the power spectrum observed by Morton et al. (2019), namely the excess in power around the four-minute oscillations and the particular distribution of the power.

Figure 22a shows the average heating in the boundary shell for the four simulations based on the power spectra shown in Fig. 3a. The simulations all use a uniform boundary shell. Two of these simulations (black lines) are driven by time series that are based on the same (observed) power spectrum. The simulation without the excess power at four minutes is represented by the red lines, and the blue lines correspond to the simulation for which the power was redistributed with a more shallow gradient (i.e. the power associated with the long-period oscillations is reduced, and the power at the short-period oscillations is increased). As the driver is always applied to the centre of the domain, the amount of energy (i.e. Poynting flux) that is injected into the system, and hence available for heating, crucially depends on the amount of energy in the driver, on the power spectrum, and the position of the waveguide. In order to isolate the effect of the spectrum, we here compare simulations for which the maximum displacement during the time series is comparable.

Fig. 22. a: average heating in boundary shell as a function of time for the simulations with four different time series, as explained in Fig. 3b. b: average heating in the boundary shell as a function of the position along the loop for the four different time series, as explained in Fig. 3b.

We find that the heating profiles for the simulations using the observed power spectrum and the simulation without the four-minute excess power are quite similar. In contrast, the simulation with a shallower spectrum consistently leads to higher heating. The heating increase in this simulation is consistent with the additional power at the higher frequencies with respect to the observed power spectrum, confirming that high-frequency oscillations dissipate more efficiently.

Figure 22b shows the time-averaged heating as a function of the coordinate z along the loop for these simulations. For three of the simulations (red and black lines), we find mostly similar heating profiles, with the heating largely uniform along the loop. The difference between the two black lines (both based on the observed power spectrum) can be explained by the different footpoint displacements (see Fig. 3c). For the loop with the larger displacement (thick black line), phase mixing is most efficient low down in the loop, as the larger footpoint displacements reduce the field-aligned consistency of the boundary shell at larger z. However, the heating profile is qualitatively different for the simulation based on a shallower spectrum. In this simulation, the combination of the enhanced power at higher frequencies and the more efficient dissipation of the high-frequency oscillations leads to a distinct heating profile, concentrated near the driven loop footpoint.

7. Discussion and conclusions

In this paper, we expand the investigation of how wave energy is transformed into plasma heating when transverse MHD waves travel along coronal loops. In particular, we based this investigation on an observed power spectrum of transverse waves as well as analysing how other parameters affect this process. We focused on (1) the initial loop structure where we consider uniform loop structures composed of an interior region and a boundary shell, or loop structures where the loop interior vanishes and the boundary shell is present only for a portion of the loop length; (2) the behaviour of the footpoints of the loop; (3) the energy available for dissipation; and (4) the power spectrum itself. The aim of this work is to show how the wave heating contribution is affected by these parameters and how the loop structure (mostly the boundary shell) evolves in these circumstances.

The evolution of the boundary shell is a crucial aspect for coronal heating models. Small spatial-scale structures where electric currents can develop and be dissipated are a universal requirement of coronal heating models, regardless of the mechanism(s) responsible for the ultimate heat deposition. Examples include wave heating models, which need thin boundary shells where the waves can quickly undergo phase mixing, or nanoflare models, which need small structures where the magnetic field is sufficiently tangled for magnetic reconnection to release the free energy. The boundary shells, where the Alfvén speed changes rapidly, are precisely those regions where we find the necessary small scales. In this paper, we show not only that waves can produce such boundary shells, but also that the observed spectrum in the solar corona produces perturbations that contributed to the formation of boundary shells. We analysed three different properties of the boundary shells: its spatial extent (both in terms of volume and in terms of surface at different loop cross-sections), its efficiency (measured by the variation of the Alfvén speed across the boundary), and its field-aligned consistency (measured by the distance along the loop over which an Alfvén speed gradient is present).

We find that the extent and efficiency of the boundary shell increase in time when transverse waves propagate along the loop, where this process accelerates after the development of uniturbulence. On the other hand, the field-aligned consistency of the boundary shell is reduced by the displacement of the loop from its initial position and increases only when the loop has time to relax to a new equilibrium. It is not clear though, whether such a new equilibrium is possible in a scenario where waves are present at all times. At the same time, we find that an initially non-uniform loop where no loop interior is present and where the density enhancement is concentrated near the lower boundary can increase the spatial extent, efficiency, and consistency of its boundary shell. In simple terms, this happens because such an initial loop structure is more able to expand in 3D and to enhance its Alfvén speed gradient with respect to a uniform loop. This result is particularly important, as it shows that transitory density enhancement (such as spicules), which is subject to the continuous propagation of waves, can have a boundary shell expanding around and above where either wave heating or other heating mechanisms can be triggered.

Figure 23 summarises the heating output in our simulations (normalised to the reference simulation) versus an estimate of the Poynting flux injected directly in the boundary shell of the loop from the lower boundary. This scatter plot shows that, depending on the density structure and the various parameters we have considered here, the heating can be enhanced by more than one order of magnitude. The key result here is that to achieve stronger heating, injecting more energy in the loop structure is not necessarily required (as the stronger driver or wider driver simulations show), but that, for example, more effective dissipation can achieve the same result. Of particular interest is the simulation of the closed loop, which shows a very small Poynting flux with respect to other simulations, but significant heating. This occurs because at the lower boundary, where the Poynting flux is calculated, energy flux, from the waves reflected at the other boundary, also leaves the domain.

Fig. 23. Scattered plot of total heating as a function of the total Poynting flux injected in the boundary shell. Both quantities are shown with respect to the reference simulation. Colours and symbols are clarified in the plot.

Overall, the heating from the dissipation of transverse MHD waves does not seem to be sufficient to compensate for the energy radiated by the plasma. We find that the wave heating can contribute only about 1% of the radiative losses. Factors that can lead to a more substantial contribution are (1) a non-uniform loop structure with a lower overall density, and thus fewer radiative losses; (2) a driver region that is larger than the loop cross-section, in order to facilitate the generation of waves, a larger energy flow into the boundary shell, and thus a larger increase in the available energy than can be converted into heating; (3) efficient reflection of the waves at the loop footpoints; (4) larger wave amplitudes; or (5) a more effective dissipation process. Although these factors could increase the heating deposition from the transverse waves, their feasibility and effectiveness needs further investigation. For example, reflection and transmission at the transition region depends on the wavelength of the waves compared to the extent and density jump of the transition region, but so does the efficiency of the wave dissipation through phase mixing. On one hand, the dissipation of high-frequency waves is more efficient, but at the same time, the energy of these waves is also more easily lost from the corona through transmission to the lower atmosphere (see e.g. Hollweg 1984a,b; Berghmans & de Bruyne 1995; De Pontieu et al. 2001; Van Damme et al. 2020). Hence, it is not, a-priori, obvious whether higher frequency waves would indeed lead to increased heating. Finally, in this work, we do not focus on the detailed development of the phase mixing and the consequent energy conversion. A possible and interesting continuation of this work would be to devise a computationally feasible simulation to investigate in detail how much our approach using the anomalous resistivity departs from a detailed description of the phase-mixing process.

A more efficient wave-energy-dissipation mechanism would imply significant damping of the wave amplitude along the loop. Can such a decay be measured? Some indications from the analysis of propagating kink waves suggest that the damping of waves can be measured in the corona, and therefore provide fundamental constraints for these theoretical studies (e.g. Tiwari et al. 2019). Finally, having driver cells larger than the loop would also increase the Poynting flux and thus the energy available. However, this can be an elusive measurement, as the size of a single magnetic thread in coronal loops is still not clear (see e.g. Brooks et al. 2016; Aschwanden & Peter 2017; Williams et al. 2020) and it is difficult to measure the flow speeds in the tenuous background corona. Moreover, it remains to be investigated how this scenario relates to the wave heating described in van Ballegooijen et al. (2011, 2017), where a horizontal motion is applied over a spatial scale smaller than the cross-section of the loop. In that granular-driven scenario, driving happens below the corona on spatial scales smaller than magnetic bright points. However, a magnetic bright point will likely be the footpoint of a number of coronal loops, as they are at least ∼100 km in size in the photosphere (Alipour & Safari 2015), and coronal loops are about ∼300 − 400 km in width (Aschwanden & Peter 2017), whereas the magnetic field expansion between the photosphere and the corona would cover a larger area. Hence, once the field expansion is taken into account, the coherence of the waves that reach the corona will be larger, and our wider driver scenario may not be in contradiction with the scenario modelled by van Ballegooijen et al. (2011). Additionally, if transverse waves are driven by p-modes, whose fronts are thought to be larger than the small-scale dynamics of the turbulent photospheric flows (Bloomfield et al. 2007; Felipe et al. 2010), the coherence would be even larger and more likely to occur. Higher amplitude waves can also lead to more effective heating, but observations seem to measure velocities on average of the order of ∼20 km s⁻¹ (Weberg et al. 2018), possibly up to ∼50 km s⁻¹ (McIntosh & De Pontieu 2012), which would not be sufficient in the scenario we present here. However, line-of-sight integration might mean that higher velocities could be present, but so far, measurements of non-thermal line widths have not given conclusive evidence of the amplitudes of coronal waves (e.g. Doschek et al. 2008; McIntosh & De Pontieu 2012; Brooks & Warren 2016).

Furthermore, our time-distance plots of the heating along the loops show that although there is a slight concentration near the footpoints, over time the heating is mostly distributed along the loop. If we focus on a specific loop location, our simulations show impulsive heating patterns. Whereas the impulsive nature of the heating appears to agree, the location of the heating, on the other hand, seems in contradiction with at least some observational studies where the heating location has been found to be concentrated away from the loop footpoint (Reale et al. 2000a,b; López Fuentes et al. 2007; Song et al. 2015). Overall, identifying whether any of the mechanisms suggested here to increase the efficiency of wave heating are viable in the solar corona will require further detailed modelling as well as higher cadence and higher resolution observations.

Our study also provides some predictions on the location of the heating along the loop for different loop and driver configurations. These results and the observable consequences are summarised in Fig. 21b and Table 2. Galsgaard et al. (1999) provided a similar summary for heating by braiding, differentiating between heating profiles depending on the time interval between two braiding events. Similarly, Reale et al. (2000b) highlighted the need to constrain the heating location with models to infer the heating mechanism and parameters. An example of this can be found in Priest et al. (2000), who tried to infer the heating mechanism from the (observed) temperature profile and concluded that heating is not localised at the footpoint. More recently, Van Doorsselaere et al. (2007) used a comparison between models and observations and suggested that heating localised near the footpoints was more likely to be caused by a resistive wave-heating mechanism (rather than viscous). Karampelas et al. (2017) confirmed this result, noting, though, that the development of Kelvin-Helmholtz instabilities can displace the heating location closer to the apex of the loop.

Finally, we point out that a key limitation of this work is that we did not adopt a realistic loop structure in terms of gravitational stratification and a proper transition region at the footpoint. In particular, gravitational stratification could imply lower densities near the apex, and hence smaller radiative losses, but at the same time might affect the efficiency of the wave heating, as the boundary-shell gradients become less steep. However, if the heating is concentrated at the footpoints, it might be able to counteract the smaller radiative losses from a stratified loop.

Finally, a key result from this paper is that transverse MHD waves can expand and sustain a boundary shell around or within loops. This is an important aspect of the debate on the role of waves as a heating mechanism, as transitory density enhancements are common features in the corona, as a consequence of, for example, spicules; plasma evaporation; coronal rain; or transverse displacements. As footpoint motions continuously drive transverse waves, this scenario considerably facilitates the formation of small-scale structures allowing the conversion of wave energy into thermal energy. Even if the mechanism is not a direct dissipation of MHD waves by phase mixing, the transverse waves clearly play a key role in the generation of small-scale structures (see also Howson et al. 2020).

Movies

Movie of Fig. 4 Access here

Movie of Fig. 8 Access here

Acknowledgments

The authors would like to thank Tom Van Doorsselaere for the helpful discussions on uniturbulence. This work has received support from the UK Science and Technology Facilities Council (Consolidated Grant ST/K000950/1), the European Union Horizon 2020 research and innovation programme (grant agreement No. 647214) and the Research Council of Norway through its Centres of Excellence scheme, project number 262622. R. J. Morton is grateful for support from the UKRI Future Leader Fellowship (RiPSAW – MR/T019891/1) and STFC (ST/T000384/1). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC@Durham equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1 and ST/R002371/1, Durham University and STFC operations grant ST/R000832/1. The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e- Infrastructure. We acknowledge the use of the open source (gitorious.org/amrvac) MPI-AMRVAC software, relying on coding efforts from C. Xia, O. Porth, R. Keppens.

References

All Tables

All Figures

	Fig. 1. Isosurfaces of density (top row) and Alfvén speed gradient (bottom row) for three different loop configurations: uniform (left-hand side), z0 = 130 Mm (centre), and z0 = 46 Mm (right-hand side).
In the text

	Fig. 2. Profile of density contrast ρc (blue curves) and interior region radius b (red curves), as a function of z for two density profiles of the loop with z0 = 46 Mm (dashed curves) and z0 = 130 Mm (continuous curves).
In the text

Fig. 3. a: spectra for transverse waves used in this study. The black curve represents the averaged observed spectrum, the red curve represents a spectrum where the power enhancement at ν0 = 4 mHz is removed, and the blue curve represents a spectrum where we used a shallower power-law index. b: velocity vx time series for the three spectra in normalised units with the same coloured legend. c: displacement associated with the velocity time series in normalised units and with the same colour legend. The black dashed curve represents the displacement derived from an additional velocity profile obtained from the same observed power spectrum.

In the text

	Fig. 4. Density and electric current density isosurfaces at two different times (t = 617 s and t = 2224 s) for the simulation with a uniform density structure of the loop. A movie of the density evolution is available online.
In the text

	Fig. 5. Density cross-sections at z = 79.6 Mm at three different times (t = 0, t = 741.4 s, and t = 2241 s). Over-plotted are velocity vectors where |v| > 5 km s−1, and contours of the gradient of the Alfvén speed ∇VA = 0.5 s−1 (green contour) and ∇VA = 1 s−1 (blue contour).
In the text

	Fig. 6. Map of expansion of boundary shell as a function of z and t for the simulation with a uniform loop (upper panel). Evolution of the boundary-shell efficiency in the same simulation (lower panel).
In the text

	Fig. 7. Cross-section of the gradient of the Alfvén speed at y = 0, at the initial condition t = 0 s (a), and at t = 1977 s (b). Evolution of the boundary-shell consistency for the simulation with uniform loop.
In the text

	Fig. 8. Isosurface of ohmic heating at H = 108 erg cm−2 s−1 at two different times (t = 617 s and t = 2224 s) in the simulation with a uniform loop. A movie of the heating evolution is available online.
In the text

	Fig. 9. Map of average heating in boundary shell along the loop as a function of time for the simulation with a uniform loop. Although the heating is generally within the range shown here, the maximum heating is 6 × 109 erg s−1 cm2.
In the text

	Fig. 10. Average heating in boundary shell as a function of time for the simulation with a uniform loop, at three different coordinates along the loop.
In the text

	Fig. 11. Map of ratio of average heating and average radiative losses in the boundary shell as functions of loop position and time, for the simulation with a uniform loop.
In the text

	Fig. 12. Isosurfaces of Alfvén speed gradient at t = 2224 s for the simulations with z0 = 130 Mm (left-hand side) and z0 = 46 Mm (right-hand side).
In the text

	Fig. 13. Map of boundary-shell expansion as a function of z and t (upper row) and evolution of the boundary-shell efficiency (lower row). for the simulations with z0 = 130 Mm (left-hand side) and z0 = 46 Mm (right-hand side).
In the text

	Fig. 14. Evolution of boundary-shell consistency for the simulation with uniform loop, with z0 = 130 Mm, and z0 = 46 Mm.
In the text

	Fig. 15. Maps of average-heating boundary-shell logarithm in the simulation with (a) z0 = 130 Mm or (b) z0 = 46 Mm normalised to the reference simulation and plotted as a function of position along the loop and time.
In the text

	Fig. 16. Maps of average-heating boundary-shell logarithm in the simulation with (a) a closed loop, or (b) a wider driver normalised to the reference simulation and plotted as a function of position along the loop and time.
In the text

	Fig. 17. Isosurface of Alfvén speed gradient for the reference simulation (left-hand side) and for the wider driver simulation (right-hand side) at t = 988 s, coloured by the modulus of the component of the Poynting flux across the isosurface.
In the text

	Fig. 18. Map of average heating boundary-shell logarithm in the simulation with a stronger driver, normalised to the reference simulation and plotted as a function of position along the loop and time.
In the text

	Fig. 19. Maps of average-heating boundary-shell logarithm in the simulation with (a) η = 1010 ηS and (b) η = 1010 ηS and j0 = 0, normalised to the reference simulation and plotted as a function of position along the loop and time.
In the text

	Fig. 20. Maps of average wave energy in the boundary shell as a function of time and position along the loop for (top) the reference simulation and (bottom) the simulation with η0 = 1010 ηS and j0 = 0.
In the text

	Fig. 21. a: average heating in boundary shell as a function of time for a number of simulations listed in the colour legend, with the average heating enhancement reported. b: average heating in the boundary shell as a function of position along the loop for the same simulations normalised to the average heating at z = 0.
In the text

	Fig. 22. a: average heating in boundary shell as a function of time for the simulations with four different time series, as explained in Fig. 3b. b: average heating in the boundary shell as a function of the position along the loop for the four different time series, as explained in Fig. 3b.
In the text

	Fig. 23. Scattered plot of total heating as a function of the total Poynting flux injected in the boundary shell. Both quantities are shown with respect to the reference simulation. Colours and symbols are clarified in the plot.
In the text