© 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 build-up 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), reflection-driven Alfvén wave turbulence (Hollweg 1986; van Ballegooijen et al. 2011; Shoda et al. 2019), turbulence triggered via the tearing mode or Kelvin-Helmholtz 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 Kelvin-Helmholtz 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 self-similar 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

(1)

where ⟨⟩ denotes the time average of either E × B/μ, the Poynting flux, or ρu²/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

(2)

hence,

(3)

and so the heating rate is given by . Therefore,

(4)

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⁻¹ s⁻¹ 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

(5)

where T is temperature and χ = 10⁻³² K^1/2 W m³ and a uniform coronal heating profile, H_c, then Priest (2014) shows that to maintain a loop with a uniform pressure, p₀, requires a heating rate of

(6)

where k_B is the Boltzmann constant and T_max is the maximum temperature. Here, T_max = 10⁶ K and p₀ = 10⁻² Pa gives H_c ≈ 10⁻⁵ W m⁻³, 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⁻¹ and in active regions of around 5 km s⁻¹. If we assume an amplitude of u₀ = 10 km s⁻¹ and a density of ρ₀ = 10⁻¹² km s⁻¹, this gives a kinetic energy density of J m⁻³/2. Therefore, the required damping rate is approximately s⁻¹.

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 closed-loop. 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 stand-alone 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,

(7)

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,

(8)

and the density,

(9)

The equations are the linearised momentum equation,

(10)

and the linearised induction equation,

(11)

where the permeability of free space is denoted with μ, the viscosity coefficient, ν, is given by

(12)

and the magnetic diffusivity coefficient, η, is given by

(13)

where ω_p is the proton cyclotron frequency given by

(14)

with e the elementary charge and m_p the mass of a proton. Here, τ_p is the proton-proton collision time, given by

(15)

where logΛ is the Coulomb logarithm ≈22 in coronal conditions (Priest 2014). Substituting T = 10⁶ K, B = 10⁻³ T and ρ = 10⁻¹² m⁻³ gives ω_pτ_p ≈ 10⁵. 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

(16)

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 Kelvin-Helmholtz 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 field-aligned 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 non-uniform 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 shorter-wavelength waves form shorter perpendicular length scales more quickly. However, if our uniform Alfvén speed represents the average Alfvén speed of a non-uniform loop, then the time-averaged heating rates will not differ significantly. The second effect is that a non-uniform Alfvén speed can give rise to reflection. If the wavelength of a wave is shorter than the length-scale 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 semi-circular 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),

(17)

where

(18)

(19)

(20)

(21)

(22)

with

(23)

where Z is the tensor with components Z_ij = ϵ_ikjb_k, where ϵ_ikj are components of the Levi-Civita symbol. The viscosity coefficients (see Braginskii 1965) are given by

(24)

assuming the strong field limit (ω_pτ_p ≫ 1). We aim to show that η₀W⁽⁰⁾, η₂W⁽²⁾, η₃W⁽³⁾, η₄W⁽⁴⁾ can be neglected in favour of η₁W⁽¹⁾. This is surprising because if ω_pτ_p ≈ 10⁵ then this means that η₁ ≈ 10⁻¹⁰η₀.

Here we assume that u and b are footpoint-driven waves on closed loops that have reached steady state such that their amplitudes do not change with time and are of the form,

(25)

(26)

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:

(27)

(28)

(29)

(30)

We expect v/v_A to be approximately in the range 10⁻²–10⁻¹ (McIntosh et al. 2011; McIntosh & De Pontieu 2012), and therefore we can neglect the heating from W⁽²⁾ in favour of W⁽⁰⁾ since we estimate η₀ ≈ 10¹⁰η₂. The goal now is first to calculate the damping rate for our wave in the case where the only dissipative contribution comes from W⁽⁰⁾. We then compare this with the damping rate for the case where the only dissipative contribution comes from W⁽¹⁾ and this can be used to estimate where W⁽¹⁾ can be neglected in favour of W⁽⁰⁾. 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⁽⁰⁾ 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⁽¹⁾ 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: η₀/ρ = 10¹⁰ m² s⁻¹, η₁/ρ = 1 m² s⁻¹, v_A = 1 Mm s⁻¹ and ∇_⊥v_A = 1 s⁻¹. It can be seen that γ_ph ≫ γ_|| for f ≤ 10⁻¹ Hz and so in this parameter space we can neglect W⁽⁰⁾ in favour of W⁽¹⁾. 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⁽⁰⁾. 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⁽⁰⁾ term decays magnetoacoustic waves more efficiently than it does Alfvén waves as it is proportional to (b/B₀)⁰ rather than (b/B₀)². However, considering the non-linear mode coupling of Alfvén waves to other wave modes is a sufficiently complex problem that we believe it is best left as a stand-alone paper. To simplify the notation and keep it in accordance with that of Priest (2014) we take the viscous force to be given by η₁∇ ⋅ W⁽¹⁾, 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

(31)

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 phase-mixed Alfvén waves perpendicular to the field lines, , can be approximated by

(32)

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,

(33)

and assume V is weakly varying in z which means ∂V/∂z ≪ k_zV. They find the solution to be

(34)

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

(35)

This gives a timescale, , given by

(36)

The goal of this section is to assess whether the phase-mixed damping rate, γ_ph, as defined in Eq. (1) can be approximated with , given by

(37)

To do this, we extend the open field solution to a closed-loop 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

(38)

where H() denotes the Heaviside step function,

(39)

(40)

(41)

Here, ⌊⌋ denote the floor function, namely the integer part to the real number. Using the trig identity,

(42)

the steady-state solution, that is the solution for t → ∞ to Eq. (38), is given by

(43)

where A is given by

(44)

and B is given by

(45)

The solution for the magnetic field perturbation, b, is very similar to that of the velocity and is,

(46)

At steady state, this simplifies to

(47)

where C is given by

(48)

and D is given by

(49)

The Poynting flux on the boundary can be simplified using Ohm’s law and neglecting the small resistive terms gives

(50)

Using Eqs. (43) and (47) the average steady-state Poynting flux, −B₀⟨ub⟩/μ, is given by,

(51)

The average steady-state kinetic wave energy density, ρ⟨u²⟩/2, can be expressed as

(52)

Therefore, the average steady-state kinetic wave energy, E_k, for a field line is

(53)

Therefore, the damping rate, γ, for a field line is

(54)

We note that for R = 1 this simplifies to

(55)

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,

(56)

In the numerical experiments the Alfvén speed, v_A, is given by

(57)

where L = 2l and v_A0 = v_A(0). This was chosen as it is the simplest Alfvén speed profile with a non-zero gradient. The background magnetic field, B₀, is uniform in the z direction. The plasma is initially static. A driver is imposed at the z = −l boundary and has the form

(58)

where ω is given by

(59)

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², where j is current density.

Plots of the average steady-state 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 steady-state 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

(60)

Fig. 4. Plots from a numerical experiment showing the rate of change of total Ohmic and viscous energy, dEη + ν/dt, total wave energy, Ewave, and their ratio. Each plot has been normalised by its respective steady-state value and t0 = L/vA0.

It also shows the wave energy, E_wave, where

(61)

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⁻¹⁵η_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 steady-state 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² s⁻¹, v_A = 1 Mm s⁻¹, ∇_⊥v_A = 1 s⁻¹, and L = 100 Mm and gives a fundamental frequency of f₁ = 5 × 10⁻³ s⁻¹. 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. Steady-state damping rate, γ, for a field line as a function of the driver frequency. The red curve shows an approximation for the damping rate. The right-hand 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

(62)

where,

(63)

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

(64)

where,

(65)

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,

(66)

Fig. 6. Energy reflection coefficient, RE, as function of frequency, f.

Expanding Eq. (64) about ξ = 0 gives, to leading order,

(67)

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

(68)

This equation can be derived by considering a partially confined wave; its amplitude, u₀, behaves as

(69)

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

(70)

Both timescales are plotted in Fig. 7, where the following parameters were used: η + ν = 1 m² s⁻¹, v_A = 1 Mm s⁻¹, ∇_⊥v_A = 1 s⁻¹, 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⁻¹ 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² s⁻¹, v_A = 1 Mm s⁻¹, ∇_⊥v_A = 1 s⁻¹, and f = 10⁻² 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,

(71)

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 phase-mixing 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⁻² 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 five-minute periods which roughly correspond to the p-mode 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 steady-state 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,

(72)

which gives,

(73)

(74)

where ω₁ is the fundamental angular frequency, ϕ_n is a random phase, α is a constant which controls the power spectrum of the driver,

(75)

(76)

(77)

(78)

This solution has the property that

(79)

which can be shown by the fact that,

(80)

for n ≤ m ∈ ℕ. This means we can easily control the power spectrum of the steady-state kinetic energy of the system with α. The average Poynting flux, −B₀⟨ub⟩/μ, is given by

(81)

The damping rate, γ, is given by

(82)

where γ_n is given by

(83)

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² s⁻¹, v_A = 1 Mm s⁻¹, ∇_⊥v_A = 1 s⁻¹, and L = 100 Mm, corresponding to a fundamental harmonic of f₁ = 5 × 10⁻³ 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

(84)

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

(85)

and for α = 1 it is approximated by

(86)

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⁻² Hz. Podesta et al. (2007) measure the power spectra in the solar wind and can measure up to 10⁻¹ 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 phase-mixed 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 non-linearities 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. Non-linearities can trigger turbulence, for example through the interaction of counter-propagating Alfvén waves (Hollweg 1986; van Ballegooijen et al. 2011; Shoda et al. 2019) or via the Kelvin-Helmholtz/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⁻⁴ s⁻¹. We believe a damping rate of the order 10⁻¹ s⁻¹ 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 cross-field viscosity becomes much stronger. In this study it is assumed that the flow remains laminar. Browning & Priest (1984) show that a phase-mixed standing Alfvén wave will trigger the Kelvin-Helmholtz 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 phase-mixing-induced Kelvin-Helmholtz 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

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

(A.1)

(A.2)

where C₁, C₂, C₃ 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

(A.3)

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₁ and the reflection coefficient, R, is given by C₂/C₃. Hence,

(A.4)

Appendix B: γ in the high leakage limit

In the high leakage limit, Eq. (38) reduces to

(B.1)

Letting t → ∞, considering only the apex of the loop (z = 0), and replacing the sine function with a complex exponential gives

(B.2)

This equation takes the form of a geometric series which can be evaluated, provided R <  1, to give

(B.3)

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

(B.4)

Therefore, the damping rate can be approximated with ω = (2n + 1)πv_A/L to give

(B.5)

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 steady-state 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, Ewave, and their ratio. Each plot has been normalised by its respective steady-state value and t0 = L/vA0.
In the text

	Fig. 5. Steady-state damping rate, γ, for a field line as a function of the driver frequency. The red curve shows an approximation for the damping rate. The right-hand 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, RE, 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