© ESO 2019

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, 2015), Parnell & De Moortel (2012), and De Moortel & Browning (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; van Ballegooijen et al. 2011) or coupling with compressive wave modes (Kudoh & Shibata 1999; 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; Thurgood & McLaughlin 2013a) 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; 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 β ≪ 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 (2013a) 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 β = 0 plasma, the amplitude of the generated density profile grows algebraically with a t² profile, but this growth is limited by the plasma pressure force if β ≠ 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 Sect. 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 Sect. 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 Sect. 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:

$$\begin{array}{c}\hfill \frac{\mathrm{D}\rho }{\mathrm{D}t}=-\rho \mathbf{\nabla }·\mathit{\upsilon },\end{array}$$(1)

$$\begin{array}{c}\hfill \rho \frac{\mathrm{D}\mathit{\upsilon }}{\mathrm{D}t}=\mathit{j}\times \mathit{B}-\mathbf{\nabla }p+{\mathit{F}}_{\nu }^{\mathrm{shock}},\end{array}$$(2)

$$\begin{array}{c}\hfill \rho \frac{\mathrm{D}ϵ}{\mathrm{D}t}=-p(\mathbf{\nabla }·\mathit{\upsilon })+j^2/\sigma +H_{\nu }^{\mathrm{shock}},\end{array}$$(3)

$$\begin{array}{c}\hfill \frac{\mathrm{D}\mathit{B}}{\mathrm{D}t}=(\mathit{B}·\mathbf{\nabla })\mathit{\upsilon }-(\mathbf{\nabla }·\mathit{\upsilon })\mathit{B}+\eta {\mathrm{\nabla }}^2\mathit{B},\end{array}$$(4)

$$\begin{array}{c}\hfill p=\frac{k_{\mathrm{B}}}{{\mu }_{\mathrm{m}}{m}_{\mathrm{p}}}\rho T,\end{array}$$(5)

where $F_{\nu }^{\text{shock}}$ and $H_{\nu }^{\text{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 μ_m = 1/2. All other variables have their usual meanings. The code uses a uniform value for the magnetic diffusivity, η, given by

$$\begin{array}{c}\hfill \eta ={10}^{-3}{\eta }_0,\end{array}$$(6)

where

$$\begin{array}{c}\hfill {\eta }_0=\frac{L_0{B}_0}{\sqrt{\mu {\rho }_0}},\end{array}$$

where B₀, L₀, and ρ₀ are normalising constants and where ρ₀ also corresponds to the initial density. The value used for η 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 η is roughly equal to 1 m² s⁻¹ (see Priest 2014, p. 79), which means that for η in the code to be physically accurate, it should be approximately

$$\begin{array}{c}\hfill \eta ={10}^{-12}{\eta }_0,\end{array}$$

if L₀ = 1 Mm, B₀ = 10⁻³ T, and ρ₀ = 10⁻¹² kg m⁻³. The effects of varying the parameter η are explored in Sect. 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 Sect. 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 ( ρ₀) and pressure (p₀). The initial magnetic field is a potential 2D x-type null point, B₀, defined as

$$\begin{array}{c}\hfill {\mathit{B}}_0=\frac{B_0}{L_0}(x,-y,0).\end{array}$$(7)

The z-direction is taken as the invariant direction, i.e. ∂/∂z ≡ 0 throughout all the experiments. Our chosen magnetic field configuration is illustrated at the bottom of Fig. 1. It was chosen to represent a simplified version of the magnetic field in a mixed polarity region. The top left image in Fig. 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, 2016), and Thurgood & McLaughlin (2013b).

Fig. 1. Top left panel: magnetogram taken from the Hinode spacecraft of a mixed-polarity region. Top right panel: a simplified diagram of the magnetic field configuration in a mixed polarity region when viewed edge on (for example when viewed on the solar limb). Bottom panel: isolates the centre of the top right panel and is the profile of the magnetic field used in the numerical experiments of this paper.

The initial uniform plasma pressure was chosen such that most of the domain is a low-beta domain. The β = 1 contour is a circle occurring at a radius of R/L₀ = 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, υ = 0 and $\widehat{\mathit{n}}·\mathrm{\nabla }=0$ for all other variables, where $\widehat{\mathit{n}}$ is a vector normal to the boundary. The computational domain is given by

$$\begin{array}{c}\hfill x_{\min }\le x\le x_{\max },y_{\min }\le y\le y_{\max },\end{array}$$

where x_max = y_max = 2L₀ and x_min = y_min = −2L₀. On the y = y_min boundary, a continuous driver of the following form is imposed:

$$\begin{array}{c}\hfill {\upsilon }_z={\upsilon }_{\mathrm{driv}}f(x)g(t),\end{array}$$(8)

where υ_driv is the driver amplitude and τ_driv is the period of the driver. The spatial profile of the driver is described by the following equation:

$$\begin{array}{c}\hfill f(x)=\{\begin{array}{cc}1\hfill & |x|\le 1.5L_0,\hfill \\ {sin}^2(\pi x/L_0)\hfill & 1.5L_0\le |x|\le 2.0L_0.\hfill \end{array}\end{array}$$(9)

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

$$\begin{array}{c}\hfill g(t)=\{\begin{array}{cc}{sin}^2(2\pi t/{\tau }_{\mathrm{driv}})\hfill & t\le \frac{1}{4}{\tau }_{\mathrm{driv}},\hfill \\ sin(2\pi t/{\tau }_{\mathrm{driv}})\hfill & \frac{1}{4}{\tau }_{\mathrm{driv}}\le t\le t_{\mathrm{end}}^{\mathrm{driv}},\hfill \\ 0\hfill & t_{\mathrm{end}}^{\mathrm{driv}}\le t\le t_{\mathrm{end}}.\hfill \end{array}\end{array}$$(10)

The period of the driver is given by

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

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 Eq. (D.6)) such that the resonance occurs on field lines which are easier to describe analytically (see Sect. 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

$$\begin{array}{c}\hfill t_{\mathrm{end}}^{\mathrm{driv}}=20{\tau }_{\mathrm{driv}},\end{array}$$(11)

$$\begin{array}{c}\hfill t_{\mathrm{end}}=25{\tau }_{\mathrm{driv}}.\end{array}$$(12)

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), and De Groof & Goossens (2002). In Sect. 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_{\mathrm{lin}}^{\mathrm{end}}=E_{\mathrm{lin}}(t_{\mathrm{end}}^{\mathrm{driv}})$ refers to the total energy input from the driver after the driving has finished. We note that E_lin is a function of υ_driv and so each experiment is normalised by a different value depending on the value of υ_driv.

3. Numerical results

Before discussing the nonlinear effects in Sect. 4, we describe the linear effects that occur in the experiments, which are obtained by using a velocity amplitude of υ_driv = 10⁻³υ_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 Fig. 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 Fig. 3, which shows the velocity component of the Alfvén wave, υ_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.

Fig. 2. Contour plots of the (normalised) Alfvén wave velocity perturbations, υz, at different times.

Fig. 3. Absolute value of the velocity, |υz|, associated with the Alfvén waves along the line y = x at different times. The line y = x is perpendicular to the field lines and this shows the variation in υz across the field lines. The figure shows that the length scale across the field lines is shortened as time progresses and so phase mixing is occurring.

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 ∇_⊥B_z gradients whereas wave steepening generates ∇_||B_z gradients, namely,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{||}=\frac{{\mathit{B}}_0·\mathrm{\nabla }}{|{\mathit{B}}_0|},{\mathrm{\nabla }}_{\perp }=\frac{\widehat{\mathit{z}}\times {\mathit{B}}_0·\mathrm{\nabla }}{|{\mathit{B}}_0|},\end{array}$$

where || refers to the component parallel to the magnetic field and ⊥ refers to the perpendicular component in the ẑ × B₀ direction. In Fig. 4, the Ohmic heating contributions from both types of gradients are plotted, where

$$\begin{array}{c}\hfill P_{||}^{\mathrm{Ohmic}}=\frac{1}{\sigma }{\left(\frac{{\mathrm{\nabla }}_{||}{B}_z}{\mu }\right)}^2,\end{array}$$

$$\begin{array}{c}\hfill P_{\perp }^{\mathrm{Ohmic}}=\frac{1}{\sigma }{\left(\frac{{\mathrm{\nabla }}_{\perp }{B}_z}{\mu }\right)}^2.\end{array}$$

Fig. 4. Ohmic heating contributions from ∇||Bz and ∇⊥Bz. Labels are provided on the right-hand side of the figure. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$ (see Sect. 2.3), which gives the average power input from the driver in an equivalent but linear and ideal set-up.

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 Fig. 2 and the final snapshots in Fig. 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 (τ_n) of each field line at the initial time are given by the following equation:

$${\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),$$(D.6)

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

$$A=\frac{B_0}{L_0}xy·$$(D.1)

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 Fig. 5 as a function of A. The period of the driver, given by Eq. (E.1), 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:

$$\frac{\mathit{xy}}{L_0^2}=\pm 4^{1-n},y\le 0,$$(E.2)

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

Fig. 5. First three harmonic periods (τn) of each field as a function of the vector potential (A) normalised by the period of the driver (τdriv).

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

$$\frac{\mathrm{d}{E}_{\mathrm{tot}}}{\mathrm{d}t}=-\frac{B_y}{\mu }{\displaystyle {\int }_{y=y_{\min }}{\upsilon }_z{B}_z\mathrm{d}x,}$$(B.6)

where E_tot is defined as

$$E_{\mathrm{tot}}={\displaystyle {\int }_S\frac{p}{\gamma -1}+\frac{B^2}{2\mu }+\frac{1}{2}\rho {\upsilon }^2\mathrm{d}S,}$$(B.3)

where S is the computational domain. The driver velocity, υ_z, is imposed on the boundary by the driver given by Eq. (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 υ_z. We introduce the dimensionless parameter K(x, t), given by

$$\begin{array}{c}\hfill \frac{B_z}{\sqrt{\mu }}=-K\sqrt{{\rho }_0}{\upsilon }_z,\end{array}$$

which defines the relationship between υ_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

$$\begin{array}{c}\hfill K_{\mathrm{driv}}=\frac{1}{\sqrt{{\rho }_0\mu }}\frac{{\int }_0^{t_{\mathrm{end}}^{\mathrm{driv}}}{B}_z{\upsilon }_z\mathrm{d}t}{{\int }_0^{t_{\mathrm{end}}^{\mathrm{driv}}}{\upsilon }_z^2\mathrm{d}t},\end{array}$$(13)

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

$$\begin{array}{c}\hfill K_{\mathrm{driv}}=\frac{-B_y/\mu {\int }_0^{t_{\mathrm{end}}^{\mathrm{driv}}}{\upsilon }_z{B}_z\mathrm{d}t}{-B_y/\mu {\int }_0^{t_{\mathrm{end}}^{\mathrm{driv}}}{\upsilon }_z(-\sqrt{{\rho }_0\mu }{\upsilon }_z)\mathrm{d}t},\end{array}$$

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}{\upsilon }_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 Fig. 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 Sect. 4). The total energy in the domain is shown in Fig. 7 where the step-like profile corresponds to the period of the driver.

Fig. 6. Driver effectiveness Kdriv (Eq. (13)) on the bottom boundary.

Fig. 7. Change in total energy (Eq. (B.3)) from its initial value (Etot0) for different driver amplitudes. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.

Fig. 8. Alfvén wave energy (EA) for different driver amplitudes. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$, which gives the total energy input from the driver in an equivalent but linear and ideal experiment.

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_{\mathrm{lin}}^{\mathrm{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

$$\begin{array}{c}\hfill e_{\mathrm{A}}=\frac{1}{2}\rho {\upsilon }_z^2+\frac{B_z^2}{2\mu }·\end{array}$$

The total Alfvén wave energy in the domain, E_A, is then given by

$$\begin{array}{c}\hfill E_{\mathrm{A}}={\displaystyle {\int }_S{e}_{\mathrm{A}}}\mathrm{d}S,\end{array}$$(14)

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

$$\begin{array}{c}\hfill E_{\mathrm{acoustic}}={\displaystyle {\int }_S\frac{1}{2}\rho ({\upsilon }_x^2+{\upsilon }_y^2)+\frac{B_x^2+B_y^2}{2\mu }}\mathrm{d}S.\end{array}$$(15)

Although the Alfvén wave energy stops growing (see Fig. 8), the total energy of the system continues to grow (see Fig. 7) until t = 20τ_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 Sect. 5 it is shown that most of the energy goes into heat and not magneto-acoustic energy. Since a steady state is reached, Fig. 6 also gives a good indication of which field lines are heated most. From Fig. 8, it can be seen that steady state is reached at about t = 5τ_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 Fig. 7) in the linear experiment.

From Fig. 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τ_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 Fig. 9. As stated in Sect. 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.

Fig. 9. Contour plot showing the change in temperature relative to the initial temperature at $t=t_{\mathrm{end}}^{\mathrm{driv}}$.

Perhaps unexpectedly, significantly more heating occurs near (x, y)/L₀ = (±2, 0) compared with (x, y)/L₀ = (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 (2013a) 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 Fig. 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, $\mathbf{\nabla }{B}_z^2/(2\mu )$, also acts to push plasma towards the antinodes, away from the nodes. The right-hand side of Fig. 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 υ_A/υ_s, where υ_s is the sound speed. Hence, if a higher plasma-beta were used then the amplitude of the density structures would be reduced.

Fig. 10. Left panel: contour of the density after the driving has finished for the experiment where a driver amplitude of υdriv = 10−2υA0 is used. Right panel: density along one of the outer most resonant field lines for all three experiments, where the same colour scheme is used as in previous plots. The value lmax is equal to twice the length of the field line and l is a variable giving the distance from the centre of the field line.

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 Fig. 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)

$$\begin{array}{c}\hfill \frac{{\partial }^2{\upsilon }_z}{\partial t^2}=\frac{B_0^2}{\mu L_0^2}\frac{1}{\rho (s)}\frac{{\partial }^2{\upsilon }_z}{\partial s^2},\end{array}$$(16)

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, ω, given by

$$\begin{array}{c}\hfill \omega =k{\upsilon }_{\mathrm{A}}\end{array}$$

then the density structures will oscillate with a frequency, ω_s, which is approximately given by

$$\begin{array}{c}\hfill {\omega }_{\mathrm{s}}=2{\upsilon }_{\mathrm{s}}k,\end{array}$$

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 υ_s ≪ υ_A as this means ω_s ≪ ω. A second simplification is made by assuming the time dependence of υ_z is given by e^iωt. By multiplying through by υ_z, replacing time derivatives with iω and using integration by parts, Eq. (16) can be written as

$$\begin{array}{c}\hfill {\omega }^2=\frac{B_0^2}{\mu L_0^2}\frac{\int {(\partial {\upsilon }_z/\partial s)}^2\mathrm{d}s}{\int \rho (s){\upsilon }_z^2\mathrm{d}s},\end{array}$$(17)

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

$$\begin{array}{c}\hfill \tau =2\pi \frac{L_0}{{\upsilon }_{\mathrm{A}0}}\sqrt{\frac{\int \rho (s){\upsilon }_z^2\mathrm{d}s}{\int {(\partial {\upsilon }_z/\partial s)}^2\mathrm{d}s}}·\end{array}$$(18)

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), Fig. 11 predicts the location of the resonance to be at about 0.6L₀ compared to the actual location of about 0.7L₀ (see Fig. 6). Given the number of approximations, this is a reasonable agreement. In addition, Fig. 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.

Fig. 11. x-coordinate of where the resonant field lines crosses the bottom boundary as a function of time (normalised by the period of the driver). The resonant location was calculated by finding the field line with a fundamental time period, given by Eq. (18), which is closest in value to the driving time period.

From Fig. 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_{\mathrm{tot}}-E_{\mathrm{tot}0})/E_{\mathrm{lin}}^{\mathrm{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 Sect. 3.4, the Ohmic heating is the dominant effect which determines the shape of the temperature profile in the experiments. In experiments where η → 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

$$\frac{\mathrm{d}{E}_{\mathrm{A}}}{\mathrm{d}t}=\frac{\mathrm{d}{E}_{\mathrm{tot}}}{\mathrm{d}t}+{\displaystyle {\int }_S\mathit{\upsilon }·\mathbf{\nabla }\left(\frac{B_z^2}{2\mu }\right)-\frac{1}{\sigma }{\left(\frac{\mathbf{\nabla }{B}_z}{\mu }\right)}^2}\mathrm{d}S,$$(C.3)

where Eq. (C.3) is derived in Appendix C. Equation (C.3) 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,

$$\begin{array}{c}\hfill P_{\mathrm{A}}^{\mathrm{pond}}(t)=-{\displaystyle {\int }_S\mathit{\upsilon }·\mathbf{\nabla }\left(\frac{B_z^2}{2\mu }\right)}\mathrm{d}S,\end{array}$$(19)

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

$$\begin{array}{c}\hfill P_{\mathrm{A}}^{\mathrm{Ohmic}}(t)=\frac{1}{\sigma }{\displaystyle {\int }_S{\left(\frac{\mathbf{\nabla }{B}_z}{\mu }\right)}^2}\mathrm{d}S,\end{array}$$(20)

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 Eq. (C.3) 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τ_driv. As stated in Sect. 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

$$\begin{array}{c}\hfill {\omega }_{\mathrm{s}}=2{\upsilon }_{\mathrm{s}}k,\end{array}$$

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τ_driv (see the black curve in Fig. 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τ_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 Sect. 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.

Fig. 12. Ohmic power, $P_{\mathrm{A}}^{\mathrm{Ohmic}}$, and the ponderomotive power, $P_{\mathrm{A}}^{\mathrm{pond}}$ as functions of time, for different driver amplitudes. Labels are provided on the right-hand side of the figure. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$, which gives the average power input from the driver in an equivalent but linear and ideal set-up.

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, σ, 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 Fig. 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 2013b,a). 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 Fig. 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 (2013b), 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.

Fig. 13. Internal and acoustic energy for different driver amplitudes. Labels are provided on the right-hand side of the figure. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.

Fig. 14. Ratio of the root mean square velocity in the bottom half of the domain to the driver amplitude.

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

$$\begin{array}{c}\hfill \langle \frac{\mathit{E}\times \mathit{B}}{\mu }\rangle ={\upsilon }_{\mathrm{Ay}}{\rho }_0{\upsilon }_{\mathrm{driv}}^2⟨K_{\mathrm{driv}}⟩,\end{array}$$(21)

for t >  τ_driv/4, where ⟨K_driv⟩ is defined as

$$\begin{array}{c}\hfill ⟨K_{\mathrm{driv}}⟩=\frac{1}{(x_{\max }-x_{\min })}{\displaystyle {\int }_{x_{\min }}^{x_{\max }}{K}_{\mathrm{driv}}(x)f(x)}\mathrm{d}x.\end{array}$$(22)

In the experiments presented in this section, ⟨K_driv⟩≈ 0.55, 0.53, 0.45 corresponding to υ_driv/υ_A0 = 10⁻³,  10⁻²,  10⁻¹, respectively. McIntosh et al. (2011) observed Alfvén waves with an average amplitude between 20 and 25 km s⁻¹ 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 ρ^−1/4, provided there is no reflection. Photospheric motions are approximately 1 − 2 km s⁻¹ (Beliën et al. 1999; 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 Fig. 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: υ_driv = 20 − 25 km s⁻¹, υ_Ay = 200 − 250 km s⁻¹, ρ₀ = (5−10) × 10⁻¹³ kg m⁻³, ⟨K_driv⟩ = 0.5, then the Poynting flux is given by

$$\begin{array}{c}\hfill {\upsilon }_{\mathrm{Ay}}{\rho }_0{\upsilon }_{\mathrm{driv}}^2⟨K_{\mathrm{driv}}⟩\approx 20-80\mathrm{W}{\mathrm{m}}^{-2}.\end{array}$$(23)

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 Sect. 2, it was shown that η is too large by about a factor of 10⁹. Figure 15 shows results from experiments where η was varied. For these figures, the driver amplitude is set equal to

$$\begin{array}{c}\hfill {\upsilon }_{\mathrm{driv}}={10}^{-3}{\upsilon }_{\mathrm{A}}^{\mathrm{norm}},\end{array}$$

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 η, but appears to converge towards a minimum as η decreases. In the η = 0 experiment, the total energy was calculated by calculating the amount of Poynting flux entering the system through the driver. For small η (and for fixed υ_driv), it appears the total energy evolution is independent of η, 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 η. 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 η for small η whereas Wright & Allan (1996) showed that the total Ohmic heating is independent of η 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 η = 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 η. The figure shows that increasing η 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.

Fig. 15. Total energy of the domain for different values of η. In the η = 0 experiment there is no energy transfer due to Ohmic heating, however, there are still energy losses through numerical dissipation. The total energy was calculated using the Poynting flux on the boundary, therefore, any numerical energy losses in the domain are accounted for. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.

Fig. 16. Driver effectiveness (Eq. (13)) for different values of η. We note that numerical dissipation occurs in all the experiments including the η = 0 experiment. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.

Even though there is no (explicit) magnetic diffusion in the η = 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 η >  10⁻⁸η₀ Ohmic heating dominates and for this reason, we only considered experiments with η >  10⁻⁸η₀.

Fig. 17. Ratio, Enumeric/EOhmic, for a range of η values. The quantity Enumeric refers to the energy in the domain which is lost through numerical diffusion and is estimated by comparing the total (volume integrated) energy in the domain with the Poynting flux through the driven boundary. The quantity EOhmic gives the total amount of energy produced by Ohmic heating.

Increasing η results in more Alfvén waves leaking across the separatrices (see Fig. 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_{\mathrm{A}}^{y>0}$, as a function of time. Figure 18 shows that increasing η results in more leakage, however when compared with Fig. 15, it can be seen that even in the higher η 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 η = 10⁻³η₀) and this tends to zero as η → 0.

Fig. 18. Alfvén wave energy in the top half of the domain, $E_{\mathrm{A}}^{y>0}$, for different values of η. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.

Our earlier assessment of the Poynting flux was based on comparing υ_rms in our experiments with the values estimated by McIntosh et al. (2011) from SDO/AIA observations and assuming that υ_driv ∼ υ_rms (see Fig. 14 for η = 10⁻³η₀). However, for smaller values of η, the ratio υ_rms/υ_driv increases as can be seen in Fig. 19. Therefore, to compare with the same, observed value of υ_rms, we would have to reduce υ_driv, resulting in a smaller Poynting flux. At the same time, for smaller values of η, 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 (υ_driv = 10⁻¹υ_A0), υ_rms/υ_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 ∼ (υ/υ_A)²t, where υ is the amplitude of a wave. Hence, we expect the effect of reducing η on υ_rms/υ_driv to be less significant than for the linear experiments shown in Fig. 19. Therefore, if observed values of υ_rms can be considered to be nonlinear, our assumption that υ_driv ∼ υ_rms is perhaps not unreasonable, even for smaller values of η.

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 η 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).

Fig. 19. Ratio of the root mean square velocity to driver amplitude for different values of magnetic diffusivity. In this case, η/η0 = [10−1, 10−2, …, 10−8] with the bottom curve corresponding to 10−1 with each successive curve corresponding to the next value of η in the list above.

Fig. 20. Maximum velocity in the domain for different values of magnetic diffusivity. In this case, η/η0 = [10−1, 10−2, …, 10−8] with the bottom curve corresponding to 10−1 with each successive curve corresponding to the next value of η in the list above.

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 (η = 10⁻³η₀), 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, υ_rms, matches observations (McIntosh et al. 2011). We show that as η decreases, the Poynting flux remains approximately constant. However, υ_rms increases and this means that for smaller η the driver amplitude must be reduced to ensure υ_rms remains fixed. We estimated from linear experiments that if a physical value of η 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 Fig. 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 (υ/υ_A)⁻². Thus, although our linear results suggest that with a realistic value of η 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; Thurgood & McLaughlin 2013a) and that the Alfvén waves generate density structures as shown in Terradas & Ofman (2004). Equation (C.3) 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.

Acknowledgments

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).

References

Appendix A: Ideal and linear analytic solution

The ideal and linear solution was calculated by solving the wave equation (Eq. (D.3)) using d’Alembert’s formula for both υ_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

$$\begin{array}{c}\hfill {\upsilon }_z={\upsilon }_{\mathrm{driv}}sin(\omega t),\end{array}$$

with no ramp-up period involving the square of a sine (as used in Eq. (10)). Once B_z and υ_z are calculated, the Poytning flux on the boundary is calculated as

$$\begin{array}{cc}\hfill -\frac{B_y}{\mu }{\upsilon }_z{B}_z& ={\upsilon }_A^y{\rho }_0{\upsilon }_{\mathrm{driv}}^{2}sin(\omega t)(sin\left[\omega (t-\frac{2ml}{{\upsilon }_{\mathrm{A}0}})\right]\hfill \\ \hfill & +2cot\left(\frac{\omega l}{{\upsilon }_{\mathrm{A}0}}\right)sin\left(\frac{m\omega l}{{\upsilon }_{\mathrm{A}0}}\right)sin\left[\omega (t-\frac{\mathit{ml}}{{\upsilon }_{\mathrm{A}0}})\right]),\hfill \end{array}$$(A.1)

for ωl/υ_A0 ≠ kπ, where k is an integer, l = l(x) is the length of the loop given by s₂ − s₁ in Appendix D, ${\upsilon }_{\mathrm{Ay}}=B_y/\sqrt{{\rho }_0\mu }$ and m is an integer given by

$$\begin{array}{c}\hfill m=\lfloor \frac{t{\upsilon }_{\mathrm{A}0}}{2l}\rfloor ·\end{array}$$

The floor brackets, ⌊ ⌋, correspond to the largest integer smaller than tv_A0/(2l). The apparent singularity at ωl/υ_A0 = kπ, can be resolved and if ωl/υ_A0 = kπ then the Poynting flux is given by

$$\begin{array}{c}\hfill -\frac{B_y}{\mu }{\upsilon }_z{B}_z={\upsilon }_A^y{\rho }_0{\upsilon }_{\mathrm{driv}}^2(2m+1){sin}^2(\omega t).\end{array}$$

This equation shows that the Poynting flux grows linearly with time along resonant field lines and so the energy grows quadratically. To calculate $E_{\mathrm{lin}}^{\mathrm{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 0 to $t_{\mathrm{end}}^{\mathrm{driv}}$.

Appendix B: Total energy evolution

Using Eqs. (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

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\frac{p}{\gamma -1}+\frac{B^2}{2\mu }+\frac{1}{2}\rho {\upsilon }^2)+\mathbf{\nabla }·\mathit{F}=0,\end{array}$$(B.1)

where

$$\begin{array}{c}\hfill \mathit{F}=\frac{\gamma p}{\gamma -1}\mathit{\upsilon }+\frac{\mathit{E}\times \mathit{B}}{\mu }+\frac{1}{2}\rho {\upsilon }^2\mathit{\upsilon }.\end{array}$$

Taking the integral over the whole domain and making use of Gauss’ divergence theorem, Eq. (B.1) can be written as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{tot}}}{\mathrm{d}t}=-{\displaystyle {\oint }_{\partial S}\left(\frac{\mathit{E}\times \mathit{B}}{\mu }\right)·\widehat{\mathit{n}}}\mathrm{d}l,\end{array}$$(B.2)

where E_tot gives the total energy in the domain S, E_tot is defined as

$$\begin{array}{c}\hfill E_{\mathrm{tot}}={\displaystyle {\int }_S\frac{p}{\gamma -1}+\frac{B^2}{2\mu }+\frac{1}{2}\rho {\upsilon }^2}\mathrm{d}S.\end{array}$$(B.3)

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

$$\begin{array}{c}\hfill \mathit{E}=\mathit{j}/\sigma -\mathit{\upsilon }\times \mathit{B}.\end{array}$$(B.4)

Therefore,

$$\begin{array}{cc}\hfill \frac{\mathit{E}\times \mathit{B}}{\mu }·\widehat{\mathit{n}}& =\eta \mathit{j}\times \mathit{B}·\widehat{\mathit{n}}-\frac{1}{\mu }(\mathit{\upsilon }\times \mathit{B})\times \mathit{B}·\widehat{\mathit{n}},\hfill \\ \hfill & =\eta \mathit{j}\times \mathit{B}·\widehat{\mathit{n}}-\frac{1}{\mu }(\mathit{\upsilon }·\mathit{B})\mathit{B}·\widehat{\mathit{n}}.\hfill \end{array}$$

This term can be simplified further because $\mathit{j}\times \mathit{B}·\widehat{\mathit{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:

$$\begin{array}{cc}\hfill \frac{\partial B_y}{\partial t}& =\mathrm{\nabla }\times (\mathit{\upsilon }\times \mathit{B})·\widehat{\mathit{y}}+\eta {\mathrm{\nabla }}^2{B}_y,\hfill \\ \hfill & =\frac{\partial }{\partial x}({\upsilon }_x{B}_y-B_x{\upsilon }_y)+\eta {\mathrm{\nabla }}^2{B}_y,\hfill \\ \hfill & =0,\hfill \end{array}$$(B.5)

where x-derivative term equals zero because υ = 0 on the boundary so it is constant along the bottom boundary; the Laplacian term equals zero because initially B_y = −B₀y/L₀ and therefore remains zero for all time. Now consider the y-component of the Lorentz force on the bottom boundary as follows:

$$\begin{array}{cc}\hfill \mathit{j}\times \mathit{B}·\widehat{\mathit{y}}& =(\frac{1}{\mu }(\mathit{B}·\mathbf{\nabla })\mathit{B}-\mathbf{\nabla }\left(\frac{B^2}{2\mu }\right))·\widehat{\mathit{y}},\hfill \\ \hfill & =\frac{1}{\mu }(B_x\frac{\partial }{\partial x}+B_y\frac{\partial }{\partial y})B_y-\frac{\partial }{\partial y}\left(\frac{B^2}{2\mu }\right),\hfill \\ \hfill & =0,\hfill \end{array}$$(B.6)

where the y-derivatives are zero because $\widehat{\mathit{n}}·\mathrm{\nabla }=0$ on the boundary and the x-derivatives are zero because B_y = B_y(y). Hence, Eq. (B.2) can be written as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{tot}}}{\mathrm{d}t}=-\frac{1}{\mu }{\displaystyle {\int }_{y=y_{\min }}{\upsilon }_z{B}_z{B}_y}\mathrm{d}x,\end{array}$$(B.7)

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

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{tot}}}{\mathrm{d}t}=-\frac{B_y}{\mu }{\displaystyle {\int }_{y=y_{\min }}{\upsilon }_z{B}_z}\mathrm{d}x.\end{array}$$(B.8)

Appendix C: Total Alfvén wave energy evolution

Using Eqs. (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

$$\begin{array}{c}\hfill \frac{\partial e_{\mathrm{A}}}{\partial t}+\mathbf{\nabla }·{\mathit{F}}_A=\mathit{\upsilon }·\mathbf{\nabla }\left(\frac{B_z^2}{2\mu }\right)-\frac{1}{\sigma }{\left(\frac{\mathbf{\nabla }{B}_z}{\mu }\right)}^2,\end{array}$$(C.1)

where F_A is the Alfvén wave energy flux and is given by

$$\begin{array}{c}\hfill {\mathit{F}}_A=\frac{1}{2}\rho {\upsilon }_z^2\mathit{\upsilon }+\mathit{E}\times {\mathit{B}}_z.\end{array}$$(C.2)

In Appendix B it was shown that $\mathit{j}\times \mathit{B}·\widehat{\mathit{n}}=0$ on the boundary, hence if υ_x = υ_y = 0 on the boundary then it can be shown that

$$\begin{array}{c}\hfill \mathit{E}\times {\mathit{B}}_z·\widehat{\mathit{n}}=\mathit{E}\times \mathit{B}·\widehat{\mathit{n}},\end{array}$$

on the boundaries of the domain. Consequently, by taking the integral of Eq. (C.1) and substituting Eq. (B.8) the following equation is obtained:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{A}}}{\mathrm{d}t}=\frac{\mathrm{d}{E}_{\mathrm{tot}}}{\mathrm{d}t}+{\displaystyle {\int }_S\mathit{\upsilon }·\mathbf{\nabla }\left(\frac{B_z^2}{2\mu }\right)-\frac{1}{\sigma }{\left(\frac{\mathbf{\nabla }{B}_z}{\mu }\right)}^2}\mathrm{d}S.\end{array}$$(C.3)

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 ≥ 0. Symmetry arguments can be used to derive the formula for the other quadrants. The change of coordinates is given by

$$\begin{array}{c}\hfill x=\sqrt{\widehat{A}}{L}_0{e}^s,y=\sqrt{\widehat{A}}{L}_0{e}^{-s},\end{array}$$

where

$$\begin{array}{c}\hfill A=A_0\widehat{A}=\frac{B_0}{L_0}xy\end{array}$$(D.1)

and A₀ = B₀L₀. The linearised ideal wave equation is given by

$$\begin{array}{c}\hfill \frac{{\partial }^2{\upsilon }_z}{\partial t^2}=\frac{{({\mathit{B}}_0·\mathbf{\nabla })}^2}{\mu {\rho }_0}{\upsilon }_z,\end{array}$$(D.2)

Using the fact that

$$\begin{array}{c}\hfill {\mathit{B}}_0·\mathbf{\nabla }=\frac{B_0}{L_0}(\frac{\mathrm{d}x}{\mathrm{d}s}\frac{\partial }{\partial x}+\frac{\mathrm{d}y}{\mathrm{d}s}\frac{\partial }{\partial y})=\frac{B_0}{L_0}\frac{\partial }{\partial s},\end{array}$$

the wave equation (Eq. (D.2)) can be rewritten as

$$\begin{array}{c}\hfill \frac{{\partial }^2{\upsilon }_z}{\partial t^2}=\frac{B_0^2}{\mu {\rho }_0{L}_0^2}\frac{{\partial }^2{\upsilon }_z}{\partial s^2}·\end{array}$$(D.3)

Thus, the harmonic periods are given by

$$\begin{array}{c}\hfill {\tau }_n=\frac{2}{n}\frac{L_0}{{\upsilon }_{\mathrm{A}0}}(s_2-s_1),\end{array}$$(D.4)

where s₁ and s₂ are the values of s at each of the footpoints. In the first quadrant, at s₁, y = y_max and at s₂, x = x_max, therefore

$$\begin{array}{c}\hfill s_1=-log\left(\frac{y_{\max }}{L_0\sqrt{\widehat{A}}}\right),s_2=log\left(\frac{x_{\max }}{L_0\sqrt{\widehat{A}}}\right)·\end{array}$$

Finally, the harmonic periods are given by

$$\begin{array}{c}\hfill {\tau }_n=\frac{2}{n}\frac{L_0}{{\upsilon }_{\mathrm{A}0}}log\left(\frac{A_0}{L_0^2}\frac{x_{\max }{y}_{\max }}{A}\right),\end{array}$$(D.5)

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

$$\begin{array}{c}\hfill {\tau }_n=\frac{2}{n}\frac{L_0}{{\upsilon }_{\mathrm{A}0}}log\left(\frac{A_0}{L_0^2}\frac{x_{\max }{y}_{\max }}{|A|}\right)·\end{array}$$(D.6)

Appendix E: Resonance locations

The driving period is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{driv}}=\frac{L_0\sqrt{\mu {\rho }_0}}{B_0}4log\left(2\right).\end{array}$$(E.1)

The driving period equals one of the harmonic periods, where τ_n = τ_driv and

$$\begin{array}{cc}\hfill \frac{|A|}{A_0}& =\frac{|xy|}{L_0^2}=\frac{x_{\max }{y}_{\max }}{L_0^2{4}^n},\hfill \\ \hfill & \Rightarrow \frac{\mathit{xy}}{L_0^2}=\pm 4^{1-n},y\le 0,\hfill \end{array}$$(E.2)

where y ≤ 0 because the driver is imposed on the bottom boundary.

All Tables

All Figures

Fig. 1. Top left panel: magnetogram taken from the Hinode spacecraft of a mixed-polarity region. Top right panel: a simplified diagram of the magnetic field configuration in a mixed polarity region when viewed edge on (for example when viewed on the solar limb). Bottom panel: isolates the centre of the top right panel and is the profile of the magnetic field used in the numerical experiments of this paper.

In the text

	Fig. 2. Contour plots of the (normalised) Alfvén wave velocity perturbations, υz, at different times.
In the text

	Fig. 3. Absolute value of the velocity, |υz|, associated with the Alfvén waves along the line y = x at different times. The line y = x is perpendicular to the field lines and this shows the variation in υz across the field lines. The figure shows that the length scale across the field lines is shortened as time progresses and so phase mixing is occurring.
In the text

	Fig. 4. Ohmic heating contributions from ∇||Bz and ∇⊥Bz. Labels are provided on the right-hand side of the figure. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$ (see Sect. 2.3), which gives the average power input from the driver in an equivalent but linear and ideal set-up.
In the text

	Fig. 5. First three harmonic periods (τn) of each field as a function of the vector potential (A) normalised by the period of the driver (τdriv).
In the text

	Fig. 6. Driver effectiveness Kdriv (Eq. (13)) on the bottom boundary.
In the text

	Fig. 7. Change in total energy (Eq. (B.3)) from its initial value (Etot0) for different driver amplitudes. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.
In the text

	Fig. 8. Alfvén wave energy (EA) for different driver amplitudes. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$, which gives the total energy input from the driver in an equivalent but linear and ideal experiment.
In the text

	Fig. 9. Contour plot showing the change in temperature relative to the initial temperature at $t=t_{\mathrm{end}}^{\mathrm{driv}}$.
In the text

Fig. 10. Left panel: contour of the density after the driving has finished for the experiment where a driver amplitude of υdriv = 10−2υA0 is used. Right panel: density along one of the outer most resonant field lines for all three experiments, where the same colour scheme is used as in previous plots. The value lmax is equal to twice the length of the field line and l is a variable giving the distance from the centre of the field line.

In the text

	Fig. 11. x-coordinate of where the resonant field lines crosses the bottom boundary as a function of time (normalised by the period of the driver). The resonant location was calculated by finding the field line with a fundamental time period, given by Eq. (18), which is closest in value to the driving time period.
In the text

Fig. 12. Ohmic power, $P_{\mathrm{A}}^{\mathrm{Ohmic}}$, and the ponderomotive power, $P_{\mathrm{A}}^{\mathrm{pond}}$ as functions of time, for different driver amplitudes. Labels are provided on the right-hand side of the figure. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}/t_{\mathrm{end}}^{\mathrm{driv}}$, which gives the average power input from the driver in an equivalent but linear and ideal set-up.

In the text

	Fig. 13. Internal and acoustic energy for different driver amplitudes. Labels are provided on the right-hand side of the figure. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.
In the text

	Fig. 14. Ratio of the root mean square velocity in the bottom half of the domain to the driver amplitude.
In the text

Fig. 15. Total energy of the domain for different values of η. In the η = 0 experiment there is no energy transfer due to Ohmic heating, however, there are still energy losses through numerical dissipation. The total energy was calculated using the Poynting flux on the boundary, therefore, any numerical energy losses in the domain are accounted for. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.

In the text

	Fig. 16. Driver effectiveness (Eq. (13)) for different values of η. We note that numerical dissipation occurs in all the experiments including the η = 0 experiment. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.
In the text

	Fig. 17. Ratio, Enumeric/EOhmic, for a range of η values. The quantity Enumeric refers to the energy in the domain which is lost through numerical diffusion and is estimated by comparing the total (volume integrated) energy in the domain with the Poynting flux through the driven boundary. The quantity EOhmic gives the total amount of energy produced by Ohmic heating.
In the text

	Fig. 18. Alfvén wave energy in the top half of the domain, $E_{\mathrm{A}}^{y>0}$, for different values of η. The plots have been normalised by $E_{\mathrm{lin}}^{\mathrm{end}}$, which gives the total energy input from the driver in an equivalent but linear and ideal set-up.
In the text

	Fig. 19. Ratio of the root mean square velocity to driver amplitude for different values of magnetic diffusivity. In this case, η/η0 = [10−1, 10−2, …, 10−8] with the bottom curve corresponding to 10−1 with each successive curve corresponding to the next value of η in the list above.
In the text

	Fig. 20. Maximum velocity in the domain for different values of magnetic diffusivity. In this case, η/η0 = [10−1, 10−2, …, 10−8] with the bottom curve corresponding to 10−1 with each successive curve corresponding to the next value of η in the list above.
In the text