Issue 
A&A
Volume 632, December 2019



Article Number  A93  
Number of page(s)  12  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201936658  
Published online  09 December 2019 
Investigating the damping rate of phasemixed Alfvén waves
School of Mathematics and Statistics, University of St Andrews, St Andrews, Fife KY16 9SS, UK
email: apkp@standrews.ac.uk
Received:
9
September
2019
Accepted:
21
October
2019
Context. This paper investigates the effectiveness of phase mixing as a coronal heating mechanism. A key quantity is the wave damping rate, γ, defined as the ratio of the heating rate to the wave energy.
Aims. We investigate whether or not laminar phasemixed Alfvén waves can have a large enough value of γ to heat the corona. We also investigate the degree to which the γ of standing Alfvén waves which have reached steadystate can be approximated with a relatively simple equation. Further foci of this study are the cause of the reduction of γ in response to leakage of waves out of a loop, the quantity of this reduction, and how increasing the number of excited harmonics affects γ.
Methods. We calculated an upper bound for γ and compared this with the γ required to heat the corona. Analytic results were verified numerically.
Results. We find that at observed frequencies γ is too small to heat the corona by approximately three orders of magnitude. Therefore, we believe that laminar phase mixing is not a viable standalone heating mechanism for coronal loops. To arrive at this conclusion, several assumptions were made. The assumptions are discussed in Sect. 2. A key assumption is that we model the waves as strictly laminar. We show that γ is largest at resonance. Equation (37) provides a good estimate for the damping rate (within approximately 10% accuracy) for resonant field lines. However, away from resonance, the equation provides a poor estimate, predicting γ to be orders of magnitude too large. We find that leakage acts to reduce γ but plays a negligible role if γ is of the order required to heat the corona. If the wave energy follows a power spectrum with slope −5/3 then γ grows logarithmically with the number of excited harmonics. If the number of excited harmonics is increased by much more than 100, then the heating is mainly caused by gradients that are parallel to the field rather than perpendicular to it. Therefore, in this case, the system is not heated mainly by phase mixing.
Key words: Sun: corona / Sun: magnetic fields / magnetohydrodynamics (MHD) / Sun: oscillations / waves
© ESO 2019
1. Introduction
The coronal heating problem relates to the question of why the temperature of the corona is over a hundred times hotter than the photosphere; 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 remain to be addressed. The corona being hotter than the chromosphere, conduction is an energy loss mechanism. Similarly, the corona is optically thin, and therefore radiation is a loss mechanism. Ohmic dissipation of electric currents and viscous dissipation of plasma motions are thought to play a major role in balancing the conductive and radiative losses in the corona (Klimchuk 2015). It is unclear whether Ohmic or viscous heating dominates the other. Moreover, the precise mechanism(s) by which these occur remains an open question. The proposed mechanisms can be split into two categories: reconnection and wave heating mechanisms. This paper focuses on phase mixing, one of the wave heating mechanisms.
Magnetohydrodynamic (MHD) waves are commonplace in the solar atmosphere and have been observed over the last two decades as a consequence of new, improved imaging and spectroscopic instruments (see, e.g. Tomczyk et al. 2007; McIntosh et al. 2011; 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 2015 and references therein). It is believed that the main mechanisms by which Alfvén waves are converted into heat are Ohmic and viscous dissipation. Both of these mechanisms are proportional to the current density and vorticity. There are several proposed mechanisms which may be able to dissipate waves at a high enough rate to heat the corona. The mechanism that this paper focuses on is phase mixing as suggested by Heyvaerts & Priest (1983). Phase mixing is the process where gradients perpendicular to the field buildup due to Alfvén waves propagating on field lines with a spatial gradient in Alfvén travel time. This process leads to neighbouring waves moving out of phase with each other, hence the name phase mixing. Other notable mechanisms are: resonant absorption (Ionson 1982), reflectiondriven Alfvén wave turbulence (Hollweg 1986; van Ballegooijen et al. 2011; Shoda et al. 2019), turbulence triggered via the tearing mode or KelvinHelmholtz instability (Heyvaerts & Priest 1983; Browning & Priest 1984; Antolin et al. 2016, 2018) and coupling with compressive modes (Kudoh & Shibata 1999; Antolin & Shibata 2010).
Phase mixing has been researched quite extensively since Heyvaerts & Priest (1983) first proposed it as a coronal heating mechanism. For example, Browning & Priest (1984) expanded on the KelvinHelmholtz analysis and argue that the phase mixing of standing Alfvén waves can trigger turbulence, which can lead to a turbulent cascade and enhanced dissipation of wave energy. Similon & Sudan (1989) and Howson et al. (2019) study phase mixing in a complex magnetic field. Parker (1991) investigated phase mixing in a coronal hole and argues that it is not valid to assume an ignorable coordinate, and therefore that Alfvén waves couple to other modes and are not subject to pure phase mixing. Hood et al. (1997, 2002) investigate phase mixing in coronal holes, and find a selfsimilar solution which enables them to analyse a more general class of solutions. These latter authors find that a single pulse decays algebraically rather than exponentially. More recently, phase mixing has been investigated in 3D (Magyar et al. 2017) and 3D coronal loops (Pagano & De Moortel 2017; Pagano et al. 2018).
Here, we aim to provide a critical analysis of laminar phase mixing in coronal loops. Through our analysis, we aim to assess whether phase mixing can provide a sufficient damping rate, γ. The damping rate, γ, is defined as
where ⟨⟩ denotes the time average of either E × B/μ, the Poynting flux, or ρu^{2}/2, the kinetic energy density. We assumed that kinetic and magnetic energy is transported into the corona mainly via Poynting flux from lower layers in the atmosphere, which is then dissipated as heat (Klimchuk 2015). This definition was chosen in part because it can be calculated relatively easily using observational data. The required heating rate has been estimated by, for example, Withbroe & Noyes (1977). The amplitude of transverse oscillations in the corona has been observed by, for example, McIntosh et al. (2011). This definition approximates the more classical definition of the damping rate. To see this, we consider a damped harmonic oscillator with amplitude, x(t), with its motion described by
hence,
and so the heating rate is given by . Therefore,
Hollweg (1984a); Hollweg (1984b) uses a very similar idea to γ except he uses a quantity called the quality factor, Q, from resonance theory, which is approximately given by Q = ω/γ, where ω is the angular frequency of a wave. We focus on using γ because it is easier to apply to a system in which multiple harmonics oscillate. Arregui (2015) also stresses the importance of the damping rate, although he does not give a precise definition. He argues, (through an order magnitude analysis) that phase mixing could take too long to reach the required length scales for the heating to become important. In this paper, we show that even if the waves are allowed to evolve to steady state, the damping rate is still too low.
A damping rate of about 10^{−1} s^{−1} is required to heat coronal loops, which is based on predictions of the required heating rate and observations of velocity fluctuations in the corona. If we assume that the radiative loss function, Λ_{rad}(T), is given by
where T is temperature and χ = 10^{−32} K^{1/2} W m^{3} and a uniform coronal heating profile, H_{c}, then Priest (2014) shows that to maintain a loop with a uniform pressure, p_{0}, requires a heating rate of
where k_{B} is the Boltzmann constant and T_{max} is the maximum temperature. Here, T_{max} = 10^{6} K and p_{0} = 10^{−2} Pa gives H_{c} ≈ 10^{−5} W m^{−3}, and this approximately agrees with Withbroe & Noyes (1977) for a 100 Mm loop. McIntosh et al. (2011) and McIntosh & De Pontieu (2012) observe amplitudes in the quiet sun of around 20 km s^{−1} and in active regions of around 5 km s^{−1}. If we assume an amplitude of u_{0} = 10 km s^{−1} and a density of ρ_{0} = 10^{−12} km s^{−1}, this gives a kinetic energy density of J m^{−3}/2. Therefore, the required damping rate is approximately s^{−1}.
This paper does not investigate the origin of the coronal waves. We generate the waves using a driver at the edge of our domain. The origin of coronal waves remains an open question. The exponential density and the steep jump in density at the transition region make it difficult for Alfvén waves generated at the photosphere to enter the corona (Cranmer & van Ballegooijen 2005). Hollweg (1984b) suggests resonances in coronal loops and spicules provide enough energy flux to the corona. Cally & Hansen (2011) and Hansen & Cally (2012) suggest that mode conversion from fast waves to Alfvén waves at the transition region enables sufficient energy flux to enter the corona. It is also possible that the corona itself generates Alfvén waves via reconnection (Cranmer 2018).
This paper is structured as follows: Section 2 shows the equations which we solve in this paper and then discusses some of the simplifications which were used in deriving these equations; in particular, how the Braginskii viscous stress tensor can be simplified. In Sect. 3 and all subsequent sections, we consider a closedloop. The goal of Sect. 3 is to calculate the damping rate of standing Alfvén waves which have reached steady state and then to show it can be approximated with a relatively simple equation. Section 4 discusses the effects of allowing waves to leak out of the loop. Finally, Sect. 5 considers the effects of exciting multiple harmonics. Results from Sect. 5 are used in the conclusions to deduce that laminar phase mixing is not a viable standalone heating mechanism in coronal loops.
2. Equations
Our model is illustrated in Fig. 1. This section presents the equations that we solve. Section 2 discusses some of the simplifications which were made to arrive at these equations. The variables we consider are the magnetic field,
Fig. 1. Diagram of our model, showing a set of closed coronal loops (vertical black lines) with a driver imposed at the bottom boundary (blue). The gradient in height of the blue arrows indicates there is a gradient in Alfvén speed. 
the velocity,
and the density,
The equations are the linearised momentum equation,
and the linearised induction equation,
where the permeability of free space is denoted with μ, the viscosity coefficient, ν, is given by
and the magnetic diffusivity coefficient, η, is given by
where ω_{p} is the proton cyclotron frequency given by
with e the elementary charge and m_{p} the mass of a proton. Here, τ_{p} is the protonproton collision time, given by
where logΛ is the Coulomb logarithm ≈22 in coronal conditions (Priest 2014). Substituting T = 10^{6} K, B = 10^{−3} T and ρ = 10^{−12} m^{−3} gives ω_{p}τ_{p} ≈ 10^{5}. We note, ω_{p}τ_{p} ≫ 1 corresponds to the strong field limit which means that conduction and viscosity are highly anisotropic. The implications of this are discussed further in Sect. 2. We have neglected spatial derivatives of the viscosity coefficient, ν, and magnetic diffusivity, η, in favour of derivatives in the velocity and magnetic field perturbations. Equations (10) and (11) can be combined to give
provided the products and squares of the viscosity and magnetic diffusivity coefficients can be neglected, where is the Alfvén speed.
Equations (7) and (8) show that we assume an invariant direction. This approximation may be valid in, for example, coronal arcades or coronal loops which could be approximately invariant along the azimuthal direction. An important consequence of keeping an invariant direction in our system is that it stabilises to the KelvinHelmholtz instability and tearing instability. To trigger an instability, a disturbance needs to have a wave vector, k, that satisfies
Our boundary conditions ensure k_{z} ≠ 0 and since ∂/∂y = 0, this ensures that k_{y} = 0. Therefore,
for all possible disturbances and so instabilities are not triggered in our system.
Equation (7) shows that we model the field lines as straight. In general, the coronal field is curved, especially in coronal arcades, but the concepts outlined here can still be applied. The wave equation describes linear Alfvén waves for a potential field in a fieldaligned coordinate system. Therefore, there exists a straight field with a density structure which will reproduce the dynamics of an arcade field.
For analytical progress, we assume the Alfvén speed is uniform along the field lines in the corona. In reality, coronal loops are stratified due to gravity and the field strength may vary as well. Using a uniform loop instead of a nonuniform loop has two key effects: the first effect is that it means the wavelength of our waves do not change as they propagate along the loop. This may affect the phase mixing as shorterwavelength waves form shorter perpendicular length scales more quickly. However, if our uniform Alfvén speed represents the average Alfvén speed of a nonuniform loop, then the timeaveraged heating rates will not differ significantly. The second effect is that a nonuniform Alfvén speed can give rise to reflection. If the wavelength of a wave is shorter than the lengthscale of the Alfvén speed variations, then to a good approximation the reflection is negligible (Leroy 1980). In the corona, we estimate the pressure scale height to be approximately 60 Mm. Therefore, we argue that for waves with a wavelength smaller than this (say < 50 Mm) the reflection within the corona itself can be ignored. We note that a semicircular loop of length 100 Mm has a vertical height of approximately 32 Mm at the apex. In this paper, we present results with a wavelength varying from 10 Mm to 200 Mm. Waves with a longer wavelength than the pressure scale height are presented as we believe their results are still informative.
Our equations to model the viscosity and resistivity are much simpler than that outlined in Braginskii (1965); the remainder of this paragraph explains how we justify our expressions. As pointed out by Russell (2019), linearising the viscosity can result in the predicted viscosity being orders of magnitude smaller than it otherwise would be. The viscous force is modelled as the divergence of the Braginskii stress tensor, σ_{brag}, given here in the same form as in MacTaggart et al. (2017),
where
with
where Z is the tensor with components Z_{ij} = ϵ_{ikj}b_{k}, where ϵ_{ikj} are components of the LeviCivita symbol. The viscosity coefficients (see Braginskii 1965) are given by
assuming the strong field limit (ω_{p}τ_{p} ≫ 1). We aim to show that η_{0}W^{(0)}, η_{2}W^{(2)}, η_{3}W^{(3)}, η_{4}W^{(4)} can be neglected in favour of η_{1}W^{(1)}. This is surprising because if ω_{p}τ_{p} ≈ 10^{5} then this means that η_{1} ≈ 10^{−10}η_{0}.
Here we assume that u and b are footpointdriven waves on closed loops that have reached steady state such that their amplitudes do not change with time and are of the form,
where k denotes the wave number. If we assume the viscosity coefficients are uniform then the viscous heating is given by σ_{brag} : ∇u. The heating contributions from each term, with u and B given by Eqs. (7) and (8) is given by:
We expect v/v_{A} to be approximately in the range 10^{−2}–10^{−1} (McIntosh et al. 2011; McIntosh & De Pontieu 2012), and therefore we can neglect the heating from W^{(2)} in favour of W^{(0)} since we estimate η_{0} ≈ 10^{10}η_{2}. The goal now is first to calculate the damping rate for our wave in the case where the only dissipative contribution comes from W^{(0)}. We then compare this with the damping rate for the case where the only dissipative contribution comes from W^{(1)} and this can be used to estimate where W^{(1)} can be neglected in favour of W^{(0)}. The average wave energy on a field line is given by
The average heating rate for the case where the only dissipative contribution comes from W^{(0)} along a field line of length L is given by
to leading order. Therefore, the damping rate is given by
to leading order, denoted with the symbol γ_{} as the heating is only dependent on gradients parallel to the magnetic field. We note that
and therefore to find the damping rate of our wave in the case where the only dissipative contribution comes from W^{(1)} we simply need to use standard phase mixing results (Heyvaerts & Priest 1983). In Sect. 3 we show that the damping rate of a phase mixed Alfvén waves can be closely approximated by
where ∇_{⊥} denotes the gradient in Alfvén speed perpendicular to the field. Both damping rates are plotted as a function of frequency, f, (where k = 2πf/v_{A}) in Fig. 2. To make this plot, the following parameters were used: η_{0}/ρ = 10^{10} m^{2} s^{−1}, η_{1}/ρ = 1 m^{2} s^{−1}, v_{A} = 1 Mm s^{−1} and ∇_{⊥}v_{A} = 1 s^{−1}. It can be seen that γ_{ph} ≫ γ_{} for f ≤ 10^{−1} Hz and so in this parameter space we can neglect W^{(0)} in favour of W^{(1)}. It is worth noting that in Eq. (10) we have neglected the possibility of mode coupling from Alfvén waves to magnetoacoustic waves and pondermotive wings (Verwichte et al. 1999) and this could enhance the importance of W^{(0)}. If
Fig. 2. Estimated damping rates for the case where the only dissipative contribution comes from W^{(0)}, denoted with γ_{}, and for the case where the only dissipative contribution comes from W^{(1)}, denoted with γ_{ph}. 
and
then
which shows that, in general, the W^{(0)} term decays magnetoacoustic waves more efficiently than it does Alfvén waves as it is proportional to (b/B_{0})^{0} rather than (b/B_{0})^{2}. However, considering the nonlinear mode coupling of Alfvén waves to other wave modes is a sufficiently complex problem that we believe it is best left as a standalone paper. To simplify the notation and keep it in accordance with that of Priest (2014) we take the viscous force to be given by η_{1}∇ ⋅ W^{(1)}, consider only the leading order term, and set
In the induction Eq. (11) we neglect the diffusion term involving derivatives parallel to the background field. We justify this because we expect the perpendicular gradients to be much stronger than the parallel gradients after the wave has phase mixed. According to Braginskii (1965) the parallel and perpendicular conductivities differ only by approximately a factor of two. We expect the gradients perpendicular to the field to be many times greater than parallel gradients as the waves phase mix. Therefore, we only consider the perpendicular gradients and set
where σ_{} is the conductivity parallel to the field.
Important conditions (Braginskii 1965) for the validity of the transport coefficients are
where L_{⊥} and L_{} are the characteristic distances in the directions perpendicular and parallel to the magnetic field, r_{T, e} is the thermal electron gyro radius and l_{T, p}, is the proton mean free path. Assuming full ionisation and a purely hydrogen corona, the electron gyro radius is given by
The proton mean free path can be written as
Priest (2014) shows that the length scale of phasemixed Alfvén waves perpendicular to the field lines, , can be approximated by
Therefore, for our parameter space, the conditions in Eq. (31) are usually satisfied if the wavelength, λ, is greater than 100 km with the first condition being easier to satisfy than the second.
3. Closed loop
Heyvaerts & Priest (1983) and Priest (2014) find the solution for an Alfvén wave on an open field line by assuming a solution of the form,
and assume V is weakly varying in z which means ∂V/∂z ≪ k_{z}V. They find the solution to be
where ∇_{⊥}v_{A} = dv_{A}/dx denotes the gradient in Aflvén speed perpendicular to the field. This solution has a characteristic damping length, l_{ph}, given by
This gives a timescale, , given by
The goal of this section is to assess whether the phasemixed damping rate, γ_{ph}, as defined in Eq. (1) can be approximated with , given by
To do this, we extend the open field solution to a closedloop of Heyvaerts & Priest (1983), as in Fig. 1. We note that Heyvaerts & Priest (1983) have already derived a solution for a closed loop, but we believe the form of our solution to be more useful for our problem. We derive our formula via a different approach to that of Heyvaerts & Priest (1983); they make use of Green’s function whereas we use a method of images approach. In Sect. 3.1 we calculate a solution analytically for the damping rate at steady state. We then verify our calculation of the Poynting flux numerically in Sect. 3.2. Finally, in Sect. 3.3 we discuss how well Eq. (37) approximates the damping rate.
3.1. Analytic solution
The goal here is to extend Eqs. (33) and (34) for a confined domain, as in Fig. 1 (but with the complex exponential replaced with a sine function). We consider a domain in which a factor, R, of the wave amplitude reflects at the transition region boundary, where the boundaries are given by
We note that R < 1 denotes the amplitude reflection coefficient which is not (in general) the same as the energy reflection coefficient, R_{E}. We solve the problem using a method of images approach, and the full solution for u is
where H() denotes the Heaviside step function,
Here, ⌊⌋ denote the floor function, namely the integer part to the real number. Using the trig identity,
the steadystate solution, that is the solution for t → ∞ to Eq. (38), is given by
where A is given by
and B is given by
The solution for the magnetic field perturbation, b, is very similar to that of the velocity and is,
At steady state, this simplifies to
where C is given by
and D is given by
The Poynting flux on the boundary can be simplified using Ohm’s law and neglecting the small resistive terms gives
Using Eqs. (43) and (47) the average steadystate Poynting flux, −B_{0}⟨ub⟩/μ, is given by,
The average steadystate kinetic wave energy density, ρ⟨u^{2}⟩/2, can be expressed as
Therefore, the average steadystate kinetic wave energy, E_{k}, for a field line is
Therefore, the damping rate, γ, for a field line is
We note that for R = 1 this simplifies to
3.2. Numerical solution
The purpose of this section is to demonstrate that Eqs. (43) and (47) do indeed give accurate solutions provided the damping is weak enough.
The numerical solution is found using the Lare2D code (Arber et al. 2001). Equations (43) and (47) are checked for R = 1, as this is the easiest boundary condition to impose numerically. The numerical domain is square and is given by,
In the numerical experiments the Alfvén speed, v_{A}, is given by
where L = 2l and v_{A0} = v_{A}(0). This was chosen as it is the simplest Alfvén speed profile with a nonzero gradient. The background magnetic field, B_{0}, is uniform in the z direction. The plasma is initially static. A driver is imposed at the z = −l boundary and has the form
where ω is given by
This driver excites resonance at only one location: the orange field line in Fig. 1, where the fundamental harmonic is excited. In Sect. 5 we investigate the effect of using a driver which excites multiple harmonics. Reflective or solid boundary conditions (Laney 1998) are otherwise imposed on all the boundaries. In other words, u = 0 and for all other variables, where is a vector normal to the boundary. The boundaries at z = ±l are designed to simulate the interaction of waves with the transition region. However, perfect reflection is only an approximation; in Sect. 4, we investigate the effect of using a partially reflective boundary. We note that the code uses isotropic incompressible viscous heating, ρν∇u : ∇u and isotropic Ohmic heating, μηj^{2}, where j is current density.
Plots of the average steadystate Poynting flux are given in Fig. 3. The Poynting flux peaks at the middle as it is the fundamental harmonic of the middle field line, which is excited. As expected, the solutions show better agreement for smaller values of η + ν (provided the resolution is high enough) as this means the weak damping assumption from Heyvaerts & Priest (1983) is then valid.
Fig. 3. Numerical and analytic steadystate average Poynting flux for each field line (given by x = a constant) for different values of η + ν. Res refers to the resolution used, i.e. the grid size. 
3.3. Damping rate
Figure 4 shows plots of the spatially integrated Ohmic and viscous heating, dE_{η + ν}/dt, given by
Fig. 4. Plots from a numerical experiment showing the rate of change of total Ohmic and viscous energy, dE_{η + ν}/dt, total wave energy, E_{wave}, and their ratio. Each plot has been normalised by its respective steadystate value and t_{0} = L/v_{A0}. 
It also shows the wave energy, E_{wave}, where
Finally, Fig. 4 shows the ratio, (dE_{η + ν}/dt)/E_{wave}; this ratio gives the damping rate, γ, for our system as a function of time. The figure was produced using data from the numerical experiment described in Sect. 3.2, with η + ν = 2^{−15}η_{norm} and a resolution of 512 × 512, where η_{norm} = 2lv_{A0}. The key result is that the damping rate increases with time and then converges towards a maximum value at steady state. This means that the damping rate at steady state represents an upper bound.
Figure 5 shows the steadystate damping rate for a field line as a function of frequency. The black curve in Fig. 5 was made using Eq. (55) with the following parameters: η + ν = 1 m^{2} s^{−1}, v_{A} = 1 Mm s^{−1}, ∇_{⊥}v_{A} = 1 s^{−1}, and L = 100 Mm and gives a fundamental frequency of f_{1} = 5 × 10^{−3} s^{−1}. A key result is that the damping rate is significantly larger at resonance and then sharply approaches zero away from resonance (even becoming negative), and therefore the resonant damping rate represents an upper bound. Figure 5 also shows (in red) our approximation for the damping rate, namely Eq. (37), which is a good approximation for the natural frequencies. However, there is a noticeable error of less than approximately 10%. The green curve in Fig. 5 was produced using Eq. (55), except the cubic exponential in Eqs. (44), (45), (48), and (49) were replaced with a linear exponential, namely
Fig. 5. Steadystate damping rate, γ, for a field line as a function of the driver frequency. The red curve shows an approximation for the damping rate. The righthand figure magnifies the black curve in a smaller frequency range to show that the curve is indeed continuous but very steep. 
The green curve shows that Eq. (37) is better at predicting the damping rate when there is a linear exponential compared with a cubic exponential. Therefore, the cubic exponential nature of phase mixing reduces the accuracy of Eq. (37) in predicting the damping rate. It is also interesting to note the fact that the damping rate approaches zero (and becomes negative) away from resonance, and this is a result of the cubic nature of phase mixing. If a linear exponential is used, the curve does not approach zero (see green curve) away from resonance.
4. Leaky loop
Here, we aim to explain why leakage decreases the damping rate and by how much. To do this, we first calculate an expression for the reflection coefficient of Alfvén waves at the transition region. After that, we analyse how leakage affects the damping rate.
Hollweg (1984b) derives an expression which gives an approximation for the energy reflection coefficient, R_{E}, for Alfvén waves at the transition region. To derive this, this latter author modelled the corona as having a uniform Alfvén speed and assumed the Alfvén speed in the chromosphere varies exponentially in the corona with pressure scale height, h. The formula he obtains is
where,
where v_{A} is the Alfvén speed in the corona, denotes the Hankel function of the second kind of zero order and is the Hankel function of the second kind of first order. We note that the energy transmission coefficient is given by T_{E} = 1 − R_{E}. The amplitude reflection coefficient is given by . We modify the Hollweg (1984b) model by including a discontinuous jump in Alfvén speed from the chromosphere to the corona, which is designed to model the transition region. The behaviour of the wave interaction with the transition region approximates that of a discontinuity provided the wavelength of the wave is greater than the width of the transition region. The Alfvén speed increases by a factor, a, at the discontinuity and we find that the inclusion of this discontinuity causes the energy reflection coefficient to be instead given by
where,
and this is derived in Appendix A. The reflection coefficient is plotted in Fig. 6 as a function of frequency. For ξ ≫ 1, the equation for the reflection coefficient reduces to the following form,
Fig. 6. Energy reflection coefficient, R_{E}, as function of frequency, f. 
Expanding Eq. (64) about ξ = 0 gives, to leading order,
The leakage timescale, τ_{leakage}, which is the timescale at which wave energy in the corona is lost through leakage into the chromosphere, is given by
This equation can be derived by considering a partially confined wave; its amplitude, u_{0}, behaves as
where ⌊⌋ denotes the floor function, which takes the integer part of the input. In other words, a factor R of the wave amplitude will be lost each time it partially reflects off either end of the loop. In Sect. 3.3 we showed that the damping rate of a resonant field line can be approximated by Eq. (37), and therefore the resistive timescale is given by
Both timescales are plotted in Fig. 7, where the following parameters were used: η + ν = 1 m^{2} s^{−1}, v_{A} = 1 Mm s^{−1}, ∇_{⊥}v_{A} = 1 s^{−1}, L = 100 Mm, and h = 150 km, with R_{E} being given by Eq. (64). It can be seen that for these parameters the leakage timescale is shorter than the resistive timescale. Therefore, leakage could play an important role in the dynamics. However, as shown in the introduction, we estimate that a viable heating mechanism needs to have a damping rate of about 10 s^{−1} which gives a timescale that is much shorter than the leakage timescale
Fig. 7. Leakage timescale, τ_{leakagae}, (see Eq. (68)) and the resistive timescale, τ_{resistive}, (see Eq. (70)). 
Leakage acts to reduce the wave energy, but it also acts to reduce the heating rate and so it is perhaps not clear how leakage will affect the damping rate. The damping rate is plotted as a function of the reflection coefficient in Fig. 8 for a resonant field line. This was produced using Eq. (54) with the following parameters: η + ν = 1 m^{2} s^{−1}, v_{A} = 1 Mm s^{−1}, ∇_{⊥}v_{A} = 1 s^{−1}, and f = 10^{−2} Hz. It can be seen that increasing the leakage, − log R, causes the damping rate, γ, to decrease. Therefore, leakage has the effect of reducing the damping rate. It can be seen that for the damping rate is largely independent of the reflection coefficient and can be approximated with Eq. (37). For the damping rate can be approximated by the blue dashed curve, which has the following equation,
Fig. 8. Damping rate, γ, as a function of the leakage, − log R. 
which is derived in Appendix B. The horizontal red dotted line in Fig. 8 was produced by replacing the cubic exponentials with linear exponentials (as in Sect. 3.3); it can be seen to be independent of leakage. This suggests that it is because of the cubic nature of phase mixing that leakage causes the damping rate to decrease. More intuitively, the leakage prevents the waves phasemixing down to very short length scales before they leak out of the loop. We have marked in magenta the value of the reflection coefficient for f = 10^{−2} Hz and h = 150 km. This shows that for this parameter space, the damping rate can be approximated by Eq. (71).
5. Multiple harmonics
As shown in the literature, for example by Morton et al. (2016, 2019), coronal loops do not oscillate at one frequency. Instead, they oscillate at a spectrum of frequencies. These latter authors show that velocity fluctuations in the corona approximately follow a power law with a bump at approximately fiveminute periods which roughly correspond to the pmode frequencies. For simplicity, we assume that our velocity fluctuations obey a single power law. Additionally, we drive with a series of sinusoidal drivers at the natural frequencies because a broadband driver drives resonances in a manner that is similar to a superposition of monochromatic drivers (Wright & Rickard 1995).
In Sect. 3.1, we show that for a driver of the form
the steadystate solution for u is given by Eq. (43). Since we are assuming linear waves, we can easily find the full solution for a driver consisting of many harmonics,
which gives,
where ω_{1} is the fundamental angular frequency, ϕ_{n} is a random phase, α is a constant which controls the power spectrum of the driver,
This solution has the property that
which can be shown by the fact that,
for n ≤ m ∈ ℕ. This means we can easily control the power spectrum of the steadystate kinetic energy of the system with α. The average Poynting flux, −B_{0}⟨ub⟩/μ, is given by
The damping rate, γ, is given by
where γ_{n} is given by
Equation (82) takes the form of a weighted average of the set {γ_{n}} with weights {n^{−α}}. The damping rate is plotted in Fig. 9 as a function of the number of excited harmonics for different values of α using the following parameter values: η + ν = 1 m^{2} s^{−1}, v_{A} = 1 Mm s^{−1}, ∇_{⊥}v_{A} = 1 s^{−1}, and L = 100 Mm, corresponding to a fundamental harmonic of f_{1} = 5 × 10^{−3} Hz. As expected, increasing the number of harmonics increases the damping rate. Using Eq. (37) as an approximation for γ_{n} and replacing the summation with an integral gives
Fig. 9. Damping rate, γ, of a resonant field line as a function of the number of excited harmonics. The approximate curves are produced by replacing the summation in Eq. (82) with an integral. 
for α ≠ 5/3 and α ≠ 1. For α = 5/3, it is approximated by
and for α = 1 it is approximated by
These approximations are also plotted as dashed curves in Fig. 9. The reason there is a noticeable error between the exact and approximate solutions in Fig. 9 is twofold: first, the approximate solution approximates γ_{n} using Eq. (37), and second, the summations in Eq. (82) have been replaced with an integral. It can be seen that for larger α, increasing N does relatively little to change the damping rate. Figure 9 uses a range of α values including 1.5 and 5/3 as these are the values predicted from MHD turbulence theory (Bruno & Carbone 2016). Morton et al. (2016) provides observations of the power spectra of velocity fluctuations in the quiet sun, the active regions, and the coronal holes. These latter authors find that the slope varies from α = 1 to α = 1.53 for higher frequencies, although they are only able to measure up to frequencies of around 10^{−2} Hz. Podesta et al. (2007) measure the power spectra in the solar wind and can measure up to 10^{−1} Hz and find the slope to be between α = 1.5 and α = 5/3. From Fig. 2, it can be seen that for the higher frequencies, the heating due to parallel gradients start to dominate. Therefore, if we use a value of N of higher than 100, we can see that parallel gradients will start to dominate the heating. Therefore, we can no longer describe the system as being heated mainly by phase mixing.
6. Conclusions
We believe an upper bound for the damping rate of laminar phasemixed Alfvén waves oscillating with a given power spectrum is given by Eq. (82) with γ_{n} approximated by (37), provided our assumptions (see Sect. 2) are valid. Here, we provide a brief argument for why we believe this to be an upper bound: In Sect. 3.3, we show that the damping rate increases with time but converges towards a maximum as the system approaches steady state (the blue curve illustrates this in Fig. 4) and our equations were derived by assuming the system had reached a steady state. This is important because (as pointed out by Arregui 2015) waves may not have a chance to reach steady state. Also, we show that the damping is largest at the natural frequencies (illustrated in Fig. 5) and our equations were derived by only exciting the resonant frequencies with the driver. In Sect. 4, we show that allowing waves to leak through the atmosphere decreases the damping rate and Eq. (37) was derived assuming perfect reflection at the transition region. We have not considered the thermodynamics, but Cargill et al. (2016) showed that the thermodynamic response due to heating reduces the density gradients, which reduces the damping rate. As stated in Sect. 2, we solve the linearised MHD equations. Prokopyszyn et al. (2019) found that provided there was no turbulence (which is ensured by the presence of an ignorable coordinate and the fact that MHD turbulence of pure Alfvén waves is a strictly 3D phenomenon (see Howes & Nielson 2013) then for u/v_{A} ⪅ 0.1 the nonlinearities have a negligible effect on the damping rate. One problem with the results from Prokopyszyn et al. (2019) is that an unphysically large value for the dissipative coefficients had to be used due to numerical constraints. Nonlinearities can trigger turbulence, for example through the interaction of counterpropagating Alfvén waves (Hollweg 1986; van Ballegooijen et al. 2011; Shoda et al. 2019) or via the KelvinHelmholtz/tearing mode instability (Heyvaerts & Priest 1983; Browning & Priest 1984). Turbulence leads to the transfer of energy into higher wavenumbers/frequencies, which causes the damping rate to increase. However, our claim is that Eq. (82) is an upper bound for laminar waves. Finally, Threlfall et al. (2011) found that Hall MHD terms produce wave dispersion and reduce the damping rate.
Provided Eq. (82) is a valid upper bound for the damping rate, this implies that phase mixed Alfvén waves at observed frequencies are unlikely to play a role in coronal heating. This can be seen from Fig. 9 which predicts a damping rate of no more than 10^{−4} s^{−1}. We believe a damping rate of the order 10^{−1} s^{−1} is required to heat the corona (as shown in the introduction). Changing our parameters will change the damping rate. The damping rate may be high enough in some locations of the coronal atmosphere, for example near null points where the crossfield viscosity becomes much stronger. In this study it is assumed that the flow remains laminar. Browning & Priest (1984) show that a phasemixed standing Alfvén wave will trigger the KelvinHelmholtz instability at its antinode, which could lead to a turbulent cascade, and these latter authors calculated that this could lead to sufficient heating in the corona. They find the instability is triggered at the antinodes of the wave as this is where the magnetic field is smallest. In theory, multiple harmonics should help to stabilise the field as the magnetic field perturbations should be more uniformly distributed across the field line if this is the case. Future work could investigate the effects of multiple harmonics on the phasemixinginduced KelvinHelmholtz instability.
Acknowledgments
We thank Andrew Wright and Ineke De Moortel for help with checking through the paper and useful discussions. This project has received funding from the Science and Technology Facilities Council (UK) through the consolidated grant ST/N000609/1.
References
 Antolin, P., & Shibata, K. 2010, ApJ, 712, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Antolin, P., De Moortel, I., Van Doorsselaere, T., & Yokoyama, T. 2016, ApJ, 830, L22 [Google Scholar]
 Antolin, P., Schmit, D., Pereira, T. M. D., De Pontieu, B., & De Moortel, I. 2018, ApJ, 856, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, J. Comput. Phys., 171, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Arregui, I. 2015, Trans. R. Soc. London Ser. A, 373, 20140261 [Google Scholar]
 Braginskii, S. 1965, Rev. Plasma Phys., 1, 205 [NASA ADS] [Google Scholar]
 Browning, P. K., & Priest, E. R. 1984, A&A, 131, 283 [NASA ADS] [Google Scholar]
 Bruno, R., & Carbone, V. 2016, in Turbulence in the Solar Wind (Berlin: Springer Verlag), Lect. Notes Phys., 928 [CrossRef] [Google Scholar]
 Cally, P. S., & Hansen, S. C. 2011, ApJ, 738, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Cargill, P. J., De Moortel, I., & Kiddie, G. 2016, ApJ, 823, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Cranmer, S. R. 2018, ApJ, 862, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265 [NASA ADS] [CrossRef] [Google Scholar]
 De Moortel, I., & Browning, P. 2015, Trans. R. Soc. London Ser. A, 373, 20140269 [NASA ADS] [CrossRef] [Google Scholar]
 De Moortel, I., & Nakariakov, V. M. 2012, Trans. R. Soc. London Ser. A, 370, 3193 [Google Scholar]
 Goossens, M., Erdélyi, R., & Ruderman, M. S. 2011, Space Sci. Rev., 158, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, S. C., & Cally, P. S. 2012, ApJ, 751, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Heyvaerts, J., & Priest, E. R. 1983, A&A, 117, 220 [NASA ADS] [Google Scholar]
 Hollweg, J. V. 1984a, Sol. Phys., 91, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Hollweg, J. V. 1984b, ApJ, 277, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Hollweg, J. V. 1986, J. Geophys. Res., 91, 4111 [Google Scholar]
 Hood, A. W., Ireland, J., & Priest, E. R. 1997, A&A, 318, 957 [NASA ADS] [Google Scholar]
 Hood, A. W., Brooks, S. J., & Wright, A. N. 2002, Trans. R. Soc. London Ser. A, 458, 2307 [NASA ADS] [CrossRef] [Google Scholar]
 Howes, G. G., & Nielson, K. D. 2013, Phys. Plasmas, 20, 072302 [NASA ADS] [CrossRef] [Google Scholar]
 Howson, T. A., De Moortel, I., Reid, J., & Hood, A. W. 2019, A&A, 629, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ionson, J. A. 1982, ApJ, 254, 318 [NASA ADS] [CrossRef] [Google Scholar]
 Klimchuk, J. A. 2006, Sol. Phys., 234, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Klimchuk, J. A. 2015, Trans. R. Soc. London Ser. A, 373, 20140256 [NASA ADS] [Google Scholar]
 Kudoh, T., & Shibata, K. 1999, ApJ, 514, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Laney, C. B. 1998, Computational Gasdynamics (New York: Cambridge University Press) [CrossRef] [Google Scholar]
 Leroy, B. 1980, A&A, 91, 136 [NASA ADS] [Google Scholar]
 MacTaggart, D., Vergori, L., & Quinn, J. 2017, J. Fluid Mech., 826, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Magyar, N., Van Doorsselaere, T., & Goossens, M. 2017, Sci. Rep., 7, 14820 [NASA ADS] [CrossRef] [Google Scholar]
 McIntosh, S. W., & De Pontieu, B. 2012, ApJ, 761, 138 [NASA ADS] [CrossRef] [Google Scholar]
 McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Morton, R. J., Tomczyk, S., & Pinto, R. F. 2016, ApJ, 828, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Morton, R. J., Weberg, M. J., & McLaughlin, J. A. 2019, Nat. Astron., 3, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Pagano, P., & De Moortel, I. 2017, A&A, 601, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pagano, P., Pascoe, D. J., & De Moortel, I. 2018, A&A, 616, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parker, E. N. 1991, ApJ, 376, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Parnell, C. E., & De Moortel, I. 2012, Trans. R. Soc. London Ser. A, 370, 3217 [Google Scholar]
 Podesta, J. J., Roberts, D. A., & Goldstein, M. L. 2007, ApJ, 664, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Priest, E. 2014, Magnetohydrodynamics of the Sun (Cambridge: Cambridge University Press) [Google Scholar]
 Prokopyszyn, A. P. K., Hood, A. W., & De Moortel, I. 2019, A&A, 624, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Russell, A. J. B. 2019, unpublished [Google Scholar]
 Shoda, M., Suzuki, T. K., AsgariTarghi, M., & Yokoyama, T. 2019, ApJ, 880, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Similon, P. L., & Sudan, R. N. 1989, ApJ, 336, 442 [NASA ADS] [CrossRef] [Google Scholar]
 Threlfall, J., McClements, K. G., & De Moortel, I. 2011, A&A, 525, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 van Ballegooijen, A. A., AsgariTarghi, M., Cranmer, S. R., & DeLuca, E. E. 2011, ApJ, 736, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Verwichte, E., Nakariakov, V. M., & Longbottom, A. W. 1999, J. Plasma Phys., 62, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, A. N., & Rickard, G. J. 1995, ApJ, 444, 458 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Reflection coefficient
Hollweg (1984b) derived an expression for the reflection coefficient (Eq. (62)). He derived this by splitting his domain into two regions, z < 0 corresponds to the chromosphere with an exponentially growing Alfvén speed with pressure scale height, h, and z > 0 corresponds to the corona with a uniform Alfvén speed. By considering a source of waves in z > 0 he calculates that in general, the waves have the form
where C_{1}, C_{2}, C_{3} are arbitrary constants and ξ is given by Eq. (63). We modify these equations slightly because we assume there is a discontinuous jump in Alfvén speed of a factor a from the chromosphere to the corona due to the jump in density. Taking v_{A} as the Alfvén speed in the corona and v_{A}/a as the speed at the top of the chromosphere then this gives the same expression for u but with b, now given by
and ξ given by (65). Despite v_{A} being discontinuous, it can be shown from Farday’s, Ampère’s, and Ohm’s law that u and b must be continuous (Hollweg 1984a). The continuity conditions can be used to eliminate C_{1} and the reflection coefficient, R, is given by C_{2}/C_{3}. Hence,
Appendix B: γ in the high leakage limit
In the high leakage limit, Eq. (38) reduces to
Letting t → ∞, considering only the apex of the loop (z = 0), and replacing the sine function with a complex exponential gives
This equation takes the form of a geometric series which can be evaluated, provided R < 1, to give
where L = 2l. The imaginary part of Eq. (B.3) has maxima if odd harmonics are excited, i.e. ω = (2n + 1)πv_{A}/L, n ∈ ℕ then exp(−iωL/v_{A}) = − 1. The gradient of Eq. (B.3) perpendicular to the field for ω = (2n + 1)πv_{A}/L is given by
Therefore, the damping rate can be approximated with ω = (2n + 1)πv_{A}/L to give
All Figures
Fig. 1. Diagram of our model, showing a set of closed coronal loops (vertical black lines) with a driver imposed at the bottom boundary (blue). The gradient in height of the blue arrows indicates there is a gradient in Alfvén speed. 

In the text 
Fig. 2. Estimated damping rates for the case where the only dissipative contribution comes from W^{(0)}, denoted with γ_{}, and for the case where the only dissipative contribution comes from W^{(1)}, denoted with γ_{ph}. 

In the text 
Fig. 3. Numerical and analytic steadystate average Poynting flux for each field line (given by x = a constant) for different values of η + ν. Res refers to the resolution used, i.e. the grid size. 

In the text 
Fig. 4. Plots from a numerical experiment showing the rate of change of total Ohmic and viscous energy, dE_{η + ν}/dt, total wave energy, E_{wave}, and their ratio. Each plot has been normalised by its respective steadystate value and t_{0} = L/v_{A0}. 

In the text 
Fig. 5. Steadystate damping rate, γ, for a field line as a function of the driver frequency. The red curve shows an approximation for the damping rate. The righthand figure magnifies the black curve in a smaller frequency range to show that the curve is indeed continuous but very steep. 

In the text 
Fig. 6. Energy reflection coefficient, R_{E}, as function of frequency, f. 

In the text 
Fig. 7. Leakage timescale, τ_{leakagae}, (see Eq. (68)) and the resistive timescale, τ_{resistive}, (see Eq. (70)). 

In the text 
Fig. 8. Damping rate, γ, as a function of the leakage, − log R. 

In the text 
Fig. 9. Damping rate, γ, of a resonant field line as a function of the number of excited harmonics. The approximate curves are produced by replacing the summation in Eq. (82) with an integral. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.