© S. J. White and E. Verwichte 2021

1. Introduction

In the last few decades, with the development of high spatial resolution telescopes in the extreme ultraviolet and x-ray regime, we have been able to resolve a variety of magnetohydrodynamic (MHD) waves in structures in the solar corona (De Moortel & Nakariakov 2012; Liu & Ofman 2014; Roberts 2019). This includes the observation of transverse loop oscillations, which were first observed in connection with eruption events by TRACE (Aschwanden et al. 1999; Nakariakov et al. 1999) and have been routinely observed since (e.g., Verwichte et al. 2009; Van Doorsselaere et al. 2008; White et al. 2012; Verwichte & Kohutova 2017). MHD wave theory reveals these oscillations to be fast magnetoacoustic in nature that exhibit periodic coherent displacements of the loop axis in a quasi-incompressible manner (Edwin & Roberts 1983; Nakariakov & Ofman 2001; Aschwanden et al. 2002; Verwichte et al. 2004). Due to their resemblance to the shear Alfvén waves on a magnetic field line, this type of oscillation has been dubbed Alfvénic. Importantly, transverse oscillations are also studied in connection with the coronal heating problem, as they naturally transfer energy from global coherent to local incoherent motions through the process of resonant absorption (Ruderman & Roberts 2002; Goossens et al. 2002; Arregui et al. 2008; Van Doorsselaere et al. 2020). Transverse waves provide a means to estimate local coronal plasma parameters such as the magnetic field and Alfvén speed (Nakariakov & Ofman 2001; Verwichte et al. 2013; Yuan & Van Doorsselaere 2016; Li et al. 2017), the study of which is known as coronal seismology (Nakariakov & Verwichte 2005).

Observations also show evidence of acoustic (slow magnetoacoustic) waves in the solar corona, that is, as standing waves in flaring loops (Wang et al. 2003a,b) and as propagating waves propagating upwards from loop points of large loops (Berghmans et al. 1999; De Moortel et al. 2000; Robbrecht et al. 2001; De Moortel 2009) which also has applications in the field of coronal seismology (see Wang et al. 2021, for a review). Acoustic waves are usually identified from their intensity perturbations (running or standing) due to periodic compressions and rarefactions, or from periodic Doppler velocity signatures aligned with the loop axis (and magnetic field) (Wang et al. 2003b). In most cases, transverse and acoustic oscillations are interpreted independently from each other. Alfvénic modes tend to be modelled in the zero plasma-β limit where there is no longitudinal velocity component and slow magnetoacoustic waves are absent. This low plasma-β uncoupled model of the MHD modes is not always justified, such as in large coronal loops (Verwichte et al. 2010) and hot, dense flaring loops. Furthermore, despite the quasi-incompressible nature of Alfvénic oscillations, observations of transverse loop oscillations often reveal accompanying variations in intensity with the same periodicity that have been attributed to the compressive nature of vertically polarised transverse loop oscillations in curved loops (Wang & Solanki 2004; Verwichte et al. 2006; Díaz 2006), line-of-sight integration effects (Verwichte et al. 2009; White & Verwichte 2012) and wave coupling at loop footpoints (Terradas et al. 2011).

Goedbloed & Halberstadt (1994) showed that MHD waves in magnetic flux tubes under line-tied conditions prohibit the existence of pure MHD waves. Line-tying conditions amount to treating the photosphere as providing complete reflection to any coronal wave. This condition is typically justified due to the large gradient in density and temperature from the photosphere to the corona. Although the propagation of waves from the lower atmosphere into the corona has received continuous interest, the leakage of coronal waves through the footpoints into the chromosphere has received less attention, though it has implications for wave damping and hence heating (Hollweg 1984; Davila 1991; Goedbloed & Halberstadt 1994; De Pontieu et al. 2001). Terradas et al. (2011) have shown, using a uniform 1d model, that line-tying conditions themselves induce coupling between fast and slow magnetoacoustic modes. A fast wave obliquely incident to the solar surface leads to the generation of slow waves with the same periodicity but with a wavelength scaled by the ratio of sound to Alfvén speeds. In particular, a standing fast mode leads to the generation of propagating slow waves from both loop footpoints with an amplitude that is proportional to the square root of the plasma-β. This coupling effect has been proposed as a possible explanation for intensity oscillations that accompanied a transverse loop oscillation with the same periodicity and was observed near the footpoint of a large loop (Verwichte et al. 2010). Large coronal loops that reach heights of over 1000 Mm are embedded in parts of the corona where the plasma-β is of the order of a few tenths of unity, and hence the Alfvén and acoustic speed values approach each other (Gary 2001; Régnier et al. 2008). The correct identification of the source of upward propagating waves into the corona, especially if that source lies in the corona itself, is important for assessing the amount of energy transport from the lower atmosphere.

In the study by Terradas et al. (2011) the corona was treated as uniform, with waves allowed to travel freely in any direction. We wish to expand upon this work by considering the effect of transverse coronal structuring. We evaluate the efficiency of linear foot point coupling for a coronal loop that supports a transverse loop oscillation.

The paper is structured as follows. In Sect. 2 an analytical model describing the linear coupling that takes into account transverse structuring is introduced. In Sect. 3 non-linear MHD simulations are performed for an isothermal 2D loop arcade model to evaluate the range of validity of the analytical findings. In Sect. 4, the numerical results are forward-modelled to produce synthetic diagnostics to reveal how the wave coupling may be observed. This is applied to the observation by Verwichte et al. (2010). In Sect. 5 we summarise our findings and provide an outlook.

2. Analytical modelling

We have chosen to model the effect of footpoint coupling of a transverse loop oscillation in a two-dimensional slab geometry. As such, the third dimension is ignored and the loop essentially represents a coronal loop arcade. This also eliminates the important mechanism of resonant absorption (Ionson 1978; Heyvaerts & Priest 1983; Goossens et al. 2011). This approach allows to study the effect of transverse structuring on footpoint coupling in a computationally feasible parametric MHD simulation study, and separate for now from the additional complexity arising from resonant absorption. The analytical modelling is performed for the same geometry to allow for a direct comparison. This analysis can be easily extended to a 3d cylindrical geometry with piece-wise continuous transverse equilibrium profiles. In this paper, the solar surface is modelled in a simplified manner as a hard line-tied boundary with zero velocities determined by the numerical domain.

Given the solar corona is at first approximation well described in the zero plasma-β limit, we chose in the analytical model to describe the transverse behaviour of the wave in this limit. The effect of a finite plasma-β is introduced as a perturbation to this so that the longitudinal dynamics follows as a perturbation driven by the transverse motion. Through comparisons with the non-linear MHD simulations, we can evaluate the efficacy of this approximation.

We introduce the coordinates x and y to be across and along the loop respectively, with the equilibrium loop axis at x = 0, and the loop footpoints situated at y = 0 and y = L, where L is the loop length. The equilibrium magnetic field B₀ is straight along the loop in the y direction and may vary in the x-direction to maintain pressure balance. The equilibrium density ρ₀ is uniform along the loop but may vary in the transverse x-direction. For this equilibrium total pressure is constant:

(1)

where p₀ and μ₀ are the equilibrium pressure and the magnetic permeability of free space respectively. For an isothermal loop, the equilibrium pressure will vary with x through the density. To maintain equilibrium, the magnetic field strength has to vary as well. For an isobaric loop, pressure and magnetic field strength are constant. In terms of the sound speed, c_S, and Alfvén speed, v_A, Eq. (1) is written as

(2)

where γ is the adiabatic index. The strength of the plasma pressure relative to the magnetic pressure is expressed by the plasma-β

(3)

where p₀ is the equilibrium plasma pressure and μ₀ the permeability of free space. In an isothermal loop the speed of sound is constant. The Alfvén speed varies as:

(4)

In the zero plasma-β limit the equilibrium magnetic field strength is constant and the Alfvén speed’s dependency on x is purely through the density. In an isobaric loop, where plasma pressure and magnetic field are constant, both sound and Alfvén speeds depend on x through the density for all values of the plasma-β.

For the equilibrium, the linearised MHD equations may be written in terms of the displacement vector ξ as:

(5)

(6)

(7)

(8)

where B, ρ and p are the perturbed magnetic field, density and plasma pressure, and with ⊥ and ∥ directions perpendicular and parallel to the loop respectively. These quantities may be eliminated in favour of the displacement.

(9)

(10)

The displacement obeys reflecting boundary conditions at the foot points, i.e. ξ(y = 0) = ξ(y = L) = 0. The displacement can thus be decomposed into sine Fourier series in y. We are primarily interested in the longitudinal dynamics in response of a transverse wave. In the zero plasma-β limit longitudinal displacements are absent. Therefore, to leading order ξ_⊥ may be approximated by the solution of the transverse wave Eq. (10) in the zero plasma-β limit:

(11)

In principle, this PDE has to be solved as an initial-value problem to model the excitation of the transverse loop oscillation. However, it is reasonable to consider ξ_⊥ to have the spatial structure of a specific trapped normal mode solution. There will still be some small transient waves, but their impact is secondary. Therefore, without much loss of generality, we take the solution of Eq. (11) to be

(12)

where ξ₀ is the displacement amplitude, k_p = pπ/L the longitudinal wave number, ω is the mode frequency and H(t) the Heaviside function. The transverse profile f(x) is the solution of

(13)

The profile f(x) tends to zero in the limit of |x| going to infinity and is symmetric around x = 0 (provided the equilibrium is also symmetric).

It is well-known that besides the linear footpoint coupling (Terradas et al. 2011), a transverse wave can also excite longitudinal dynamics through the non-linear ponderomotive force (e.g., Hollweg 1971; Verwichte et al. 1999; Terradas & Ofman 2004). Since we ignore to leading order the effect of the longitudinal dynamics on the transverse dynamics, we can add the ponderomotive acceleration term

(14)

to the right-hand side of Eq. (9):

(15)

The two right-hand terms of Eq. (15) represent transverse wave driving terms due to acoustic coupling and ponderomotive force, respectively. As the PDE is linear for ξ_∥ we construct the solution as the superposition of the response to each term separately, i.e.

(16)

2.1. Acoustic coupling

We consider the first of the three right-hand driving terms and solve Eq. (15) as an initial-value problem:

(17)

We decompose ξ_∥, a and cos(k_py) into a sine Fourier series

(18)

with

(19)

Equation (17) is for each Fourier component:

(20)

which can be straight-forwardly solved using the Laplace transform as

(21)

where is the convolution operator in time. The displacement ξ_∥, s takes the form of a convolution of the transverse wave time evolution of two counter-propagating sound waves that in total obey the boundary conditions:

(22)

where

(23)

We may take two approaches towards calculating the convolution in expression (22). Firstly, we solve the convolution for each Fourier mode separately. This results in the general solution

(24)

where

(25)

(26)

The two counter-propagating acoustic waves are again apparent, but their spatial structure, for example wavelength, is obscured by the series expansion. The third term represents a standing wave that is the longitudinal counterpart of the transverse oscillation. More physical insight can be gained from calculating the convolution using ϕ(s) in terms of the absolute value of a sine. Without loss of insight, we restrict the solution to p = 1 and times up to one acoustic loop crossing time, i.e. t ≤ L/c_S. We split the convolution integral into four regions in time-space separated by world lines t − y/c_S = 0 originating from y = 0 and t − (L − y)/c_S = 0 originating from y = L. We may then write Eq. (22) as

(27)

where G(s) = 2H(s)−1. The amplitude of the solution has the typical form for the polarisation between slow and fast modes. The first two terms of Eq. (27) confirm two counter-propagating acoustic sinusoidal waves travelling in from the foot points, with frequency ω and wavelength λ_S = 2πc_S/ω = λ_F(k₁c_S/ω) that is proportional to the ratio of the sound speed over transverse wave phase speed. This is consistent with the finding of Terradas et al. (2011). The third term is again the longitudinal counterpart of the transverse oscillation. The final two terms are acoustic transients.

We note that the dependency of the solution on the transverse coordinate x is passively through df/dx, independent of the transverse temperature structuring. At each location, sound waves propagate from the footpoints at the local sound speed. Since for a transverse oscillation df/dx is an anti-symmetric with respect to the loop axis, the sound waves will have the same transverse profile. However, there is an important difference between an isothermal and isobaric loop. In the former, the sound waves propagate in the y-direction at the same speed across all transverse locations in x. For a non-isothermal loop, the sound waves propagate at varying speeds and with time phase-mixing will occur (De Moortel et al. 2004; Voitenko et al. 2005). In that case, dissipative effects such as viscosity need to be included which introduces additional terms including −ν ∂³ξ_∥/∂x²∂t to the left-hand side of Eq. (15).

2.2. Ponderomotive acceleration

The third driving term in Eq. (15) leads to the problem

(28)

Following the same procedure to solve the initial-value problem we find

(29)

Again, solving the convolution leads to:

(30)

Note that in the zero plasma-β limit the solution has an unphysical, secular solution of the form (Verwichte et al. 1999):

(31)

The ponderomotive solution (30) is a standing mode with half the wavelength and represents the acoustic response to plasma being pushed towards the anti-nodes of the transverse velocity wave field. The transverse structure is symmetric with respect to x = 0. An order of magnitude ratio of Eqs. (24) and (30) is proportional to β^3/2ξ₀/L. It shows that the ponderomotive force needs to be taken into account when the wave amplitude is non-linear and the plasma-β is low.

3. Numerical modelling

For the numerical modelling, we use Lare2d, a Lagrangian remap code that solves the non-linear MHD equations, here under the assumption of a fully ionised, conductive plasma (Arber et al. 2001). We choose an average ion mass appropriate for the solar corona, where m_p is the proton mass. In the two-dimensional domain we model a loop as a slab with a magnetic field aligned along the loop in the y-direction. On the loop axis, we choose a magnetic field strength B₀ = 2.0 × 10⁻³T and mass density ρ₀ = 5.0 × 10⁻¹² kg m⁻³. Between simulation runs we vary the temperature and loop length, L, but they are generally of the order of millions of degrees K and hundreds of Mm, respectively.

We first perform simulations with uniform transverse density ρ(x, y) = ρ₀ and temperature profiles T(x, y) = T₀ so that we can compare our results with those established by Terradas et al. (2011). Transverse structuring is then added to the model. For this, we choose a symmetric Epstein density profile to model a continuous profile from internal to external loop media (Nakariakov & Roberts 1995):

(32)

where ρ_i and ρ_e are the mass densities on the loop axis (internal) and in the limit of |x|→∞ (external), respectively. We choose a loop half-width a = 0.05L and a density contrast ρ_i/ρ_e = 10. We vary the plasma-β of the simulation by varying the plasma temperature while keeping the magnetic pressure constant, this implies that the Alfvén speed remains unchanged and hence there is little variation in the oscillation frequency between simulation runs. The simulation is initialised with a perturbation in velocity of the form

(33)

with velocity amplitude, v₀, longitudinal wave number, k_p = π/L and symmetric transverse perturbation profile, f(x). This choice favours the excitation of the fundamental kink mode. In the linear limit, the velocity components v_x and v_y correspond to v_⊥ and v_∥, respectively. This allow us to directly compare the perturbation components between the numerical and analytical modelling. However, when the perturbation amplitude becomes non-linear, these two types of components will deviate from each other.

3.1. Wave coupling in a homogeneous medium

In order to compare with the linear modelling done by Terradas et al. (2011), we perform simulations with a uniform transverse density ρ(x, y) = ρ₀ across the whole simulation domain. The transverse velocity profile is of the form f(x) = sin(k_⊥x), where k_⊥ = 10/L is the perpendicular wave number and v₀ = 0.001v_A. The domain is set to fit a complete standing wave along the loop in y and a whole number of wavelengths across in x, i.e., 0 Mm ≤ y ≤ 240 Mm and −2π/k_⊥ ≤ x ≤ 2π/k_⊥. We introduce closed boundaries at both loop footpoints in y and periodic boundaries across the loop in x. The number of numerical grid points used is 280 × 352. We vary the sound speed c_S through the temperature. We perform simulations with the following sound to Alfvén speed ratios: c_S/v_A= 0.190, 0.379 and 0.602. From Eq. (27), we see that the amplitude of each propagating slow wave in the uniform media is of the form

(34)

which matches the expression found by Terradas et al. (2011).

To allow for a direct comparison we choose for the analytical model the same fast wave frequency equal as found in the simulations. Figure 1 shows the velocity profiles for a set of sound speeds we obtained both numerically and analytically. There is good agreement between the numerical and analytical profiles. The discrepancies in phase and amplitude may be explained due to the slight mismatch in the fast wave frequency ω. It can be seen that the line tying boundary conditions are sufficient to bring about linear wave coupling in which the fast mode drives propagating slow modes from both loop footpoints. After time τ_a, leg = L/2c_S the propagating slow modes meet and interfere forming a quasi-standing wave. In coronal loops multiple hundreds of Mm in length we expect these waves to experience significant damping before reaching the loop top (e.g., De Moortel et al. 2002). We evaluate the parallel and perpendicular velocities at their maxima, at x = 0 and x = π/2k_⊥, respectively. The maximum amplitudes of the slow modes agree with those predicted by Eq. (34), with a maximum error of 3.3% due to the finite spatial and temporal resolutions.

Fig. 1. Profiles of the parallel and perpendicular velocity as a function of y, for three times. Velocities are normalised with respect to the initial velocity amplitude. Time is normalised by the Alfvén travel time τA = 2L/vA. The solid line represents the perpendicular velocity, which is approximately the same for varying sound speeds. The transverse velocity is taken at x = πL/20 where maximum and the parallel velocity at x = 0 where also maximum. The dotted dashed, dashed and dotted lines are the parallel velocity profiles for cS/vA = 0.190, 0.379, 0.602, respectively. The velocity has been multiplied by a factor four for clarity. Left panels: results from the Lare2d numerical simulation. Right panels: analytical solution Eq. (24).

Through the wave coupling, the slow wave will have the same frequency as the fast wave, but a shorter wavelength. We make use of the approximate dispersion relations ω ≃ k_Fv_A and ω ≃ k_∥, Sc_S, where k_F and k_∥, S are the total wave number for the fast mode and the longitudinal wave number for the slow mode, respectively. We may then find that the waves’ wavelengths are related as . This is evident when comparing between the simulations with ratios c_S/v_A = 0.190 and 0.379. As the ratio doubles, so does the wavelength of the slow mode. For larger ratios, the approximate dispersion relations require additional finite-β contributions. These initial tests encourage us that our model of slow mode coupling is valid in the linear regime we can therefore explore this wave coupling mechanism further by the inclusion of transverse structuring and non-linearity.

3.2. Wave coupling in an isothermal loop with transverse structuring

To better understand the fundamentals of this mode coupling in our 2D loop simulation it is useful to study a completely isothermal loop, so that there is no change in the sound speed across the loop.

The equilibrium magnetic field strength increases away from the loop axis to maintain total pressure balance. The Alfvén speed varies across the loop as well.

For the transverse loop structuring via a density enhancement of the form given by Eq.(32) it has been shown by Nakariakov & Roberts (1995) that in the low-β limit, the transverse fast wave profile is described by the ODE,

(35)

where v_A(x) is defined by Eq. (4) in the zero plasma-β limit, using the equilibrium density profile (32). Equation (35) has for a trapped kink mode the solution,

(36)

where v_A∞ is the Alfvén speed in the limit x → ∞, v_ph = ω/k_∥ is the longitudinal phase speed which satisfies v_A0 < v_ph < v_A∞ and which satisfies the dispersion relation,

(37)

with v_A0 the Alfvén speed at x = 0. We solve the dispersion relation (37) numerically for the phase speed and subsequently deduce α. The simulation is initialised with a velocity perturbation (33) with the form for f(x) given by Eq. (36). This will excite the fundamental fast kink mode, even in a finite-β plasma. Though, we expect some additional small transients in that case.

The numerical domain in the x-direction spans −2L ≤ x ≤ 2L, which is large enough to cover the exponential decay of the wave solution across the loop and minimise issues with the boundaries. The simulations are run with a numerical grid size of 4800 × 1200. We first neglect non-linear effects by taking an initial perturbation amplitude v₀ = 0.001v_A. Figure 2 shows the excited transverse fast oscillation with a period of around 4min, which matches the period found from Eq. (37) to within 20 seconds. We also see that the fast oscillation drives slow mode propagation from the footpoints similarly as in the uniform case. The significant difference with the uniform case occurs in the slow wave profile in the x-direction. Equation (24) shows that for the Epstein equilibrium density profile this profile is equal to:

(38)

Fig. 2. Results of the Lare2d simulation for a model with cS/vA = 0.190, temperature 1 MK and plasma-β = 0.043. Velocities are normalised by the internal Alfén speed vA0. Time is normalised by the Alfvén travel time τA = 2L/vA. The two horizontal dashed lines indicate positions |x|=a. Top panels: transverse velocity vx/vA0 in the spatial numerical domain for four times. Bottom panels: longitudinal velocity vy/vA0 in the spatial numerical domain for the same four times.

which is anti-symmetric and evanescent with x. We therefore see two propagating slow waves on each edge of the loop and in anti-phase with each other. Also, the waves propagating from each foot point have opposite phases. This transverse wave profile has important consequences for their observability, as will be discussed later.

3.3. Analysis of the slow wave density perturbations

It is useful to perform a parametric analysis of the slow mode density perturbations as a function of plasma-β for several different initial velocity amplitudes. This allows us to both, determine the theoretical maximum possible density variations due to wave coupling in a loop, and assess the limits of validity of our linear-analytical model. For these simulations a loop length of L = 800 Mm is chosen, to draw out the spatial structure. We measure the maximum density perturbation at the loop footpoint where the longitudinal counterpart of the fast mode is absent and we can decompose the density perturbation into the components due to acoustic coupling and the ponderomotive force. These components can be distinguished by their different frequencies, i.e. ω for the slow wave and 2k_pc_S and 2ω for the ponderomotive response. We apply a spatial high-pass filter to the time series to isolate the acoustic coupling components.

As a benchmark, we repeat the same bandpass filtering for the analytically obtained density perturbation, which is found after substituting the displacement components Eqs. (27) and (30) into the continuity Eq. (6). We note that the period used in the analytical expression is chosen to be equal to those found in the corresponding simulation, to take into account the decrease of the oscillation period with the plasma-β.

Figure 3a shows the dependency of the maximum slow-mode density amplitude on plasma-β for different initial velocity amplitudes. The analytical model is in excellent agreement with the numerical data for velocity amplitudes in the linear regime where v₀ ≪ v_A0. As predicted from Eq. (27), the maximum slow mode density perturbation has a β^1/2 dependency. For larger values of the plasma-β it becomes more difficult to disentangle the contribution of acoustic coupling from that from the ponderomotive force, which becomes more pronounced at larger velocity amplitudes due to the dependency of the ponderomotive force. This may therefore explain the discrepancies seen in curves with v₀ > 0.1. For values of the plasma-β greater than 0.1, the density perturbation amplitude has reached a plateau and does not vary much. This implies that the measurement of the slow wave amplitude is not an accurate measure of the plasma-β in a loop.

Fig. 3. (a) Relative maximum slow mode density perturbation as a function of plasma-β for four different initial velocity amplitude v0. Lines represent the predicted analytical solution and crosses represent the numerical values. (b) Relative density perturbation amplitude due to the pondermotive force as a function of v0 for four different plasma-β values, taken at the point x = 0 and y = 0 where maximum.

3.4. Analysis of the ponderomotive density perturbations

We consider the contribution of the ponderomotive force to the density perturbation. The ponderomotive force pushes plasma towards anti-nodes in the transverse wave field. These variations are observed at time scales of L/c_S as a decrease in density at the footpoint and an increase at the loop top. We measure this contribution to the density perturbation at the loop footpoint x = 0, where both the acoustic wave coupling and fast mode contributions are zero. and where the ponderomotive effect has an extremum. We measure the absolute maximum of the density perturbation as a function of initial velocity amplitude v₀ for a number of different plasma-β. Figure 3a shows that for small perturbations there is excellent agreement between the analytical and numerical modelling, with both showing the predicted dependency in the quasi-linear regime. For v₀ ≳ 0.1v_A and more pronounced for smaller plasma-β, there is a clear discrepancy between the analytical and numerical results, because of the stronger non-linearity. In the zero plasma-β limit, the ponderomotive force can cause the density to become locally unphysically infinite (Verwichte et al. 1999). For v₀ ≳ 0.1v_A the density variations due to the ponderomotive force dominate over the variations due to the acoustic wave coupling. Given a plasma-β of 0.16 and v₀ = 0.35 we find a ponderomotive density variation of 18% which translates into a variation in intensity of 39%, which should be readily observable but because of its long wavelength may become obscured by any initial flows related to the event leading to the loop oscillations.

4. Synthetic diagnostics

4.1. Implementation

A possible observational signature is afforded by the intensity perturbation taken along the line-of-sight (LOS), which for the optically thin corona is typically given by

(39)

where s is the spatial coordinate along the LOS, G(λ, T) is the contribution function for a specific resonant emission line or instrumental bandpass, and n_e is the electron number density. When we consider isothermal equilibria and quasi-linear perturbations, we can neglect the contribution function when calculating lines of sight. Lines of sight are chosen across the straight slab simulations such that they are physically feasible if translated to a curved loop. Again, we choose a loop of length L = 800 Mm and examine the behaviour during the initial time interval before the propagating slow waves from both footpoints meet and start to interfere at τ_a, leg. The analytical modelling allows for the differentiation of the various contributions to the LOS intensity perturbation in the quasi-linear regime as

(40)

where I_a, I_pmf and I_f are Lagrangian intensity variations from acoustic wave coupling, ponderomotive effect and the fast mode density perturbation. Note that there is important contribution from advection (Cooper et al. 2003). These contributions are derived by substituting the displacements found in Eqs. (12), (27) and (30) into the continuity equation Eq. (6). Figure 4 shows the square of the density perturbation solely due to the slow waves. The anti-symmetric nature of the perturbation in the x-direction is again apparent. Therefore, the choice of a LOS ray angle, θ, is important for the visibility of the wave in intensity. Figure 5 shows for different plasma-β the relative slow mode intensity amplitude as a function of θ. We count θ from the y-direction (along the loop axis) to the x-direction (perpendicular to the loop). For θ = ±90° the LOS integrated intensity perturbation will be greatly reduced as in the linear regime the contributions from the compression and rarefaction cancel each other exactly. Also, for θ ≲ arctan(4a/L) = 11°, the LOS ray passes through at least half the loop. This will lead to a reduced intensity perturbation as the many compressions and rarefactions along the LOS tend to cancel each other. There exists an optimal LOS angle, θ_max, which passes through an extrema in emission on both sides of the loop. This angle depends on the wavelength of the generated slow modes and hence on the plasma-β. Using geometrical considerations, it is approximately equal to

(41)

Fig. 4. Snapshots of slow mode emission measure perturbation (ρa/ρf, max)2 obtained analytically for plasma-β = 0.04 and velocity amplitude v0 = 0.001vA. The half-loop travel time τa, leg = L/2cS = 2750 s. Note we have normalised the emission measure by the maximum fast mode intensity perturbation. The solid line passing through the fourth image is an example of an optimum LOS ray passing through two positive anti-nodes.

where we place the slow wave field extrema at approximately |x|= ± 0.75a. We then find for the appropriate values of frequency, optimal angles θ_max = 44° ,35° and 27° for plasma-β values 0.04, 0.08 and 0.16, respectively. This matches up well with the angle of the maximum in Fig. 5. he value of θ_max decreases with plasma-β and is explained given the increase in acoustic wavelength with plasma-β so that a ray needs to move further along the loop to reach the next slow mode anti-node on the other side of the loop axis. The amplitude of these maxima follows the expected β^1/2 dependency in the linear regime. This analysis shows that the observation of the slow waves strongly depends on the three-dimensional orientation of the loop and the line-of-sight of the observer.

Fig. 5. Relative slow mode intensity as a function of the LOS angle taken with respect to the loop axis. The slow mode intensity is obtained analytically for plasma-β values of 0.04, 0.08 and 0.16. The relative intensity is divided by the normalised initial perturbation amplitude v0/vA0.

We analyse the total intensity perturbation found in the numerical simulations for several suitable LOS rays as a function of time. We choose for each particular plasma-β two LOS angles: θ_max and θ = 80°. For each angle, we have two rays, one taken at the loop footpoint and another at approximately a quarter of the loop length. They are shown in Fig. 6.

Fig. 6. LOS rays taken through the loop from which we obtain intensity time series. A ray is taken at the footpoint for an angle of 80° given by the solid line and rays are taken which pass through two maxima for the specific plasma-β.

In Fig. 7 we compare the numerical intensity perturbations as a function of time with those predicted by the analytical model for the chosen line-of-sights shown in Fig. 6. We have divided each relative intensity perturbation by the initial velocity amplitude to make results comparable across perturbation scales. There is a good agreement between the numerical and analytical results. There is some disparity at the footpoint for small values of the plasma-β, which may be explained by disparities in the ponderomotive contribution which scale as β⁻¹ so that errors are amplified at lower plasma-β. Through the analytical model, we are able to characterise directly the contribution to the intensity from the coupling to running acoustic waves. For rays taken near the footpoint at 80° with respect to the loop axis, the intensity is dominated by the running acoustic waves. Its intensity amplitude increases with β and follows the expected β^1/2 dependency.

Fig. 7. Relative intensity variations as a function of normalised time, τS = L/cS, integrated along line-of-sights shown in Fig. 6. Each panel corresponds to a specific plasma-β and LOS angle with respect to the loop axis. Thick (thin) red lines correspond to the intensity derived analytically that includes (excludes) the acoustic coupling term. Line-style and colour match those in Fig. 6. The relative intensity is divided by the initial perturbation factor v0/vA = 0.001.

This is to be expected since for such an angle LOS confusion is minimal. The intensity variations are for typical transverse oscillations expected to have relative amplitudes of a few percentages, which is within the realms of observation.

We consider rays taken almost a quarter of the way along the loop and chosen with an angle such as to pass through two consecutive maxima in the slow mode wave field, i.e. at the angle of maximum intensity amplitude in Fig. 5. We observe initially no slow mode perturbation until the propagating slow waves pass through the LOS ray. This is reflected by the deviation of the analytical profile without these waves from the numerical result. At this position, the intensity is dominated by the density perturbation counterpart of the fast wave. The slow mode perturbation again increases with β.

4.2. Application to observation on September 15 2001

Terradas et al. (2011) presented a possible explanation for the intensity variations reported by Verwichte et al. (2010) in terms of linear wave coupling at line-tied loop footpoints. However, as their model did not include transverse structuring, the question of observability of this phenomena could not be fully addressed. Verwichte et al. (2010) found in observations from TRACE and EIT/SoHo a set of long period transverse loop oscillations accompanied by variations in intensity with the same period. For a particular loop oscillation also intensity variations at the loop footpoint were observed with the same period but in anti-phase with respect to the intensity variations in the loop leg. We simulate these observations using our model with equilibrium values chosen to closely resemble those of the loop studied by Verwichte et al. (2010). We take an isothermal loop of length L = 690 Mm, loop width a = 20 Mm, an electron number density on the loop axis of n_e = 2.5 × 10¹⁵ m⁻³, density contrast ρ_e/ρ_i = 0.5 and equilibrium magnetic field along the loop of 1 × 10⁻³ T. The sound speed is equal to c_S = 141 km s⁻¹. The Alfvén speed on the loop axis is v_A0 = 398 km s⁻¹. The loop axis plasma-β = 0.15. For this equilibrium, we solve dispersion relation (37) and retrieve the transverse profile, f(x), phase velocity v_ph = 580 km s⁻¹ and period P = 40 min. We add the caveat that in a slab geometry the phase speed of the transverse oscillations is close to the external Alfvén speed rather than the kink speed. Therefore, the chosen equilibrium values differ from those for a cylindrical loop model in order to match the observed periodicity.

The simulation is initialised with the fundamental transverse oscillation with a velocity amplitude v₀ = ξ₀ω = 0.270v_A, where the displacement amplitude ξ₀ = 44 Mm is taken from Verwichte et al. (2010). The density amplitude analysis presented here gives for a plasma-β > 0.1 density variations due to acoustic coupling with amplitudes of around 2−5%. This translates into intensity variations with amplitude of 4 − 10%. This matches well with the reported intensity amplitudes of the order of 10% by Verwichte et al. (2010).

We wish to re-evaluate the observation of intensity variations at the loop footpoint made by Verwichte et al. (2010) and shown in their Fig. 12, here reproduced in Fig. 8. We therefore map the simulated physical quantities from a straight slab onto a semi-circular grid to match observations more closer when taking LOS intensities. Figure 9 shows the curved loop along with the LOS rays taken along the loop up to 22 Mm to simulate intensity observation made at the loop footpoint by Verwichte et al. (2010). The oscillation starts at 11:19 UT and we artificially add damping with a damping time τ = 5190 s to more closely replicate the observational intensity variations.

Fig. 8. Intensity variations at the footpoint of a large loop which supports a transverse oscillation, taken from TRACE observation on September 15 2001. Top: relative intensity of the TRACE 2d loop region, averaged over the longitudinal coordinate from the foot point up to 22 Mm, as a function of time and transverse distance. The vertical dashed line indicates the reference time t0. Middle: intensity of the same region, averaged over transverse distance, as a function of time. The dashed line is a smooth trend that is taken to be the background intensity. Bottom: relative intensity fitted with a damped oscillator between times t0 and t1, is shown as a dashed curve. Reproduced from Verwichte et al. (2010).

Fig. 9. Straight loop simulation mapped to a semi-circular grid. Only one half of the numerical domain is shown. Lines of sight are taken at the foot point from s = 0 to 22 Mm height, indicated by the solid blue lines.

The intensity variations along LOS rays at the footpoint are shown in Fig. 10. The anti-symmetric wave pattern between the two sides of the loop is still visible, though confused by the LOS integration. This pattern corresponds to the 40 min slow/fast wave oscillation. Figure 8 also shows a twin-track wave pattern, though clear evidence of an anti-symmetric pattern is absent. We conclude that this pattern may be evidence of slow waves on two sides of the loop, or alternatively of two separate loops.

Fig. 10. (a): Relative intensity variations averaged along a set of LOS rays as shown in Fig. 9, as a function of time, obtained numerically simulated by lare2d for an oscillation beginning at 11:19 UT. (b): Relative intensity averaged over the transverse coordinate as a function of time, taken from the same region, as a function of time. The dashed line is a smooth slow trend fitted using a Savitzky-Golay filter (Savitzky & Golay 1964). (c): Relative intensity as a function of time found after subtracting the slow trend. The dashed line is a damped sinusoidal fit. (d) Relative intensity as a function of time found after subtracting the slow trend from a time series where the relative intensity in panel a has been averaged only over the transverse interval for x between −228 and −185 Mm in the numerical domain. The dashed line is a damped sinusoidal fit.

There is also an overall pattern present with a period of about 79 min, common across transverse coordinate x. In order to obtain a clearer time profile we average the intensity over the transverse coordinate x (see the middle panel of Fig. 10). Here, the 79 min oscillation becomes apparent. This corresponds to the oscillations driven by the ponderomotive force, which have predicted periods of L/c_S = 81.6 min and P/2 = 80 min. These perturbations dominate the overall signal with an intensity amplitude of 16.2%. With this knowledge we re-examine Fig. 8. Their second panel shows the relative intensity with time. The slow wave periodicity is superimposed on a trend that consists of a large dip that corresponds to the density depletion associated with the CME that excited the oscillation, and a small-amplitude oscillation with a period of approximately 80 min that is most apparent in the recovery phase of the dip. The analysis here leads us to interpret this oscillation as the response to the ponderomotive force. A similar but less apparent 80 min periodicity is present in the intensity perturbation in the loop leg.

Subtracting a smooth fit from this signal leaves the remaining slow wave oscillations (bottom panel of Fig. 10). We fitted these with a damped sinusoidal curve and retrieved the period P = 40 min as expected. The fit ignores the effect of slow mode superposition which occurs after τ_a, leg = 40.8 min. We observe a maximum slow mode perturbation of around 3%. This is smaller than the 15% reported by Verwichte et al. (2010) although it is important to note due to the anti-symmetric slow mode profile much of the signal is lost in averaging across the loop and we study a much thicker loop for computational ease. If we focus only on one transverse half of the loop to limit the destructive superposition, we find a maximum amplitude of 4.5%. This is also still smaller than the observed amplitude. But we note that the reported value of this amplitude is based on a model extrapolation as only one loop leg was visible in the TRACE field of view. Thus, there is some uncertainty on the amplitude value. Also, the lack of detailed modelling of the lower atmosphere in our numerical model may be the cause for this difference.

5. Discussion

Terradas et al. (2011) presented a possible explanation for observations of intensity variations that accompany Alfvénic loop oscillations, in terms of linear coupling between Alfvénic and slow magnetoacoustic waves at the loop footpoints. They had modelled this in a uniform plasma using linear wave theory. In this work, we have extended the modelling of this coupling by the use of fully non-linear MHD simulations and the inclusion of the effect of transverse structuring which is a common feature for Alfvénic loop oscillations. We have chosen to model the loop as a two-dimensional slab arcade for numerical convenience and also to clearly separate this physical mechanism from the resonant absorption that should take place in the loop as well. We have performed both analytical and numerical modelling and find overall good agreement between the two.

We find that the generation of propagating slow waves through wave coupling at the line-tied boundaries persists when transverse structuring is included. The important difference with the uniform case is that for a trapped mode the transverse wave profile is fixed by the dispersion relation and evanescent on both sides of the loop. This leads to an anti-symmetric pattern for the propagating slow wave field. The slow waves have the same period as the transverse (fast, Alfvénic) oscillation but a shorter wavelength given by λ_S = c_SP. The strength of the wave coupling is proportional to the square root of the plasma-β and to the initial perturbation amplitude, in agreement with the finding of Terradas et al. (2011). The numerical study allowed us to show that this relation persists for a wide range of sub-Alfvénic velocity amplitudes. This is also reflected by the good agreement between analytical and numerical modelling.

We have modelled the density and subsequent intensity perturbations arising from the propagating slow waves. The transverse loop structuring leads to a trapped Alfvénic oscillation and the propagating slow waves exhibit an anti-symmetric wave pattern in the direction of transverse wave polarisation. We have shown that this has an important consequence for the ability to observe these waves. Due to destructive superposition, the slow waves will not be easily visible along line-of-sights in the directions parallel with the wave polarisation or along the loop axis. We have shown that there is an optimal angle of the line-of-sight when we observe through a maximum in the slow wave field at each side of the loop. Therefore, when looking for observational evidence of these slow waves, the three-dimensional loop geometry with respect to the observer needs to be carefully considered. The analysis also shows that in order for the intensity amplitudes to be sufficiently large for observations, transverse oscillations with velocity amplitudes of at least 10% the local Alfvén speed and a plasma-β > 0.05 are required. The studies by Verwichte et al. (2010) and Terradas et al. (2011) have established that this wave coupling mechanism can be observed. Future observational studies that focus on the footpoints of a large transversely oscillating loop in chromospheric, transition region and coronal temperatures at high spatial resolution, using observations from for example SDO, IRIS and DKIST, should be able to establish further detailed evidence.

We have shown that in a discussion of wave coupling from a transverse to slow wave applied to finite-amplitude oscillations in coronal loops, the impact of the ponderomotive force needs to be considered. Both analytical and numerical results show that the ponderomotive force leads to a loss of density at the loop footpoint as matter is pushed towards the loop top. In terms of intensity this effect acts on a longer acoustic time scale with the square of the wave amplitude. For a loop oscillating with a velocity amplitude of more than about 10% of the Alfvén speed the ponderomotive response should be readily observable. However, this signal may be masked by flows associated with a typical eruption event which generated the loop oscillation. Evidence of the ponderomotively driven acoustic oscillations may provide an additional useful coronal diagnostic tool alongside transverse oscillations. Using the modelling presented here, we are able to succesfully identify the pondermotive response as an oscillation in intensity of about 80 min in the observation reported by Verwichte et al. (2010). This period is consistent within 3% with the theoretically expected value of L/c_S with the loop length estimated from motional stereoscopy and sound speed corresponding to the peak temperature of the instrument bandpass. Therefore, the ponderomotive response provided an independent determination of the accuracy of the performed wave analysis.

The presented model can be improved in several ways, which are subject of a currently undertaken follow-up study. First, the choice of an isothermal loop is arbitrary. In reality, we may expect the loop temperature to differ from the external corona. The slow wave coupling will still occur at the footpoints. Because the sound speed now varies across the loop, the slow waves, in their propagation up the loop legs, will have their wave phase inside and outside the loop increasingly deviate from each other. This will lead to phase-mixing and enhanced dissipation (Voitenko et al. 2005). Second, in our modelling, the mechanism of wave coupling at the loop footpoints relies on line-tied boundary conditions imposed at the boundary of the domain. Instead, a more realistic lower atmosphere which includes a transition region and chromosphere can be included. Terradas et al. (2011) has already demonstrated using linear modelling in a transversely uniform medium that wave coupling will occur in such a scenario. We therefore expect wave coupling to occur at the transition region as well as the photospheric boundary. However transverse structuring may lead to transversal wave leakage in the chromosphere and we wish to assess its impact on the slow wave flux into the corona.

The analytical modelling presented here can easily be extended to a three-dimensional loop model where a coronal loop is modelled in cylindrical geometry, provided the transverse Alfvén speed profile is piece-wise continuous. Using the same perturbation ordering as in this study, we may express the parallel dynamics with the ODE

(42)

with the dilation in the plane perpendicular to the magnetic field

(43)

which for a standing kink (Alfvénic) wave that is linearly polarised in the x-direction has the form

(44)

where

(45)

with J₁(s) the Bessel function of the first kind of order one and K₁(s) the modified Bessel function of the second kind of order one. κ_i, e are the internal and external radial wave numbers, which are determined by solving the relevant dispersion relation for trapped modes. We may apply the same analysis for the acoustic wave coupling by substituting the term df/dx in Eq. (24) by g(r) cos θ/a. In this geometry, the propagating slow waves will also have an anti-symmetric nature, with the waves localised towards the edge of the loop with opposite phases at opposite locations in x, in the direction of wave polarisation. In cylindrical geometry the frequency of the transverse mode is approximately proportional to the kink speed, whilst in a slab geometry it is approximately proportional to the higher external Alfvén speed. Because the acoustic coupling mechanism is inversely proportional to the factor , we expect the coupling to be enhanced in cylindrical geometry. Note that the resonant absorption mechanism, which is present when the radial equilibrium profile varies continuously, has an azimuthal structure is concentrated at the edges in y. Ultimately, the effect of slow wave coupling should be modelled in the presence of resonant absorption.

Acknowledgments

S.J.W. acknowledges the support of the UK STFC Ph.D. studentship.

References

All Figures

Fig. 1. Profiles of the parallel and perpendicular velocity as a function of y, for three times. Velocities are normalised with respect to the initial velocity amplitude. Time is normalised by the Alfvén travel time τA = 2L/vA. The solid line represents the perpendicular velocity, which is approximately the same for varying sound speeds. The transverse velocity is taken at x = πL/20 where maximum and the parallel velocity at x = 0 where also maximum. The dotted dashed, dashed and dotted lines are the parallel velocity profiles for cS/vA = 0.190, 0.379, 0.602, respectively. The velocity has been multiplied by a factor four for clarity. Left panels: results from the Lare2d numerical simulation. Right panels: analytical solution Eq. (24).

In the text

Fig. 2. Results of the Lare2d simulation for a model with cS/vA = 0.190, temperature 1 MK and plasma-β = 0.043. Velocities are normalised by the internal Alfén speed vA0. Time is normalised by the Alfvén travel time τA = 2L/vA. The two horizontal dashed lines indicate positions |x|=a. Top panels: transverse velocity vx/vA0 in the spatial numerical domain for four times. Bottom panels: longitudinal velocity vy/vA0 in the spatial numerical domain for the same four times.

In the text

Fig. 3. (a) Relative maximum slow mode density perturbation as a function of plasma-β for four different initial velocity amplitude v0. Lines represent the predicted analytical solution and crosses represent the numerical values. (b) Relative density perturbation amplitude due to the pondermotive force as a function of v0 for four different plasma-β values, taken at the point x = 0 and y = 0 where maximum.

In the text

Fig. 4. Snapshots of slow mode emission measure perturbation (ρa/ρf, max)2 obtained analytically for plasma-β = 0.04 and velocity amplitude v0 = 0.001vA. The half-loop travel time τa, leg = L/2cS = 2750 s. Note we have normalised the emission measure by the maximum fast mode intensity perturbation. The solid line passing through the fourth image is an example of an optimum LOS ray passing through two positive anti-nodes.

In the text

	Fig. 5. Relative slow mode intensity as a function of the LOS angle taken with respect to the loop axis. The slow mode intensity is obtained analytically for plasma-β values of 0.04, 0.08 and 0.16. The relative intensity is divided by the normalised initial perturbation amplitude v0/vA0.
In the text

	Fig. 6. LOS rays taken through the loop from which we obtain intensity time series. A ray is taken at the footpoint for an angle of 80° given by the solid line and rays are taken which pass through two maxima for the specific plasma-β.
In the text

Fig. 7. Relative intensity variations as a function of normalised time, τS = L/cS, integrated along line-of-sights shown in Fig. 6. Each panel corresponds to a specific plasma-β and LOS angle with respect to the loop axis. Thick (thin) red lines correspond to the intensity derived analytically that includes (excludes) the acoustic coupling term. Line-style and colour match those in Fig. 6. The relative intensity is divided by the initial perturbation factor v0/vA = 0.001.

In the text

Fig. 8. Intensity variations at the footpoint of a large loop which supports a transverse oscillation, taken from TRACE observation on September 15 2001. Top: relative intensity of the TRACE 2d loop region, averaged over the longitudinal coordinate from the foot point up to 22 Mm, as a function of time and transverse distance. The vertical dashed line indicates the reference time t0. Middle: intensity of the same region, averaged over transverse distance, as a function of time. The dashed line is a smooth trend that is taken to be the background intensity. Bottom: relative intensity fitted with a damped oscillator between times t0 and t1, is shown as a dashed curve. Reproduced from Verwichte et al. (2010).

In the text

	Fig. 9. Straight loop simulation mapped to a semi-circular grid. Only one half of the numerical domain is shown. Lines of sight are taken at the foot point from s = 0 to 22 Mm height, indicated by the solid blue lines.
In the text

Fig. 10. (a): Relative intensity variations averaged along a set of LOS rays as shown in Fig. 9, as a function of time, obtained numerically simulated by lare2d for an oscillation beginning at 11:19 UT. (b): Relative intensity averaged over the transverse coordinate as a function of time, taken from the same region, as a function of time. The dashed line is a smooth slow trend fitted using a Savitzky-Golay filter (Savitzky & Golay 1964). (c): Relative intensity as a function of time found after subtracting the slow trend. The dashed line is a damped sinusoidal fit. (d) Relative intensity as a function of time found after subtracting the slow trend from a time series where the relative intensity in panel a has been averaged only over the transverse interval for x between −228 and −185 Mm in the numerical domain. The dashed line is a damped sinusoidal fit.

In the text