Phase mixing of nonlinear Alfven waves

TL;DR

This study uses 2.5D numerical simulations to investigate how nonlinear Alfvén waves undergo phase mixing, affecting energy transfer and dissipation in the solar corona, with implications for coronal heating.

Contribution

It introduces a detailed numerical model of nonlinear Alfvén wave phase mixing near null points, highlighting the effects of nonlinearities on energy transfer and resonant locations.

Findings

Nonlinearities reduce energy entering the domain.

Density structures alter natural frequencies and shift resonant locations.

Poynting flux is comparable to solar coronal energy requirements.

Abstract

Aims: This paper presents 2.5D numerical experiments of Alfv\'en wave phase mixing and aims to assess the effects of nonlinearities on wave behaviour and dissipation. In addition, this paper aims to quantify how effective the model presented in this work is at providing energy to the coronal volume. Methods: The model is presented and explored through the use of several numerical experiments which were carried out using the Lare2D code. The experiments study footpoint driven Alfv\'en waves in the neighbourhood of a two-dimensional x-type null point with initially uniform density and plasma pressure. A continuous sinusoidal driver with a constant frequency is used. Each experiment uses different driver amplitudes to compare weakly nonlinear experiments with linear experiments. Results: We find that the wave trains phase-mix owing to variations in the length of each field line and variations in the field strength. The nonlinearities reduce the amount of energy entering the domain, as they reduce the effectiveness of the driver, but they have relatively little effect on the damping rate (for the range of amplitudes studied). The nonlinearities produce density structures which change the natural frequencies of the field lines and hence cause the resonant locations to move. The shifting of the resonant location causes the Poynting flux associated with the driver to decrease. Reducing the magnetic diffusivity increases the energy build-up on the resonant field lines, however, it has little effect on the total amount of energy entering the system. From an order of magnitude estimate, we show that the Poynting flux in our experiments is comparable to the energy requirements of the quiet Sun corona. However a (possibly unphysically) large amount of magnetic diffusion was used however and it remains unclear if the model is able to provide enough energy under actual coronal conditions.

Tables

Table 1. Table 1: Total coronal energy losses from conduction, radiation, and the solar wind in different regions of the corona, based on Withbroe & Noyes ( 1977 ) .

	Quiet	Coronal	Active
	Sun	hole	region
Total coronal
energy losses	$3\times {10}^2$	$8\times {10}^2$	${10}^4$
(W ${\text{m}}^{-2}$)

Equations

$$\frac{\text{D}\rho }{\text{D}t}=-\rho \text{@vec}\nabla \cdot \text{@vec}v,$$

$$\rho \frac{\text{D}\text{@vec}v}{\text{D}t}=\text{@vec}j\times \text{@vec}B-\text{@vec}\nabla p+\text{@vec}{F}_{\nu }^{shock},$$

$$\rho \frac{\text{D}ϵ}{\text{D}t}=-p(\text{@vec}\nabla \cdot \text{@vec}v)+j^2/\sigma +H_{\nu }^{shock},$$

$$\frac{\text{D}\text{@vec}B}{\text{D}t}=\left(\text{@vec}B\cdot \text{@vec}\nabla \right)\text{@vec}v-\left(\text{@vec}\nabla \cdot \text{@vec}v\right)\text{@vec}B+\eta {\nabla }^2\text{@vec}B,$$

$$p=\frac{k_B}{{\mu }_m{m}_p}\rho T,$$

$$\eta =10^{-3}{\eta }_0,$$

$${\eta }_0=\frac{L_0{B}_0}{\sqrt{\mu {\rho }_0}},$$

$$\eta =10^{-12}{\eta }_0,$$

$$\text{@vec}{B}_0=\frac{B_0}{L_0}(x,-y,0).$$

$$x_{min}\le x\le x_{max},y_{min}\le y\le y_{max},$$

$$v_z=v_{driv}f(x)g(t),$$

$$f(x)=\cases{1}&|x|\leq 1.5L_{0},\\ \sin^{2}(\pi x/L_{0})1.5L_{0}\leq|x|\leq 2.0L_{0}.{}$$

$$g(t)=\cases{\sin}^{2}(2\pi t/\tau_{driv})&t\leq\frac{1}{4}\tau_{driv},\\ \sin(2\pi t/\tau_{driv})\frac{1}{4}\tau_{driv}\leq t\leq t_{end}^{driv},\\ 0t_{end}^{driv}\leq t\leq t_{end}.\\ {}$$

$${\tau }_{driv}=\frac{L_0\sqrt{\mu {\rho }_0}}{B_0}4\log \left(2\right),$$

$$t_{end}^{driv}=20{\tau }_{driv},$$

$$t_{end}=25{\tau }_{driv}.$$

$${\nabla }_{||}=\frac{\text{@vec}{B}_0\cdot \nabla }{|\text{@vec}{B}_0|},{\nabla }_{\perp }=\frac{\text{@vec}\hat{z}\times \text{@vec}{B}_0\cdot \nabla }{|\text{@vec}{B}_0|},$$

$$P_{||}^{Ohmic}=\frac{1}{\sigma }{\left(\frac{{\nabla }_{||}{B}_z}{\mu }\right)}^2,$$

$$P_{\perp }^{Ohmic}=\frac{1}{\sigma }{\left(\frac{{\nabla }_{\perp }{B}_z}{\mu }\right)}^2.$$

$${\tau }_n=\frac{L_0\sqrt{\mu {\rho }_0}}{B_0}\frac{2}{n}\log \left(\frac{A_0}{L_0^2}\frac{x_{max}{y}_{max}}{|A|}\right),$$

$$A=\frac{B_0}{L_0}xy.$$

$$\frac{xy}{L_0^2}=\pm 4^{1-n},y\le 0,$$

$$\frac{d{E}_{tot}}{dt}=-\frac{B_y}{\mu }{\int }_{y=y_{min}}{v}_z{B}_{z}dx,$$

$$E_{tot}={\int }_S\frac{p}{\gamma -1}+\frac{B^2}{2\mu }+\frac{1}{2}\rho v^{2}dS,$$

$$\frac{B_z}{\sqrt{\mu }}=-K\sqrt{{\rho }_0}{v}_z,$$

$$K_{driv}=\frac{1}{\sqrt{{\rho }_0\mu }}\frac{{\int }_0^{t_{end}^{driv}}{B}_z{v}_{z}dt}{{\int }_0^{t_{end}^{driv}}{v}_z^{2}dt},$$

$$K_{driv}=\frac{-B_y/\mu {\int }_0^{t_{end}^{driv}}{v}_z{B}_{z}dt}{-B_y/\mu {\int }_0^{t_{end}^{driv}}{v}_z(-\sqrt{{\rho }_0\mu }{v}_z)dt},$$

$$e_A=\frac{1}{2}\rho v_z^2+\frac{B_z^2}{2\mu }.$$

$$E_A={\int }_S{e}_{A}dS,$$

$$E_{acoustic}={\int }_S\frac{1}{2}\rho \left(v_x^2+v_y^2\right)+\frac{B_x^2+B_y^2}{2\mu }dS.$$

Full text

*Methods. *The model is presented and explored through the use of several numerical experiments which were carried out using the Lare2D code. The experiments study footpoint driven Alfvén waves in the neighbourhood of a two-dimensional x-type null point with initially uniform density and plasma pressure. A continuous sinusoidal driver with a constant frequency is used. Each experiment uses different driver amplitudes to compare weakly nonlinear experiments with linear experiments.

*Results. *We find that the wave trains phase-mix owing to variations in the length of each field line and variations in the field strength. The nonlinearities reduce the amount of energy entering the domain, as they reduce the effectiveness of the driver, but they have relatively little effect on the damping rate (for the range of amplitudes studied). The nonlinearities produce density structures which change the natural frequencies of the field lines and hence cause the resonant locations to move. The shifting of the resonant location causes the Poynting flux associated with the driver to decrease. Reducing the magnetic diffusivity increases the energy build-up on the resonant field lines, however, it has little effect on the total amount of energy entering the system. From an order of magnitude estimate, we show that the Poynting flux in our experiments is comparable to the energy requirements of the quiet Sun corona. However a (possibly unphysically) large amount of magnetic diffusion was used however and it remains unclear if the model is able to provide enough energy under actual coronal conditions.

1 Introduction

It has been hypothesised as early as Schatzman (1949) that the solar corona could be heated largely by the dissipation of magnetohydrodynamic (MHD) waves generated in the lower layers of the Sun. Heating by MHD waves is still one of the mechanisms under consideration for heating the corona; see for example Klimchuk (2006), Parnell & De Moortel (2012), De Moortel & Browning (2015) and Klimchuk (2015) for an overview of the coronal heating problem and the open questions that still need to be addressed. Magnetohydrodynamic waves are commonplace in the solar atmosphere and have been observed over the last two decades as a consequence of improved imaging and spectroscopic instruments (see e.g. De Moortel & Nakariakov (2012)).

A review of the linear behaviour of MHD waves can be found in, for example Goossens et al. (2011). The dissipation of Alfvén waves has been the basis of many coronal heating models (see review by Arregui et al. (2008), Arregui (2015) and references therein). The main mechanisms for converting the energy associated with Alfvén waves into heat are Ohmic and viscous dissipation. Both of these heating mechanisms are proportional to gradients in either magnetic field or velocity. The steeper the gradients, the more efficient the wave energy is converted into thermal energy. Two main mechanisms for generating steep gradients have been proposed, namely, phase mixing (Heyvaerts & Priest, 1983) and resonant absorption (Ionson, 1978). Phase mixing occurs when Alfvén waves propagating on magnetic surfaces move out of phase with neighbouring waves on nearby surfaces and this means that steep cross field gradients form. Resonant absorption occurs when driven standing Alfvén waves resonate while their neighbours on different magnetic surfaces are not resonating. Other mechanisms such as a turbulent cascade of wave energy (e.g. Hollweg (1986), Cranmer et al. (2007) and van Ballegooijen et al. (2011)) or coupling with compressive wave modes (Kudoh & Shibata (1999) and Antolin & Shibata (2010)) have also been proposed to increase the damping rate of Alfvén waves.

Phase mixing has been investigated as a mechanism for heating in many parts of the atmosphere, including coronal holes (e.g. Hood et al., 2002) and inside flux tubes (e.g. Pagano et al., 2018). In this paper, we evaluate the effect of nonlinearities on the amount of energy which can be provided to a coronal domain by Alfvén waves to account for the energy lost through optically thin radiation and thermal conduction. The energy required to keep a loop at coronal temperatures has been researched extensively (see for example Rosner et al. (1978), Martens (2010), Priest (2014) and references therein).

Nonlinear effects have been studied in a variety of settings. The magnetic tension force associated with an Alfvén wave is a linear force whereas the associated magnetic pressure force is a nonlinear force, often called the ponderomotive force (e.g. Verwichte et al. (1999).

In this paper, our model combines several nonlinear effects, the most important of which is the generation of density structures (e.g. Terradas & Ofman (2004)). Other nonlinear effects considered in our model are as follows: nonlinear coupling from Alfvén waves to magnetoacoustic waves (e.g. Verwichte et al. (1999), Tsiklauri et al. (2001) and Thurgood & McLaughlin (2013b) and Alfvén wave steepening due to the Alfvén speed being dependent on the perturbed magnetic energy associated with the Alfvén wave (e.g. Verwichte et al. (1999) and Tsiklauri et al. (2002)). However, these latter effects appear to play a less significant role.

The generation of density structures and the coupling to magnetoacoustic waves are generated by the ponderomotive force. Verwichte et al. (1999) showed that every time two Alfvén pulses superimpose each other they generate slow magnetoacoustic waves (in $\beta\ll 1$ plasma) due to a force they call the cross-ponderomotive force, which is a subset of the nonlinear magnetic pressure force. The same is true whenever Alfvén waves reflect off a solid boundary. Thurgood & McLaughlin (2013b) showed that if there are gradients in the Alfvén speed perpendicular to the magnetic field then fast waves are generated. Terradas & Ofman (2004) showed that a line-tied standing Alfvén wave pushes plasma towards the antinodes and away from the nodes due to the nonlinear magnetic pressure force and this creates a loop aligned density profile. They showed that in a $\beta=0$ plasma, the amplitude of the generated density profile grows algebraically with a $t^{2}$ profile, but this growth is limited by the plasma pressure force if $\beta\neq 0$.

In this paper, we focus on the nonlinear aspects of Alfvén wave propagation and dissipation in a simplified version of a coronal arcade system. The paper is structured as follows: In Section 2 the phase mixing model is presented and a linear and ideal solution to the model is calculated, which is used to compare with the nonlinear experiments. In Section 3 the linear results are presented. Section 4 assesses how the nonlinearities affect the heating in the model. Section 5 is a discussion on quantifying how effective the system is at converting wave energy into heat for typical coronal values. Finally, in Section 6 conclusions are presented.

2 Method and set-up

2.1 Equations

All the numerical experiments presented in this paper are performed using the MHD code Lare2D (Arber et al., 2001). The code solves the following set of MHD equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\@vec{F}_{\nu}^{shock}$ and $H_{\nu}^{shock}$ are terms related to the shock viscosity of the code, which is based on the edge viscosity formulation in Caramana et al. (1998). The Boltzmann constant is denoted with $k_{B}$, the mass of a proton is denoted with $m_{p}$, and the mass fraction of the ions in proton masses is denoted with $\mu_{m}=1/2$. All other variables have their usual meanings. The code uses a uniform value for the magnetic diffusivity, $\eta$, given by

[TABLE]

where

[TABLE]

where $B_{0}$, $L_{0}$, and $\rho_{0}$ are normalising constants and where $\rho_{0}$ also corresponds to the initial density. The value used for $\eta$ is unphysically large, but was chosen to be as small as possible without the effects of numerical diffusion and dispersion becoming too large. In the corona, the value of $\eta$ is roughly equal to 1$\,\text{m}^{2}\,s^{-1}$ (see Priest, 2014, p. 79), which means that for $\eta$ in the code to be physically accurate, it should be approximately

[TABLE]

if $L_{0}=1$Mm, $B_{0}=10^{-3}$T, and $\rho_{0}=10^{-12}$kg $\text{m}^{-3}$. The effects of varying the parameter $\eta$ are explored in Section 5. In our experiments, the coefficient of kinematic viscosity is set to zero. Hence, there is no heating due to this term. However, shock viscosities are included, along with any heating associated with shocks. This was chosen partly because according to Van Doorsselaere et al. (2007) observational evidence favours a resistive (wave) heating mechanism for coronal loops over viscous dissipation. In Section 3.4, some of the consequences of setting the coefficent of kinematic viscosity to zero are discussed.

2.2 Initial conditions

Phase mixing and resonant absorption require gradients in the Alfvén speed and this is often achieved through the use of an imposed density profile. A recent paper by Cargill et al. (2016) demonstrated that wave heating cannot sustain the assumed density structure. For this reason, we do not assume any density structures in this paper and instead, a uniform initial density profile is used. Our experiment relies on a gradient in magnetic field strength and also a variation in field line length to provide the conditions necessary for phase mixing.

An initial static equilibrium is set up with uniform density ($\rho_{0}$) and pressure ($p_{0}$). The initial magnetic field is a potential 2D x-type null point, $\@vec{B}_{0}$, defined as

[TABLE]

The $z$-direction is taken as the invariant direction, i.e. $\partial/\partial z\equiv 0$ throughout all the experiments. Our chosen magnetic field configuration is illustrated at the bottom of Figure 1. It was chosen to represent a simplified version of the magnetic field in a mixed polarity region. The top left image in Figure 1 shows a magnetogram of a (generic) mixed polarity region and was taken from the Hinode spacecraft. A mixed polarity field was used because it contains strong variations in field strength and field line length. Two simplifications were made: the first is that the field is approximated using a 2D model and the second is that an x-type null point configuration was used. The x-point field is qualitatively similar to the top right field, particularly close to the null point. The top right field is a simplified diagram of the magnetic field in a mixed polarity region if viewed on the limb. Waves near an x-point field have been studied extensively; see for example McLaughlin et al. (2011), McLaughlin (2013), Thurgood & McLaughlin (2013a) and McLaughlin (2016).

The initial uniform plasma pressure was chosen such that most of the domain is a low-beta domain. The $\beta=1$ contour is a circle occurring at a radius of $R/L_{0}=0.1$ about the null point. Since the magnetic field strength increases linearly with radius and the density is initially constant, the Alfvén speed increases linearly with radius.

2.3 Boundary conditions

To simulate the steep jump in density and temperature between the chromosphere and the corona, reflective boundary conditions are used (see Laney, 1998, p. 434). In other words, $\@vec{v}=0$ and $\@vec{\hat{n}}\cdot\nabla=0$ for all other variables, where $\@vec{\hat{n}}$ is a vector normal to the boundary. The computational domain is given by

[TABLE]

where $x_{max}=y_{max}=2L_{0}$ and $x_{min}=y_{min}=-2L_{0}$. On the $y=y_{min}$ boundary, a continuous driver of the following form is imposed:

[TABLE]

where $v_{driv}$ is the driver amplitude and $\tau_{driv}$ is the period of the driver. The spatial profile of the driver is described by the following equation:

[TABLE]

Equation (9) implies that the spatial profile of the driver is constant along most of the boundary and smoothly goes to zero at the ends. To ensure the driver smoothly generates a wave train, the time profile of the driver is given by

[TABLE]

The period of the driver is given by

[TABLE]

where the function $\log(2)$ was chosen for convenience rather than any physical reason. In particular, it was chosen to be of the same form as the resonant field line locations (see Equation (38)) such that the resonance occurs on field lines which are easier to describe analytically (see Section 3.2). After 20 driving time periods have elapsed, the driver is switched off and the experiments continue to run for another 5 driving time periods, such that

[TABLE]

[TABLE]

One key simplification which is made is that the frequency spectrum of the driver is discrete, while the frequency spectrum of the driver which generates waves in the corona likely resembles that of a broadband spectrum. For reference see for example Wright & Rickard (1995), De Groof et al. (2002), De Groof & Goossens (2002). In Section 5, the effects of a more random driver are briefly discussed.

To enable us to assess how nonlinearities affect the energy evolution in the experiments, we analytically calculate the energy evolution in a similar set-up which is ideal and linear. The details of this calculation are given in Appendix A. When referring to energy values from the linear reference calculation, the following notation is used: $E_{lin}$ refers to the total energy input from the driver in the analytical ideal and linear set-up and $E_{lin}^{end}=E_{lin}(t_{end}^{driv})$ refers to the total energy input from the driver after the driving has finished. We note that $E_{lin}$ is a function of $v_{driv}$ and so each experiment is normalised by a different value depending on the value of $v_{driv}$.

3 Numerical results

Before discussing the nonlinear effects in Section 4, we describe the linear effects that occur in the experiments, which are obtained by using a velocity amplitude of $v_{driv}=10^{-3}v_{A0}$ in the numerical experiments. In order to reduce the number of Figures in the paper the graphs present linear and nonlinear results. The linear results are represented by solid blue lines and this section focusses on these results.

3.1 Phase Mixing

The left-hand side of Figure 2 shows that the wave-front of the Alfvén wave is initially parallel to the $x$-axis. The right-hand side of the figure shows that at a later time, the Alfvén waves are out of phase with their neighbours. The phase mixing can also be seen in Figure 3, which shows the velocity component of the Alfvén wave, $v_{z}$, along the line $y=x$ at multiple times. It can be seen that the length scale across the field lines has shortened owing to phase mixing. In this particular set-up, there are two reasons why the phase mixing occurs. The first reason is that there is a gradient in Alfvén speed across the field lines due to variations in magnetic field strength. The second reason is that there is a variation in the length of each field line and so different waves reflect at different times relative to their neighbours, again leading to out-of-phase waves on neighbouring field lines.

The remainder of this subsection shows that the dominant reason for the formation of magnetic field gradients is indeed because of phase mixing and not other effects such as wave steepening as the waves approach the null. Phase mixing generates gradients in the magnetic field across different magnetic field lines, whereas wave steepening generates gradients parallel to the magnetic field lines. In other words, phase mixing generates $\nabla_{\perp}B_{z}$ gradients whereas wave steepening generates $\nabla_{||}B_{z}$ gradients, namely,

[TABLE]

where $||$ refers to the component parallel to the magnetic field and $\perp$ refers to the perpendicular component in the $\hat{z}\times\@vec{B}_{0}$ direction. In Figure 4, the Ohmic heating contributions from both types of gradients are plotted, where

[TABLE]

[TABLE]

Figure 4 clearly shows that substantially more heating occurs from gradients perpendicular to the magnetic field rather than parallel gradients, which confirms that phase mixing is indeed the dominant mechanism for generating the heating in this experiment.

3.2 Resonance

On the right-hand side of Figure 2 and the final snapshots in Figure 3, the signs of resonance occurring on discrete field lines can be seen. This section focusses on explaining why and where the resonance occurs. The harmonic time periods ($\tau_{n}$) of each field line at the initial time are given by the following equation:

[TABLE]

which is derived in Appendix D. In this case, $\@vec{A}=A\@vec{\hat{z}}$ is the vector potential of the magnetic field lines with flux function $A$, given by

[TABLE]

In 2D, lines of constant $A$ define the field lines and $A=0$ corresponds to the separatrices passing through the null point. The periods are plotted in Figure 5 as a function of $A$. The period of the driver, given by Equation (39), is overplotted as the red horizontal line. Resonance occurs on field lines where the period of the driver equals one of the harmonic periods. It can be seen that changing the period of the driver merely changes the location of the resonance but does not remove the resonance. In Appendix E, the locations at which resonance occurs for these waves are shown to lie on field lines described by the following equation:

[TABLE]

where the integer $n$ is the harmonic number. The above formula provides the resonance locations for linear waves in an ideal plasma. However, in Section 4, we show that some of the resonance locations shift outwards in the nonlinear experiments.

3.3 Energy evolution and driver effectiveness

With the exception of the driver, the velocity components are set to zero on all the boundaries. This implies that the Poynting flux is the only term responsible for changes in the total energy of the system (see Appendix B for proof). The change in total energy, $E_{tot}$, of the system is given by

[TABLE]

where $E_{tot}$ is defined as

[TABLE]

where $S$ is the computational domain. The driver velocity, $v_{z}$, is imposed on the boundary by the driver given by Equation (8). However, $B_{z}$ is free to adjust its value, which means the Poynting flux can be positive as well as negative. The value $B_{y}$ is positive and does not change on the driver boundary. Hence, the Poynting flux is determined by the relationship between $B_{z}$ and $v_{z}$. We introduce the dimensionless parameter $K(x,t)$, given by

[TABLE]

which defines the relationship between $v_{z}$ and $B_{z}$. The equation $K=1$ corresponds to a single upward propagating linear Alfvén wave. We define the driver effectiveness, $K_{driv}(x)$, on the bottom boundary as

[TABLE]

where $K_{driv}$ is a weighted time average of $K$. $K_{driv}$ gives a measure of how effectively the driver provides energy to a given field line. To see this more clearly, the above equation is equivalent to

[TABLE]

where the denominator corresponds to the Poynting flux through the boundary for an upward propagating linear Alfvén wave.

For a field line which has only upward propagating linear Alfvén waves at the driven boundary, $B_{z}=-\sqrt{\mu\rho_{0}}v_{z}$, hence, $K_{driv}=1$. A field line may only have upward propagating waves at the driven boundary because the waves are efficiently damped before they can reflect or the field is open. If the waves are reflected then the relationship is more complicated as there are now upward propagating and downward propagating waves at the driven boundary. In most cases, reflection acts to reduce driver effectiveness as reflection opens up the possibility for destructive interference to occur. However, if the frequency of the driver is such that it sets up a resonance, then this has the effect of increasing the driver effectiveness. This reflects the fact that resonant field lines have the property that a relatively small driver amplitude produces a large growth in energy.

The driver effectiveness on the bottom boundary is plotted as a function of $x$ in Figure 6. The plot clearly shows the locations of the resonating field lines, where most of the Poynting flux is concentrated. Figure 6 confirms that near the resonating field lines the driver effectiveness is greater than unity and away from the resonant lines it is much less than unity. The position of the resonant field line varies as the amplitude of the driver increases, due to nonlinear effects (see Section 4). The total energy in the domain is shown in Figure 7 where the step-like profile corresponds to the period of the driver.

Figure 8 shows that the total energy associated with the Alfvén waves, $E_{A}$, grows to a maximum, after which the energy oscillates about its time-averaged value (about 0.12$E_{lin}^{end}$). The Alfvén waves oscillate in the $z$-direction and so the energy density of an Alfvén wave, $e_{A}$, at a point in space is given by

[TABLE]

The total Alfvén wave energy in the domain, $E_{A}$, is then given by

[TABLE]

where $S$ is the computational domain. The magneto-acoustic energy, $E_{acoustic}$, is defined to be the magnetic and kinetic energy associated with all perturbations which are not Alfvén waves and is defined as

[TABLE]

Although the Alfvén wave energy stops growing (see Figure 8), the total energy of the system continues to grow (see Figure 7) until $t=20\tau_{driv}$ when the driver is switched off. Therefore, a steady state is reached, where all the energy generated by the time average of the Poynting flux goes into thermal energy and magneto-acoustic energy. In Section 5 it is shown that most of the energy goes into heat and not magneto-acoustic energy. Since a steady state is reached, Figure 6 also gives a good indication of which field lines are heated most. From Figure 8, it can be seen that steady state is reached at about $t=5\tau_{driv}$ in the linear experiment (blue curve). It is interesting to note that the transition from a transient state to a steady state does not have a noticeable impact on the total energy evolution (see Figure 7) in the linear experiment.

From Figure 7, it can be seen that the growth of the total energy is linear. This is not surprising during steady state, as the amplitudes of the waves have stopped growing and so all the Poynting flux is transferred into heat at a constant rate. However, even during the transient phase ($t<5\tau_{driv}$), the energy growth is linear. The linear growth during the transient phase occurs because the amplitude of the resonant wave grows quadratically with time, whereas the width of the resonant region decreases linearly with time (in accordance with resonant absorption theory (Ionson, 1978)). Hence, the total energy growth is linear.

3.4 Location of the heating

Most of the heating occurs at the nodes of the resonating standing field lines and this can be seen in Figure 9. As stated in Section 2.1, there is no viscous dissipation (besides shock viscosity) and therefore most of the heating occurs from Ohmic heating. Since most of the heating is generated by gradients in $B_{z}$ perpendicular to the magnetic field this means that most of the heating occurs where $B_{z}$ is largest. For standing Alfvén waves, the magnetic field component of the Alfvén wave is largest at the nodes of the standing wave and so most of the heating occurs at the nodes of the standing waves. Conversely, viscous dissipation acts on gradients in velocity and this would lead to heating occurring at the antinodes. However, Van Doorsselaere et al. (2007) showed from observational evidence that heating occurs mainly at loop footpoints and from this, they inferred that resistive heating dominates over viscous heating for wave heating mechanisms.

Perhaps unexpectedly, significantly more heating occurs near $(x,y)/L_{0}=(\pm 2,0)$ compared with $(x,y)/L_{0}=(0,-2)$. A similar phenomenon has been shown in, for example, McLaughlin (2013). The author showed that the wavefronts of Alfvén waves which are not reflective remain planar as they approach the null and the current builds up exponentially with time, causing most of the heating to occur at the horizontal separatrices and furthest from the null.

4 Nonlinear aspects

The effects of nonlinear Alfvén waves have been researched extensively; see for example Verwichte et al. (1999) for details on Alfvén waves in 1D and Thurgood & McLaughlin (2013b) for details in 2D and Terradas & Ofman (2004) for details on standing Alfvén waves. We analyse the effects of nonlinear density structures and investigate how the damping rate is affected by the nonlinearities.

4.1 Density structures

The nonlinearities generate density structures which can be seen in Figure 10, where red/yellow indicates density enhancement and green/blue indicates density reduction in the contour plots. Density structures are most pronounced along the resonating field lines where standing waves have been established. The density is largest at the antinodes of the standing waves and smallest at the nodes for two reasons. Firstly, as we only consider resistive heating and do not include viscosity, most of the heating takes place at the nodes. This causes the plasma pressure to increase at the nodes and hence a plasma pressure force is set up which pushes plasma away from the nodes towards the antinodes. If viscosity was included, it is likely some heating would also occur near the loop apex, reducing this effect. The second reason is that $B_{z}$ has its maximum amplitude at the nodes and smallest amplitude at the antinodes. This means that the nonlinear magnetic pressure force/ponderomotive force, $\@vec{\nabla}B_{z}^{2}/(2\mu)$, also acts to push plasma towards the antinodes, away from the nodes. The right-hand side of Figure 10 shows the density along one of the resonating field lines which are oscillating at the fundamental harmonic. It can be seen that plasma has been concentrated towards the apex/antinode of the field line and moved away from the footpoints/nodes of the field lines. For a more detailed analysis of the density structures formed by standing Alfvén waves, see Terradas & Ofman (2004). These authors show that the amplitude of the density structures is proportional to $v_{A}/v_{s}$, where $v_{s}$ is the sound speed. Hence, if a higher plasma-beta were used then the amplitude of the density structures would be reduced.

One key effect these density structures have is that they cause the natural periods of the field lines to change, which in turn changes the location of the resonance. Indeed, in Figure 6, it can be seen that the peaks have shifted away from the origin for higher amplitude drivers. The density structures result in the natural periods of the field lines increasing. As the Alfvén speed is now a function of position along the field lines, even in a field-aligned coordinate system, the period is not simply given by the wavelength divided by the wave speed. The density structures have an enhancement in density at the antinodes and as the amplitude of the plasma velocity is highest at the antinodes, changes to the Alfvén speed at the antinodes affect the period the most. Since the density is enhanced at the antinodes, this results in a decrease in Alfvén speed at the antinodes and thus an increase in the period.

More rigorously, this result can be derived by considering the wave equation (see Appendix D)

[TABLE]

where $s$ is related to the distance along a field line. For simplicity, the density is assumed to be a function of only space and not time. In Terradas & Ofman (2004), it was shown that for a standing Alfvén wave with angular frequency, $\omega$, given by

[TABLE]

then the density structures will oscillate with a frequency, $\omega_{s}$, which is approximately given by

[TABLE]

where $k$ is the wave number of the standing Alfvén wave. Hence, the simplification that the density is constant in time can be made if $v_{s}\ll v_{A}$ as this means $\omega_{s}\ll\omega$. A second simplification is made by assuming the time dependence of $v_{z}$ is given by $e^{i\omega t}$. By multiplying through by $v_{z}$, replacing time derivatives with $i\omega$ and using integration by parts, Equation (16) can be written as

[TABLE]

provided the integrals are not taken at a time when the denominator is equal to zero. Hence,

[TABLE]

Equation (18) confirms that the value of the density near where the velocity is largest (the antinode) is the most important. Figure 11 shows an approximation of the location of the field line (in terms of the $x$-coordinate of where the field line crosses the bottom boundary) with a fundamental harmonic period which is closest in value to the period of the driver. In other words, the figure shows the approximate location of the fundamental harmonic resonant field line. We note that only the coordinates in the $x>0$ side of the domain are shown. For the nonlinear experiment (black line), Figure 11 predicts the location of the resonance to be at about 0.6$L_{0}$ compared to the actual location of about $0.7L_{0}$ (see Figure 6). Given the number of approximations, this is a reasonable agreement. In addition, Figure 11 gives an indication of how quickly and in which direction the resonance moves. To derive the resonance location more rigorously, the wave equation would need to be solved as a Sturm-Liouville problem and the eigenvalues which satisfy the boundary conditions would give the resonant frequencies.

From Figure 7, it can be seen that the nonlinearities cause the system to be less effective at extracting energy from the driver. This is demonstrated by the fact that the final value for $(E_{tot}-E_{tot0})/E_{lin}^{end}$ is lower when a higher amplitude driver is used. One of the reasons the nonlinearities reduce the driver effectiveness is because they cause the resonance to shift location. Resonant field lines are effective at extracting energy from the driver because they generate a build-up in $B_{z}$, which increases the Poynting flux. Shifting the resonance location results in energy at the previous resonance location being lost owing to destructive interference. Moreover, the nonlinearities reduce the density at the footpoints, reducing the energy associated with the Alfvén waves generated by the boundary driver, resulting in less energy entering the system.

In addition to generating density structures, the ponderomotive force also acts to create a similarly shaped temperature profile. However, as shown in Section 3.4, the Ohmic heating is the dominant effect which determines the shape of the temperature profile in the experiments. In experiments where $\eta\rightarrow 0$, there is no longer any Ohmic heating and the ponderomotive force now enhances the temperature at the antinodes owing to plasma compression.

4.2 Damping rate

The term damping rate refers to the rate at which Alfvén wave energy is converted into other forms of energy. The change in total Alfvén wave energy is given by

[TABLE]

where Equation (32) is derived in Appendix C. Equation (32) shows that in the absence of a driver there are two terms which affect the energy evolution of the Alfvén wave energy. The first term,

[TABLE]

is referred to as the ponderomotive power, which is related to the work done by the magnetic pressure force/ponderomotive force. The second term,

[TABLE]

is referred to as the Alfvén Ohmic power and is related to the Ohmic heating. The ponderomotive power can both increase or decrease the Alfvén wave energy whereas the direct effect of the Ohmic heating only ever acts to reduce the Alfvén wave energy. It is worth noting that Equation 32 shows that it is possible for flows perpendicular to the invariant direction to increase the amplitude of the Alfvén waves, however, it is impossible for the flows to change the amplitude if an Alfvén wave initially has zero amplitude.

Figure 12 shows both the ponderomotive and Ohmic powers as functions of time. It can be seen that the net effect of the ponderomotive power is to increase the damping rate, however, its contribution is small compared to the Ohmic power. The ponderomotive power only has a large contribution at $t=2.5\tau_{driv}$. As stated in Section 4.1, Terradas & Ofman (2004) showed that for a 1D ideal standing Alfvén wave which is not driven, the amplitude of the density structures oscillates between a maximum and minimum at an angular frequency given by

[TABLE]

where $k$ is the wave number of the Alfvén wave. These authors showed that the longitudinal velocity also oscillates with this frequency. Hence, by energy conservation, if the energy of the longitudinal flows oscillates, the energy of the Alfvén waves must also oscillate and this energy change is carried out by the ponderomotive power. Therefore, the initial increase and then decrease in the ponderomotive power up to $t=3.5\tau_{driv}$ (see the black curve in Figure 12) can be understood as the magnetic pressure force generating longitudinal flows and then the restoring plasma pressure force acting to return the density structures to their equilibrium position. After the initial peak in ponderomotive power, the power becomes negative around $t=4\tau_{driv}$, reflecting the fact that the plasma pressure force is acting to push the density structures back to their equilibrium position. However, by this point, the system has reached steady state and so now there is sufficient Ohmic heating for the associated pressure forces to dominate the motion of the longitudinal flows. As stated in Section 4.1, the pressure forces associated with the Ohmic heating act with the magnetic pressure force to push plasma away from the nodes. Hence, once steady state is reached, the ponderomotive power rarely becomes negative and oscillates with a period which is half that of the driving period.

Figure 12 shows that the Ohmic power dominates over the ponderomotive power. One reason this property holds is that the ponderomotive power is a fourth-order nonlinear term because the velocity perturbations in the plane are second order, whereas the Ohmic power is a second-order nonlinear term. The Ohmic heating depends on the conductivity, $\sigma$, and so in a highly conductive plasma, it may at first seem that the ponderomotive power should dominate. However, if the plasma is highly conducting this often results in very short length scales forming and the Ohmic heating is inversely proportional to the square of the length scales while the ponderomotive power is only inversely to proportional to the length scale. Hence, we suggest that this property is also likely to hold in similar configurations.

Figure 8 shows that in all the experiments, the Alfvén wave energy decays to zero at nearly the same rate when the driver is switched off. This is somewhat unexpected given that, in Figure 10, it can be seen that the nonlinearities generate large density structures and density structuring is usually associated with the enhancement of phase mixing and thus increased wave dissipation. It would seem that because the density is only redistributed along the field lines in this figure, this does little to change the time taken for Alfvén waves to travel from one footpoint to another and so does little to enhance the phase mixing. In addition, the nonlinearities cause there to be an increase in coupling to compressive modes (Thurgood & McLaughlin (2013a) and Thurgood & McLaughlin (2013b)). As stated in the introduction, coupling to compressive modes has been proposed as an efficient mechanism for damping Alfvén waves. However, its contribution is relatively small; the coupling is a nonlinear effect and very quickly becomes negligible as the amplitude of the Alfvén wave decreases to zero. Evidence for this can be seen in Figure 12, where the ponderomotive power reaches zero much more quickly than the Ohmic power. For a more detailed analysis of the evolution of the compressive modes see Thurgood & McLaughlin (2013a), who study a nonlinear Alfvén pulse near an x-point with a similar magnetic field to that presented in this work.

5 Discussion

This section aims to assess whether the model and driver presented in this work can provide sufficient heat to balance conductive and radiative losses in the quiet Sun. Let us assume that 100% of the net energy provided by the driver (the Poynting flux) goes into heating the domain. The viability of this assumption is discussed at the end of this section.

The average Poynting flux at a point along the driver boundary is given by

[TABLE]

for $t>\tau_{driv}/4$, where $\langle K_{driv}\rangle$ is defined as

[TABLE]

In the experiments presented in this section, $\langle K_{driv}\rangle\approx$ 0.55, 0.53, 0.45 corresponding to $v_{driv}/v_{A0}=10^{-3},\ 10^{-2},\ 10^{-1}$, respectively. McIntosh et al. (2011) observed Alfvén waves with an average amplitude between $20-25$ km $\text{s}^{-1}$ at a height of 15 Mm in the quiet region of the solar corona. This velocity seems plausible, as the amplitude of an Alfvén wave scales with $\rho^{-1/4}$, provided there is no reflection. Photospheric motions are approximately $1-2$ km/s (Beliën et al. (1999) and Moriyasu & Shibata (2004)), giving a velocity amplitude around 100 times larger at the top of the chromosphere. There is some reflection at the transition region, and therefore the value at the top of the transition region observed by McIntosh et al. (2011) is not unreasonable. From Figure 14 it can be seen that at steady state, the ratio of the average amplitude to driver amplitude is approximately unity. Therefore, if the following values (taken from McIntosh et al. (2011)) are used: $v_{driv}=20-25$ km $\text{s}^{-1}$, $v_{Ay}=200-250$ km $\text{s}^{-1}$, $\rho_{0}=(5-10)\times 10^{-13}$ kg $\text{m}^{-3}$, $\langle K_{driv}\rangle=0.5$, then the Poynting flux is given by

[TABLE]

At steady state, the wave energy stops growing and so 100% Poynting flux provided by the driver goes into heat. The Poynting flux obtained above is of the order of the required flux to balance energy losses in the quiet Sun (see Table 1), suggesting that phase mixing of Alfvén waves, as in this model, may indeed play a significant role in the heating of the corona. However, the model presented made many simplifications and the remainder of this section discusses the potential consequences of some of these simplifications.

In Section 2, it was shown that $\eta$ is too large by about a factor of $10^{9}$. Figure 15 shows results from experiments where $\eta$ was varied. For these figures, the driver amplitude is set equal to

[TABLE]

and hence, the results are mostly linear. Also, the shock viscosity was switched off in these experiments such that the experiment can be considered close to ideal, although some amount of numerical diffusion is still present. We see that the total energy increases with $\eta$, but appears to converge towards a minimum as $\eta$ decreases. In the $\eta=0$ experiment, the total energy was calculated by calculating the amount of Poynting flux entering the system through the driver. For small $\eta$ (and for fixed $v_{driv}$), it appears the total energy evolution is independent of $\eta$, A similar phenomenon has been reported in the literature. For example, Wright & Allan (1996) measured the total Ohmic heating in a similar phase mixing experiment and proved analytically that once steady state is reached, the total spatially integrated Ohmic dissipation is independent of $\eta$. There is a subtle difference between the result derived in Wright & Allan (1996) and the result obtained in this work. The total energy evolution becomes independent of $\eta$ for small $\eta$ whereas Wright & Allan (1996) showed that the total Ohmic heating is independent of $\eta$ once steady state is reached. The proof in Wright & Allan (1996) assumes that the simulation has run sufficiently long for steady state to be reached, however, in this paper, for the $\eta=0$ experiment, steady state is not reached during the simulation time. Figure 16 helps to illustrate why the total energy converges to a limit for smaller $\eta$. The figure shows that increasing $\eta$ acts to increase the width of the resonating region, however, it also reduces the height and so there is little change to the total energy.

Even though there is no (explicit) magnetic diffusion in the $\eta=0$ experiment, wave energy is still dissipated through numerical diffusion. Figure 17 is presented such that the importance of numerical diffusion can be assessed. The figure shows the ratio of $E_{numeric}$ to $E_{Ohmic}$, where $E_{numeric}$ gives the energy in the domain which is lost through numerical diffusion and $E_{Ohmic}$ gives the total amount of energy produced by Ohmic heating. It can be seen that for $\eta>10^{-8}\eta_{0}$ Ohmic heating dominates and for this reason, we only considered experiments with $\eta>10^{-8}\eta_{0}$.

Increasing $\eta$ results in more Alfvén waves leaking across the separatrices (see Figure 18). None of the field lines in the bottom half of the domain enter the top half of the domain, hence, only a small fraction of the Alfvén waves leak into the top half of the domain. The Alfvén waves that travel into the top half of the domain do so via magnetic diffusion, which enables Alfvén waves to leak onto neighbouring field lines. Figure 18 shows the amount of Alfvén wave energy in the top half of the domain, denoted by $E_{A}^{y>0}$, as a function of time. Figure 18 shows that increasing $\eta$ results in more leakage, however when compared with Figure 15, it can be seen that even in the higher $\eta$ experiments the amount of leakage is still negligible when compared with the total energy increase of the system. Although nonlinearities appear to increase the amount of wave leakage, even in our most nonlinear experiments, the Alfvén wave energy which leaks across the separatrices is no more than 0.1% (for $\eta=10^{-3}\eta_{0}$) and this tends to zero as $\eta\rightarrow 0$.

Our earlier assessment of the Poynting flux was based on comparing $v_{rms}$ in our experiments with the values estimated by McIntosh et al. (2011) from SDO/AIA observations and assuming that $v_{driv}\sim v_{rms}$ (see Fig. 14 for $\eta=10^{-3}\eta_{0}$). However, for smaller values of $\eta$, the ratio $v_{rms}/v_{driv}$ increases as can be seen in Fig. 19. Therefore, to compare with the same, observed value of $v_{rms}$, we would have to reduce $v_{driv}$, resulting in a smaller Poynting flux. At the same time, for smaller values of $\eta$, nonlinear effects become more important as the maximum velocities in the domain grow substantially (see Fig. 20). McIntosh et al. (2011) argued that the observed amplitudes in the quiet Sun are of the order of 10% of the local Alfvén speed. From Fig. 14, we can see that for amplitudes of this order ($v_{driv}=10^{-1}v_{A0}$), $v_{rms}/v_{driv}$ is smaller than in the corresponding linear experiments. This can be explained by the fact that according to Verwichte et al. (1999) nonlinear damping mechanisms grow with $\sim(v/v_{A})^{2}t$, where $v$ is the amplitude of a wave. Hence, we expect the effect of reducing $\eta$ on $v_{rms}/v_{driv}$ to be less significant than for the linear experiments shown in Fig. 19. Therefore, if observed values of $v_{rms}$ can be considered to be nonlinear, our assumption that $v_{driv}\sim v_{rms}$ is perhaps not unreasonable, even for smaller values of $\eta$.

Finally, we point out that the assumption that 100% of the net energy provided by the driver goes into heating the domain depends on the frequency of the driver. For a constant frequency driver, the system eventually reaches a steady state where 100% of the Poynting flux from the driver goes into heat. However, for a non-constant driver frequency, the system may be unable to reach steady state before the frequency profile of the driver changes, particularly, if a smaller value of $\eta$ were used. Using a random driver could reduce the driver effectiveness as changes to the driver frequency could lead to destructive interference with pre-existing waves. The effects of using a random driver imposed on a zero-dimensional harmonic oscillator are reviewed in Masoliver & Porra (1993) and Gitterman (2013). A random driver in a 2D phase mixing experiment is investigated in Wright & Rickard (1995). The effects of a random driver imposed on a 3D coronal loop are studied in De Groof et al. (2002) and De Groof & Goossens (2002).

6 Conclusions

In this paper, we demonstrate that phase mixing can occur because of variations in field strength and field line length without the need for variations in density. The model deliberately does not impose any initial density structures as demonstrated by Cargill et al. (2016) that heating from phase mixing cannot sustain the density structures self-consistently.

We find that the nonlinearities reduce the driver effectiveness, which results in the total amount of heating being reduced. For the range of driver amplitudes studied in this paper, the reduction in effectiveness is found to be about 20%. Density structures are generated by the ponderomotive force and by pressure forces associated with the heating; this causes the resonance location to shift, which means energy build up is smaller than it would be otherwise. In addition to this, since the density at the boundary is reduced, the energy associated with the Alfvén waves entering the system is also reduced and so less energy enters the domain. The nonlinearities have a comparatively small effect on the damping rate (for the range of amplitudes studied in this work), where the damping rate is related to the rate at which the energy associated with the Alfvén waves is converted into other forms of energy.

We calculated an order of magnitude estimate of the Poynting flux to determine if the model presented in this work provides enough energy to balance conductive and radiative losses in the coronal region. We find that the Poytning flux provided in the model, with a large magnetic diffusion ($\eta=10^{-}3\eta_{0}$), indeed provides energy of the order necessary to balance conductive and radiative losses in the quiet Sun corona (but not active regions). We did not consider coronal holes because they are typically composed of open magnetic field lines and our model addresses closed loops. The order of magnitude estimate was constrained by ensuring that the steady-state, root-mean-square velocity, $v_{rms}$, matches observations (McIntosh et al., 2011). We show that as $\eta$ decreases, the Poynting flux remains approximately constant. However, $v_{rms}$ increases and this means that for smaller $\eta$ the driver amplitude must be reduced to ensure $v_{rms}$ remains fixed. We estimated from linear experiments that if a physical value of $\eta$ were used, the driver amplitude would have to be reduced by approximately a factor of 10, resulting in a decrease in Poynting flux by a factor of 100. From Figure 14 it can be seen that nonlinearities reduce the root-mean-square velocity. One possible mechanism for this could be the nonlinear self-modification of Alfvén waves which, as shown by Verwichte et al. (1999), results in the formation of strong currents and hence strong Ohmic dissipation in a time that is proportional to $(v/v_{A})^{-2}$. Thus, although our linear results suggest that with a realistic value of $\eta$ the model does not produce enough heat to balance losses in the corona, it is still plausible that there might be enough heat in a nonlinear model.

It has long been known that Alfvén waves can mode convert to magnetoacoustic waves (as described in Verwichte et al. (1999) and Thurgood & McLaughlin (2013b)) and that the Alfvén waves generate density structures as shown in Terradas & Ofman (2004). Equation (32) shows that the Alfvén waves can transfer energy into flows perpendicular to the invariant direction and vice versa by doing work through the magnetic pressure force (ponderomotive force). Therefore, if a large velocity is imposed in the longitudinal direction, this could result in large changes to the energy of the Alfvén wave.

Acknowledgements. This research has received funding from the Science and Technology Facilities Council (UK) through the consolidated grant ST/N000609/1 and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 647214).

Appendix A Ideal and linear analytic solution

The ideal and linear solution was calculated by solving the wave equation (Equation 35) using d’Alembert’s formula for both $v_{z}$ and $B_{z}$. We note that the solution is a superposition of Heaviside functions which enter the domain owing to wave reflection. A simplification is made by assuming the driver is given by

[TABLE]

with no ramp-up period involving the square of a sine (as used in Equation (10)). Once $B_{z}$ and $v_{z}$ are calculated, the Poytning flux on the boundary is calculated as

[TABLE]

for $\omega l/v_{A0}\neq k\pi$, where $k$ is an integer, $l=l(x)$ is the length of the loop given by $s_{2}-s_{1}$ in Appendix D, $v_{Ay}=B_{y}/\sqrt{\rho_{0}\mu}$ and $m$ is an integer given by

[TABLE]

The floor brackets, $\lfloor\;\rfloor$, correspond to the largest integer smaller than $tv_{A0}/(2l)$. The apparent singularity at $\omega l/v_{A0}=k\pi$, can be resolved and if $\omega l/v_{A0}=k\pi$ then the Poynting flux is given by

[TABLE]

This equation shows that the Poynting flux grows linearly with time along resonant field lines and so the energy grows quadratically. To calculate $E_{lin}^{end}$, the Poytning flux was then integrated along the bottom boundary in space and time, where $x$ goes from $x_{min}$ to $x_{max}$ and $t$ goes from [math] to $t_{end}^{driv}$.

Appendix B Total energy evolution

Using Equations (1) - (4) as well as Faraday’s law, it can be shown that the rate of change of energy at a point in space in the domain is given by

[TABLE]

where

[TABLE]

Taking the integral over the whole domain and making use of Gauss’ divergence theorem, Equation (24) can be written as

[TABLE]

where $E_{tot}$ gives the total energy in the domain $S$, $E_{tot}$ is defined as

[TABLE]

Most of the terms in $\@vec{F}$ can be eliminated because $\@vec{v}\cdot\@vec{\hat{n}}=0$ on every boundary. The Poynting flux term can be simplified further by making use of Ohm’s law as follows:

[TABLE]

Therefore,

[TABLE]

This term can be simplified further because $\@vec{j}\times\@vec{B}\cdot\@vec{\hat{n}}=0$ on the boundary. To demonstrate this, we first show that $B_{y}$ does not change on the bottom boundary. Consider the $y$-component of the induction equation as follows:

[TABLE]

where $x$-derivative term equals zero because $\@vec{v}=0$ on the boundary so it is constant along the bottom boundary; the Laplacian term equals zero because initially $B_{y}=-B_{0}y/L_{0}$ and therefore remains zero for all time. Now consider the $y$-component of the Lorentz force on the bottom boundary as follows:

[TABLE]

where the $y$-derivatives are zero because $\@vec{\hat{n}}\cdot\nabla=0$ on the boundary and the $x$-derivatives are zero because $B_{y}=B_{y}(y)$. Hence, Equation (25) can be written as

[TABLE]

where the integral is taken only over the bottom boundary as this is where the driver is located. Since $B_{y}$ does not depend on $x$ it can be taken out of the integral to give

[TABLE]

Appendix C Total Alfvén wave energy evolution

Using Equations (1) - (4) as well as Faraday’s law it can be shown that the rate of change of Alfvén wave energy density is given by

[TABLE]

where $\@vec{F}_{A}$ is the Alfvén wave energy flux and is given by

[TABLE]

In Appendix B it was shown that $\@vec{j}\times\@vec{B}\cdot\@vec{\hat{n}}=0$ on the boundary, hence if $v_{x}=v_{y}=0$ on the boundary then it can be shown that

[TABLE]

on the boundaries of the domain. Consequently, by taking the integral of Equation (30) and substituting Equation (29) the following equation is obtained:

[TABLE]

Appendix D Harmonic periods

The goal of this appendix is to derive an expression for the harmonic series associated with each field line as a function of the vector potential ($A$) at the initial time step. To do this, a change of coordinates is used to rewrite the wave equation in a form such that the wave speed is constant. For now, an expression is derived for the first quadrant where $A\geq 0$. Symmetry arguments can be used to derive the formula for the other quadrants. The change of coordinates is given by

[TABLE]

where

[TABLE]

and $A_{0}=B_{0}L_{0}$. The linearised ideal wave equation is given by

[TABLE]

Using the fact that

[TABLE]

the wave equation (Equation (34)) can be rewritten as

[TABLE]

Thus, the harmonic periods are given by

[TABLE]

where $s_{1}$ and $s_{2}$ are the values of $s$ at each of the footpoints. In the first quadrant, at $s_{1}$, $y=y_{max}$ and at $s_{2}$, $x=x_{max}$ , therefore

[TABLE]

Finally, the harmonic periods are given by

[TABLE]

using symmetry arguments it can be shown that the formula for all quadrants is given by

[TABLE]

Appendix E Resonance locations

The driving period is given by

[TABLE]

The driving period equals one of the harmonic periods, where $\tau_{n}=\tau_{driv}$ and

[TABLE]

[TABLE]

where $y\leq 0$ because the driver is imposed on the bottom boundary.