Heating and Cooling in Transversely Oscillating Coronal Loops Powered by Broadband, Multi-Directional Wave Drivers

Abstract

Recent studies have identified the potential for coronal wave heating to balance radiative losses in a transversely oscillating low-density loop undergoing resonant absorption, phase mixing and the Kelvin–Helmholtz instability. This result relied on a continuous, resonant oscillatory driver acting on one of the loop footpoints and similar setups with non-resonant driving produce insufficient heating. Here, we consider broadband and multi-directional drivers with power in both resonant and non-resonant frequencies. Using three-dimensional magnetohydrodynamic simulations, we impose transverse, continuous velocity drivers at the footpoints of a coronal loop, which is dense in comparison to the background plasma. We include the effects of optically thin radiation and a uniform background heating term that maintains the temperature of the external plasma but is insufficient to balance energy losses within the loop. For both broadband and multi-directional drivers, we find that the energy dissipation rates are sufficient to balance the average energy losses throughout the simulation volume. Resonant components of the wave driver efficiently inject energy into the system and these frequencies dominate the energetics. Although the mean radiative losses are balanced, the loop core cools in all cases as the wave heating rates are locally insufficient, despite the relatively low density considered here.

1. Introduction

Contemporary observations reveal the proliferation of oscillatory behaviour throughout the coronal volume (e.g., [1,2,3,4,5,6,7,8] and reviews [9,10,11,12]). These waves can be associated with a significant energy flux (e.g., [13]) and the dissipation of this energy has been widely proposed as a mechanism for maintaining the high temperatures of the solar corona (e.g., [14,15,16]. However, in typical corona conditions, the rate of energy dissipation through resistive or viscous effects is expected to be low. Therefore, in order to enhance these rates, energy needs to be transferred from large-scale waves to much smaller spatial scales, where energy is dissipated more readily.

Over the previous decades, a variety of physical processes have been proposed for generating the required small scales from the observed large-scale oscillations. These include resonant absorption/mode coupling [17] and phase mixing [18], or both effects combined (e.g., [19,20]. These mechanisms develop as a result of cross-field gradients in the local Alfvén speed, which is typically assumed to be associated with density gradients in the boundary of coronal loops (between the dense interior and tenuous external plasma). In non-linear regimes, the combination of these processes can result in the oscillating plasma becoming unstable to the Kelvin–Helmholtz instability (KHI). The general paradigm is as follows.

1.

Energy in large-scale transverse standing waves is transferred to azimuthal Alfvén waves on resonant field lines in the boundary of the oscillating loop.

2.

Cross-field gradients in the local Alfvén speed and therefore the natural Alfvén frequency allows phase mixing to progress in the loop boundary, generating increasingly small length scales.

3.

Together, these processes induce radial gradients in the azimuthal velocity profile. This velocity shear can be unstable to the magnetic KHI.

4.

In unstable regimes, the KHI initially develops by generating characteristic vortices which deform the density profile of the oscillating loop.

5.

These vortices are also unstable, driving the development of yet smaller vortices, and inducing a cascade of energy to small scales.

6.

In non-ideal plasmas, the dissipation length scale is inevitably reached and energy is converted to heat. Thus, the instability development can enhance the rate of wave heating in non-ideal plasmas.

This process is reviewed in Ref. [21] and has been studied in numerical simulations (e.g., [22,23,24,25,26,27,28,29,30,31,32,33,34]) and in analytical settings (e.g., [35,36,37]). The development of the instability is robust for a variety of setups, including in the presence of gravitational stratification [38], optically-thin radiation [39,40,41], for purely torsional waves [42] and in multi-stranded structures [43]. Despite this, the instability growth rate can be reduced by a variety of factors, including twisted or sheared magnetic field [27,31,36,44], short wavelengths [44] and high dissipation and/or low numerical resolution [26]. Wave driven turbulence can also develop for propagating waves in inhomogeneous media, again leading to enhanced energy dissipation rates [45,46,47,48,49].

Using observed wave amplitudes, the authors of Ref. [50] argued that the instability is unlikely to generate sufficient heating to balance expected radiative losses; there is simply insufficient wave energy. Indeed, more generally, numerical studies of resonant absorption and phase mixing have shown that the resultant wave heating is typically insufficient to balance radiative losses (e.g., [51]). This is even the case with continuous multi-frequency and/or multi-directional driving, which better reflect wave-excitation mechanisms in the solar corona [52,53,54]. Despite this, in Ref. [40], it is shown that for a continuously driven system, this process was able to balance energy losses from optically thin radiative losses within a low-density oscillating coronal loop. However, there are several important points which still need to be addressed before this mechanism is accepted as having an important role in coronal heating. For example, the frequency of the wave driver used in the study exactly matched the natural kink frequency (e.g., [10,55]) of the system. As such, the oscillatory forcing was resonant and energy was injected into the loop as efficiently as possible. Secondly, the density used was fairly low and whilst perhaps reflecting quiet Sun conditions, it is not representative of the much denser coronal loops observed in active regions. The density is very important in this context as radiative effects (and hence energy losses) scale with the density squared. In a subsequent study, presented in Ref. [41], it is determined that for either non-resonant drivers or for denser loops, the KHI was unable to balance energy losses. In addition, neither of these studies considered the effects of thermal conduction, which would further enhance energy losses and reduce the likelihood of this mechanism being able to maintain coronal conditions.

In this paper, we investigate the effects of the driver frequency and polarisation in more detail. Whilst the non-resonant driving did not provide sufficient energy flux to balance radiative losses in Ref. [41], it is also unrepresentative of the broadband nature of coronal waves (e.g., [56,57,58]). Although most of the power in oscillatory drivers will inevitably be associated with non-resonant frequencies, even a small component of wave power at frequencies that coincide with the natural harmonics of the system will allow resonances to be excited and thus enhance energy injection rates. As such, to ascertain the suitability of this heating process in more realistic conditions, we investigate the effects of broadband and multi-directional drivers on energy injection, current formation and dissipation rates. In particular, we investigate whether broadband and multi-directional drivers power wave heating rates that can balance optically thin radiative losses within a coronal loop.

2. Methods and Basic Equations

The simulations presented within this paper, employed the numerical code, Lare3d [59], which advances the fully-resistive, magnetohydrodynamics (MHD) equations in normalised form. The equations are given by:

$$\frac{\mathrm{D}\rho }{\mathrm{D}t}=-\rho \nabla ·\mathbf{v},$$

(1)

$$\rho \frac{\mathrm{D}\mathbf{v}}{\mathrm{D}t}=\mathbf{j}\times \mathbf{B}-\nabla P-\rho \mathbf{g}+{\mathbf{F}}_{\mathrm{visc}},$$

(2)

$$\rho \frac{\mathrm{D}ϵ}{\mathrm{D}t}=-P(\nabla ·\mathbf{v})-{\rho }^2\mathsf{\Lambda }\left(T\right)+\eta {\left|\mathbf{j}\right|}^2+Q_{\mathrm{visc}}+Q_{\mathrm{bg}},$$

(3)

$$\frac{\mathrm{D}\mathbf{B}}{\mathrm{D}t}=\left(\mathbf{B}·\nabla \right)\mathbf{v}-\left(\nabla ·\mathbf{v}\right)\mathbf{B}-\nabla \times \left(\eta \nabla \times \mathbf{B}\right),$$

(4)

$$P=2k_{B}nT.$$

(5)

where t denotes the time, $\rho$ is the plasma density, $\mathbf{v}$ is the velocity, $\mathbf{j}$ is the current density, $\mathbf{B}$ is the magnetic field, P is the gas pressure, $ϵ$ is the specific internal energy, n is the number density, T is the temperature, $\eta$ is the magnetic diffusivity, and $k_B$ is the Boltzmann constant. Here, $\mathrm{D}/\mathrm{D}t$ refers to the total time derivative.The effects of optically thin radiation, $\mathsf{\Lambda }\left(T\right)$, and a uniform background heating, $Q_{\mathrm{bg}}$, are included in the energy Equation (3). The optically thin radiative loss function, $\mathsf{\Lambda }\left(T\right)$, is described in Ref. [60] and the heating term is described in detail below.

Equations (1)–(5) above include the effects of the magnetic diffusivity and shock viscosity terms. The magnetic diffusivity is approximately $2\times {10}^8{\mathrm{m}}^2{\mathrm{s}}^{-1}$. Due to numerical constraints, this is necessarily orders of magnitude larger than the expected value within a coronal plasma. The shock viscosities appear in the momentum Equation (2) in the form of a force, ${\mathbf{F}}_{\mathrm{visc}}$, and in the energy Equation (3) as a heating term, $Q_{\mathrm{visc}}$. These viscosity terms help ensure numerical stability and are described in Ref. [61]. Let us note that as the KHI is sensitive to the transport coefficients [26], the wave dynamics and heating rates are strongly affected by these values.

2.1. Initial Conditions

Following the setups of Refs. [40,41], our model represents a dense, three-dimensional coronal loop embedded in a magnetic field aligned with the z axis. The initial density is defined by

$$\rho (x,y)={\rho }_e+\left(\frac{{\rho }_i+{\rho }_e}{2}\right)\left\{1+tanh\left(\frac{r-a}{b}\right)\right\},$$

(6)

where $r=\sqrt{x^2+y^2}$ and ${\rho }_e$ and ${\rho }_i=3{\rho }_e$ represent the external and internal density, respectively. In this setup, ${\rho }_e=1.67\times {10}^{-13}\mathrm{kg}{\mathrm{m}}^{-3}$. This matches the initial density profile imposed in Ref. [40] and the low-density case in Ref. [41]. The constants a and b define the loop radius and the width of the boundary between the internal and external regions, respectively. The constants a and b are fixed such that the radius is 1 Mm and the thickness of the boundary is approximately 0.4 Mm. The loop has length 200 Mm and is initially invariant along its length (aligned with the z direction). The temperature throughout the domain is initially uniform and is set to be 1 MK. Therefore, in order to ensure an initial total pressure balance (for an equilibrium state), the magnetic field strength is reduced within the density-enhanced loop. As the plasma-$\beta \ll 1$, this dip in the magnetic field strength is only approximately $0.7\%$ of the background field strength (30 G). Figure 1 shows the initial density (solid line, left-hand axis) and the initial $B_z$ profile (dashed line, right-hand axis).

Due to the inclusion of optically thin radiation, the plasma will cool during the simulation (unless heated at a sufficiently high rate). Therefore, a background heating term, $Q_{\mathrm{bg}}$, is included in the energy Equation (3). This is chosen to exactly balance the initial radiative losses in the external plasma. However, even though this term is included uniformly throughout the domain, it will not balance the energy losses in the dense loop. This is because the optically thin losses scale with the square of the plasma density. As such, the cooling is nine times larger inside the loop than in the external plasma and the background heating cannot balance this enhanced radiation. This approach allows us to determine whether wave heating can balance the majority of the energy losses in the flux tube and was adopted in Refs. [40,41].

The numerical volume had dimensions: $-4\le x,y\le 4$ Mm and $0\le z\le 200$ Mm, and consists of 512 × 512 × 200 grid points in the x, y, and z directions, respectively. In this setup, the effects of gravitational stratification, loop curvature and thermal conduction are neglected.

2.2. Boundary Conditions

The fundamental kink frequency, ${\omega }_k$, of the system (e.g., [55,62,63]) was calculated using Equation (8) in Ref. [10]. This gives a resonant period of approximately 86 s. In other words, a transverse wave driver with this period will resonantly excite the fundamental kink oscillation. For these simulations, we consider drivers of the form,

$$v_x\left(t\right)=\sum _{i=1}^N{v}_{0}sin\left({\omega }_{i}t+{\psi }_i\right)cos{\varphi }_i,$$

(7)

$$v_y\left(t\right)=\sum _{i=1}^N{v}_{0}sin\left({\omega }_{i}t+{\psi }_i\right)sin{\varphi }_i.$$

(8)

Here, the sum is taken over N single harmonic wave components to generate a multi-frequency driver. The component frequencies, ${\omega }_i$, are defined by

$${\omega }_i=\frac{{\omega }_k}{2}+\frac{i-1}{N-1}\left(\frac{3{\omega }_k}{2}\right).$$

(9)

The component frequencies, ${\omega }_i$, vary uniformly between ${\omega }_k/2$ and $2{\omega }_k$. This means we excite wave modes with periods between 43 s and 172 s. For each component, ${\psi }_i$ represents a random phase shift, such that the wave components are not all in phase and ${\varphi }_i$ represents the polarisation angle of the component. Two sets of simulations are considered: ${\varphi }_i=\pi /2$ (for all i), and ${\varphi }_i$ is a randomly selected from the interval $[0,2\pi ]$. In the former case, the wave driver oscillates in the y direction ($v_x=0$). For the latter case, the wave polarisation shifts during the simulation. The amplitude of each component is fixed ($v_0$), as such, the power spectrum for these waves is flat (akin to white noise) over the range of frequencies being considered. Whilst this does not necessarily represent waves in the solar corona (e.g., Figure 6 in Ref. [57]), it permits an initial analysis of the effects of power in a range of frequencies on this system. For these simulations we choose $N=20$ frequency bins. We also compare the effects of these wave drivers with the single-polarisation, resonantly-driven case presented in Ref. [41]. The fixed amplitude, $v_0$, is selected such that the mean of the magnitude of the velocity is the same across all simulations, i.e., the drivers all have the same average kinetic energy.

On the lower z (driven) boundary, $v_z=0$ and all other variables have zero-gradients. The upper z boundary stipulates $\mathbf{v}=0$, and uses a zero-gradient condition for all other variables. As such, it is a reflecting boundary. The x and y boundaries are periodic.

Since no flow is permitted into or out of the computational domain (through the upper and lower z boundaries), the only non-zero energy flux through the boundaries is the Poynting flux. This is sensitive to both the transverse driver and how the magnetic field responds during the simulations. As the upper z boundary acts as a perfect reflector, continuous sinusoidal driving can excite MHD wave resonances if the driver frequency coincides with a natural harmonic of the system. As such, the interaction between the boundary driver and the reflected waves can significantly change the energy injection rates for small changes in the frequency of the driver (e.g., from resonant to non-resonant driving). An upper, perfectly reflecting boundary is frequently used for modelling a closed coronal loop, however, is not applicable for open field lines that require a careful treatment for the upper open boundary condition (e.g., [64,65]). In such cases, reflections can still be generated due to density gradients along the field, however, these tend to be much less efficient (e.g., [66,67]). The specific effects of our imposed drivers on the injected Poynting flux are discussed in Section 3.1.

3. Results

Four simulations that investigate the effects of different imposed wave drivers are considered here, namely:

1.

ResUni: resonantly driven $(\omega ={\omega }_k)$, single polarisation $(\varphi =\pi /2)$;

2.

BbUni: broadband (multi-frequency, as defined above), single polarisation $(\varphi =\pi /2)$;

3.

ResMulti: resonantly driven $(\omega ={\omega }_k)$, changing polarisation;

4.

BbMulti: multi-frequency and changing polarisation.

For brevity, in what follows, the acronyms listed above are used.

Figure 2 shows the form of $v_y$ as a function of time for the ResUni (black curve) and BbUni (red curve) simulations.

Let us start by considering the evolution of the density profile in each of these cases. As the instability develops, it increasingly deforms the loop density profile, initially by forming characteristic Kelvin–Helmholtz (K-H) vortices, and then progressively driving turbulent-like flows in the loop cross-section. The rate and extent of this deformation is indicative of the instability growth rate in a particular system. Figure 3 shows the evolution of the density cross-section at the loop-apex ($z=100$ Mm) for each of the simulations. Each row corresponds to a different simulation and each column to a different time. Let us note that only a limited section of the horizontal plane is shown to display the density deformation in greater detail.

For the broadband and multi-directional drivers, the growth rate of the instability is reduced. This is shown by the later formation of K-H vortices and the reduced deformation in the core of the loop towards the end of the simulation (final column). The reduced growth rate occurs for several reasons, including reduced wave amplitudes (due to less power in the resonant frequencies), less coherent azimuthal Alfvén waves (particularly for multi-directional drivers), the presence of propagating modes (due to non-harmonic frequencies) and smaller wavelengths (due to the higher frequency component of the driving). These factors (partially) stabilise the system by reducing the radial velocity shear in the loop boundary and/or enhancing the stabilising magnetic tension force that opposes the bending of field lines and hence the formation of K-H vortices. The shorter wavelengths can also increase the azimuthal wavenumber of the vortices that form (e.g., [44]). This is apparent by comparing the size of the emerging vortices in the first two panels in the top row (ResUni) with those in the second panel of the second row (BbUni) and the last two panels in the final row (BbMulti).

In the simulations with a single oscillation direction (ResUni and BbUni), the loop is only displaced in the y direction. For the ResMulti and BbMulti simulations, on the other hand, the loop apex can also move in the x direction (see, e.g., Figure 3, bottom, right). This has important consequences for the evolution of the azimuthal Alfvén wave in the loop boundary (excited by resonant absorption/mode coupling) and hence the location of the instability. For the unidirectional drivers (ResUni and BbUni), the azimuthal Alfvén wave forms tangentially to the constant direction of the wave driver, with the largest amplitudes developing where the lines $x=\pm 1$ Mm intersect the loop apex. As such, the greatest velocity shear forms on the left and right sides of the loop (as viewed in Figure 3) and hence the KHI forms here initially. For multi-directional drivers, however, this is no longer a preferred location and vortices can form all around the loop perimeter. To see this in Figure 3, compare the vortex locations in the upper two rows with those in the lower two rows.

Several of the factors that reduce the instability growth rate can be observed by tracking the evolution of the velocity profile in each simulation. To this end, in Figure 4, we show $v_y$ for the ResUni and BbUni simulations. The vertical cuts are taken parallel to the x axis through the centre of the oscillating flux tubes and the y component of the velocity shows flows perpendicular to this plane. Figure 4 shows a time ($t=300$ s) when both flux tubes are approximately straight so this plane cuts through the central axis of each loop. This is after several periods of wave driving but prior to the development of the instability.

In Figure 4, upper, one can see that a fundamental standing mode is excited in the ResUni loop with nodes at the two z boundaries and an antinode at $z=100$ Mm. The positive velocity (red) in the core of the flux tube represents the transverse (kink) oscillation and the two blue bands (along $x\approx \pm 1$ Mm) are associated with the azimuthal Alfvén wave. One can see that these are out-of-phase with the velocity in the loop core. Consequently, there is a large velocity shear across the cross-section of the flux tube (e.g., along the line $z=100$ Mm). Whilst the amplitude of this shear changes along the length of the flux tube (due to lower velocities away from the wave antinode), this is coherent along the length of the flux tube. As the transverse and Alfvén modes continue to oscillate out-of phase, the velocity shear also oscillates but maintains its spatial structure. As such, conditions at the loop apex are favourable for the instability to develop.

In Figure 4, lower, one sees no evidence of a fundamental mode at this time. Instead, the reduced wavelength within the loop core shows the existence of higher frequency waves. Whilst there is still out-of-phase behaviour between the loop exterior and interior, this is less coherent along the loop length, and the Alfvén wave bands have also not formed consistently along the loop length. As no standing mode has been excited (or in any case is masked by the magnitude of the propagating mode), conditions at the loop apex are now much less favourable (there is no consistent velocity shear) for the onset of the KHI. As a result, it takes longer for the instability to develop. The other two simulations (ResMulti and BbMulti) also exhibit similar incoherence in the radial velocity shear and hence lower growth rates.

As the magnetic field is essentially frozen-in to the coronal plasma, the K-H vortices deform and twist the embedded magnetic field. As such, gradients in the magnetic field develop, enhancing Ohmic heating rates and driving magnetic reconnection [44]. At the same time, the velocity shear is associated with large velocity gradients, enhancing the dissipation of energy through viscous effects. As a result, the instability is interesting in the context of the coronal heating problem as it can significantly enhance wave heating rates. For the current study, it is of interest whether these enhanced heating rates can balance the optically-thin radiative losses emanating from the dense loop. Figure 5 shows how the KHI-driven heating affects the apex temperature profile in each of the four simulations. One immediately sees the close relationship between the temperature profiles and the density cuts in the final column of Figure 3. In all cases, a modest increase in temperature in the loop boundary (from approximately 1 MK to 1.5 MK) and a decrease (to approximately 0.5 MK) in the loop core is visible. There are two pertinent factors which cause this behaviour. Firstly, in all simulations, the effects of the KHI are greatest in the loop boundary and thus we expect the highest heating rates here. Secondly, the loop boundary has lower densities than the loop core and thus the radiative cooling is significantly lower. Together, these result in enhanced temperatures forming around the K-H vortices. Let us also note that simulations with greater density deformation (e.g., ResUni case in Figure 5, upper, left), produce heating over a greater area, reflecting both enhanced heating and lower cooling rates across the loop cross-section.

In these simulations, adiabatic effects (due to the compression and rarefaction of plasma) will also modify the temperature. In order to verify that the temperature change observed in Figure 5 is not dominated by these reversible effects, we compare the magnitude of the adiabatic term, $-P(\nabla ·\mathbf{v})$, with the left hand side of the energy Equation (3). Figure 6 shows the result of this comparison along the $y=0$ Mm, $z=100$ Mm line in the ResUni simulation. The cumulative contribution of the adiabatic term (blue line) and the $\rho \mathrm{D}ϵ/\mathrm{D}t$ term (red line) over the duration of the simulation are shown. One can see that the contribution of the adiabatic effects is small (in comparison to the red curve) and so can conclude that the temperature change in the simulation is dominated by irreversible energy dissipation and radiation (not by adiabatic effects). Let us note that the peaks in the red curve correspond to heating in the boundary of the loop and the central dip demonstrates radiative cooling in the dense loop core.

The heating in the simulations considered is driven by resistivity acting on gradients in the magnetic field (Ohmic heating) and viscosity acting on gradients in the velocity field (viscous heating). The relative sizes of these heating components depend on the specific transport coefficients and this, in turn, determines the spatial distribution of the heating. For a fundamental standing mode, the perturbed velocity field is largest at the loop apex and the perturbed magnetic field is largest at the loop foot points (e.g., [68]). As such, if viscous heating dominates, the temperature increase will be greatest at the loop apex, and conversely, if Ohmic heating dominates, the largest temperature increases will be located near the magnetic foot points. However, when there are additional frequencies (and hence additional wavelengths) present, this picture is less clear-cut.

Figure 7 shows how the spatial distribution of gradients in the magnetic and velocity fields differ between the ResUni (Figure 7, upper) and BbMulti (Figure 7, lower) simulations. In particular, we consider the current density, $\mathbf{j}$, and the vorticity, $\mathbf{\omega }=\nabla \times \mathbf{v}$, to show gradients in the magnetic and velocity fields, respectively. In order to include both of these values within the same panel, the $\left|\mathbf{\omega }\right|-\left|\mathbf{j}\right|$ in the $y=0$ Mm plane is shown. As such, positive values (shown in red/pink) indicate regions where the vorticity is larger than the current, and negative values (blue) show the opposite.

Figure 7, upper (ResUni), shows the classic features of a fundamental standing kink/Alfvén mode, with the largest currents at the foot points and the largest vorticities at the loop apex. One can also note that as the velocity and magnetic field perturbations are out-of-phase with each other, the vorticity and currents are also out-of-phase. For example, values of x that exhibit large currents (e.g., $x=\pm 0.75$ Mm) show relatively small vorticities, and vice versa. The narrow banding, which is particularly apparent in the high-vorticity region, is indicative of phase-mixing in the loop boundary. On the other hand, for Figure 7, lower (BbMulti), evidence of the fundamental harmonic is much less pronounced. For example, the largest vorticities are displaced from the apex towards the $z=0$ Mm foot point. This is indicative of the propagating modes associated with the non-harmonic frequencies in the wave driver. Despite this, one can see that the largest gradients in the velocity field still form around the boundary of the flux tube, where the azimuthal Alfvén wave is excited. Additionally, the general prevalence of large currents near the foot points and large vorticities near the apex remains. This is because the fundamental mode is still being resonantly driven by one component of the imposed driver, and thus its characteristics are still apparent.

Figure 8 quantifies how the current density (Figure 8, left) and vorticity (Figure 8, right) is distributed along the z direction. In particular, the mean current density and vorticity are shown as a function of z. The quantities are averaged over x, y, and time, t, and for each plot, the two extreme cases are shown: ResUni (black curves) and BbMulti (blue curves). The curves are normalised by their corresponding maxima. By comparing the blue and black curves, one can see that the distribution of both the current density and the vorticity are modified. Although currents and vorticities are always larger at the loop foot points and loop apex, respectively, the higher frequencies in the BbMulti simulation result in the spatial gradients being distributed more evenly through the simulation volume. As a result, the spatial distribution of heating is modified between these simulations.

3.1. Energetics

The imposed wave drivers inject energy through the Poynting flux. For an initially uniform magnetic field, $\mathbf{B}=B_0\widehat{\mathbf{z}}$, driven by a simple harmonic oscillating driver, $\mathbf{v}=v_{0}sin\omega t\widehat{\mathbf{y}}$, acting at magnetic foot points (at $z=0$), the Poynting flux, S, can be expressed as

$$S=\frac{-1}{{\mu }_0}{\int }_A\mathbf{E}\times \mathbf{B}·\mathrm{d}\mathbf{S}=\frac{1}{{\mu }_0}{B}_0{\int }_A{v}_y{B}_y\mathrm{d}x\mathrm{d}y.$$

(10)

Here, the ideal plasma ($\eta =0$) is assumed and the integral is calculated over the area on which the driver acts. If no wave reflections interact with the driven boundary (e.g., for an open field or before any waves have returned following reflection from a distant foot point), the velocity and magnetic fields satisfy:

$$v_y=v_{0}sin\left\{\frac{\omega }{v_A}(v_{A}t-z)\right\},$$

(11)

$$B_y=\frac{v_0{B}_0}{v_A}sin\left\{\frac{\omega }{v_A}(z-v_{A}t)\right\}.$$

(12)

Here, $v_A$ is the Alfvén speed. By integrating Equation (10) over time, one finds that the total energy injected through the $z=0$ plane between $t=0$ and a later time, $t=\tau$, is given by

$$\frac{B_0^2{v}_0^{2}A}{v_A}{\int }_{t=0}^{t=\tau }{sin}^2\omega t\mathrm{d}t=\frac{B_0^2{v}_0^{2}A}{v_A}\left(\frac{\tau }{2}-\frac{sin2\omega \tau }{4\omega }\right),$$

(13)

where A is the area, over which the driver acts (e.g., the lower z boundary in our simulations). One can see that the cumulative energy injection is given by a linear growth term and an oscillatory modification. As the linear term is not a function of the frequency, over long times (such that the sinusoidal term is negligible), the injected energy for a purely propagating wave mode does not depend on the driver frequency.

However, in the case considered here, reflected waves from the upper z boundary soon interact with the driven boundary. In particular, whilst we still control the imposed velocity driver ($v_y$ in Equation (10)), the magnetic field (e.g., $b_y$ in Equation (10)) at the boundary is dependent on the evolution of the system. For resonant wave frequencies, the reflected waves guarantee an increase in the Poynting flux injected into the system. However, for non-resonant frequencies, over long times, the wave driver will remove as much energy from the system as it injects. As such, for broadband oscillatory drivers, it is the resonant frequencies (for the fundamental mode and higher harmonics) that will dominate the energetics of the system.

For the BbUni simulation, one can decompose the imposed velocity driver into its individual frequency components (see Equation (8)). From this decomposition, one can then consider how much each frequency contributes to the overall energy injection rate (e.g., using Equation (10)). Figure 9, left, shows the contribution to the total energy injection rate (over the duration of the simulation) from each component of the wave driver. The curve is normalised to its maximum. One can see that the energy injection peaks at the kink frequency (${\omega }_k$, the resonant fundamental mode for the flux tube) and has a secondary peak close to 1.4 ${\omega }_k$ (corresponding to the fundamental frequency of standing Alfvén waves in the external plasma). As the kink mode is evanescent outside the flux tube (e.g., [10,55]), at large distances from the loop, standing Alfvén modes are resonantly excited instead. One can also note that there is another peak close to 2 ${\omega }_k$, suggesting that the highest frequency modes in the driving are exciting the second harmonic.

Figure 9, right, shows the spatial distribution of the total Poynting flux (over the simulation duration) associated with the $\omega ={\omega }_k$ component of the driving. In agreement with Figure 10 of Ref. [41] (which showed a similar plot for the ResUni case), one can see that kink mode energy is preferentially injected in the core of the loop and extracted in the loop boundaries. This is because the excited azimuthal Alfvén waves are out-of-phase with the oscillation in the loop core (despite having the same frequency). As such, the velocity driver consistently opposes the Alfvén wave oscillation and thus removes energy from the system.

In Ref. [41], we showed that non-resonant, single frequency drivers were not able to inject and dissipate sufficient energy to balance the optically thin radiative losses. However, for the our broadband cases here, whilst much of the driver power is at non-resonant frequencies, there is still power at resonant frequencies (e.g., the fundamental kink mode, $\omega ={\omega }_k$, and the second harmonic, $\omega =2{\omega }_k$). In order to examine the energy injected by the drivers considered, in Figure 10, the evolution of the volume integrated wave energy (Figure 10, left) and the change in the specific internal energy (Figure 10, right) are shown. For the wave energy, we calculate the total kinetic energy plus the total magnetic energy minus the initial magnetic energy. For both plots, the change in energy is shown relative to the initial total energy (kinetic, magnetic and thermal). Let us note that, in both cases, the energy change is modest (≪0.1%), in comparison with the initial energy in the system.

There are several important factors that affect the amount of wave energy in the simulations at any given time. For the ResUni case, although the energy injection rate is initially efficient, it does not exhibit the largest wave energies. This is because:

• the other wave drivers also excite resonant modes, ensuring they inject energy efficiently;
• in the ResUni case, the KHI quickly (within a few wave periods) deforms the loop density profile and hence detunes the system from this single frequency driver, reducing subsequent energy injection rates;
• the earlier formation of the KHI results in wave energy dissipating sooner in the ResUni case;
• other wave drivers can be resonant for field lines that are distant from the loop, where the local Alfvén frequency is different from the kink frequency; as such, they can inject more energy into the external plasma;
• the multi-directional drivers excite waves in the x direction too; this means that the system has a greater capacity for wave energy than the single direction cases.

This final point is important in establishing why the ResMulti (orange curve) and BbMulti (blue curve) simulations exhibit more wave energy than the uni-directional driver cases (red and black curves).

Figure 10, right, shows how the volume-integrated thermal energy changes during the four simulations. In all cases, the final thermal energy is similar or greater than the initial thermal energy (with just a slight decrease in the BbUni simulation). This suggests that generally these wave drivers are able to balance the radiative losses from the flux tube. Despite the greater wave energies in the ResMulti and BbMulti cases, for most of the duration, the ResUni simulation has the largest thermal energy. There are two main reasons for this, both associated with the KHI. First, the relatively early onset in the ResUni case, means that it exhibits the greatest heating rates across all simulations from $t=500$ s for approximately 1000 s (until heating rates are enhanced by the instability in other cases). Second, the instability reduces the volume integral of ${\rho }^2$ in the computational domain. As such, the radiative losses decrease, reducing the mean cooling rate.

These two effects are well seen in Figure 11. Figure 11, left, shows the decrease in the volume integral of ${\rho }^2$ across the four simulations. All curves are shown normalised to the corresponding initial value. One can see that all simulations show a drop in the integral as density is redistributed throughout the domain. We note that density is conserved in all cases and these curves reflect plasma moving from high-density regions to low-density regions, particularly as the instability develops. The drop is greatest and begins soonest in the ResUni simulation, reducing the radiative cooling in comparison to other simulations. Although the drop is reasonably small (∼2%), this integral is computed over the entire domain, and the reduction in the radiative losses within the dense flux tube is larger.

Figure 11, right, shows the cumulative heating (total viscous and Ohmic heating) that has occurred up to a given time within each simulation. The time derivative gives the instantaneous heating rate. Although the heating rate is initially largest in the ResUni simulation, subsequently, the approximately equal gradients show that the heating rate is similar across all simulations. At late times, the BbMulti simulation exhibits enhanced heating rates as the KHI begins to dissipate the enhanced wave energies (see Figure 10) effectively. It is only the combined effects of both quantities, shown in Figure 11, which yield the change in thermal energy discussed above.

3.2. Synthetic Emission

In order to compare how each of these simulations would manifest in real observations, we generated synthetic intensities using the FoMo code [69]. We consider the 171 Å channel as observed by the Atmospheric Imaging Assembly (AIA) [70] aboard the Solar Dynamics Observatory, which detects emission from plasma with temperatures of approximately 1 MK. Figure 12 shows synthetic intensities from each of the four simulations. The line-of-sight (LOS) is in the y direction, which is parallel to the wave driver velocity in the ResUni and BbUni simulations (Figure 12, two upper plots).

Despite the enhanced density in the centre of each flux tube (e.g., see Figure 3), the radiation from these regions does not lead to brighter emission in the centre of the loops when viewed in this channel. This is because the cooler plasma (see Figure 5) only emits weakly in this channel. This is consistent with Figure 12 from Ref. [41]. Indeed, a low emission core is detectable between bright bands in the loop boundary in the second panel (BbUni) of Figure 12, top second plot (BbUni). In all plots, the synthetic emission shows apparent substructure within the coronal loop. These are reminiscent of coronal strands but, as has been reported in similar studies (e.g., [71]), these are manifestations of just the K-H vortices that develop along the length of the structure. In the unidirectional cases (Figure 12, two upper plots), these substructures form parallel to the loop axis and can be relatively homogeneous along the length of the loop (e.g., the bright strands along $x=\pm 0.5$ Mm in Figure 12, top second plot. This, however, is not the case in the simulations where the imposed driving changes directions (Figure 12, two lower plots). Here, the more complex motions lead to emission, which appears to show evidence of magnetic strands with an element of twist around the flux tube axis (particularly in the BbMulti case). More generally, the displacement of the loop across the plane-of-the-sky (POS; in the x direction here) is visible here. Let us note that different viewing angles (not parallel to the wave driver) would show similar displacements for the ResUni and BbUni cases.

One can also track the observable signatures of the evolution of the KHI in these systems. Figure 13 shows a time–distance dependence of the 171 Å emission emanating from the loop apex. In particular, we stack emission from the loop apex ($z=100$ Mm line, see Figure 12) in time. Figure 13, two upper plots, shows emission from the ResUni simulation along two different LOS. In Figure 13, top, the LOS is parallel to the y axis (as in Figure 12), and in Figure 13, top second plot, the LOS is rotated to be parallel to the x axis. As such, Figure 13, top, is parallel to the driver polarisation, and in Figure 13, top second, the wave appears as a transverse motion in the POS. Figure 13, two lower plots, shows similar plots for the ResMulti and BbMulti simulations, respectively. In both of these cases, the LOS is parallel to the y axis. Since the driver polarisation changes during the simulation, this angle is no longer parallel to the direction of the oscillation and, hence, both LOS and POS velocities develop. As such, characteristics of both of Figure 13, two upper plots, are important here.

In Figure 13, top, the formation of the K-H vortices is eminently visible after $t\approx 700$ s. However, we note that this may require much higher spatial resolution beyond the observing power of AIA. Whilst the kink mode is not visible in the intensity in this LOS, Doppler velocities would reveal the nature of the oscillation. In Figure 13, top second plot, however, the transverse oscillation is visible in the intensity and identifying the wave period is simple enough. Let us note that in this continuously driven system, we do not see any decay of the oscillation amplitude, and the disintegration of the loop emission is again indicative of the formation of the KHI. Leaving aside observational considerations, such results are not directly applicable to solar observations, as wave excitation is more likely to be similar to the ResMulti and BbMulti cases (Figure 13, two lower plots). Here, one can see that the transverse oscillation is poorly defined and signatures of the KHI are more subtle. However, for the given spatial resolution, the change in the structure from $t\approx 2500$ s, is indicative of the density deformation associated with the KHI.

4. Discussion

In this paper, we have presented the results of three-dimensional MHD simulations of wave heating in transversely oscillating low density coronal loops. In these models, the wave dynamics drive the formation of the KHI, which, in turn, deforms the loop density profile and leads to large gradients in both the perturbed magnetic and velocity fields. Through this carefully-studied process, wave heating rates are significantly enhanced, increasing the possibility of waves being important in maintaining coronal temperatures. In an extension to previous results, we have shown that heating powered by broadband wave drivers can be sufficient to balance optically thin radiative losses, provided that there is power in resonant frequencies (e.g., fundamental and higher harmonic kink and Alfvén frequencies). Indeed, by considering the power injected by different frequency components of the broadband driving, we show that these resonant frequencies dominate the energetics of the system. We also found that multi-directional drivers produced more energetic systems than their uni-directional counterparts. Although the more complex, multi-directional drivers typically result in delayed instability onset (due to less coherent velocity shears), once the KHI forms, they can generate greater heating rates.

Despite this positive result for wave heating, a number of limitations in this model still remain. Firstly, the heating rates in these simulations remain relatively low, and radiative losses can only be balanced because of the low densities considered here. For higher densities more typical of active region loops, the radiative losses can be 100 times greater (assuming plasma densities are an order of magnitude larger than in the current setup). As such, the wave heating rates found in these simulations will not be sufficient to balance radiative losses in denser loops [41]. Furthermore, we currently neglect cooling due to thermal conduction, which transfers heat to lower atmospheric layers. This effect can be significant and is often larger than the radiative cooling in the corona, particularly for shorter loops. Even in the long loops (200 Mm) considered here, energy losses from the corona due to conduction can be comparable to radiative losses. As such, the heating requirements are larger (perhaps significantly) than those in the simulations considered here. Additionally, we note that, due to computational constraints, the dissipation coefficients used here are many orders of magnitude larger than typical coronal values and thus heating rates are likely much higher than should be expected (e.g., [72]). However, if the instability induces true MHD turbulence, the energy dissipation rates may be largely independent of the dissipation coefficients (e.g., [48]). Consequently, whether the corona is truly turbulent is an important open question, which has significant consequences for the energetics of the system [73].

As with many wave heating models, the current study relies on an existing coronal loop, which is assumed to be dense in relation to its environment. Although the exact nature of the density profile may be unimportant, a cross-field gradient in the natural Alfvén frequency of field lines must be maintained for several wave periods to allow the instability to develop. Hitherto, no wave heating models have been able to self-consistently generate or maintain such a density profile. Indeed, in Ref. [74], it is argued that wave heating cannot sustain a dense loop core against energy losses. As a result, the dense, cooling plasma would quickly drain downwards to the lower atmosphere. Furthermore, in a related numerical setting [75], it was shown that evaporative upflows of chromospheric and transition region plasma driven by wave heating may be small. Even in the case here (which does not permit the draining of plasma), we see the dense loop core cooling as it experiences relatively weak heating due to low velocity and magnetic field gradients developing in this region. However, this may be less of an issue in simulations with larger wave amplitudes and/or higher Reynolds numbers as, in such cases, the entire loop cross-section is disrupted (e.g., [30]). This leads to lower average densities in the loop core and heating over a wider region.

The complex wave drivers considered here (particularly the BbMulti simulation) are likely to be more representative of the continuous wave excitation mechanisms in the real Sun. For example, it is unlikely that an oscillatory velocity in the photosphere maintains a single polarisation angle for thousands of seconds (as in the ResUni case). Despite this, the white noise (power constant over all frequencies) driver implemented in these simulations is unlikely to be representative of atmospheric wave drivers (e.g., there may be preferred frequencies which reflect a transfer of energy from p-modes to atmospheric oscillations [57]). Regardless of the exact nature of the oscillatory driving, it is likely that resonant frequencies play an important role in the energetics. However, the growth of resonances relies on the excitation mechanism containing power at the correct frequency for many wave periods. In particular, this component of the driving must have a constant phase shift over this time. Whether this is realistic assumption for the complex solar atmosphere remains unclear. Furthermore, if the background medium is evolving (e.g., field aligned plasma flows generated in response to heating and cooling), it is more likely that transient resonances develop. Then, it is unclear whether the KHI and significant wave heating can be driven in these more dynamic configurations. As such, more study is required to understand the effects of realistic and self-consistent wave driving on oscillating loops to ascertain the propensity for coronal wave heating.