An Analytical Model of the Kelvin-Helmholtz Instability of Transverse Coronal Loop Oscillations

TL;DR

This paper develops the first analytical model of Kelvin-Helmholtz instability in transverse coronal loop oscillations, revealing that magnetic shear cannot prevent the instability and introducing the concept of loop σ-stability.

Contribution

It presents an analytical framework using Mathieu's equation to analyze KH instability in coronal loops, including effects of magnetic shear and loop twist.

Findings

The interface is always unstable to KH instability.

Magnetic shear can reduce but not eliminate the instability growth rate.

Weakly twisted loops can be σ-stable, with growth times longer than damping times.

Abstract

Recent numerical simulations have demonstrated that transverse coronal loop oscillations are susceptible to the Kelvin-Helmholtz (KH) instability due to the counter-streaming motions at the loop boundary. We present the first analytical model of this phenomenon. The region at the loop boundary where the shearing motions are greatest is treated as a straight interface separating time-periodic counter-streaming flows. In order to consider a twisted tube, the magnetic field at one side of the interface is inclined. We show that the evolution of the displacement at the interface is governed by Mathieu's equation and we use this equation to study the stability of the interface. We prove that the interface is always unstable, and that, under certain conditions, the magnetic shear may reduce the instability growth rate. The result, that the magnetic shear cannot stabilise the interface, explains the numerically found fact that the magnetic twist does not prevent the onset of the KH instability at the boundary of an oscillating magnetic tube. We also introduce the notion of the loop $\sigma$-stability. We say that a transversally oscillating loop is $\sigma$-stable if the KH instability growth time is larger than the damping time of the kink oscillation. We show that even relatively weakly twisted loops are $\sigma$-stable.

Equations

$${\mathbf{B}}_{\mathbf{i}}=(B_i\sin \theta ,0,B_i\cos \theta ),$$

$${\mathbf{B}}_{\mathbf{e}}=(0,0,B_e),$$

$${\mathbf{U}}_{\mathbf{i}}=(U\cos (\Omega t)\cos \theta ,0,-U\cos (\Omega t)\sin \theta ),$$

$${\mathbf{U}}_{\mathbf{e}}=(-U\cos (\Omega t),0,0),$$

$$\begin{array}{rl}\frac{\mathrm{D}\mathbf{v}}{\mathrm{D}t}& =-\frac{1}{{\rho }_{i,e}}\nabla p_T+\frac{1}{{\mu }_0{\rho }_{i,e}}({\mathbf{B}}_{\mathbf{i},\mathbf{e}}\cdot \nabla )\mathbf{b},\\ \frac{\mathrm{D}\mathbf{b}}{\mathrm{D}t}& =({\mathbf{B}}_{\mathbf{i},\mathbf{e}}\cdot \nabla )\mathbf{v},\\ \nabla \cdot \mathbf{v}& =0,\\ \nabla \cdot \mathbf{b}& =0,\end{array}$$

$$\frac{\mathrm{D}}{\mathrm{D}t}=\{\begin{array}{l}\begin{array}{r}\frac{\partial }{\partial t}+U\cos (\Omega t)\cos \theta \frac{\partial }{\partial x}\\ -U\cos (\Omega t)\sin \theta \frac{\partial }{\partial z},y<0,\end{array}\\ \frac{\partial }{\partial t}-U\cos (\Omega t)\frac{\partial }{\partial x},y>0.\end{array}$$

$$\frac{{\mathrm{D}}^2\mathbf{\xi }}{\mathrm{D}{t}^2}-\frac{1}{{\mu }_0{\rho }_{i,e}}({\mathbf{B}}_{\mathbf{i},\mathbf{e}}\cdot \nabla {)}^2\mathbf{\xi }=-\frac{1}{{\rho }_{i,e}}\nabla p_T.$$

$${\nabla }^2{p}_T=0.$$

$$\hat{p}_{T}(y)=p_{0}\left\{\begin{array}[]{cc}\mathrm{e}^{ky},&y<0,\\ \mathrm{e}^{-ky},&y>0,\end{array}\right.$$

$${\left(\frac{\partial }{\partial t}+i{k}_{x}U\cos (\Omega t)\cos \theta -i{k}_{z}U\cos (\Omega t)\sin \theta \right)}^2{\hat{\xi }}_y$$

$$+v_{Ai}^2{\left(k_x\sin \theta +k_z\cos \theta \right)}^2{\hat{\xi }}_y=-\frac{1}{{\rho }_i}\frac{\partial p_T}{\partial y},$$

$${\left(\frac{\partial }{\partial t}-i{k}_{x}U\cos (\Omega t)\right)}^2{\hat{\xi }}_y+v_{Ae}^2{k}_z^2{\hat{\xi }}_y=-\frac{1}{{\rho }_e}\frac{\partial p_T}{\partial y},$$

$$\displaystyle\begin{split}&\bigg{\{}\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}+2iA\cos(\Omega t)\frac{\mathrm{d}}{\mathrm{d}t}-i\Omega A\sin(\Omega t)\\ &\hphantom{xx}-B\cos^{2}(\Omega t)+C\bigg{\}}\hat{\xi}_{y}=0,\\ &A=\frac{U\big{[}\rho_{i}(k_{x}\cos\theta-k_{z}\sin\theta)-\rho_{e}k_{x}\big{]}}{\rho_{i}+\rho_{e}},\\ &B=\frac{U^{2}\big{[}\rho_{i}\left(k_{x}\cos\theta-k_{z}\sin\theta\right)^{2}+\rho_{e}k_{x}^{2}\big{]}}{\rho_{i}+\rho_{e}},\\ &C=\frac{\rho_{i}v_{Ai}^{2}\left(k_{x}\sin\theta+k_{z}\cos\theta\right)^{2}+\rho_{e}v_{Ae}^{2}k_{z}^{2}}{\rho_{i}+\rho_{e}},\end{split}$$

$$k_x=k\cos \varphi ,k_z=k\sin \varphi .$$

$$g(t)=\exp \left\{-\frac{iA}{\Omega }\sin (\Omega t)\right\},$$

$$\frac{{\mathrm{d}}^2\eta }{\mathrm{d}{\tau }^2}+[a-2q\cos (2\tau )]\eta =0,$$

$$\begin{array}{rl}q& =\frac{r{\kappa }^2{M}_A^2[\cos (\theta +\varphi )+\cos \varphi {]}^2}{4(1+r{)}^2},\\ \alpha & =\frac{{\kappa }^2[{\sin }^2(\theta +\varphi )+r{\bar{v}}_A^2{\sin }^2\varphi ]}{1+r},\\ a& =\alpha -2q,\end{array}$$

$$\begin{array}{rl}({\rho }_i+& {\rho }_e){\omega }^2-2Uk[{\rho }_i\cos (\theta +\varphi )-{\rho }_e\cos \varphi ]\omega \\ +& U^2{k}^2[{\rho }_i{\cos }^2(\theta +\varphi )+{\rho }_e{\cos }^2\varphi ]\\ -& {\rho }_i{v}_{Ai}^2{k}^2{\sin }^2(\theta +\varphi )-{\rho }_e{v}_{Ae}^2{k}^2{\sin }^2\varphi =0,\end{array}$$

$$M_{A0}^2=\frac{(1+r)[{\sin }^2(\theta +\varphi )+r{\bar{v}}_A^2{\sin }^2\varphi ]}{r[\cos (\theta +\varphi )+\cos \varphi {]}^2}.$$

$${\varphi }_0=-\arctan \left(\frac{\sin \theta }{\cos \theta +r{\bar{v}}_A^2}\right).$$

$$\min \{M_{A0}^2\}=\frac{{\bar{v}}_A^2(1+r){\tan }^2(\theta /2)}{1+r{\bar{v}}_A^2}.$$

$${\eta }_{+}(\tau )={\mathrm{e}}^{\mu \tau }P(a,q,\tau ),$$

$${\mathrm{e}}^{\mu \tau }=\exp (\Re (\mu )\Omega t)\exp (i\Im (\mu )\Omega t),$$

$$a=Kq,K=\frac{4M_{A0}^2}{M_A^2}-2.$$

$$\frac{{\mathrm{d}}^2\eta }{\mathrm{d}{\tilde{\tau }}^2}+[\tilde{a}-2\tilde{q}\cos (2\tilde{\tau }/\kappa )]\eta =0.$$

$$2\tilde{q}\cos (2\tilde{\tau }/\kappa )-\tilde{a}\ge 4h^2\tilde{q},$$

$$\frac{{\mathrm{d}}^2\eta }{\mathrm{d}{\tilde{\tau }}^2}-4h^2\tilde{q}\eta =0,$$

$$\eta_{1}=\eta_{0}\exp(2h\tilde{q}^{1/2}\tilde{\tau})=\eta_{0}\exp\big{(}\tau\sqrt{q(1-K/2)}\big{)},$$

$${\eta }_1={\eta }_0,\frac{\mathrm{d}{\eta }_1}{\mathrm{d}\tilde{\tau }}=2h{\eta }_0{\tilde{q}}^{1/2}\text{mbox}at\tilde{\tau }=0.$$

$${\eta }_2={\eta }_0,\frac{\mathrm{d}{\eta }_2}{\mathrm{d}\tau }=2h{\eta }_0{\kappa }^{-1}{\tilde{q}}^{1/2}\text{mbox}at\tilde{\tau }=0.$$

$$\bar{\eta }=1,\frac{\mathrm{d}\bar{\eta }}{\mathrm{d}\tau }=0\text{mbox}at\tau =0.$$

$$\cosh (\pi \mu )=\bar{\eta }(\pi ).$$

Full text

Sun: corona — Sun: oscillations — Sun: magnetic fields — instabilities — plasmas — magnetohydrodynamics (MHD)

1 Introduction

Transverse oscillations of coronal loops have been a subject of extensive study since their original observation on 14 July 1998 by the Transition Region and Coronal Explorer (TRACE) (Aschwanden et al., 1999; Nakariakov et al., 1999). For a review of the theory of these oscillations see Ruderman & Erdélyi (2009).

In particular, the damping mechanism of transverse loop oscillations has received much attention (e.g. Ruderman & Roberts, 2002; Goossens et al., 2002; Van Doorsselaere et al., 2004; Dymova & Ruderman, 2006; Williamson & Erdélyi, 2014), with the caveat that many studies have relied on the assumption that the oscillations are in the linear regime. The nonlinear damping of transverse coronal loop oscillations has also been studied, both analytically (Ruderman et al., 2010; Ruderman & Goossens, 2014; Ruderman, 2017), as well as numerically (e.g. Terradas & Ofman, 2004; Magyar & Van Doorsselaere, 2016a). The numerical studies revealed important effects, such as that of the ponderomotive force, and the presence of the Kelvin-Helmholtz instability (KHI) at the loop boundaries. More recently, Goddard & Nakariakov (2016) carried out a statistical study of observations of the damping of coronal loop kink oscillations.

Terradas et al. (2008) suggested that a kink oscillation may render a flux tube unstable due to the shear motions at the boundaries. The authors found that, for a smooth transition layer, the instability developed rapidly where the difference between the internal and external flow amplitudes was the greatest. However, increasing the thickness of the transitional layer significantly decreased the growth rate of the instability. It is worth noting that the KHI in smooth transition layers via other mechanisms (e.g. phase mixing, resonant absorption) had also received attention previously (see, for example, Heyvaerts & Priest, 1983; Ofman et al., 1994; Poedts et al., 1997). For a recent review on modelling the KHI see, e.g. Zhelyazkov (2015).

The topic of the transverse wave induced Kelvin-Helmholtz (TWIKH) instability was subsequently investigated by Antolin et al. (2014), who suggested that this phenomenon may be responsible for the fine strand-like structure observed in some coronal loops. In their numerical modelling these authors found that this structure is formed near the loop boundary even when the oscillation amplitude is very small, about 3 km/s. This result implies that the TWIKH instability develops even for very small oscillation amplitudes. The TWIKH instability has since been studied by Antolin et al. (2016); Magyar & Van Doorsselaere (2016a, b); Antolin et al. (2017); Karampelas et al. (2017); Howson et al. (2017a, b); Karampelas & Van Doorsselaere (2018), who considered various aspects of the instability onset, growth rate and observational properties.

The configuration of the equilibrium magnetic field is an important aspect of TWIKH instabilities. It was suggested by Terradas et al. (2008) that a twisted magnetic field may suppress the instability. The effect of twist on the stability of transverse loop transverse oscillations was studied numerically by Howson et al. (2017b) who investigated the energetics of the instability of a magnetically twisted coronal loop and found that its evolution is affected by the strength of the azimuthal component of the magnetic field. The authors also found that, when magnetic twist is present, the KHI leads to greater Ohmic dissipation as a result of the production of larger currents. Furthermore, Terradas et al. (2018) studied the evolution of the instability and found that the magnetic twist increases the instability growth time.

Numerical simulations have provided some insight into the development of the KHI, but have not thoroughly established what the conditions are needed for its onset. In this paper, we find these requirements analytically by modelling the boundary of the flux tube where the shearing is greatest as a single interface separating regions of different densities and magnetic fields, and performing a local stability analysis. We emulate the effect of the transverse oscillation by subjecting each region to temporally periodic counter-streaming flows.

Although this work is the first local analysis of the TWIKH instability with oscillating flows, the KHI in the presence of transverse shear and twisted magnetic fields has previously been studied by Soler et al. (2010) and Zaqarashvili et al. (2015). The aforementioned studies, however, consider steady flows in a cylindrical geometry, while this paper is concerned with the analysis of temporally periodic flows in a Cartesian geometry.

The paper is organised as follows: in Section 2, we introduce a Cartesian model of the boundary of a twisted flux tube, and derive the governing equation for the displacement. The stability of the flow is analysed in Section 3, followed by applications to transverse coronal loop oscillations in Section 4. Section 5 contains the summary of the obtained results and our conclusions.

2 The Governing Equation

It is well established that a magnetic flux tube undergoing transverse oscillation is prone to the Kelvin-Helmholtz instability due to the shearing motions at the boundaries (Terradas et al., 2008). Considering only the fundamental mode of oscillation, we wish to obtain the TWIKH instability criterion. We start by considering a magnetically twisted flux tube of length $L$. For mathematical simplicity, we consider the boundary of the tube to be a tangential interface, meaning there is no smooth boundary layer connecting the interior with the exterior. The amplitude of a fundamental transverse oscillation is greatest at the half-length of the tube, $L/2$, where the shearing is the greatest. We consider a plane $\Pi$ orthogonal to the tube axis and crossing it at its half-length. The intersection of this plane with the tube boundary is a circle. We also assume that the kink oscillation of the magnetic tube is linearly polarised and introduce the angle $\varphi$ in the plane $\Pi$, measured from the direction of the oscillation velocity in the counter-clockwise direction. Then, the shear velocity at the tube boundary takes its maximum at $\varphi=\pi/2$ and $\varphi=3\pi/2$, i.e. at the two points where it is parallel to the oscillation velocity (see Figure 1).

In order to study the effect of the shearing motions around this region, we model it as a single interface separating temporally periodic counter-streaming flows. We introduce the Cartesian coordinate system $x$, $y$, $z$ with the $x$-axis parallel to the direction of the polarisation of the kink oscillation, and the $z$-axis parallel to the tube axis. The interior and exterior of the tube are represented by the regions $y<0$ and $y>0$, respectively. The equilibrium quantities in these regions are denoted by the subscripts $i$ and $e$, respectively.

We assume that the equilibrium magnetic field is in the $xz$-plane. Since we wish to obtain the stability criteria both for straight and twisted tubes, we assume that the equilibrium magnetic field is parallel to the $z$-axis in the region $y>0$, and makes an angle $\theta$ with respect to the $z$-axis in the region $y<0$. Here, $\theta$ corresponds to the degree of twist (Figure 2a), which should be small since highly twisted magnetic flux tubes are prone to other types of instabilities, such as the kink instability, with which we are not concerned in the present study (e.g. Shafranov, 1958; Kruskal et al., 1958; Hood & Priest, 1979). In the case of a non-twisted flux tube, $\theta=0$.

In the present model, the background flows are similar to the velocity field at the boundary of a cylindrical flux tube undergoing a transverse oscillation. In transverse oscillations of coronal loops, the displacement of the flux tube boundary is almost perpendicular to the background magnetic field in the low-beta plasma approximation (see, e.g. Ruderman, 2007), therefore, we consider unperturbed magnetic fields and flow velocities of the form

[TABLE]

as illustrated in Figure 2b. Here, the period of the oscillatory flow, $2\pi/\Omega$, corresponds to the period of oscillation of the flux tube.

The kink oscillation of a coronal loop creates not only the oscillating velocity, but also the oscillating magnetic field orthogonal to the background field $\bf B$. However, in our model we carry out a local analysis of the stability of the region near the middle of the loop where the amplitude of oscillating velocity takes maximum. Since the oscillating magnetic field has a node at the middle of the loop, that is its amplitude is zero there, we do not take this oscillating magnetic field into account in our model.

It is worth noting that the problem of oscillatory counter-streaming flows has been previously studied by, e.g. Kelly (1965) and Roberts (1973). Our model is an improvement since we do not only consider parallel flows. Furthermore, our model differs from that of Roberts (1973) since we consider magnetic fields perpendicular to the flows on each side of the interface.

We study the dynamics of the outlined problem in the framework of linear ideal MHD. In the thin flux tube approximation, typically valid for transverse loop oscillations, the effects of compressibility are not significant. As such, we may use the approximation of incompressible plasma, which greatly simplifies the analysis. Thus, the set of governing equations is

[TABLE]

where $\mathbf{v},\mathbf{b}$ and $p_{T}$ are the perturbations of the velocity, magnetic field, and total pressure (magnetic plus plasma), $\rho_{i,e}$ are the background internal and external densities, and $\mu_{0}$ is the magnetic permeability of free space. $\mathrm{D}/\mathrm{D}t$ is the material derivative defined by

[TABLE]

Equation (1) must be supplemented with the conditions that $p_{T}$ and $\xi_{y}$ are continuous at $y=0$.

We now introduce the Lagrangian displacement $\bm{\xi}=\bm{\xi}(\mathbf{x},t)$, which is related to the velocity perturbation by $\mathbf{v}(\mathbf{x},t)=\mathrm{D}\bm{\xi}/\mathrm{D}t$. Combining the momentum and induction equations and substituting the expression for $\mathbf{v}$ in terms of the displacement yields

[TABLE]

Taking the divergence of this equation, and using $\nabla\cdot\bm{\xi}=0$, we obtain Laplace’s equation for the total pressure

[TABLE]

We Fourier-decompose all variables and write them in the form $f=\hat{f}\exp[i(k_{x}x+k_{z}z)]$. We immediately obtain that the solution to Equation (3) satisfying the condition that it is continuous at $y=0$ is

[TABLE]

where $p_{0}$ is an arbitrary constant, $\mathbf{k}=(k_{x},0,k_{z})$ is the wave vector, and $k=\sqrt{k_{x}^{2}+k_{z}^{2}}$.

The Fourier-decomposed $y$-component of Equation (2) reads

[TABLE]

for $y<0$, and

[TABLE]

for $y>0$. Here, $v_{Ai,e}^{2}=B_{i,e}^{2}/\mu_{0}\rho_{i,e}$ are the Alfvén speeds on either side of the interface. We substitute Equation (4) into Equations (2) and (6), take $y=0$, and eliminate the constant $p_{0}$ from the obtained equations. As a result, we arrive at the equation for the displacement of the boundary,

[TABLE]

where $\hat{\xi}_{y}$ is calculated at $y=0$.

It is now convenient to introduce the magnitude of the wave vector $k$, and the angle between the wave vector and the $x$-axis, $\phi$. We may, then, write

[TABLE]

Now, making the variable substitution $\hat{\xi}_{y}(t)=g(t)\eta(t)$, where

[TABLE]

we reduce Equation (LABEL:eq7) to

[TABLE]

where

[TABLE]

$\tau=\Omega t$, $r=\rho_{e}/\rho_{i}$ is the density ratio, $M_{A}=U/v_{Ai}$ is the Alfvén Mach number, $\bar{v}_{A}=v_{Ae}/v_{Ai}$ is the ratio of Alfvén speeds, and $\kappa=kv_{Ai}/\Omega$ is the dimensionless wavenumber.

It is important to note that, since $|g(t)|=1$, the variable substitution does not affect the stability analysis. Hence, unstable perturbations of the boundary correspond to unstable solutions of Equation (10). Equation (10) is known as Mathieu’s equation (McLachlan, 1946). It is interesting that Mathieu’s equation also arises in quite a different kind of MHD problem. Namely, it describes the amplification of MHD waves by periodic external forcing (e.g. Zaqarashvili, 2000; Zaqarashvili et al., 2002, 2005), and the Rayleigh-Taylor instability of a magnetic interface in the presence of oscillating gravity (Ruderman, 2018).

3 Investigation of Stability

In this section, we use Equation (10) to study the stability of the tangential discontinuity with an oscillating shear velocity. For comparison, we first briefly outline the well-known results related to the stability of a tangential discontinuity separating steady flows. To the best of our knowledge, these results were first obtained by Syrovatskii (1957) (see also Chandrasekhar, 1961).

3.1 Stability of Steady Flows

Before analysing the fully time dependent governing Equation (10), we return to Equation (LABEL:eq7) and set $\Omega=0$, in order to perform the analysis of the configuration in the presence of steady flows. Since the coefficients in Equation (LABEL:eq7) are now independent of $t$, we can look for the solution to this equation proportional to $\mathrm{e}^{-i\omega t}$ and obtain the dispersion equation

[TABLE]

where $\omega$ is the angular frequency of the perturbation.

We note that if the roots to Equation (12) are real, then $\hat{\xi}_{y}(t)$ is oscillatory and the system is neutrally stable. However, if complex conjugate roots exist, one of the roots has a positive imaginary part, meaning that $|\mathrm{e}^{-i\omega t}|\to\infty$ as $t\to\infty$, and the equilibrium configuration is unstable. Equation (12) has complex roots when its discriminant is negative, which occurs when $M_{A}>M_{A0}$, where

[TABLE]

The right-hand side of Equation (13) is singular for $\theta=(2n+1)\pi$ and $\theta+2\phi=(2n+1)\pi$, where $n$ is any integer number. The interface is stable for any value of $U$, for $\theta$ and $\phi$ satisfying either of the singularity conditions. We can see that for $\theta=(2n+1)\pi$, the velocity has the same magnitude and direction on both sides of the interface, meaning that there is no velocity jump across the interface. Hence, the equilibrium is static in the reference frame moving with the speed $U$ in the positive $x$-direction and, consequently, the presence of flow does not cause instability. In the second case, the interface is stable with respect to perturbations having wave vectors defined by $\phi=-\frac{1}{2}\theta+\left(n+\frac{1}{2}\right)\pi$. The projection of the velocity on these wave vectors is the same on both sides of the interface, that is, there is no jump in the velocity projection across the interface. Hence, these perturbations are stable for any value of $U$.

The Alfvén Mach number, $M_{A0}$, takes its minimum value with respect to $\phi$ at $\phi=\phi_{0}$, where

[TABLE]

Substituting Equation (14) into Equation (13), we obtain

[TABLE]

It follows that the system is stable for any value of $M_{A}$ below $\min\{M_{A0}\}$, while there are always unstable perturbations when $M_{A}>\min\{M_{A0}\}$. Equation (15) suggests there are no stable perturbations for $\theta=0$, and is in angreement with Syrovatskii (1957).

We note that the instability growth rate is proportional to $k$, which implies that the growth rate tends to infinity as $k\to\infty$. Since the growth rate is unbounded, we say that the initial value problem describing the evolution of the surface of discontinuity is ill-posed. This behaviour is further studied in Section 3.3.

3.2 Stability of Oscillating Flows

We now use Equation (10) to study the stability for arbitrary values of the equilibrium quantities. Floquet’s theorem states that Equation (10) has a solution of the form

[TABLE]

where $\mu=\mu(a,q)$ is the characteristic exponent, and $P(a,q,\tau)$ is a periodic function in $\tau$, with period $\pi$ (see, e.g., McLachlan, 1946; Abramowitz & Stegun, 1965). Since Equation (10) is invariant with respect to the substitution $-\tau\to\tau$ it follows that $\eta_{-}(\tau)=\mathrm{e}^{-\mu\tau}P(a,q,-\tau)$ is also a solution to this equation. Then, the general solution to Equation (10) is the linear combination of $\eta_{+}(\tau)$ and $\eta_{-}(\tau)$ unless $i\mu$ is an integer number.

The parameter $\mu$ determines the nature of solutions to Mathieu’s equation. We may always assume that $\Re(\mu)>0$, unless $\mu$ is purely imaginary, where $\Re$ indicates the real part of a quantity. Since we may write

[TABLE]

where $\Im$ indicates the imaginary part of a quantity, it follows that purely imaginary values of $\mu$ correspond to neutrally stable solutions, while real and complex values correspond to unstable solutions. Hence, $\Re(\mu)>0$ corresponds to an unstable perturbation. Unfortunately, $\mu$ cannot be easily computed analytically, and, for this reason, we perform a numerical analysis to gain further insight.

Following McLachlan (1946), we plot the stability diagram of Equation (10) in the $qa$-plane (Figure 3a). In accordance with the definition of $q$ in Equation (11), we only consider $q>0$. The white and hatched regions correspond to purely imaginary and real/complex values of $\mu$, respectively, and thus, to stable and unstable solutions to Equation (10). The contours bounding the regions are defined by the condition that $i\mu$ is an integer number, so that Equation (10) has either $\pi$ or $2\pi$-periodic solutions when the point $(q,a)$ is on one of these contours. These contours are called the characteristic curves, and are defined by the equations $a=a_{j}(q)$ and $a=b_{j}(q)$. These functions satisfy the inequalities $a_{j}<b_{j+1}<a_{j+1}$, where $j=0,1,2,\dots$ The curves $a_{j}(q)$ and $b_{j}(q)$ are shown by solid and dotted lines, respectively, in Figure 3a. The asymptotic behaviour of $a_{j}(q)$ and $b_{j+1}(q)$ for large $q$ is given by $a_{j}(q)\sim b_{j+1}(q)\sim-2q$ (Abramowitz & Stegun, 1965).

Complementary to the above, Figure 3b shows the values of the characteristic exponent $\mu$. Purely imaginary solutions are plotted in white, and are separated from real/complex solutions by the characteristic curves, while the real part of $\mu$ is plotted in contours in the unstable regions.

The coefficients in Equation (10) depend on six dimensionless parameters. Four of these parameters, $r$, $\theta$, $M_{A}$, and $\bar{v}_{A}$, are only dependent on the equilibrium quantities, while the other two, $\kappa$, and $\phi$, are related to particular perturbations, and are thus arbitrary. Hence, we must study the behaviour of solutions to Equation (10) for all possible values of these two parameters. It is also straightforward to see that $q$ and $a$ are invariant with respect to the substitution $\phi+\pi\to\phi$. This enables us to only consider values of $\phi$ in the interval $[-\pi/2,\pi/2]$.

We now wish to study the behaviour of solutions to Equation (10) for arbitrary $\kappa$. We begin by noting that, when $\phi$ is fixed and $\kappa$ varies from 0 to $\infty$ we obtain a straight line in the $qa$-plane. Using Equations (11), the equation of this line may be written as

[TABLE]

From Equations (13) and (16), we note that $K>-2$ for any $\theta\neq 0$ and any values of the other parameters. Considering the asymptotic behaviours of the characteristic curves, it follows that the line $a=Kq$ always intersects all curves $a=a_{j}(q)$ and $a=b_{j+1}(q)$, for $j=0,1,\dots$ Hence, there always exist some values of $\kappa$ and $\phi$ for which perturbations are unstable, regardless of the values of the other parameters. This implies that the tangential discontinuity separating oscillating flows is unstable for any value of $M_{A}$, which is qualitatively different from the discontinuity separating steady flows considered in Subsection 3.1. In the case of no magnetic shearing, when $\theta=0$, perturbations with $\phi=0$ and any $\kappa$ are unstable since the line $a=Kq$ will always be under the curve $a_{0}(q)$. This is illustrated by the red line in Figure 3a. The straight lines in Figure 3a are further discussed in Subsection 4.1.

3.3 The Initial Value Problem

We now consider the initial value problem for Equation (10). We fix $\phi$ and study how the properties of the initial value problem depend on $M_{A}$. First, we consider $M_{A}>M_{A0}(\phi)$, which, implies that $K<2$ due to Equation (16), and we prove that, in this case, the instability growth rate is unbounded. Let us introduce the scaled variables $\tilde{a}=\kappa^{-2}a$, $\tilde{q}=\kappa^{-2}q$, and $\tilde{\tau}=\kappa\tau$, and rewrite Equation (10) as

[TABLE]

It is important to note that $\tilde{a}$ and $\tilde{q}$ are independent of $\kappa$, and $\tilde{a}=K\tilde{q}$. We consider this equation on the interval $\tilde{\tau}\in[0,\tilde{\tau}_{0}]$, where $\tilde{\tau}_{0}=\kappa\arcsin h$ and $h=\frac{1}{2}\sqrt{1-K/2}$. Since $K>-2$, it follows that

[TABLE]

for $\tilde{\tau}\in[0,\tilde{\tau}_{0}]$.

We now consider equation

[TABLE]

and a solution to this equation

[TABLE]

where $\eta_{0}$ is an arbitrary constant. This solution satisfies the initial conditions

[TABLE]

We also consider a solution $\eta_{2}$ to Equation (17) satisfying the same initial conditions. Then, it follows from Equation (18) and the comparison theorem (e.g. Coddington & Levinson, 1955) that $\eta_{2}\geq\eta_{1}$ for $\tilde{\tau}\in[0,\tilde{\tau}_{0}]$. The initial condition for $\eta_{2}$ can be rewritten as

[TABLE]

This result implies that $\eta_{2}$ and $\mathrm{d}\eta_{2}/\mathrm{d}\tau$ are bounded at $\tau=0$ for $\kappa\in(0,\infty)$. Then, it follows from the inequality $\eta_{2}\geq\eta_{1}$ and Equation (20) that, for any $\tau_{0}\in(0,\arcsin h)$, there is such a solution to Equation (17) that it is bounded together with its first derivative at $\tau=0$ for any value of $\kappa$, but it is unbounded at $\tau=\tau_{0}$ as $\kappa\to\infty$. Hence, the instability growth rate is unbounded. This result implies that the initial value problem describing the evolution of the perturbed discontinuity is ill-posed when $M_{A}>\min\{M_{A0}\}$.

Now, we assume that $M_{A}<M_{A0}(\phi)$, so that, in accordance with Equation (16), $K>2$ and $a>2q$. We calculate the instability increment for $\kappa\gg 1$. Let $\bar{\eta}(\tau)$ be the solution to Equation (10), satisfying the initial conditions

[TABLE]

Then, the characteristic exponent is defined by the equation (Abramowitz & Stegun, 1965)

[TABLE]

We use the WKB method and look for a solution to Equation (10) in the form $\eta_{+}=\mathrm{e}^{\kappa\Theta}$. Substituting this expression into Equation (10) we obtain

[TABLE]

We impose the condition $\Theta=0$ at $\tau=0$. Then, we look for the solution to this equation in the form of expansion

[TABLE]

Substituting this expansion into Equation (25) and collecting terms of the order of unity we obtain

[TABLE]

The solution to this equation satisfying the condition $\Theta_{1}=0$ at $\tau=0$ is

[TABLE]

where we chose the plus sign at the square root.

In the next order approximation we collect terms of the order of $\kappa^{-1}$ in Equation (25) to obtain

[TABLE]

Using Equation (28) we find that the solution to this equation satisfying the condition $\Theta_{2}=0$ at $\tau=0$ is

[TABLE]

Recall that $\eta_{-}(\tau)=\eta_{+}(-\tau)$ is also a solution to Equation (10). Then, since $\Theta_{1}(\tau)$ is an odd function and $\Theta_{2}(\tau)$ is an even function, it follows that

[TABLE]

Introducing the notation $\chi=\Theta_{1}(\pi)$ and $\gamma=\Theta_{2}(\pi)$ we transform Equation (10) to

[TABLE]

When the absolute value of the right-hand side of this equation does not exceed unity the two values of $\mu$ satisfying this equation are purely imaginary and the corresponding wave mode is neutrally stable. When the absolute value of the right-hand side is larger than unity one of the two values of $\mu$ satisfying this equation has positive real part and the corresponding wave mode grows exponentially. However, we can observe that the right-hand side of Equation (32) is bounded for any $\kappa$. This implies that the real part of $\mu$ is also bounded, and the same is true for the growth rate. We made this conclusion for a particular value of $\phi$ and $M_{A}<M_{A0}(\phi)$. If we now assume that $M_{A}<\min\{M_{A0}\}$, then the growth rate of any wave mode is bounded. This means that the initial value problem describing the evolution of the discontinuity is well-posed when $M_{A}<\min\{M_{A0}\}$. From Equation (15) we see that this condition may be written in the approximate form as

[TABLE]

since, typically, $\theta\ll 1$.

4 Application to Transverse Coronal Loop Oscillations

The aim of this section is twofold. First, we further elaborate the analysis of Section 3 by considering the $\sigma$-stability of Equation (10). Afterwards, we apply some of the results obtained in Subsections 3.2 and 3.3 to the stability of coronal loop oscillations.

4.1 The $\sigma$-stability

We now use the concept of $\sigma$-stability, first introduced by Goedbloed & Sakanaka (1974) and Sakanaka & Goedbloed (1974). This concept is used in studies of thermonuclear plasma confinement where it is necessary that perturbation amplitudes remain sufficiently small on some relevant time scale. An equilibrium is $\sigma$-stable if the amplitudes of unstable perturbations grow at most like $\exp(\sigma t)$.

We apply the concept of $\sigma$-stability to the analysis of the KH instability induced by transverse oscillations of solar coronal loops. We say that a transverse coronal loop oscillation is $\sigma$-stable if the growth time of the KH instability exceeds the damping time due to resonant absorption. It is important to note that, in this paper, we only consider the KH instability due to the transverse oscillation of coronal loops without a transitional layer. If a transitional layer is present, the KH instability may still occur in coronal loops after the transverse oscillation is damped (Terradas et al., 2018) as a result of increased shearing motions due to resonant absorption (Heyvaerts & Priest, 1983; Browning & Priest, 1984).

Let $t_{D}=\alpha P$ be the damping time, where $P=2\pi/\Omega$ is the oscillation period, and $\alpha$ varies from 1 to 5 (see, e.g., Goddard & Nakariakov, 2016). It follows from our definition that $\sigma=1/\Omega t_{D}$, or

[TABLE]

When $\alpha$ varies from 1 to 5, $\sigma$ decreases from approximately 0.16 to 0.03. We see that, in any case, the interface cannot be $\sigma$-stable if the maximum growth rate exceeds 0.16, which implies that if the interface is $\sigma$-stable then the increment is much less than unity. It is shown in Appendix A that, in this case, the maximum growth rate for fixed $\phi$ is approximately equal to $1/2K$. Then, the maximum growth rate for all values of $\phi$ is $1/2K_{m}$, where $K_{m}=\min_{\phi}K$. Hence, the $\sigma$-stability condition reads

[TABLE]

To estimate $K_{m}$ we take as typical values $r=1/3$ and $\bar{v}_{A}^{2}=3$. Then, using Equations (15) and (35), and taking into account the fact that, typically, $\theta\ll 1$, we reduce the $\sigma$-stability criterion to

[TABLE]

The typical displacement of a kink-oscillating coronal loop is of the order of the loop radius. Then, the ratio of the velocity to $v_{Ai}$ is of the order of the loop radius and length. Hence, the typical value is $M_{A}=0.01$. It follows from Equation (36) that the interface is $\sigma$-stable if $\theta\gtrsim 1^{\circ}$ for $\alpha=1$, and $\sigma$-stable if $\theta\gtrsim 2^{\circ}$ for $\alpha=5$. Similar to Terradas et al. (2018) we define the number of turns of a magnetic field as

[TABLE]

where $B_{\phi}$ and $B_{z}$ are the azimuthal and axial components of the magnetic field in cylindrical coordinates with the $z$-axis coinciding with the loop axis, and $R$ is the radius of the loop cross-section. Now we use the relation $B_{\phi}/B_{z}=\theta$ valid for small $\theta$ and $R/L=100$ as a typical value for coronal loops. We obtain that even the maximum value $\theta=2^{\circ}$ corresponds to only about a half-turn of magnetic field lines from one loop footpoint to the other. Hence, the loop boundary is $\sigma$-stable even for a very moderate magnetic twist.

In Figure 4, we present the values of $\mu$ associated with the three straight lines in Figure 3. We assumed that $r=1/3$, $\bar{v}_{A}^{2}=3$, $M_{A}=0.01$, and $\phi=\phi_{0}$ so that $K=K_{m}$. For $\theta=0$, $\mu$ is a monotonically increasing function of $\kappa$, and perturbations with any $q$ are unstable. The green curve corresponds to $\theta=0.5^{\circ}$, and is unbounded as $\kappa\to\infty$ since $\min M_{A0}\approx 0.0062<M_{A}$. Finally, the blue curve, which corresponds to $\theta=1^{\circ}$, is bounded for $\kappa\in(0,\infty)$ since $\min M_{A0}\approx 0.0123>M_{A}$. The equation of the dashed line is $\mu=0.16$, and we see that the loop with $\theta=1^{\circ}$ is $\sigma$-stable for $\sigma$ defined in Equation (34) with $\alpha=1$.

We note that if a magnetic loop is $\sigma$-stable, then the initial value problem describing the evolution of its boundary perturbation is well-posed. However, the converse is not always true. The initial value problem is well-posed if the growth rate is bounded, but it may still be very large. On the other hand, a magnetic loop is $\sigma$-stable when the maximum growth rate is below a definite and, usually, sufficiently small number.

4.2 The $\sigma$-stability in Numerical Models

We compare our results with those of Howson et al. (2017b) and Terradas et al. (2018), who studied numerical models of the TWIKH instability in twisted magnetic flux tubes. Both models consider flux tubes with a finite-width transitional layer, where the density decreases from a high value in the core region of the flux tube to a low value in the surrounding plasma. The presence of the transitional layer results in damping of kink oscillations due to resonant absorption, such that the concept of $\sigma$-stability is applicable. Since we do not consider the effects of resonant absorption in the present work, we may only make a qualitative comparison between results.

Howson et al. (2017b) considered both twisted and untwisted tubes, subject to a transverse oscillation with a period of the fundamental mode of 280 s. Both the oscillation period and damping time were practically unaffected by the magnetic twist. Using the dependence of the oscillation amplitude on time presented in Howson et al. (2017b), we estimate that the damping time of the transverse oscillation was approximately 1000 s. We also estimate that the instability growth time increases from approximately 600 s in the case of the untwisted tube to approximately 700 s when the twist is maximal, which signals a relatively weak dependence of damping time on the degree of twist. The increase in growth time with increase in twist qualitatively agrees with the results obtained in the present work.

We have shown in the previous subsections that, in a tube with a sharp boundary (i.e. no transitional layer), the instability growth time is zero. Therefore, it is clear that the presence of a transitional layer strongly reduces the instability growth rate, and, in the model studied by Howson et al. (2017b), the effect of the transitional layer on the instability increment is stronger than the effect of twist. Since the damping time was larger than the instability growth time, the oscillations studied by Howson et al. (2017b) were $\sigma$-unstable for all values of twist.

Terradas et al. (2018) also studied kink oscillations of twisted tubes with transitional layer of thickness $l$. They considered three values of the transitional layer thickness, $l/R=0.3,\,1$, and 2, where $R$ is the tube radius. They also considered several values of the magnetic twist, with the turn of magnetic lines varying from 0 (no twist) to 1.65 turns.

Similarly to Howson et al. (2017b), Terradas et al. (2018) obtained that the damping time is practically independent of the twist. It was approximately equal to $4P$ for $l/R=0.3$, where $P$ is the oscillation period. They did not give the value of damping time for other values of the transitional layer thickness. However, since Terradas et al. (2018) obtained that the numerically calculated values of damping time agree very well with those given by the analytical expression, we can use the fact that the damping time is inversely proportional to $l/R$. We obtain the estimates that the damping time is about $1.2P$ for $l/R=1$ and $0.6P$ for $l/R=2$. Even if we underestimated the damping time, then the first time is definitely less than $2P$, and the second one is less than $P$.

The authors also estimated the instability growth time. They obtained that it strongly depends on the degree of twist. For $l/R=0.3$ it increases from about $1.5P$ to about $3P$ when the turn of magnetic field lines varies from 0 to 1.65. Hence, it is always smaller than the damping time meaning that the oscillations are $\sigma$-unstable. When $l/R=1$, the instability growth time increases from about $2.5P$ to about $7.5P$. Finally, when $l/R=2$ the instability growth time is about $5P$ when there is no twist, and quickly becomes larger than $10P$ when the twist increases. Hence, the oscillations are always $\sigma$-stable when $l/R=1$ and $l/R=2$. Since they are $\sigma$-stable even when there is no twist, it is obvious that there is a substantial contribution of the transitional layer in the reduction of the instability increment. However, it is also obvious that the twist substantially contributes in this reduction.

4.3 Coronal Loop Parameters

The model that we outlined in the previous sections can be only applied for the local analysis of the stability of the boundary of an oscillating magnetic tube. In this analysis, we can consider oscillations with the characteristic scale in the azimuthal direction that is much smaller than the tube radius $R$, and the characteristic scale in the axial direction that is much smaller than the tube length $L$. Hence, we take

[TABLE]

where $m$ and $n$ are sufficiently large integer numbers. Using Equations (8) and (37) we obtain

[TABLE]

We assume that $n\lesssim|m|$. Since in coronal magnetic loops $R\ll L$, it follows that we may use the approximate expressions

[TABLE]

Throughout this section we assume that $\bar{v}_{A}^{2}=r^{-1}$. This assumption holds if the magnitudes of the interior and exterior magnetic fields are equal, which is typically true for coronal loops. We also assume that $\theta\ll 1$. Then, we obtain the approximate expressions

[TABLE]

[TABLE]

The condition $M_{A}^{2}<\min\{M_{A0}^{2}\}$ gives

[TABLE]

If we take $r=1/3$, the right-hand side of this inequality is approximately equal to $M_{A}$, that is it is of the order of 0.01. Hence, the inequality (42) can be satisfied even for quite moderated twist. If the inequality is satisfied, then the IVP describing the evolution of the tube boundary is well-posed and the growth rate of perturbations is bounded.

In Figures 5, we show the dependence of the growth rate on $m$ for $n=1$ (left) and $n=4$ (right), $M_{A}=0.01$, $r=1/3$, $\bar{v}_{A}^{2}=3$, $R/L=200$, and $\theta=0^{\circ}$ (red), $\theta=0.5^{\circ}$ (green) and $\theta=1^{\circ}$ (blue). We note that, obviously, $n=1$ does not satisfy the condition that $n$ is large, so we considered $n=1$ only for comparison. While, for $n=1$, the points in the $qa$-plane corresponding to $\theta=0^{\circ}$ are virtually unchanged as compared to the line in Figure 3, for $n=4$ they are shifted upwards considerably. This is also the case for $\theta=0.5^{\circ}$. We see that for $n=1$ there are some modes which are unstable in the range selected, for $n=4$ there are no such modes. There may be unstable modes for $\theta=0.5^{\circ}$ and $n=4$, but only for very large $m$. In terms of the IVP, for $\theta=1^{\circ}$, corresponding to a well-posed solution, no value of $m$ corresponds to an unstable solution in the $qa$-plane. In general, well-posed solutions seem to be unstable only for very large $m$. These results are significant since they suggest that very localised longitudinal perturbations of the flux tube are generally more stable.

5 Summary and Discussion

In this work, we performed the first local stability analysis of the transverse wave induced Kelvin-Helmholtz instability of twisted solar coronal loops. We modelled the region on the loop boundary where the shear flows are the greatest as a tangential discontinuity separating time-periodic counter-streaming flows. To model the magnetic twist in coronal loops we assumed that the equilibrium magnetic fields on either side of the discontinuity are not parallel. The flow velocities at the two sides of the discontinuity have opposite directions and equal magnitudes oscillating harmonically. For the sake of mathematical simplicity, we assumed that the plasma on both sides of the interface is incompressible. Using the linearised set of ideal MHD equations, we derived the governing equation describing the evolution of the shape of the tangential discontinuity, known as Mathieu’s equation.

We employed Mathieu’s equation to study the stability of the discontinuity. For comparison, we first presented the results of the stability analysis in the case of steady flows, which we obtained by setting the flow oscillation frequency to zero. In this case, the stability of the discontinuity is determined by the Alfvén Mach number, which is defined as the ratio of the background velocity magnitude to the Alfvén speed at one side of the interface. The discontinuity is unstable when the Alfvén Mach number exceeds a critical value, and the instability growth rate is proportional to the wavenumber, and thus unbounded. This implies that the initial value problem describing the evolution of the perturbed discontinuity is ill-posed. We note that the critical Alfvén number is zero when there is no magnetic shear.

In contrast to the interface separating steady flows, the tilted magnetic field cannot stabilise the discontinuity if the flows oscillate. A similar result was obtained by Roberts (1973) in the case of MHD tangential discontinuity with the magnetic field having the same direction at both sides and the flow velocity parallel to the magnetic field.

Even though the interface is always unstable, the critical Alfvén Mach number still plays an important role in the stability properties. We showed that the growth rate of the instability is unbounded when the Alfvén Mach number exceeds the instability threshold, and thus the initial value problem is ill-posed. Hence, in this case the stability properties are qualitatively the same as in the case of steady flows. On the other hand, when the the Alfvén Mach number is below its critical value, the instability increment is bounded, and the initial value problem is well-posed.

In Section 4.1, we introduced the definition of $\sigma$-stability for kink oscillating coronal loops, which states that the loop is $\sigma$-stable if the growth time of the instability exceeds the resonant damping time of the transverse oscillation. We obtained the criterion for the $\sigma$-stability and showed that, for parameters typical for transverse coronal loop oscillations, even moderate magnetic twist makes the loop boundary $\sigma$-stable.

In Section 4.3, we used our model to perform a local stability analysis of the sections of the loop boundary where the amplitudes of the shear flows are the greatest (see Figures 1 and 2). The local analysis is only valid for perturbations with the azimuthal wavelength much smaller than the radius of the loop cross-section $R$, and the axial wavelength much smaller than the loop length $L$. In accordance with these latter assumptions, we took $k_{x}=m/R$ and $k_{z}=\pi n/L$, where $k_{x}$ is the component of the wave vector in the azimuthal direction, and $k_{z}$ is the component of the wave vector in the axial direction, and $|m|$ and $n$ are positive integer numbers. We note that, while $n$ is positive, $m$ can be either positive or negative. We found that the nature of solutions is changed by this new definition of the parameters. While, previously, all solutions were unstable regardless of the background parameters, the discretisation of the parameter space has introduced the possibility that unstable solutions exist only for sufficiently large values of $|m|$.

It is worth noting that our study does not include the effects of strong shear induced by resonant absorption, which may be significant in the generation of the KHI, as suggested by Antolin et al. (2014). The numerical studies by Howson et al. (2017b) and Terradas et al. (2018) showed that the presence of the transitional layer leads to an increase in the instability growth time. This suggests that the main driver of the KH instability is the shear motion at the magnetic tube boundary due to the transverse oscillation, as opposed to the shearing caused by resonant absorption.

Our model may be expanded such that more accurate quantitative results about transverse loop oscillations are obtained. A transitional layer, where the oscillating velocity continuously varies from one side to the other, may be included. A further extension may consider a continuous variation of density from one side to the other, such that the effects of resonant absorption are also considered. Both of these generalizations are likely to be mathematically complicated.

A different possible application of the present model relates to prominence oscillations (Arregui et al., 2012). Assuming that the magnetic field has the same magnitude inside and outside the structure, for a typical density contrast of $r=100$, Equation (16) yields that $K\ll 1$, for $\theta\neq 0$. This suggests that, unless the magnetic fields inside and outside the prominence are perfectly aligned, the growth time of perturbations is very small.

Finally, we make the following comment. Usually it is written in papers dealing with the numerical study of the KH instability of oscillating coronal magnetic loops that this instability occurs in the nonlinear regime. However, in our paper the background state is given by the linear solution describing the kink oscillation. The stability analysis is also based on the use of the linear MHD. This clearly shows that the KH instability of oscillating coronal loops is not related to the nonlinearity at all.

Acknowledgments: MB, MSR and RE are grateful to the Science and Technology Facilities Council (STFC, grant number ST/M000826/1) UK and the Royal Society (UK) for the support received. TVD was supported by GOA-2015-014 (KU Leuven). This work was based on discussions at the ISSI (Bern, Switzerland, March 2017). This project has also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 724326).

Appendix A The Maximum Growth Rate

As we have already stated before, the characteristic exponent, $\mu$, is determined by the equation

[TABLE]

where $\bar{\eta}(\tau)$ is the solution to the initial value problem to Equation (10) with

[TABLE]

(Abramowitz & Stegun, 1965). When a perturbation is unstable, its growth rate is given by $\gamma=\Re(\mu)$. In the context of the $\sigma$-stability analysis, we assume that the growth time of the instability is much larger than the oscillation period. In terms of dimensionless quantities, this condition is written as $\gamma\ll 1$. The numerical investigation shows that this condition is only satisfied for all values of $q$ when $K\gg 1$. In accordance with this, we introduce the small parameter $\epsilon=1/K$. Figure 3 shows that $a$ is close to $j^{2}$ on parts of the line $a=Kq$ corresponding to unstable perturbations when $K\gg 1$, where $j=1,2,\dots$ We obtain $a=j^{2}$ taking $q=j^{2}\epsilon$, which implies that $q={\cal O}(\epsilon)$.

First we study the case with $j=1$. Using the expansion valid for small $q$ (Abramowitz & Stegun, 1965),

[TABLE]

we obtain that the line $a=Kq$ in Figure 3 intersects the curves $a=b_{1}(q)$ and $a=a_{1}(q)$ at $q\approx\epsilon-{\cal O}(\epsilon^{2})$ and $q\approx\epsilon+{\cal O}(\epsilon^{2})$, respectively. Then $q=\epsilon+\bar{q}\epsilon^{2}$ on the part of the curve $a=Kq$ between the intersection points, where $\bar{q}$ is a free parameter varying from approximately $-1$ to approximately 1. It follows that $q=\epsilon+\bar{q}\epsilon^{2}$ on the line $a=Kq$ between the intersection points, where $\bar{q}$ is a free parameter. The equation of the curve $a=Kq$ is now rewritten as $a=1+\bar{q}\epsilon$, and Equation (10) becomes

[TABLE]

To calculate the increment we need to find the solution $\bar{\eta}(\tau)$ to this equation satisfying the initial conditions Equation (A2). To do this we use the regular perturbation method with

[TABLE]

Substituting Equation (A5) into Equations (A2) and (A4), and collecting the terms of the order of unity, we obtain

[TABLE]

and the associated initial conditions

[TABLE]

The solution to this initial value problem is

[TABLE]

Collecting term of the order of $\epsilon$ yields

[TABLE]

[TABLE]

After straightforward calculation we obtain

[TABLE]

Finally we collect terms of the order of $\epsilon^{2}$ to obtain

[TABLE]

[TABLE]

The solution to this initial value problem is given by

[TABLE]

Using Equations (A8), (A11), and (A14) we obtain

[TABLE]

It follows from this equation that

[TABLE]

This result implies that

[TABLE]

where $\gamma_{m}$ is the maximum value of the instability increment when the point $(a,q)$ is on the part of line $a=Kq$ that is between the curves $a=b_{1}(q)$ and $a=a_{1}(q)$.

Now we consider the part of line $a=Kq$ that is between the curves $a=b_{j}(q)$ and $a=a_{j}(q)$, $j=2,3,\dots$ For $q\ll 1$ we have $b_{1}(q)=n^{2}+{\cal O}(q^{2})$ and $a_{1}(q)=n^{2}+{\cal O}(q^{2})$, where $n$ is a natural number (Abramowitz & Stegun, 1965). Since $K=\epsilon^{-1}$, it follows that $q=n^{2}\epsilon(1+\bar{q}\epsilon^{2})$ and $a=n^{2}(1+\bar{q}\epsilon^{2})$, where $\bar{q}$ is again a free parameter. Substituting these expressions in Eq. (A1) we transform it to

[TABLE]

Then we again look for the solution in the form of the expansion given by Eq. (A5). Substituting this expansion in Equations (10) and (A2), and collecting terms of the order of unity we obtain

[TABLE]

[TABLE]

The solution to this initial value problem is

[TABLE]

Collecting terms of the order of $\epsilon$ yields

[TABLE]

[TABLE]

After straightforward calculation we obtain

[TABLE]

for $j=2$, and

[TABLE]

for $j>2$. Collecting terms of the order of $\epsilon^{2}$ we obtain

[TABLE]

[TABLE]

The solution to this initial value problem is given by

[TABLE]

for $j=2$, and by

[TABLE]

for $j>2$. Using Eqs. (A21), (A24), (A25), (A28), and (A29), we obtain

[TABLE]

It follows from this equation that $\mu={\cal O}(\epsilon^{3})$ for even $j$ and $\mu=i+{\cal O}(\epsilon^{3})$ for odd $j$, and thus $\gamma={\cal O}(\epsilon^{3/2})$, that is $\gamma\ll\gamma_{m}$. Hence, $\gamma_{m}=1/2K$ is the maximum value of the instability increment with respect to $q$ when $K=\epsilon^{-1}$.