© ESO 2021

1. Introduction

Although it has been known for a long time that the corona is three orders of magnitude hotter than the photosphere, a definite explanation for this unexpected feature has continued to elude researchers. One possible explanation that has been advanced over the years, and which is supported by both observations and theoretical models, is the deposition of energy of magnetohydrodynamic (MHD) waves into the corona (Parnell & De Moortel 2012; Arregui 2015; Van Doorsselaere et al. 2020). There are several mechanisms that allow for the conversion of wave energy into heating of the surrounding plasma. Possibilities include dissipation by Ohmic resistivity, values of which were calculated by Kovitya & Cram (1983) for the photosphere and used by Chen et al. (2018, 2021) and Geeraerts et al. (2020) to infer its importance in damping slow waves in a photospheric pore. Mode coupling (Pascoe et al. 2010, 2012; Hollweg et al. 2013; De Moortel et al. 2016) and resonant absorption, both in the Alfvén (Hollweg & Yang 1988; Hollweg et al. 1990; Goossens et al. 1992, 2002; Soler et al. 2013) and cusp (Cadez et al. 1997; Erdélyi et al. 2001; Soler et al. 2009; Yu et al. 2017) continua, have also been studied for their potential to damp waves by transferring their energy to local oscillations. Another mechanism for damping waves in atmospheric structures and which has received considerable attention recently is the Kelvin-Helmholtz instability (KHI). Indeed, the transition of wave energy to turbulence and plasma heating on smaller scales in the coronal plasma is known to be facilitated by this instability (Heyvaerts & Priest 1983; Ofman et al. 1994; Karpen et al. 1994; Karampelas et al. 2017; Afanasyev et al. 2019; Hillier et al. 2020; Shi et al. 2021). The KHI has been observed in the solar atmosphere, for example on coronal mass ejection flanks (Ofman & Thompson 2011; Foullon et al. 2011), and, more recently, on coronal loops (Samanta et al. 2019).

Slow magnetosonic waves were observed in solar coronal loops more than a decade ago, both as propagating waves and standing waves. Upwardly propagating disturbances along coronal loops have been observed, for example, by Berghmans & Clette (1999), Nightingale et al. (1999), and De Moortel et al. (2000). These have been interpreted as propagating slow modes through a theoretical model by Nakariakov et al. (2000). Other observations of disturbances in coronal loops by the SUMER spectrometer have been analyzed and interpreted by Wang et al. (2002, 2003a,b, 2007) as standing slow modes, whereas Kumar et al. (2013) and Mandal et al. (2016) reported the observation of a reflecting, also referred to as sloshing, slow wave in a coronal loop. Wang (2011) provides a review of observations and modeling of standing slow waves in coronal loops. These modes have also been studied through numerical simulations, for example, by De Moortel & Hood (2003), who found that thermal conduction is an important damping mechanism for propagating slow waves in coronal conditions, and by Mandal et al. (2016), who reported about the frequency dependent damping of propagating slow waves in a coronal loop. Structures in the lower solar atmosphere, such as sunspots and pores in the photosphere, have also been observed to harbor slow magnetosonic waves (Dorotovič et al. 2008, 2014; Morton et al. 2011; Grant et al. 2015; Moreels et al. 2015; Freij et al. 2016), which were then classified as either surface or body waves (Moreels et al. 2013; Keys et al. 2018).

The KHI and its growth rate at the interface between two aligned sheared stationary flows have been known for a long time (Chandrasekhar 1961). Although in the purely hydrodynamics (HD) model the instability always develops in the presence of a shear flow, in the MHD model an aligned magnetic field can prevent its triggering. A natural question that then arises in the context of wave heating is under which conditions the KHI would develop in the presence of an oscillating shear flow. Indeed, it is possible that waves propagating or standing in solar coronal structures, such as loops, can trigger the KHI on the boundary and hereby convey some of their energy to the plasma background. Zaqarashvili et al. (2015), for example, studied the KHI in rotating and twisted magnetized jets in the solar atmosphere in the presence of a stationary shear flow and used their derivations to discuss the stability of standing kink and torsional Alfvén waves at the velocity antinodes. They found that the standing waves are always unstable whereas the propagating waves are stable. Transverse oscillations are known to be unstable to the KHI in coronal loops and numerical studies on this topic include Terradas et al. (2008), Antolin et al. (2014, 2015, 2017), Magyar et al. (2015), Guo et al. (2019), Karampelas et al. (2019) and Pascoe et al. (2020).

There have also been several analytical studies regarding the stability of the interface between oscillating sheared flows. Kelly (1965) looked at the HD case of two parallel sheared flows aligned with the interface, whereas Roberts (1973) studied the same setup in the MHD model. More recently, Barbulescu et al. (2019) and Hillier et al. (2019) investigated the stability of transverse oscillations in coronal loops by modeling them locally at the loop boundary as a cartesian interface between sheared background flows, with the background velocity perpendicular to the background magnetic field. Each of these studies revealed that the interface between oscillating sheared flows is always unstable, in contrast to the constant sheared flows case. All of these studies relying on the simplifying assumption that the fluid is incompressible, it is worth asking whether their conclusions remain unchanged when compression is included.

The goal of this paper is therefore twofold. Firstly, it aims at extending the known incompressible model of Roberts (1973) for a plasma with an oscillating velocity field aligned with both the magnetic field and the interface to a compressible version. The focus lies on finding expressions for the eigenfunctions and, in particular, to derive their evolution over time. The interest in doing this lies in identifying the shortcomings made by the approximation of the incompressible model, at the cost of considerably more involved analytical derivations. This is the subject of Sects. 2 and 3. Secondly, it aims at expressing more general instability conditions as compared to the incompressible model, by taking into account the subtleties arising due to the inclusion of compression. In particular, it will allow us to assess the stability of certain solar atmospheric structures harboring a standing slow wave. We will do so by using this model as a local approximation of the structure’s boundary, at the position where the velocity shear is the greatest (for instance in a cusp resonance). We will also compare our findings for the local stability of slow waves to those of Barbulescu et al. (2019) and Hillier et al. (2019) for the local stability of fast kink waves. This is the subject of Sect. 4.

2. Model

We derived an analytical model for a compressible plasma at the boundary of solar atmospheric structures which can, in a first rough approximation, be modeled as a cylinder with a discontinuous boundary separating two regions of homogeneous plasma with different properties. Such structures include coronal loops, prominence threads, sunspots and photospheric pores. The physical setup here is the same as in Roberts (1973), except we included compression. We point out that, in a more realistic setup, a smooth transition layer would have to be included at the interface. This would result in the possibility of resonance occuring between the main oscillation of the structure and local slow mode oscillations in the boundary layer, as studied theoretically by Yu et al. (2017). The analytical derivations of Goossens et al. (2021) show that, when slow waves are resonantly absorbed in the cusp continuum, both the azimuthal component of vorticity and the parallel component of the plasma displacement are large. The huge amount of vorticity could indicate the possibility of the KHI developing in those conditions. However, the sharp spatial variation of the parallel displacement and the truly discontinuous interface separating two plasma regions are absent from the present model. This should be kept in mind when drawing conclusions regarding shear flows in resonances.

In order to be able to make progress in the analytical derivations, the model uses a Cartesian coordinate system (x, y, z), where x = 0 is the interface and the z-direction is the longitudinal direction (i.e., along the cylinder’s axis). The region x < 0 represents the interior of the structure, whereas the region x > 0 represents the surrounding plasma. This is a model for the local stability at the boundary, in the region of the structure where the shear in longitudinal velocity is the greatest, that is to say, at an antinode of the longitudinal component of the velocity eigenfunction. For the longitudinally fundamental slow mode, this would be at the middle of the structure (see Fig. 1), whereas for the first longitudinal overtone, for example, this would be at a quarter of the structure’s length measured from either end. The time variable is denoted by t.

Fig. 1. Sketch of a longitudinal cut along the axis of a coronal loop harboring a fundamental standing slow body sausage mode, at the velocity’s maximal amplitude (left) and the local cartesian model at the boundary (right). The arrows on the left represent the velocity field and their lengths are to scale with the relative local magnitude of the field for a slow body sausage mode at a given time. The lengths and directions of the magnetic field arrows on the right figure are consistent for that same slow mode, whereas for the velocity arrows only the direction is consistent (the length of the exterior velocity arrow having been increased for visual clarification).

After linearizing the ideal MHD equations, the equilibrium quantities are denoted by the subscript 0 and the perturbed quantities are denoted by the subscript 1. The regions on each side of the interface are two homogeneous but different plasmas. This means each region has its own values for the background quantities, which are assumed spatially constant. The background magnetic field is assumed to be a straight and constant axial field along the z-coordinate: B₀ = B_0z1_z. Furthermore, we assume that the background flow is oscillating with a certain frequency ω₀, which represents the frequency of the standing slow wave: v₀ = V₀cos(ω₀t)1_z. It would, of course, be more accurate to also include a background oscillation for the other quantities, such as magnetic field, pressure, and density. In the context of solar atmospheric structures, the magnetic field oscillations that occur because of this external forcing in the background can, however, be neglected in a first approximation with respect to the strong longitudinal magnetic field. As for the density and pressure background oscillations, they can be neglected in this model because they are in antiphase with respect to the velocity in their longitudinal profile and thus have a node where the logitudinal component of velocity has an antinode. For slow modes, the longitudinal component of the velocity is typically much larger than its other components (see left of Fig. 1), which can thus be neglected at the former’s antinode as well.

The perturbed density, thermal pressure, velocity and magnetic field are denoted with ρ₁, p₁, v₁ and B₁, respectively. In what follows the perturbed eulerian total pressure will be used as well, which is given by $P_1=p_1+\frac{{\mathit{B}}_0·{\mathit{B}}_1}{{\mathrm{\mu }}_0}$. Although gravity certainly has a role in solar atmospheric wave dynamics, it is neglected in this model in order to make some analytical progress.

The linearized compressible ideal MHD equations can, under these assumptions, be written as follows:

$$\begin{array}{cc}& {\displaystyle \frac{\mathrm{D}{\rho }_1}{\mathrm{D}t}+{\rho }_0(\mathrm{\nabla }·{\mathit{v}}_1)=0,}\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\rho }_1\frac{\partial {\mathit{v}}_0}{\partial t}+{\rho }_0\frac{\mathrm{D}{\mathit{v}}_1}{\mathrm{D}t}=-\mathrm{\nabla }{P}_1+\frac{1}{{\mathrm{\mu }}_0}({\mathit{B}}_0·\mathrm{\nabla }){\mathit{B}}_1,\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\mathrm{D}{\mathit{B}}_1}{\mathrm{D}t}=-{\mathit{B}}_0(\mathrm{\nabla }·{\mathit{v}}_1)+({\mathit{B}}_0·\mathrm{\nabla }){\mathit{v}}_1,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\mathrm{D}{p}_1}{\mathrm{D}t}+{\rho }_0{v}_{\mathrm{s}}^2(\mathrm{\nabla }·{\mathit{v}}_1)=0,\hfill \end{array}$$(4)

where $\frac{\mathrm{D}f}{\mathrm{D}t}=\frac{\partial f}{\partial t}+({\mathit{v}}_0·\mathrm{\nabla })f$ is the Lagrangian derivative of a quantity f, ${\mathit{v}}_{\mathrm{s}}=\sqrt{\frac{\gamma p_0}{{\rho }_0}}$ is the speed of sound and μ₀ the magnetic permeability of free space. The first equation, Eq. (1), has this simpler form because ∇ρ₀ = 0 in each region (i.e., both inside and outside the cylinder). Since the background quantities depend only on x and t, the perturbed quantities can be Fourier-analyzed in the y and z coordinates and are thus assumed to have the following form: $f_1={\overline{f}}_1(x,t)exp\left\{i(k_{\mathit{y}}\mathit{y}+k_{z}z)\right\}$, for each of the perturbed quantities p₁, ρ₁, v_1x, v_1y, v_1z, B_1x, B_1y and B_1z.

3. Expressions for the perturbed quantities

In this section, we derive the governing equations for the evolution of linear MHD perturbations in a compressible plasma for the model described in the previous section. In the next section, we then try to use the obtained information to describe the stability of the interface in compressible plasma conditions.

3.1. The central quantity ∇ ⋅ ξ

In what follows, an expression is derived for the compression term ∇ ⋅ ξ. This term is central to finding expressions for all other physical quantities, in the sense that they can all be derived solely from it.

By using the Lagrangian displacement ξ defined by $\frac{\mathrm{D}\mathit{\xi }}{\mathrm{D}t}={\mathit{v}}_1$, it can be shown (see Appendix A) that the compression term ∇ ⋅ ξ must satisfy the following partial differential equation:

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\mathrm{D}}^4(\mathrm{\nabla }·\mathit{\xi })}{\mathrm{D}{t}^4}}& +i{k}_z{\omega }_0{V}_0\frac{{\mathrm{D}}^2(sin\left({\omega }_{0}t\right)(\mathrm{\nabla }·\mathit{\xi }))}{\mathrm{D}{t}^2}\hfill \\ \hfill & -(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)\frac{{\mathrm{D}}^2}{D{t}^2}\left(\frac{{\partial }^2(\mathrm{\nabla }·\mathit{\xi })}{\partial x^2}\right)+k^2(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)\frac{{\mathrm{D}}^2(\mathrm{\nabla }·\mathit{\xi })}{\mathrm{D}{t}^2}\hfill \\ \hfill & -i{k}_z{V}_0{\omega }_{0}sin\left({\omega }_{0}t\right)v_{\mathrm{A}}^2\frac{{\partial }^2(\mathrm{\nabla }·\mathit{\xi })}{\partial x^2}\hfill \\ \hfill & +i{k}_z{k}^2{V}_0{\omega }_{0}sin\left({\omega }_{0}t\right)v_{\mathrm{A}}^2(\mathrm{\nabla }·\mathit{\xi })-k_z^2{v}_{\mathrm{A}}^2{v}_{\mathrm{s}}^2\frac{{\partial }^2(\mathrm{\nabla }·\mathit{\xi })}{\partial x^2}\hfill \\ \hfill & +k_z^2{k}^2{v}_{\mathrm{A}}^2{v}_{\mathrm{s}}^2(\mathrm{\nabla }·\mathit{\xi })=0,\hfill \end{array}$$(5)

where $k=\sqrt{k_{\mathit{y}}^2+k_z^2}$ and ${\mathit{v}}_{\mathrm{A}}=B_{0z}/\sqrt{{\mathrm{\mu }}_0{\rho }_0}$ is the Alfvén speed. This partial differential equation is of fourth order in t and second order in x. Note that if V₀ = 0, such that the time dependence of the background quantities disappears, we retrieve an ordinary differential equation (ODE) of second order in x which is the governing equation for the compression term ∇ ⋅ ξ in a plasma without a background flow. This is the equation used for example in Edwin & Roberts (1983) to study wave modes in a magnetic cylinder in the framework of ideal MHD in a plasma without background flow.

By drawing inspiration from the derivations of Barbulescu et al. (2019), one finds that Eq. (5) can take a simpler form if ∇ ⋅ ξ is written as

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{\xi }=f(x,t)g(t)exp\left\{i(k_{y}y+k_{z}z)\right\},\end{array}$$(6)

where $g(t)=exp\left\{\frac{-i{k}_z{V}_{0}sin\left({\omega }_{0}t\right)}{{\omega }_0}\right\}$. Since g(t) has modulus equal to 1 for all t ∈ ℝ, the stability of ∇ ⋅ ξ over time is entirely determined by f. Inserting expression (6) into Eq. (5) simplifies this equation considerably, now taking the following form:

$$\begin{array}{cc}\hfill \frac{{\partial }^{4}f(x,t)}{\partial t^4}& -({v_{\mathrm{A}}}^2+{v_{\mathrm{s}}}^2)\frac{{\partial }^{4}f(x,t)}{\partial x^2\partial t^2}\hfill \\ \hfill & +[i{\omega }_0{k}_z{V}_{0}sin\left({\omega }_{0}t\right)+k^2({v_{\mathrm{A}}}^2+{v_{\mathrm{s}}}^2)]\frac{{\partial }^{2}f(x,t)}{\partial t^2}\hfill \\ \hfill & -[i{\omega }_0{k}_z{V}_0{v_{\mathrm{A}}}^{2}sin\left({\omega }_{0}t\right)+{k_z}^2{v_{\mathrm{A}}}^2{v_{\mathrm{s}}}^2]\frac{{\partial }^{2}f(x,t)}{\partial x^2}\hfill \\ \hfill & +2i{\omega }_0^2{k}_z{V}_{0}cos\left({\omega }_{0}t\right)\frac{\partial f(x,t)}{\partial t}\hfill \\ \hfill & +[i{\omega }_0{k}_z{V}_{0}sin\left({\omega }_{0}t\right)(k^2{v}_{\mathrm{A}}^2-{\omega }_0^2)\hfill \\ \hfill & +k^2{k_z}^2{v_{\mathrm{A}}}^2{v_{\mathrm{s}}}^2]f(x,t)=0.\hfill \end{array}$$(7)

Equation (7) has a solution in the form of F(x)G(t), where F and G satisfy the following ODEs (with m a constant):

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^{2}F(x)}{\mathrm{d}{x}^2}+m^{2}F(x)=0}\end{array}$$(8)

and

$$\begin{array}{cc}& {\displaystyle \frac{{\mathrm{d}}^{4}G(\tau )}{\mathrm{d}{\tau }^4}+(a_1+q_{1}sin\left(\tau \right))\frac{{\mathrm{d}}^{2}G(\tau )}{\mathrm{d}{\tau }^2}+q_{3}cos\left(\tau \right)\frac{\mathrm{d}G(\tau )}{\mathrm{d}\tau }}\hfill \\ \hfill & +(a_2+q_{2}sin\left(\tau \right))G(\tau )=0,\hfill \end{array}$$(9)

with τ = ω₀t, and where

$$\begin{array}{cc}& a_1=\frac{(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)(m^2+k^2)}{{\omega }_0^2},q_1=\frac{i{k}_z{V}_0}{{\omega }_0},\hfill \\ \hfill & {\displaystyle q_3=\frac{2i{k}_z{V}_0}{{\omega }_0},}\hfill \\ \hfill & a_2=\frac{k_z^2{v}_{\mathrm{A}}^2{v}_{\mathrm{s}}^2(m^2+k^2)}{{\omega }_0^4},q_2=\frac{i{k}_z{V}_0[(m^2+k^2)v_{\mathrm{A}}^2-{\omega }_0^2]}{{\omega }_0^3}·\hfill \end{array}$$

These two equations are related by the constant m, which occurs in both of them.

3.1.1. The spatial function F

Equation (8) has a simple analytical solution, namely

$$\begin{array}{c}\hfill F(x)=C_{1}exp\left(imx\right)+C_{2}exp(-imx),\end{array}$$(10)

with C₁ and C₂ arbitrary constants. From this it can be inferred that m plays the role of the wavenumber along x, which is the direction normal to the interface. The focus of this paper goes to standing slow waves in solar atmospheric structures that can be modeled as a straight cylinder with a circular base and a discontinuous boundary separating two homogeneous but different plasma regions. Therefore, each region has its own value for m (namely m_i for the interior and m_e for the exterior), each being able to only take on specific values depending on the boundary conditions at the interface. Finding expressions for m_i and m_e is, however, not straightforward. This will be discussed in Sect. 3.1.2. We note that for the purpose of studying the local stability of the interface, m will be taken purely imaginary in order to have evanescent spatial profiles in the normal direction.

3.1.2. The temporal function G

Equation (9) does not have a simple analytical solution. This equation governs the stability of the compression term ∇ ⋅ ξ over time, its parameters depending on the background quantities v_A, v_s, V₀ and ω₀. The function h(t) = G(t)g(t), where G is the solution to Eq. (9), actually describes the general time evolution of ∇ ⋅ ξ in a compressible homogeneous plasma of inifinite extent with an oscillating background velocity which is parallel to the straight and constant equilibrium magnetic field. Although there is no closed-form analytical solution, the boundedness of ∇ ⋅ ξ over time can be derived from the properties of Eq. (9). Indeed, this linear fourth-order ODE has periodic coefficients with the same period and hence it falls in the category of equations which obey Floquet theory (Chicone 2008).

Defining τ = ω₀t, it can be rewritten as a 4 × 4 system of first-order ODEs of the form x′(τ) = A(τ)x(τ), where the coefficient matrix A(τ) is periodic with a period of 2π. One can find four linearly independent fundamental solution vectors x₁(τ), x₂(τ), x₃(τ) and x₄(τ), which depend on their respective initial conditions. The matrix obtained by putting these four vectors into its columns is called a fundamental solution matrix. If we take as initial condition matrix the identity matrix, Floquet theory states that the corresponding fundamental solution matrix (which we denote by X(τ)) evaluated at one period (i.e., at τ = 2π in our case) is intimately linked to the boundedness of the solutions of the ODE. Indeed, denoting the eigenvalues of X(2π) by λ₁, λ₂, λ₃ and λ₄, Floquet’s theorem states that for each distinct eigenvalue λ, if we write λ = e^2πμ where μ ∈ ℂ is called the Floquet exponent, there exists an independent solution of the system which has the form

$$\begin{array}{c}\hfill \mathit{x}(\tau )=e^{\mathrm{\mu }\tau }\mathbf{\Phi }(\tau ),\end{array}$$(11)

where Φ contains 2π-periodic functions. This means that a solution to Eq. (9) is unstable if and only if |λ| > 1 for at least one eigenvalue λ of X(2π). If there are less than four distinct eigenvalues of X(2π), it is possible there are less than four independent Floquet solutions of the form (11). If the eigenspace of X(2π) has dimension 4, there are still four independent eigenvectors Φ_i and thus four independent Floquet solutions. If the eigenspace of X(2π) has dimension less than 4, then there are less than four independent Floquet solutions. In that case, extra independent solutions have to be found by using the Jordan normal form of X(2π) (see e.g., Cesari 1963 for more details). These extra independent Jordan solutions are always unstable. The stability of ∇ ⋅ ξ is thus entirely determined by the Floquet exponents μ.

A solution to Eq. (9) can be written in the form of a series. Knowing that each independent Floquet solution has the form e^μτΦ(τ) with Φ a 2π-periodic function, one can assume that such a basic independent solution can be written as

$$\begin{array}{c}\hfill {\displaystyle G(\tau )=e^{\mathrm{\mu }\tau }{\displaystyle \sum _{j=-\infty }^{\infty }{\phi }_j{e}^{ij\tau },}}\end{array}$$(12)

where the φ_j are unknown coefficients. Writing μ = iν for convenience, the following recurence relation between the coefficients φ_j can then be derived by inserting expression (12) into Eq. (9):

$$\begin{array}{c}\hfill -{\epsilon }_j{\phi }_{j-1}+{\phi }_j+{\epsilon }_j{\phi }_{j+1}=0,\end{array}$$(13)

where

$$\begin{array}{c}\hfill {\displaystyle {\epsilon }_j=\frac{\frac{1}{2}{k}_z{V}_0{\omega }_0[{(j+\nu )}^2{\omega }_0^2-K^2{v}_{\mathrm{A}}^2]}{{(j+\nu )}^4{\omega }_0^4-{(j+\nu )}^2{\omega }_0^2{K}^2(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)+k_z^2{K}^2{v}_{\mathrm{A}}^2{v}_{\mathrm{s}}^2}}\end{array}$$(14)

with $K=\sqrt{k^2+m^2}$. Equation (13) are nontrivially satisfied if and only if the following infinite determinant vanishes:

$$\begin{array}{c}\hfill \mathrm{\Delta }=\left|\begin{array}{cccc}\ddots & \ddots & & 0\\ \ddots & 1& {\epsilon }_{-1}& \\ & -{\epsilon }_0& 1& {\epsilon }_0& \\ & -{\epsilon }_1& 1& \ddots \\ 0& & \ddots & \ddots \end{array}\right|.\end{array}$$(15)

This is a Fredholm determinant. Denoting with A the operator defined by the corresponding infinite matrix, this Fredholm determinant is well-defined if the operator A − I (where I is the identity operator) is a trace class operator on ℓ²(ℤ), the Hilbert space of square-summable sequences of complex numbers with entire index. It can be shown that this is the case here (see for example Sträng 2005 for the method to follow), except if the denominator of one of the ε_j vanishes. This happens if the perturbation, which is a normal mode with wave vector K = (m, k_y, k_z), satisfies the equation

$$\begin{array}{c}\hfill {\displaystyle {(j+\nu )}^2{\omega }_0^2=\frac{K^2(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)}{2}\{1\pm {\left[{\displaystyle 1-\frac{4k_z^2{v}_{\mathrm{A}}^2{v}_{\mathrm{s}}^2}{K^2{(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)}^2}}\right]}^{1/2}\},}\end{array}$$(16)

for a j ∈ ℤ. This represents a resonance between the background oscillator with frequency ω₀ and a magnetosonic wave, in a homogeneous plasma of infinite extent. In this case m is a free parameter (like k_y and k_z) and can potentially take any real value. In the model we describe in this paper, with an interface separating two such homogeneous plasmas, only certain specific values of m (different in both regions) are physically possible. Because surface waves on the interface are actually of interest in this case, m should be imaginary. In Sect. 4 we study the resonance of these surface waves with the background oscillator.

Considering Δ as a function of the two unknowns ν and m, it is thus well-defined for every (ν, m) except in the poles of the ε_j. It can also easily be checked with Eqs. (14) and (15) that Δ(ν, m) = Δ(ν + 1, m) and Δ(ν, m) = Δ(−ν, m). As previously stated, there are in general four independent Floquet solutions of the form (12) to Eq. (9). Hence, all solutions for ν of Δ(ν, m) = 0 (as functions of m) which correspond to distinct solutions of the differential equation lie on the strip −0.5 < Re[ν] ≤ 0.5. These distinct solutions relate as follows: ν₂ = −ν₁ and ν₄ = −ν₃. The four independent Floquet solutions in the most general case are thus of the form e^μ₁τΦ₁(τ), e^−μ₁τΦ₂(τ), e^μ₃τΦ₃(τ) and e^−μ₃τΦ₄(τ), where Φ₁, Φ₂, Φ₃ and Φ₄ are 2π-periodic functions. Recalling the earlier discussion in this section, we note that it is possible that there are less than four independent Floquet solutions if at least two eigenvalues λ_j = e^2πiν_j of X(2π) are equal. From the properties of Δ, we can see that this happens in two situations. One possibility is if we have Re[ν_j] = n/2 for some n ∈ ℤ and Im[ν_j] = 0, for one of the solutions ν_j. Indeed, the eigenvalues e^2πiν_j and e^−2πiν_j are equal in that case. Another possibility is if, for two solutions ν_j and ν_l, we have Re[ν_j] = Re[ν_l]+n for some n ∈ ℤ and Im[ν_j] = Im[ν_l]. In this case e^2πiν_j and e^2πiν_l will also be equal.

We find the following formula for Δ, derived in Appendix B:

$$\begin{array}{c}\hfill {\displaystyle \mathrm{\Delta }=1+\sum _{n=1}^{\infty }\left[{\displaystyle \sum _{j_1=-\infty }^{\infty }\sum _{j_2=j_1+2}^{\infty }\dots \sum _{j_n=j_{n-1}+2}^{\infty }\left({\displaystyle \prod _{l=1}^n{\epsilon }_{j_l}{\epsilon }_{j_l+1}}\right)}\right].}\end{array}$$(17)

It is clear from Eq. (14) that ${\epsilon }_j\approx \frac{k_z{V}_0}{2j^2{\omega }_0}$ as |j| → ∞. We can then try to use the fact that ε_j drops off as 1/j² to approximate the cumbersome formula (17). If we can assume that |Im[ν]| ≪ |Re[ν]|, then since −0.5 < Re[ν] ≤ 0.5 the quantity ν becomes negligible with respect to j already for quite small |j| in that case. In Sect. 4, we see that this is a good assumption in both coronal and photospheric conditions, if the Alfvén Mach number M_A = (V_0i − V_0e)/v_Ai is small enough.

One could then try considering, as a first approximation for Δ, only the first few terms on the right-hand side of Eq. (17). Taking n = 1 and j₁ ∈ { − 1, 0, 1}, this would yield:

$$\begin{array}{c}\hfill \mathrm{\Delta }\approx 1+{\epsilon }_{-1}{\epsilon }_0+{\epsilon }_0{\epsilon }_1.\end{array}$$(18)

Equation (18) corresponds to Eq. (15) with the infinite determinant in the right-hand side truncated to a 3 × 3 determinant centered on the row with ε₀. If one starts from the incompressible versions of Eqs. (1)–(4), one can derive that m = ik in an incompressible plasma. The approximation of truncating the series in Eq. (17) is only valid for large enough |j| and away from the poles of the ε_j. Under the assumption that |K|v_A ≪ ω₀ and |K|v_s ≪ ω₀, this condition is fulfilled. We note that in the solar atmosphere, v_A and v_s are rougly of the same order of magnitude. Since K = k² − |m|² for surface waves, we have K → 0 in the incompressible limit. Equation (18) thus gives us an approximation for Δ at least in a weakly compressible plasma, but could maybe even be correct more generally.

Analytical solutions for m_i and m_e in function of ν can then be derived from Eq. (18). The obtained expressions for m_i(ν) and m_e(ν) are very complicated but, introducing numerical values relevant for the physical conditions of the solar structure being considered, we can use them together with a third relation involving ν, m_i and m_e. This other relation can theoretically be the one derived in the next subsection, although in practice one of those we derive in Sect. 4 under approximating circumstances will probably be preferred. We note that ν is the same on both sides of the interface, as will be explained in the next section. In contrast to this, m is in general not identical on both sides. It can also be seen that Δ(ν, m) = Δ(ν, −m). Hence, if m is a solution then so is −m. This is reflected in the form of the spatial function F in Eq. (10).

3.2. Other perturbed physical quantities

With h(t) = G(t)g(t), the compression term can be written as ∇ ⋅ ξ = F(x)h(t)exp{i(k_yy + k_zz)}. From this and Eqs. (1)–(4), the following expressions for the perturbed quantities can then easily be derived:

$$\begin{array}{cc}& {\displaystyle {\xi }_x=-\frac{\mathit{im}}{k^2+m^2}\stackrel{\sim }{F}(x)h_1(t)exp\left\{i(k_{y}y+k_{z}z)\right\},}\hfill \end{array}$$(19)

$$\begin{array}{cc}& {\displaystyle {\xi }_y=-\frac{i{k}_y}{k^2+m^2}F(x)h_2(t)exp\left\{i(k_{y}y+k_{z}z)\right\},}\hfill \end{array}$$(20)

$$\begin{array}{cc}& {\displaystyle {\xi }_z=-\frac{i{k}_z}{k^2+m^2}F(x)h_3(t)exp\left\{i(k_{y}y+k_{z}z)\right\},}\hfill \end{array}$$(21)

$$\begin{array}{cc}& {\displaystyle B_{1x}=B_{0z}\frac{k_{z}m}{k^2+m^2}\stackrel{\sim }{F}(x)h_1(t)exp\left\{i(k_{y}y+k_{z}z)\right\},}\hfill \end{array}$$(22)

$$\begin{array}{cc}& {\displaystyle B_{1y}=B_{0z}\frac{k_z{k}_y}{k^2+m^2}F(x)h_2(t)exp\left\{i(k_{y}y+k_{z}z)\right\},}\hfill \end{array}$$(23)

$$\begin{array}{cc}& B_{1z}=B_{0z}F(x)\left({\displaystyle \frac{k_z^2}{k^2+m^2}{h}_3(t)-h(t)}\right)exp\left\{i(k_{y}y+k_{z}z)\right\},\hfill \end{array}$$(24)

$$\begin{array}{cc}& p_1=-{\rho }_0{v}_{\mathrm{s}}^{2}F(x)h(t)exp\left\{i(k_{y}y+k_{z}z)\right\},\hfill \end{array}$$(25)

$$\begin{array}{cc}& {\rho }_1=-{\rho }_{0}F(x)h(t)exp\left\{i(k_{y}y+k_{z}z)\right\},\hfill \end{array}$$(26)

$$\begin{array}{cc}& P_1={\rho }_0\left({\displaystyle \frac{k_z^2{v}_{\mathrm{A}}^2}{k^2+m^2}{h}_3(t)-(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)h(t)}\right)\hfill \\ \hfill & F(x)exp\left\{i(k_{y}y+k_{z}z)\right\},\hfill \end{array}$$(27)

with $\stackrel{\sim }{F}(x)=C_1{e}^{\mathit{imx}}-C_2{e}^{-imx}$ and where h₁, h₂, h₃ are still unknown but have to satisfy

$$\begin{array}{c}\hfill {\displaystyle \frac{m^2}{k^2+m^2}{h}_1(t)+\frac{k_y^2}{k^2+m^2}{h}_2(t)+\frac{k_z^2}{k^2+m^2}{h}_3(t)=h(t),}\end{array}$$(28)

by the definition ∇ ⋅ ξ = ∂ξ_x/∂x + ∂ξ_y/∂y + ∂ξ_z/∂z. Each of the expressions (19)–(27) is different for each region on both sides of the interface. The functions $\stackrel{\sim }{F}$, F, h, h₁, h₂ and h₃ in particular will have a different version in both regions. For the spatial functions $\stackrel{\sim }{F}$ and F, only one of the two coefficients C₁ and C₂ must be retained in each region. We make the choice that C₁ = 0 in the x < 0 region, and C₂ = 0 in the x > 0 region.

Writing h₁(t) = G₁(t)g(t), h₂(t) = G₂(t)g(t) and h₃(t) = G₃(t)g(t), the following expressions for G₁, G₂ and G₃ can be derived from Eqs. (1)–(3), (19)–(21) and (27) (see Appendix C):

$$\begin{array}{cc}\hfill G_1(\tau )& =G_2(\tau )\hfill \\ \hfill & {\displaystyle =K^2{v}_{\mathrm{s}}^2{e}^{\mathrm{\mu }\tau }\sum _{j=-\infty }^{\infty }\frac{1}{{(j+\nu )}^2{\omega }_0^2-K^2{v}_{\mathrm{A}}^2}{\phi }_j{e}^{ij\tau },}\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill G_3(\tau )& {\displaystyle =\frac{K^2}{k_z^2}{e}^{\mathrm{\mu }\tau }\sum _{j=-\infty }^{\infty }(1-\frac{(k_y^2+m^2)v_{\mathrm{s}}^2}{{(j+\nu )}^2{\omega }_0^2-K^2{v}_{\mathrm{A}}^2}){\phi }_j{e}^{ij\tau }.}\hfill \end{array}$$(30)

It can also be checked that Eq. (28) is indeed fulfilled with these expressions. Now, the following two boundary conditions have to be fulfilled at the interface between the two regions (i.e., at x = 0) for physical reasons:

$$\begin{array}{cc}& [{\xi }_x]=0,\hfill \end{array}$$(31)

$$\begin{array}{cc}& [P_1]=0,\hfill \end{array}$$(32)

where [f] = lim_x ↓ 0f(x)−lim_x ↑ 0f(x) denotes the jump in a quantity f across the interface. We note that, similarly as for the frequency of a normal mode in the case without background flow, ν is the same inside and outside. It has proven to be too difficult to show mathematically, but it can be explained as follows. When a perturbation is unstable, its growth rate is determined solely by ν: the growth of every perturbed quantity is namely expressed by the factor e^−Im[ν]t. Therefore, since quantities such as ξ_x and P₁ have to be continuous at the interface x = 0, ν must be the same on both sides of the interface.

From the two Eqs. (31) and (32) linking the interior and exterior solutions, a relation can be derived which has to be satisfied in order for the system determined by these equations to have nontrivial solutions:

$$\begin{array}{c}\hfill {\displaystyle \sum _{j=-\infty }^{\infty }({\rho }_{0\mathrm{i}}{m}_{\mathrm{e}}{\zeta }_{A,j}+{\rho }_{0\mathrm{e}}{m}_{\mathrm{i}}{\zeta }_{B,j})e^{ij\tau }=0,}\end{array}$$(33)

where

$$\begin{array}{cc}& {\displaystyle {\zeta }_{A,j}=\sum _{l=-\infty }^{\infty }\frac{2{\omega }_0^2{v}_{\mathrm{si}}^2{v}_{\mathrm{se}}^2[{(l+\nu )}^2{\omega }_0^2-k_z^2{v}_{\mathrm{Ai}}^2]{\phi }_{l,\mathrm{i}}{\phi }_{j-l,\mathrm{e}}}{[{(l+\nu )}^2{\omega }_0^2-K_{\mathrm{i}}^2{v}_{\mathrm{Ai}}^2][{(j-l+\nu )}^2{\omega }_0^2-K_{\mathrm{e}}^2{v}_{\mathrm{Ae}}^2]},}\hfill \\ \hfill & {\displaystyle {\zeta }_{B,j}=\sum _{l=-\infty }^{\infty }\frac{2{\omega }_0^2{v}_{\mathrm{si}}^2{v}_{\mathrm{se}}^2[{(l+\nu )}^2{\omega }_0^2-k_z^2{v}_{\mathrm{Ae}}^2]{\phi }_{l,\mathrm{e}}{\phi }_{j-l,\mathrm{i}}}{[{(l+\nu )}^2{\omega }_0^2-K_{\mathrm{e}}^2{v}_{\mathrm{Ae}}^2][{(j-l+\nu )}^2{\omega }_0^2-K_{\mathrm{i}}^2{v}_{\mathrm{Ai}}^2]}.}\hfill \end{array}$$

Equation (33), arising from the boundary conditions at the interface, is usually called the dispersion relation when there is no oscillating background flow. Both m_i and m_e being determined by Eq. (15) (by inserting respectively interior and exterior values for the background quantities), Eq. (33) determines the only remaining unknown, ν.

From Eq. (33), the following has to hold for every j ∈ ℤ:

$$\begin{array}{c}\hfill {\displaystyle \frac{m_{\mathrm{e}}}{m_{\mathrm{i}}}=-\frac{{\rho }_{0\mathrm{e}}{\zeta }_{B,j}}{{\rho }_{0\mathrm{i}}{\zeta }_{A,j}}·}\end{array}$$(34)

For j = 0 this gives us the relation

$$\begin{array}{c}\hfill {\displaystyle \frac{m_{\mathrm{e}}}{m_{\mathrm{i}}}=-\frac{{\rho }_{0\mathrm{e}}}{{\rho }_{0\mathrm{i}}}\frac{{\displaystyle \sum _{l=-\infty }^{\infty }\frac{[{(l+\nu )}^2{\omega }_0^2-k_z^2{v}_{\mathrm{Ae}}^2]{\phi }_{l,\mathrm{e}}{\phi }_{-l,\mathrm{i}}}{[{(l+\nu )}^2{\omega }_0^2-K_{\mathrm{e}}^2{v}_{\mathrm{Ae}}^2][{(l-\nu )}^2{\omega }_0^2-K_{\mathrm{i}}^2{v}_{\mathrm{Ai}}^2]}}}{{\displaystyle \sum _{l=-\infty }^{\infty }\frac{[{(l-\nu )}^2{\omega }_0^2-k_z^2{v}_{\mathrm{Ai}}^2]{\phi }_{l,\mathrm{e}}{\phi }_{-l,\mathrm{i}}}{[{(l+\nu )}^2{\omega }_0^2-K_{\mathrm{e}}^2{v}_{\mathrm{Ae}}^2][{(l-\nu )}^2{\omega }_0^2-K_{\mathrm{i}}^2{v}_{\mathrm{Ai}}^2]}}}·}\end{array}$$(35)

We note that this is not an explicit solution for m_e/m_i, because both m_i and m_e appear on the right-hand side, through K_i and K_e respectively. While Eq. (35) seems difficult to use directly, it gives us some information. Indeed, we learn from this equation that, whereas in the incompressible case we have m_e/m_i = 1, in the compressible case m_e/m_i seems to depend on k_z as well as k_y, the latter appearing only in K_i and K_e. In the next section, we see that the only perturbation quantities the stability of the interface depends on are k_z and m_e/m_i. This means that, while the abstraction is made that only the longitudinal wavenumber k_z is important for the stability of the interface in an incompressible plasma, the perpendicular wavenumber k_y seems to have some influence as well in the more realistic case of a compressible plasma.

4. Stability of the interface

4.1. Governing equation

The results derived in the preceding sections permit us to say something about the stability of the interface at the structure’s boundary. This stability is determined by the temporal evolution of the normal displacement ξ_x, which is governed by the following equation:

$$\begin{array}{c}\hfill {\displaystyle \frac{D^2{\xi }_x}{\mathrm{D}{t}^2}+k_z^2{v}_{\mathrm{A}}^2{\xi }_x=\frac{-1}{{\rho }_0}\frac{\partial P_1}{\partial x}·}\end{array}$$(36)

There is again a different version of Eq. (36) for each region, namely for x < 0 and for x > 0. For each of the two regions, the respective versions of expression (27) for P₁ can then be inserted in the respective versions of Eq. (36). These can then be put together by expressing C₁ in function of C₂ thanks to Eq. (32), and yield the following equation governing the displacement of the interface over time:

$$\begin{array}{cc}& \left({\displaystyle {\rho }_{0\mathrm{e}}{m}_{\mathrm{i}}\frac{D_{\mathrm{e}}^2\left({\xi }_x{|}_{x=0}\right)}{\mathrm{D}{t}^2}+{\rho }_{0\mathrm{i}}{m}_{\mathrm{e}}\frac{D_{\mathrm{i}}^2\left({\xi }_x{|}_{x=0}\right)}{\mathrm{D}{t}^2}}\right)\hfill \\ \hfill & +k_z^2({\rho }_{0\mathrm{e}}{v}_{\mathrm{Ae}}^2{m}_{\mathrm{i}}{\xi }_x{|}_{x=0}+{\rho }_{0\mathrm{i}}{v}_{\mathrm{Ai}}^2{m}_{\mathrm{e}}{\xi }_x{|}_{x=0})=0,\hfill \end{array}$$(37)

where D_i/Dt = ∂/∂t + ik_zV_0isin(ω₀t), D_e/Dt = ∂/∂t + ik_zV_0esin(ω₀t), and ξ_x|_x = 0 is ξ_x evaluated at x = 0.

Using a similar trick as was used for ∇ ⋅ ξ before and which was also used by Barbulescu et al. (2019) in a different setup, we write ξ_x|_x = 0 as

$$\begin{array}{c}\hfill {\xi }_x{|}_{x=0}(t)=\eta (t)exp\left\{\frac{-iAsin\left({\omega }_{0}t\right)}{{\omega }_0}\right\},\end{array}$$(38)

with $A=k_z\frac{V_{0\mathrm{e}}{m}_{\mathrm{i}}{\rho }_{0\mathrm{e}}+V_{0\mathrm{i}}{m}_{\mathrm{e}}{\rho }_{0\mathrm{i}}}{m_{\mathrm{i}}{\rho }_{0\mathrm{e}}+m_{\mathrm{e}}{\rho }_{0\mathrm{i}}}$. This assumption changes nothing to the stability of ξ_x|_x = 0(t) since the moduli of ξ_x|_x = 0(t) and η(t) are the same for every t ∈ ℝ⁺, and therefore ξ_x|_x = 0 is unbounded if and only if η is unbounded. Inserting Eq. (38) in Eq. (37) greatly simplifies the equation, which can be worked out to yield the following ODE for η:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2\eta (\tau )}{\mathrm{d}{\tau }^2}+(a-2qcos\left(2\tau \right))\eta (\tau )=0,}\end{array}$$(39)

with τ = ω₀t and

$$\begin{array}{cc}& {\displaystyle a=\frac{k_z^2}{{\omega }_0^2}\frac{{\rho }_{0\mathrm{i}}{v}_{\mathrm{Ai}}^2{m}_{\mathrm{e}}+{\rho }_{0\mathrm{e}}{v}_{\mathrm{Ae}}^2{m}_{\mathrm{i}}}{{\rho }_{0\mathrm{i}}{m}_{\mathrm{e}}+{\rho }_{0\mathrm{e}}{m}_{\mathrm{i}}}-\frac{k_z^2{\rho }_{0\mathrm{i}}{\rho }_{0\mathrm{e}}{m}_{\mathrm{i}}{m}_{\mathrm{e}}{(V_{0\mathrm{i}}-V_{0\mathrm{e}})}^2}{2{\omega }_0^2{({\rho }_{0\mathrm{i}}{m}_{\mathrm{e}}+{\rho }_{0\mathrm{e}}{m}_{\mathrm{i}})}^2},}\hfill \end{array}$$(40)

$$\begin{array}{cc}& {\displaystyle q=\frac{k_z^2{\rho }_{0\mathrm{i}}{\rho }_{0\mathrm{e}}{m}_{\mathrm{i}}{m}_{\mathrm{e}}{(V_{0\mathrm{i}}-V_{0\mathrm{e}})}^2}{4{\omega }_0^2{({\rho }_{0\mathrm{i}}{m}_{\mathrm{e}}+{\rho }_{0\mathrm{e}}{m}_{\mathrm{i}})}^2}·}\hfill \end{array}$$(41)

The parameters a and q can be rewritten in normalized form as follows:

$$\begin{array}{cc}& a={\kappa }_z^2\left[{\displaystyle \frac{\overline{m}+r{\overline{v}}_{\mathrm{A}}^2}{\overline{m}+r}-\frac{r\overline{m}{M}_{\mathrm{A}}^2}{2{(\overline{m}+r)}^2}}\right],\hfill \end{array}$$(42)

$$\begin{array}{cc}& {\displaystyle q={\kappa }_z^2\frac{r\overline{m}{M}_{\mathrm{A}}^2}{4{(\overline{m}+r)}^2},}\hfill \end{array}$$(43)

where we introduced the normalized quantities κ_z = k_zv_Ai/ω₀, ${\overline{\mathit{v}}}_{\mathrm{A}}={\mathit{v}}_{\mathrm{Ae}}/{\mathit{v}}_{\mathrm{Ai}}$, r = ρ_0e/ρ_0i, $\overline{m}=m_{\mathrm{e}}/m_{\mathrm{i}}$ and M_A = (V_0i − V_0e)/v_Ai.

We notice that the first term of a consists of an expression that resembles the squared kink frequency and which we call the pseudo squared kink expression. The second term of a resembles a Doppler-shifting correction and equals −2q. In an incompressible plasma, the first term would be the squared kink frequency (since then m_i, m_e → k). We also note that if there is no background flow (i.e., V_0i = V_0e = 0), we are left with only the pseudo squared kink expression in the factor in front of η(τ) in Eq. (39). The solution is then a surface wave with frequency ν determined by the dispersion relation

$$\begin{array}{c}\hfill {\displaystyle {\nu }^2={\kappa }_z^2\frac{\overline{m}+r{\overline{v}}_{\mathrm{A}}^2}{\overline{m}+r},}\end{array}$$(44)

where m_i and m_e are determined in function of ν through their respective version of Eq. (17). Equation (44) is a transcendental equation in ν with more than one solution, which means that there are several different modes in a compressible plasma. This is in contrast with an incompressible plasma, where there is only one mode, for which the expression of the frequency is readily readable from Eq. (44) and is equal to the kink frequency. These results about modes at an interface in a compressible plasma without background flow are already known (see for example Priest 2014).

4.2. Instability conditions

Equation (39) is the Mathieu equation, which also pertains to the class of equations obeying Floquet theory and which has been extensively studied for example by McLachlan (1947). It does not have an analytical solution but the stability of its solution in function of its parameters is well-known. The stability of the solutions of the Mathieu equation (39) in function of the parameters a and q is often represented in the so-called stability diagram (see Fig. 2). The white zones in the figure are regions where the solutions are stable, whereas the gray zones are region where the solutions are unstable.

Fig. 2. Stability diagram of the Mathieu equation. The white zones represent regions of the parameter space where the solution is stable, whereas the gray zones represent regions of the parameter space where the solution is unstable.

4.2.1. Constant shear flow

In order to gain some insight on the nature of the potential instabilities, we follow Hillier et al. (2019), who studied the stability of an interface in a similar model. We start by looking at the case with a constant background shear flow (i.e., we assume ω₀ = 0). We have to modify the procedure above a little bit in order to be mathematically correct: instead of Eq. (38), we assume ξ_x|_x = 0(t) = η(t)exp{−iAt} in this particular case. The obtained equation is then

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{d}}^2\eta (\tau )}{\mathrm{d}{\tau }^2}+{\kappa }_z^2\left[{\displaystyle \frac{\overline{m}+r{\overline{v}}_{\mathrm{A}}^2}{\overline{m}+r}-\frac{r\overline{m}{M}_{\mathrm{A}}^2}{{(\overline{m}+r)}^2}}\right]\eta (\tau )=0.}\end{array}$$(45)

It has normal mode solutions with normalized frequency ν defined by the dispersion relation

$$\begin{array}{c}\hfill {\nu }^2={\kappa }_z^2\left[{\displaystyle \frac{\overline{m}+r{\overline{v}}_{\mathrm{A}}^2}{\overline{m}+r}-\frac{r\overline{m}{M}_{\mathrm{A}}^2}{{(\overline{m}+r)}^2}}\right]·\end{array}$$(46)

This transcendental equation in ν has multiple solutions as well, and thus we see that there are also different modes in the case of a constant shear flow in a compressible plasma. A particular mode will be unstable to the KHI if

$$\begin{array}{c}\hfill {\displaystyle M_{\mathrm{A}}^2>\frac{(\overline{m}+r)(\overline{m}+r{\overline{v}}_{\mathrm{A}}^2)}{r\overline{m}}}\end{array}$$(47)

is satisfied. Different modes will thus have different instability conditions in a compressible plasma. If we take the incompressible limit, the above derivations match the already well-known results of a constant shear flow in an incompressible plasma (see for example Chandrasekhar 1961).

4.2.2. Oscillating shear flow

In the case of an oscillatory background flow (ω₀ ≠ 0), Eq. (39) governs the evolution of the normal displacement at the interface over time. In the same way as for Eq. (9), Floquet theory allows us to write an independent solution to the Mathieu equation as follows (McLachlan 1947):

$$\begin{array}{c}\hfill {\displaystyle \eta (\tau )=e^{i\nu \tau }\sum _{j=-\infty }^{\infty }{\varphi }_j{e}^{2ij\tau }.}\end{array}$$(48)

This resembles a normal mode with frequency ν, but with a periodic function replacing the constant. One can then find the following recursion relation by inserting Eq. (48) into Eq. (39):

$$\begin{array}{c}\hfill {ϵ}_j{\varphi }_{j-1}+{\varphi }_j+{ϵ}_j{\varphi }_{j+1}=0,\end{array}$$(49)

where

$$\begin{array}{c}\hfill {\displaystyle {ϵ}_j=\frac{q}{{(2j+\nu )}^2-a}·}\end{array}$$(50)

If the denominator of one of the ϵ_j vanishes, there can be a resonance between the background oscillator and the induced surface waves.

Now some analytical progress can be made in the limiting case q → 0, corresponding to a small squared shear rate $k_z^2{(V_{\text{0i}}-V_{\text{0e}})}^2$ with respect to the squared frequency ${\omega }_0^2$. Indeed, for q ≪ 1, we see from Eq. (49) that in order to have a nontrivial solution (48) we need the denominator of at least one of the ϵ_j to vanish. For that, we need to have

$$\begin{array}{c}\hfill {(2j+\nu )}^2-a\approx 0\end{array}$$(51)

for a j ∈ ℤ. In case ν² ≈ a, that is to say, if

$$\begin{array}{c}\hfill {\nu }^2\approx {\kappa }_z^2\left[{\displaystyle \frac{\overline{m}+r{\overline{v}}_{\mathrm{A}}^2}{\overline{m}+r}-\frac{r\overline{m}{M}_{\mathrm{A}}^2}{2{(\overline{m}+r)}^2}}\right],\end{array}$$(52)

then Eq. (16) is fulfilled for j = 0. Since Eq. (52) is again a transcendental equation in ν with multiple solutions, there are multiple modes possible. These modes are surface waves with frequency ν. Such a mode will be unstable if

$$\begin{array}{c}\hfill {\displaystyle M_{\mathrm{A}}^2>2\frac{(\overline{m}+r)(\overline{m}+r{\overline{v}}_{\mathrm{A}}^2)}{r\overline{m}}}\end{array}$$(53)

is satisfied. This renders the right-hand side in Eq. (52), which is a, negative. Hence, for q ≪ 1, the region of the parameters space which corresponds to the KHI is a < 0. It can be seen on Fig. 2 that this region is indeed uniformly unstable. We note that the right-hand side in Eq. (53) is a factor of 2 larger than in Eq. (47). In a first approximation, we could thus evaluate the minimum value of M_A for the KHI to develop in the case of an oscillatory shear flow with q ≪ 1 to be $\sqrt{2}$ times higher than in the case of a constant shear flow with q ≪ 1. This is assuming that the corresponding values of m_i and m_e are equal in both cases. In an incompressible plasma, this is obviously correct. In a compressible plasma, however, there is no reason to suggest a priori that this is the case in general. Nevertheless, the values might still be close, such as in a weakly compressible plasma for example.

If on top of Eq. (51) being satistfied for j = 0 it is also satisfied for another j ∈ ℤ, then there is a resonance between a mode satisfying Eq. (52) and the background oscillator. From ν² = a and Eq. (16) with j ≠ 0, this means that we must have

$$\begin{array}{c}\hfill a=j^2.\end{array}$$(54)

In Fig. 2 we see the regions of resonance instability around a = j², with j ∈ ℤ₀. Hence, for q ≪ 1, the KHI and the resonance instability are clearly distinct features and are mutually exclusive.

As was also mentionned by Hillier et al. (2019) (themselves based on Bender & Orszag 1999), the dominant resonance in this limit is the one with j = 1. Its instability reagion is borded by the curves

$$\begin{array}{c}\hfill a=j^2\pm q+O(q^2)\mathrm{as}q\to 0\end{array}$$(55)

and its maximum growth rate is

$$\begin{array}{c}\hfill {\displaystyle \mathrm{Im}[\nu ]=\frac{q}{2}+O(q^2)\mathrm{as}q\to 0.}\end{array}$$(56)

We could thus express this growth rate through the equation $\mathrm{Im}[\nu ]=\frac{q}{2}$ as an approximation. This is an implicit equation in ν, because q depends on it through $\overline{m}$.

4.3. Solar applications: Relevant region of the parameter space

From Eqs. (42) and (43), we see the parameters a and q are related as follows:

$$\begin{array}{c}\hfill a=\left\{{\displaystyle \frac{4(\overline{m}+r)(\overline{m}+r{\overline{v}}_{\mathrm{A}}^2)}{r\overline{m}{M}_{\mathrm{A}}^2}-2}\right\}q.\end{array}$$(57)

In the incompressible limit, $\overline{m}\to 1$ and we find the same relation as Barbulescu et al. (2019) for their model. For fixed values of the background quantities and for a fixed value of $\overline{m}$, Eq. (57) represents a line in the aq-plane. We note, however, that changing the value of k_z and/or k_y a priori changes the value of $\overline{m}$ as well and that hence the modes do not lie all on the same line. This is unlike the case of an incompressible plasma, where one stays on the same line when changing k_z, as was already found for example by Barbulescu et al. (2019) in a similar model.

We find that the slope in Eq. (57) is minimized for $\overline{m}=r{\overline{\mathit{v}}}_{\mathrm{A}}$, in which case the line’s equation is

$$\begin{array}{c}\hfill a=\{\frac{4{(1+{\overline{v}}_{\mathrm{A}})}^2}{M_{\mathrm{A}}^2}-2\}q.\end{array}$$(58)

For each set of fixed background quantities, the only physically possible values of a and q thus lie in the region of the stability diagram between the positive a-axis and the line defined by Eq. (58). We see that, in this model, it is not physically possible to have values of (q, a) in the region below the line through the origin and with slope −2.

4.3.1. Coronal loops

A first application to a solar atmospheric structure would be a coronal loop, where slow waves occur near footpoints (De Moortel et al. 2002) or during flares (Wang et al. 2002; Kumar et al. 2013). For a standing slow wave under realistic coronal loop conditions, for which we took (v_Ae, v_si, v_se) = (5/2, 1/2, 1/4)v_Ai, the line of minimum slope, Eq. (58), lies almost on the positive a-axis. Indeed, with ${\overline{\mathit{v}}}_{\mathrm{A}}=5/2$ and for M_A = 0.066 for example, we find the line of minimum slope to be given by a = 11415q. If we assume incompressibility and fix $\overline{m}=1$, Eq. (57) defines a line as well and we find a = 12486q. For the less realistic value of M_A = 1, we find that the line of minimum slope is given by a = 47q whereas in the incompressible limit the line (57) becomes a = 52q. This is summarized in Table 1. As we can see from Fig. 2, the region of the aq-diagram relevant for coronal loop conditions is thus stable with respect to the KHI. It could be unstable with respect to resonance; however, since in reality a finite cylindrical structure is closed both azimuthally and longitudinally, the respective wavenumbers will take on discrete values. In the aq-diagram, the possible values will thus be a set of discrete points (a_n, q_n). Seeing how the vast majority of the area is a stable region in the relevant region between the positive a-axis and the minimum-slope line (58), it is likely that this parametric resonance instability is avoided as well.

We conclude that the oscillating shear velocity of a standing slow wave in a coronal loop would not trigger the KHI on its own, without the involvement of other physical processes. This is in clear contrast with the fast kink waves. Indeed, according to the incompressible models of Barbulescu et al. (2019) and Hillier et al. (2019), which are similar to this model but where the velocity and magnetic fields are perpendicular instead of parallel, fast kink waves are predicted to excite the KHI on the loop boundary. This has also been confirmed by numerical simulations.

4.3.2. Photospheric pores

As was mentioned in the introduction, propagating slow waves have also been observed in photospheric pores. One would be curious whether standing versions of these modes could trigger the KHI under these conditions. The model can be used to examine the stability of slow waves in these structures as well. However, the qualitative conclusion is the same as in the case of coronal loops: with photospheric pore conditions set to (v_Ae, v_si, v_se) = (1/4, 1/2, 3/4)v_Ai, we have ${\overline{\mathit{v}}}_{\mathrm{A}}=1/4$, and, taking M_A = 0.066, the minimum slope is 1450 whereas for the slope of the line defined by Eq. (57) is 1615 in the incompressible limit. As a comparison, even for the unrealistic value of M_A = 1, the minimum slope is 4.25 and the slope in Eq. (57) in the incompressible limit is 5. This is again summarized in Table 1. These values are an order of magnitude lower then for their coronal counterpart, but looking at Fig. 2 we see that the same conclusion remains: every possible pair of parameters (a, q) lie in a region that does not permit the development of the KHI. Similarly as for a coronal loop in the previous subsection, it is likely that the parametric instability arising from resonance with the driver is also avoided and that the pore therefore remains stable.

Additionally, one could investigate the stability of resonant slow waves, which can occur in the presence of a transition layer at the pore’s boundary. Indeed, slow surface waves are known to be resonantly absorbed in photospheric conditions. It is worth wondering whether the great longitudinal velocity shear that arises around the resonant point could be enough to trigger the KHI. In ideal MHD, the longitudinal velocity of resonant slow modes displays a hyperbolic singularity (Sakurai et al. 1991), but in the presence of finite resistivity this singularity disappears to make place for a sharp but continuous profile. Kovitya & Cram (1983) reported values of electrical resistivity in sunspots, with a minimum resistivity corresponding to a magnetic Reynolds number of 10⁷. Erdelyi (1997) derived an expression for the longitudinal Lagrangian displacement in the dissipative layer around the cusp resonance in the presence of finite electrical resistivity (see Eq. (30) therein). Assuming a sinusoidal transition profile with width l = 0.1R (where R is the radius of the pore) in the squared cusp and sound speeds, a ratio of external to internal magnetic field of B_0ze/B_0zi = 0.33, a cusp resonant position at r_C = 0.955R (based on the numerical computations of Goossens et al. 2021), a value for the perturbed total pressure at the resonant position equal to its value on the interface for the corresponding surface mode in the absence of the transition layer, a real part of the frequency equal to the frequency of the corresponding surface mode in the absence of the transition layer, and a magnetic Reynolds number of 10⁷ to have a lowerbound on resistivity and thus an upperbound on realistic values for the longitudinal velocity at the resonance, we found a value of M_A = 0.36 around the cusp resonance point for a sausage mode with k_zR = 2. This value of k_zR is within the range of validity for the longitudinal wavenumber of the observed slow surface sausage mode in a photospheric pore by Grant et al. (2015). The value found for M_A implies a minimum slope of 46 in Eq. (58) and a slope of 52 for Eq. (57) in the incompressible limit. These values are still well above the values for the slopes found for M_A = 1 in photospheric conditions. Through the model derived in this paper, we conclude from this that even resonant slow waves do not display a large enough velocity shear for an instability to develop.

4.4. Limitations of the model

We point out that the local model discussed in this paper is only valid under certain conditions. Firstly, the azimuthal wavenumber of the perturbation must be large, such that the azimuthal direction can be approximated by the cartesian direction of y. This means that we must have k_yR ≫ 1, with R the structure’s radius. Secondly, for the amplitude of the oscillating background velocity to be assumed constant, the longitudinal wavelength of the perturbation must be small with respect to the structure’s length L. We must thus have k_zL ≫ 1. Lastly, we note that our model might not be suitable in the presence of a smooth transition layer between the interior and the exterior of the structure. Indeed, resonant absorption of the standing slow wave in the background could occur in that case (Yu et al. 2017), which leads to a steep and continuous variation in the longitudinal component of the velocity along the x direction. Although one can use the model to estimate the effect of the logitudinal velocity shear in resonant absorption on the stability of the interface, as we did in last subsection, one must keep in mind that certain assumptions on which the model is based are not fulfilled in the presence of a transition layer. Firstly, there is no true discontinuous interface and, secondly, the assumption of a constant amplitude profile along x for the oscillating background velocity is violated.

5. Conclusion

In this paper, we developed an analytical model for the local stability of a cylindrical solar atmospheric structure harboring a standing slow wave. We assumed that the structure is a straight cylinder with a circular base, of which the boundary is an interface discontinuously separating two homogeneous and compressible plasmas. The magnetic fields on both sides were assumed to be constant, straight and aligned with the interface. The velocity fields were assumed to be oscillating in time but spatially constant, in order to model the standing slow wave in the background.

We used linearized MHD to derive an equation for the compression term, ∇ ⋅ ξ. The expressions for all the other perturbed quantities could then be expressed in terms of the solution of ∇ ⋅ ξ. The spatial part of the solution along the direction normal to the interface could be analytically found to be a normal mode, with wavenumber m. In contrast to the incompressible approximation, the value of this wavenumber m in a compressible plasma is different on both sides of the interface. We saw that m_e/m_i is dependent on both k_z and k_y, which entails that the stability of the interface depends on k_y as well as k_z. This is in contrast with the incompressible version of this model, in which only k_z has an influence on said stability.

We then found that the governing equation for the displacement component of the perturbation normal to the interface is a Mathieu equation. The stability of the solution of this equation in function of the different involved parameters being known from the litterature, we were able to describe two kinds of instabilities that could theoretically arise in the presence of an oscillating shear background flow. As an application to the solar atmosphere, we found that the physically relevant region of the parameter space corresponds to an interface that is locally stable with respect to the KHI, both in a coronal loop and in a photospheric pore. Although the interface can be unstable to a parametric resonance between the background slow wave and the induced perturbations on the interface, we concluded that it is unlikely to happen in reality. Even in the case of resonance in the cusp continuum, we concluded from our model that the longitudinal velocity shear of the resonant slow waves is not enough to trigger an instability.

We ended by noting that this model is only applicable under specific conditions. It is indeed a local model, in the sense that the spatial scale of azimuthal variations must be small compared to the structure’s radius and the spatial scale of longitudinal variations must be small with respect to the loop’s length. Furthermore, for a structure with a smooth transition layer the assumptions of a discontinuous transition at the boundary and a constant amplitude for the background oscillating velocity fields would be violated. This needs to be kept in mind when using this model to infer stability conditions around resonances.

Acknowledgments

M.G. was supported by the C1 Grant TRACEspace of Internal Funds KU Leuven (number C14/19/089). TVD was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 724326) and the C1 grant TRACEspace of Internal Funds KU Leuven.

References

Appendix A: Deriving the governing equation for ∇ ⋅ ξ

In this appendix, we derive Eq. (5) from Eqs. (1)–(4). Equations (1), (3) and (4) can be rewritten using $\mathit{v}=\frac{\mathrm{D}\mathit{\xi }}{\mathrm{D}t}$ to yield the following expressions for ρ₁, B₁ and p₁:

$$\begin{array}{cc}& {\rho }_1=-{\rho }_0\mathrm{\nabla }·\mathit{\xi }\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& {\mathit{B}}_1=-{\mathit{B}}_0(\mathrm{\nabla }·\mathit{\xi })+({\mathit{B}}_0·\mathrm{\nabla })\mathit{\xi }\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& p_1=-{\rho }_0{v}_{\mathrm{s}}^2\mathrm{\nabla }·\mathit{\xi }.\hfill \end{array}$$(A.3)

Using Eqs. (A.1), (A.3), and (2) can be rewritten as follows:

$$\begin{array}{cc}\hfill {\displaystyle {\omega }_0{V}_{0}sin\left({\omega }_{0}t\right)(\mathrm{\nabla }·\mathit{\xi }){\mathbf{1}}_z+\frac{{\mathrm{D}}^2\mathit{\xi }}{\mathrm{D}{t}^2}=}& \frac{-1}{{\rho }_0}\mathrm{\nabla }{P}_1\hfill \\ \hfill & +i{k}_z{v}_{\mathrm{A}}^2[-{\mathbf{1}}_z(\mathrm{\nabla }·\mathit{\xi })+i{k}_z\mathit{\xi }],\hfill \end{array}$$(A.4)

where ${\mathit{v}}_{\mathrm{A}}=B_{0z}/\sqrt{{\mathrm{\mu }}_0{\rho }_0}$ is the Alfvén speed. By now taking the divergence on both sides of Eq. (A.4) we get the following:

$$\begin{array}{c}\hfill {\displaystyle i{k}_z{\omega }_0{V}_{0}sin\left({\omega }_{0}t\right)(\mathrm{\nabla }·\mathit{\xi })+\frac{{\mathrm{D}}^2(\mathrm{\nabla }·\mathit{\xi })}{\mathrm{D}{t}^2}=\frac{-1}{{\rho }_0}{\mathrm{\nabla }}^2{P}_1,}\end{array}$$(A.5)

where ∇² denotes the Laplace operator.

Next, from the definition of total pressure P₁ = p₁ + B₀ ⋅ B₁/μ₀ we obtain the following equation with the use of Eqs. (A.2) and (A.3):

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{{\rho }_0}{P}_1+(v_{\mathrm{A}}^2+v_{\mathrm{s}}^2)(\mathrm{\nabla }·\mathit{\xi })-i{k}_z{v}_{\mathrm{A}}^2{\xi }_z=0.}\end{array}$$(A.6)

Taking the z-component of Eq. (A.4) and using Eq. (A.6) yields

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{D}}^2{\xi }_z}{\mathrm{D}{t}^2}=-{\omega }_0{V}_{0}sin\left({\omega }_{0}t\right)(\mathrm{\nabla }·\mathit{\xi })+i{k}_z{v}_{\mathrm{s}}^2(\mathrm{\nabla }·\mathit{\xi }),}\end{array}$$(A.7)

while taking twice the Lagrangian derivative on both sides of Eq. (A.6) and rearanging the terms yields

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{D}}^2{\xi }_z}{\mathrm{D}{t}^2}=\frac{1}{i{k}_z{v}_{\mathrm{A}}^2{\rho }_0}\frac{{\mathrm{D}}^2{P}_1}{\mathrm{D}{t}^2}+\frac{v_{\mathrm{A}}^2+v_{\mathrm{s}}^2}{i{k}_z{v}_{\mathrm{A}}^2}\frac{{\mathrm{D}}^2(\mathrm{\nabla }·\mathit{\xi })}{\mathrm{D}{t}^2}·}\end{array}$$(A.8)

We can then combine Eqs. (A.7) and (A.8) into the following equation:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{i{k}_z{v}_{\mathrm{A}}^2{\rho }_0}\frac{{\mathrm{D}}^2{P}_1}{\mathrm{D}{t}^2}+\frac{v_{\mathrm{A}}^2+v_{\mathrm{s}}^2}{i{k}_z{v}_{\mathrm{A}}^2}\frac{{\mathrm{D}}^2(\mathrm{\nabla }·\mathit{\xi })}{\mathrm{D}{t}^2}=}& -{\omega }_0{V}_{0}sin\left({\omega }_{0}t\right)(\mathrm{\nabla }·\mathit{\xi })\hfill \\ \hfill & +i{k}_z{v}_{\mathrm{s}}^2(\mathrm{\nabla }·\mathit{\xi }).\hfill \end{array}$$(A.9)