Propagation of Leaky MHD Waves at Discontinuities with Tilted Magnetic Field

We investigate the characteristics of magneto-acoustic surface waves propagating at a single density interface, in the presence of an inclined magnetic field. For linear wave propagation, the dispersion relation is obtained and analytical solutions are derived for small inclination angle. The inclination of the field renders the frequency of the waves complex, where the imaginary part describes wave attenuation, due to lateral energy leakage.


Introduction
The problem of wave propagation in solar and space plasmas is one of the most important aspects of plasma dynamics. Waves can carry energy across different layers of the solar atmosphere and can dissipate their energy so they contribute to the process of plasma heating in the upper atmosphere (see e.g. Einaudi, Chiuderi, and Califano, 1993;Erdélyi and Ballai, 2007;Arregui, 2015, etc.). During their progression, waves also carry the imprint of the environment in which they propagate. Therefore, their seismological study can reveal physical parameters that cannot be measured directly or indirectly, such as magnetic field strengths in the tenuous corona, magnitude of various transport coefficients, intrinsic self-organisation of the plasma, etc. (see e.g. Antia, 1986;Gough et al., 1996;Gizon and Birch, 2005;Andries et al., 2009;De Moortel and Nakariakov, 2012;Mathioudakis, Jess, and Erdélyi, 2013;Arregui, Oliver, and Ballester, 2018, etc.).
Often waves are propagating along magnetic field lines, which therefore act as tracers, while the structuring of the magnetic fields guide waves along them. While plasma structuring is often concurrent with the structuring of magnetic fields, there are several examples where waves propagate along discontinuities present in plasmas, where the magnetic field B E. Vickers Figure 1 The plasma is structured into two infinite regions of constant density and gas pressure, with a sharp interface at z = 0. A constant magnetic field crosses the interface and is inclined at an angle θ to the interface.
to the derivation of the general governing equation for the plasmas either side of the interface, which is given for an arbitrary inclination of the magnetic field, by Equation 8. In order to make analytical progress, in Section 4 we introduce the small inclination angle approximation of the magnetic field, with respect to the interface and solutions to the governing equations for either side of the interface are found, using a perturbation technique. The dispersion relation corresponding to this approximation is obtained in Section 5 and is given in dimensional variables by Equation 29 and in dimensionless form by Equation 30. In Section 6, solutions to the dispersion relation are determined, both analytically and numerically and finally, our results are summarised and discussed in Section 7.

Initial Equilibrium
We aim to investigate the propagation of linear and compressible MHD waves along a contact discontinuity in density and temperature, aligned with z = 0, with density and temperature constant either side of this discontinuity. A constant magnetic field permeates the plasma at an angle, θ , to the discontinuity in the (x, z) plane and has the form B 0 = B 0 (cos θ, 0, sin θ). In our working model, with a tilted magnetic field, the restoring force will be the tangential component of the Lorentz force, with magnetic tension acting on any displacement transversal to the field and magnetic pressure acting on any displacement that changes the magnetic field strength. The effect of gravity is neglected, and we restrict ourselves to study the dynamics of two-dimensional perturbations in a 3D Cartesian geometry. A schematic representation of the equilibrium configuration is shown in Figure 1.
The dynamics of waves can be described within the framework of ideal magnetohydrodynamics (MHD). We assume that the perturbations of physical quantities are small compared to their equilibrium values, therefore we can use the linearised version of the MHD system of equations. All physical quantities (except velocity) will be written as a sum of their constant equilibrium values and their perturbations, so that the density, pressure, magnetic field and velocity are given by ρ 0 + ρ, p 0 + p, B 0 + b and v, respectively.
The system of ideal linearised MHD equations for a homogeneous equilibrium, therefore, can be written as where c 0 = (γp 0 /ρ 0 ) 1/2 is the adiabatic sound speed, γ is the adiabatic index, and μ is the magnetic permeability of free space. In addition, the connection between the magnetic field and the plasma fluid is ensured by the linearised, ideal induction equation, Finally, the system of equations has to be supplemented with the solenoidal condition, ∇ · b = 0, as well as with the ideal gas law. These equations will be used in the next section to derive the dispersion relation for waves propagating along the interface.

Governing Equation
In general the propagation of waves can be studied with the help of their dispersion relation, i.e. the relation between the frequency of waves and the wavenumber, in terms of characteristic speeds and quantities specific to the medium in which they propagate. First, general solutions are found for the plasma regions either side of the interface, the dispersion relation may be obtained after matching the solutions in the two regions, at the interface z = 0. The mathematical procedure we employ in the present study would classify our task as an eigenvalue problem. However, when comparing to an initial value problem investigation, there are differences between a standard eigenvalue problem and the problem concerning leaky waves. As pointed out by Ruderman and Roberts (2006), while standard eigenvalue solutions correspond to the asymptotic behaviour of the time-dependent solutions, the leaky mode solutions are instead intermediate asymptotics, where the time-scale is greater than the period of the wave, but less than the attenuation time.
Let us begin by writing the momentum equation in terms of the inclination angle of the equilibrium magnetic field, θ , as where B x = B 0 cos θ , B z = B 0 sin θ . After differentiating the x and z components of this equation with respect to time and using the energy and induction equations, we obtain two equations for the horizontal and vertical velocity perturbations, where the Alfvén speed is given by v A = B 0 /(μρ 0 ) 1/2 . All of the perturbed quantities will oscillate with frequency, ω, and real wavenumber, k. Therefore, we take any perturbations to be of the form f =f (z) exp[i(kx − ωt)], where f is the amplitude of perturbations that depend on z. Next, we Fourier analyse the above system of equations to arrive at the system of coupled differential equations, Since all coefficients in Equations 6 and 7 are constants, we can eliminate the component of velocity parallel to the interface,v x to obtain a single fourth-order differential equation, for the component of velocity perpendicular to the interface,v z , If the inclination of the magnetic field is omitted (i.e. θ = 0), we recover the governing equation for compressional waves derived by Roberts (1981). Since we are looking for wavelike behaviour, the solution of the above equation will be of the formv z ∼ exp [ z]. Given that the governing equation (Equation 8) is a fourth-order differential equation, the general solution will be of the form where C i are constants and i are the roots of the characteristic equation The values of will be used in Section 4 to determine the dispersion relation for the waves propagating along the interface. We note that a similar equation is derived for the plasma regimes both above and below the interface, with appropriate characteristic speeds. We will show by contradiction that in order to obtain propagating solutions, must be complex. Let us assume that is purely real, then the left-hand side of Equation 9 is a complex analytical function and may therefore be split into real and imaginary parts as where u and v are the real functions, In order for Equation 10 to be satisfied, we require that both functions u(ω, k) and v(ω, k), are equal to zero simultaneously. Setting v(ω, k) = 0 gives the solutions = 0, ±k. The solutions = 0 is not a solution to u(ω, k) = 0. Substituting = ±k into Equation 8 gives ω = 0, so this is not a propagating solution. This proves that, in order for the wave to propagate, must be a complex quantity. The imaginary component of represents an oscillatory component to the variation ofv z with respect to the transverse coordinate. This will introduce energy flow into the system, either towards or away from the interface, for certain solutions. Since no energy source is specified, this only makes physical sense if the energy flow is away from the interface. These leaky wave solutions correspond to the case where the group speed is positive above the interface and negative below and the effect of lateral energy leakage is an attenuation of the waves.

Solving the Governing Equation -a Perturbation Technique
In order to make analytical progress, we assume that the angle between the magnetic field lines and the interface is small and so the inclination induces only a small change to the waves' properties in each homogeneous semi-infinite volume, compared to the case with parallel magnetic field. This allows us to make approximations in θ , letting cos θ ≈ 1 and sin θ ≈ θ . The governing equation of the quantity , for small values of inclination angle, is given by In order to physically account for the transition of perturbed quantities from one side of the interface to the other one, we consider a thin boundary layer (embracing the interface), in which the transition takes place, with width less than 2θ .
Despite the small angle, even the first term of Equation 11 is comparable to the other terms, since it is multiplied by the highest derivative ofv z , which can be large. Let us now apply the method of dominant balance to find the roots of the fourth-order polynomial. First, since some terms are smaller than others, they could not possibly be part of the dominant balance and may be ignored. We are going to seek solutions in the form of an asymptotic power series In leading order (i.e. terms proportional to O(1)), we obtain which means that where Interestingly, this quantity coincides with the effective wavenumber determined for magnetoacoustic modes obtained by Roberts (1981) in the case of tangential discontinuity, if we assume that m is real when ω is real. However, in the case of leaky modes, ω is complex and so too is m, introducing an oscillatory behaviour in the z-direction. A similar expression can be derived for the two plasma regions. For simplicity, we will introduce the subscripts + and −, to refer to equivalent parameters in the plasmas above and below the interface, respectively. Hence, in the z < 0 region we will use m − , whereas in the upper region (z > 0), we will use m + . Here, the signs of m − and m + should be chosen in such a way that the real parts of m − and m + are positive. We choose the signs of 0 above and below the interface, such that these match the solutions for the case of the parallel magnetic field, i.e. physical solutions are given by Equation 11 is a fourth-order polynomial and the remaining two roots must still be found. This can be done by rescaling the problem. For some unknown exponent Q (to be determined), let us set in Equation 11, where y is bounded and bounded away from zero as θ → 0. That is why Equation 11 becomes (after dividing by c 2 We will now find the correct value of Q by using the principle of dominant balance, so that the rescaled equation is consistent as θ → 0, if at least two terms correspond to the same power of θ (this is called balance). In addition, the balance is dominant in the sense that every term not involved in the balance corresponds to a higher power of θ , and therefore must be smaller than the balancing terms.
It can be shown that, when balancing the first three terms, the balance will occur for Q = −1 and the balancing terms are proportional to O(θ −2 ) while the other terms are ∼ O(1). With this value of Q we obtain Multiplying the above equation by μ = θ 2 we have Now, we write y also in the form of an asymptotic series in terms of μ as y = y 0 + y 1 μ + y 2 μ 2 + · · · In the leading order, we obtain which has four roots for y 0 , i.e.
The y 0 = 0 solutions contradict our assumptions that y is bounded away from zero as θ → 0 and must be disregarded, so the remaining two solutions will be Returning now to the original variables, the four roots of the polynomial in are (in the leading order) Higher-order terms in θ are neglected, as they are small compared to the leading-order terms. Let us briefly discuss the form and sign of the last two roots. The key ingredient in both expressions is y 0(C,D) = − i c T (kc T ∓ ω). Since we are dealing with an inclined interface, as mentioned earlier, we expect that modes will be leaky and the frequency of waves can be written as ω = ω r + iω i . Introducing this expression into the form of y 0(C,D) we obtain Since attenuation of waves due to leakage corresponds to ω i < 0, Equation 21 shows that (y 0C ) > 0, while (y 0D ) < 0. To ensure that the group speed is positive above the interface and negative below it, the physically acceptable solutions for y 0C and y 0D must have the corresponding sign. Hence, for the z < 0 region we are going to use y 0D , while in the region above the interface we will need to use a similar root as y 0C but written for the corresponding characteristic speeds.
Once again, we use the subscripts + and − to refer to the upper and lower plasma regions and so we write y 0D , below the interface, as y − and y 0C above the interface as y + . Therefore the expression ofv z in the lower (z < 0) region iŝ The corresponding expression for v z , in the upper (z > 0) region iŝ v z+ = F + e −m+z + G + e θ −1 y+z , where the critical speeds for the relevant plasma are used to determine m ± and y ± in that region.
It is clear that, as θ → 0, the expressions forv z± do not approach the variables we would obtain for the tangential discontinuity (for their expressions see, e.g. Roberts, 1981), due to the terms proportional to G ± . As previously mentioned, the amplitude of the transversal component of velocity varies over z in an oscillatory fashion, due to the complex nature of the exponents. In fact, the amplitude of the z-component of the velocity increases as the inclination angle decreases. Equally, the laterally oscillating leaking wave will have its wavelength increasing proportionally to the inclination angle.
In a recent paper Ruderman et al. (2018) it was shown that, in the case of a tilted magnetic field, the perturbed quantities do not show a simple transition from their values at the contact discontinuity to the tangential discontinuity, as θ → 0, however, their averaged values tend towards those corresponding to the tangential discontinuity. This averaging technique may also be applied to the solutions found here for v z , confirming the above statement.

Dispersion Relation of Waves Along Discontinuities
Since four constants are involved in the two expressions ofv z , four boundary conditions will be needed to find the values of those constants and determine the dispersion relation. As the magnetic field intersects the interface, this is a contact discontinuity and so we require continuity of both components of velocity,v x andv z , the kinetic pressure,p, and both components of the magnetic field,b x andb z , respectively. In this equilibrium, continuity of b z is implied by the continuity of v, so we have the correct number of boundary conditions to find all unknowns.
First, we must findv x . Using the components of the momentum equation, we can find thatv x may also be written in terms of the same exponentials asv z , i.e.
where the expressions for f ± and g ± are given by The condition that the vertical component of velocity,v z , is continuous across the boundary gives Continuity of the parallel component of velocity,v x , results in The condition that the parallel component of the magnetic field, b x , is continuous, gives to first order of θ . Finally, pressure balance across the interface gives Due to the particular choice of discontinuity, ρ − /ρ + = c 2 0+ /c 2 0− , so multipliers cancel in this equation.
These conditions may be written together as a matrix equation, We seek eigen-mode solutions, by solving the dispersion relation, which is given by the equation (29)

Dispersion Relation in Dimensionless Form
The two key parameters that can affect the characteristics of waves are the ratio between the densities of the two plasma regions (d = ρ − /ρ + ) and the plasma-β, defined as β = 2c 2 0 /γ v 2 A . That is why we determine the dispersion relation in terms of these two parameters and we write the equations in their dimensionless form. We introduce the phase speed, relative to the Alfvén speed of the lower plasma, asc ph = c ph /v A− = ω/kv A− . Due to the particular form of the magnetic field, we see that d = c 2 0+ /c 2 0− = v 2 A+ /v 2 A− and hence β − = β + = β. The effective wavenumbers in dimensionless form are The ratios between tangential and transversal components of velocity are given in dimensionless form by In this new, dimensionless, form, the dispersion relation is now given by the equation

Numerical Solutions and Discussion of Results
In what follows, we solve the dispersion relation (Equation 30) numerically for varying density ratio and plasma-β. For simplicity, we choose to plot only forward propagating waves.
We should note here that in the present configuration, the symmetry between forward and backward propagating waves is broken and this is due to the magnetic force that is perpendicular to the equilibrium magnetic field, affecting the forward and backward propagating waves in a different way. However, the differences between the dimensionless phase speed of forward and backward propagating waves are small due to the chosen approximation of small field inclination. It is very likely that significant differences may be obtained for arbitrary angles, but this would require a full numerical investigation. First, we plot the dispersion curves of forward propagating waves for a fixed density ratio. In Figure 2, we set d = 0.1 (upper panel) and d = 10 (lower panel) and allow the plasma-β to vary over two orders of magnitude covering the spectrum of possible values in the solar atmospheric plasma. When the plasma above the interface is heavier (here d = 0.1) the dispersion relation allows the propagation of three modes. Mode (A) is present for larger values of plasma-β (β > 0.6), with phase speed close to the sound speed of the upper plasma, c e . The imaginary part of the solution increases in magnitude for higher values of plasma-β, so that, for β = 10, the period is approximately half of the attenuation time. Since the phase speed is higher than the Alfvén speed in the upper region, this mode must be a highly attenuated fast mode. Its propagation speed increases with plasma-β.
Mode (B) is present for the entire range of plasma-β and the phase speed is, for most of the range, close to the cusp speed of the plasma in the upper region, c T e , increasing to a speed between the two cusp speeds for high plasma-β. For β → 0, we find that mode (B) also tends to zero, meaning that we are dealing with a slow wave. The imaginary part of the solution stays fairly steady, between 0 and −0.1, so these waves show a rather weak attenuation.
Mode (C) is only present for β > 0.5, since for lower values of plasma-β solutions do not satisfy the conditions set for the imaginary part of the frequency (they do not "leakout"). For β > 1, the phase speed of this mode decreases with increasing plasma-β, while the imaginary part increases in magnitude greatly, so that, for β > 1, the period of the wave is greater than the attenuation time, meaning that the expected life time of these modes is rather short. This mode should be labelled a slow mode, as it is slower than either cusp speeds, however, it shows a rather peculiar behaviour. Its phase speed does not increase with plasma-β, a feature characteristic for slow waves. Due to their high attenuation, modes (A) and (C) are unlikely to be observable, due to having an attenuation time much shorter than the wave period, especially in the high-β regime.
In Figure 2b, we show solutions to the dispersion relation for d = 10 and we can see two modes of propagation. Mode (A) has a phase speed between the sound and cusp speeds in the upper plasma region, though is only present for β < 3, since for higher plasma-β the mode is no longer a leaky mode. For low values of plasma-β, mode (A) tends to zero, so this is clearly a slow mode. In the region corresponding to β < 1, the attenuation of this mode is very small, however, the attenuation rate increases with the value of plasma-β. At a value of β = 2 the attenuation of the mode decreases again. Figure 2a, has phase speed close to the lowest value of the cusp speeds for lower plasma-β and its phase speed decreases to zero around β = 1. Again, based on its propagation speed we label this mode as a typical slow mode, however, it has the same peculiar behaviour with respect to plasma-β. This mode, again, is highly attenuated, especially for β > 1. Figure 3 shows the variation of the solutions of the dispersion relation with density ratio, for low and high values of plasma-β. In Figure 3a, we set β = 0.1 (a typical solar upper solar atmospheric condition) and the numerical investigation of the dispersion relation reveals the existence of three modes, all with relatively small imaginary part, meaning these modes are weakly attenuated. Mode (A) has phase speed below the cusp and sound speeds of the plasma in the lower region and a small imaginary part, which decreases in magnitude towards d = 1, at which point the imaginary part becomes zero. Mode (B) is a physical solution for d > 0.8 and has phase speed very close to the cusp speed of the lower plasma region. Although the dimensionless phase speed of these two modes are very similar in the region where d = 1, they have a rather distinctive imaginary part. Mode (C) exists for d > 1 and has phase speed close to both the sound and the cusp speeds of the upper plasma and this mode shows the largest attenuation among all possible modes. The d = 1 value corresponds to the situation when the difference between the two regions disappear and there is no interface. This situation was earlier studied by Cally and Schunker (2006) and the dispersion relation for MHD waves in this context, for small inclination angle between the wave vector and the magnetic field, may be easily solved to give c ph ≈ c 0 or v A . This agrees with the value of mode (C) at this point and explains the decrease in attenuation of modes (A) and (C), towards d = 1.

Mode (B), similar to the mode (C) in
In Figure 3b, we plot the variation of the dimensionless phase speed of waves for a given value of plasma-β (here taken to be 10) and we let the density ratio of the two regions vary. This regime of parameters is more relevant to photospheric conditions. We can see, the interface enables the propagation of two surface modes. Comparing these solutions to the ones obtained for low-coronal conditions (i.e. β < 1) it is obvious that under photospheric conditions these modes have a much stronger attenuation, the leakage of waves is more accentuated in plasmas with β > 1. Mode (A) has phase speed very close to the sound speed of the upper plasma region and an imaginary part to the solutions with fairly large amplitude, which decreases towards zero, with increasing density ratio. Given its propagation speed, this mode is a fast MHD mode. This mode does not exist when the plasmas in the upper region becomes heavier. In this case, mode (A) does not satisfy the condition imposed on the imaginary part of its frequency. Mode (B) shows a rather interesting characteristics as its propagation speed is sub-Alfvénic for d < 1 (the Alfvén speed in the lower region is taken as reference), and it becomes super-Alfvénic for d > 1. For the entire domain of its definition, the propagation speed of waves stays close to the sound speed c 0 . Given the lower phase speed and the smaller degree of attenuation, we suggest that mode (B) is a slow MHD mode.

Conclusions
Our study considers the problem of waves propagating in a two-dimensional configuration along an interface in the presence of an inclined magnetic field. For a small-angle approximation, the dispersion relation is derived analytically, by employing the effective wavenumbers and using a method of dominant balance. The effective wavenumbers are shown to be complex, so wave amplitude decays laterally in an oscillatory pattern away from the interface. The complex effective wavenumbers, in turn, give rise to complex solutions for the frequency of the waves, which would give amplification or attenuation. However, since no outside energy source exists in this situation, the only physical solutions are leaky waves.
Solutions to the dispersion relation, for varying density ratio and plasma-β, were found numerically. The imaginary components to these solutions imply attenuation of the wave, which corresponds to energy leaking away from the interface, towards |z| → ∞. We have, hence, shown that the introduction of magnetic field inclination introduces energy flow to the system, compared to the case with parallel magnetic field and so even a small angle between the interface and the magnetic field produces a qualitative change in the modes which may propagate. Thus, a contact discontinuity in density and temperature in the presence of an inclined constant magnetic field may support the propagation of surface leaky waves. However, quantities averaged over a thin boundary layer do tend towards the solutions for a tangential discontinuity, as the inclination angle tends towards zero, so the contact and tangential discontinuity solutions are qualitatively comparable.
These results may also have considerable applications to the study of waves in the solar atmosphere. In particular, the penumbra of sunspots have been shown to have highly inclined magnetic field lines and, at high resolution, running penumbral waves were detected (see e.g. Giovanelli, 1972;Zirin and Stein, 1972). Understanding that the outer edge of the sunspot may support leaky waves, more readily than trapped waves, could give a different explanation to any observed damping of these running penumbral waves. While, inside the sunspot, the sharp density variation is horizontal, in the outer edge of the penumbra and the outer canopy, the sharp density gradient is vertical. It has been shown (in e.g. Jess et al., 2013) that field inclination (from the horizontal) in the penumbra is 0 -60 • , but also that the field inclination in the magnetic canopy is 0 -35 • . Hence, small-angle approximation may have relevance to a much wider range of applications, especially across the solar transition region, where the density variation is sudden. This study may also be of particular use for the investigation of transition region quakes (TRQs), since these waves have large vertical length scales, so the transition region may be viewed as a single interface. When considering TRQs as surface waves propagating along a density discontinuity, the fact that interfaces with inclined fields can only support leaky waves may help to explain how energy is transferred into the solar corona, through wave leakage.