1. Introduction
Vortices are frequently observed in the atmosphere of the earth and other planets such as Jupiter and Saturn. On the earth they sometimes appear as intense vortices such as tornadoes and cyclones, causing disasters. On Jupiter long-lived vortices like the Great Red Spot are one of the long-standing mysteries studied by a number of researchers. Sets of multiple vortices, such as vortex pairs or arrays, are sometimes formed. For example, an array of counter-rotating vortices like a von Kármán vortex street is often observed in the wake of an isolated island (Potylitsin & Peltier Reference Potylitsin and Peltier1998). On Jupiter anticyclones and cyclones formed a von Kármán vortex street for about 50 years (Youssef & Marcus Reference Youssef and Marcus2003). On Saturn numerical simulations suggest that a vortex street is responsible for the stable polygonal shape of a jet observed in the polar region (Morales-Juberías et al. Reference Morales-Juberías, Sayanagi, Dowling and Ingersoll2011); Saturn's great storm was observed to create multiple vortices (Sayanagi et al. Reference Sayanagi, Dyudina, Ewald, Fischer, Ingersoll, Kurth, Muro, Porco and West2013). These arrays of vortices can be generated by instabilities of a jet flow and a shear flow (the Kelvin–Helmholtz instability) and by other mechanisms.
One of the important features of arrays of vortices is that they possess hyperbolic stagnation points in general. When there are hyperbolic stagnation points, the flow is subject to hyperbolic instability (Friedlander & Vishik Reference Friedlander and Vishik1991; Lifschitz & Hameiri Reference Lifschitz and Hameiri1991; Sipp & Jacquin Reference Sipp and Jacquin1998; Pralits, Giannetti & Brandt Reference Pralits, Giannetti and Brandt2013); in this paper, we call this the pure hyperbolic instability in order to differentiate it from the strato-hyperbolic instability introduced below. Our recent work showed that a new type of instability called strato-hyperbolic instability exists in stably stratified vortices having hyperbolic stagnation points (Suzuki, Hirota & Hattori Reference Suzuki, Hirota and Hattori2018). Since the arrays of vortices in the atmosphere are affected by density stratification and planetary rotation, the role of the strato-hyperbolic instability as well as the pure hyperbolic instability is of much interest in understanding and predicting their motion. These instabilities can destabilize the arrays of vortices in the atmosphere; they can survive with a certain level of turbulence or break down completely depending on whether the nonlinear evolution which follows the linear instability saturates or not.
There are also other instabilities that can arise in a set of vortices or an array of vortices in stratified and/or rotating fluids: the elliptic instability (Miyazaki & Fukumoto Reference Miyazaki and Fukumoto1992; Leblanc & Cambon Reference Leblanc and Cambon1998; Leweke & Williamson Reference Leweke and Williamson1998; Miyazaki & Adachi Reference Miyazaki and Adachi1998; Guimbard et al. Reference Guimbard, Le Dizès, Le Bars, Le Gal and Leblanc2000; Otheguy, Billant & Chomaz Reference Otheguy, Billant and Chomaz2006a; Aspden & Vanneste Reference Aspden and Vanneste2009), the centrifugal instability (Leblanc & Cambon Reference Leblanc and Cambon1998; Potylitsin & Peltier Reference Potylitsin and Peltier1998, Reference Potylitsin and Peltier1999), the zigzag instability (Billant Reference Billant2000; Billant & Chomaz Reference Billant and Chomaz2000a,Reference Billant and Chomazb,Reference Billant and Chomazc; Otheguy, Billant & Chomaz Reference Otheguy, Billant and Chomaz2006b; Deloncle, Billant & Chomaz Reference Deloncle, Billant and Chomaz2008; Waite & Smolarkiewicz Reference Waite and Smolarkiewicz2008; Billant et al. Reference Billant, Deloncle, Chomaz and Otheguy2010), the radiative instability (Le Dizès & Billant Reference Le Dizès and Billant2009) and the transient growth (Arratia, Caulfield & Chomaz Reference Arratia, Caulfield and Chomaz2013; Gau & Hattori Reference Gau and Hattori2014). Leblanc & Cambon (Reference Leblanc and Cambon1998) investigated the effects of rotation on the linear stability of the Stuart vortices in non-stratified fluids by numerical analysis; the centrifugal, elliptic and pure hyperbolic instabilities were found, although the results are rather limited owing to low resolution. Leblanc & Godeferd (Reference Leblanc and Godeferd1999) showed the structures of a mode of the pure hyperbolic instability in the two-dimensional (2-D) Taylor–Green vortices by direct numerical simulation (DNS); the Reynolds number was rather low at 
$Re=400$. Potylitsin & Peltier (Reference Potylitsin and Peltier1998) studied the effects of density stratification on the stability of periodic vortices on the 
$f$-plane by numerical analysis; the base flow is a quasi-steady state obtained by relaxation at low Reynolds numbers. According to Potylitsin & Peltier (Reference Potylitsin and Peltier1998), anticyclonic vortices are strongly destabilized by weak rotation but stabilized by strong rotation; they also claimed that strong stratification stabilizes the vortices. These results were obtained from numerical simulations with limited resolution at low Reynolds numbers. Potylitsin & Peltier (Reference Potylitsin and Peltier1999) investigated the stability of the Stuart vortices in rotating non-stratified fluids by numerical analysis. They found three types of instability: the elliptic, the centrifugal and the (pure) hyperbolic instabilities.
Although several important aspects of the instabilities of arrays of vortices in stratified and/or rotating fluids have been elucidated, our understanding is far from complete. In particular, the vertical wavenumber is often bounded from below due to geometric constraints, as the vertical scale is much smaller than the horizontal vortices; for example, a large typhoon or cyclone extends about 1000 km horizontally, while the vertical scale is about 10 km at most, which is 
$1/100$ times the horizontal scale. This implies that the stability properties at high vertical wavenumbers (when based on the horizontal scale) are important. Strong stratification also makes the characteristic length scale in the vertical direction small (Billant & Chomaz Reference Billant and Chomaz2001). In addition, results obtained with low numerical resolution or at low Reynolds numbers may not apply to actual large-scale vortices at high Reynolds numbers. Thus, stability properties of arrays of vortices in stratified and/or rotating fluids should be further explored with higher resolution.
Suzuki et al. (Reference Suzuki, Hirota and Hattori2018) studied the stability of stably stratified vortices. Three base flows which possess hyperbolic stagnation points were considered: the 2-D Taylor–Green vortices, the Stuart vortices and the Lamb–Chaplygin dipole. In addition to the elliptic and the pure hyperbolic instabilities, the strato-hyperbolic instability, which is caused by hyperbolic instability near the hyperbolic stagnation points and phase shift by the internal gravity waves, was found. Which instability is dominant depends on the base flow and the stratification. However, these results are mostly based on the local stability analysis. Although the local analysis gives important and useful information at high wavenumbers, global stability analysis should be performed not only to confirm the local results but also to reveal mode structures and to clarify how the dominant instability changes depending on the base flow and the parameters.
In this paper, we study the stability of arrays of stably stratified vortices by global stability analysis. Our aim is to reveal the detailed stability properties of arrays of vortices in stratified fluids. The effects of rotation are not included in the present paper; the study of these is our next step after clarifying the effects of stratification. The characteristics of the strato-hyperbolic instability as well as the pure hyperbolic instability are also investigated in detail.
It is useful at this point to summarize the physical mechanism of the elliptic, the pure hyperbolic and the strato-hyperbolic instabilities in interpreting the results presented in this paper. It is vortex stretching due to straining flow that essentially causes all three instabilities. The difference is the mechanism which maintains the stretching effects, causing them to give rise to instabilities. In the elliptic instability the vorticity of disturbances changes its direction periodically as they rotate around an elliptic stagnation point; instability occurs when the frequency of the direction change is tuned to favour the stretching direction and avoid the compression direction, so that the stretching effects are dominant (Waleffe Reference Waleffe1990). In the pure hyperbolic instability, disturbances which are advected to a hyperbolic stagnation point are stretched for a long time (infinitely long in the short-wave limit) (Friedlander & Vishik Reference Friedlander and Vishik1991; Lifschitz & Hameiri Reference Lifschitz and Hameiri1991); an instability arises if this growth is not cancelled out by other mechanisms. In the strato-hyperbolic instability the internal gravity waves rotate the vorticity of the disturbances to change its angle or, in other words, to shift the phase; this phase shift causes the stretching effects near the hyperbolic stagnation points to be maintained when the phase shift satisfies a resonance condition (Suzuki et al. Reference Suzuki, Hirota and Hattori2018). These three instabilities are short-wave instabilities, whose growth rates converge to finite values in the short-wave and inviscid limits. This is contrasted to the zigzag instability, which is essentially a long-wave instability; it resembles the Crow instability of a vortex pair in non-stratified flows (Billant Reference Billant2000), although the wavelength can be short for strong stratification. The three instabilities above should be distinguished from the secondary instabilities of the Kelvin–Helmholtz billows (Caulfield & Peltier Reference Caulfield and Peltier2000), in which the directions of gravity and vorticity of the base flow are perpendicular, while they are parallel in the present paper.
This paper is organized as follows. In § 2 details of the stability analysis are explained. Two base flows are chosen: the 2-D Taylor–Green vortices and the Stuart vortices. In §§ 3 and 4, we present the results on the 2-D Taylor–Green vortices and the Stuart vortices, respectively. We conclude in § 5.
2. Methods
2.1. Set-up of the problem
We consider the linear stability of a stably stratified flow under the Boussinesq approximation. Two base flows are considered: the 2-D Taylor–Green vortices (figure 1) and the Stuart vortices (figure 2). The 2-D Taylor–Green vortices are an array of vortices doubly periodic in horizontal directions. Each vortex sits in a rectangular cell whose vertices are hyperbolic points. The sign of vorticity in each cell is opposite to that in the neighbouring cells; this forms a staggered lattice of vortices. On the other hand, the Stuart vortices are a one-dimensional array of periodic vortices. The sign of vorticity is the same for all vortices. There is a hyperbolic point between each pair of adjacent vortices. The vorticity is parallel to the vertical direction. Both the 2-D Taylor–Green vortices and the Stuart vortices are steady in the absence of diffusion.
Figure 1. Streamlines of 2-D Taylor–Green vortices; 
$\varepsilon _e=0$. The solid (red) and dashed (blue) lines correspond to positive and negative values of the stream function, respectively.
Figure 2. Stuart vortices; 
$C=1.2$. (a) Streamlines: closed streamlines (red), separatrix (green), open streamlines (blue). (b) Contours of vorticity distribution.
Viscosity is taken into account, while diffusion of density is neglected except in § 3.4, where the effects of diffusion are checked. We consider high-Reynolds-number flows throughout the paper. The base flow is assumed steady, since the growth of instabilities is much faster than the time evolution of the base flow due to viscous diffusion. The velocity, pressure and density fields are decomposed as
\begin{gather} \boldsymbol{u}=\boldsymbol{u}_b+\boldsymbol{u}', \end{gather}
\begin{gather}p=p_b+p', \end{gather}
\begin{gather}\rho=\rho_{0}+{\alpha}z+\rho', \end{gather}
where 
$(\boldsymbol {u}_b, p_b, \rho _b=\rho _0 +\alpha z)$ and 
$(\boldsymbol {u}', p', \rho ')=(u', v', w', p', \rho ')$ are the base flow and the disturbance, the direction of the gravity force is taken as 
$-\boldsymbol {e}_z$ and 
$\alpha ={\partial \rho _b}/\partial {z}<0$ is a constant. The magnitude of the disturbance is infinitesimally small. The governing equations in non-dimensionalized form read
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'=0, \end{gather}
\begin{gather}\frac{\partial{\boldsymbol{u}'}}{\partial t}+(\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u}_b+(\boldsymbol{u}_b\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}'= -\boldsymbol{\nabla}{p'}-\rho'\boldsymbol{e}_z+\frac{1}{Re}\nabla^{2} \boldsymbol{u}', \end{gather}
\begin{gather}\frac{\partial{\rho'}}{\partial{t}}+(\boldsymbol{u}_b\boldsymbol{\cdot}\boldsymbol{\nabla}){\rho'}- \frac{1}{F_h^{2}}w' =0, \end{gather}
where 
$Re =U_0L_0/\nu$ is the Reynolds number, 
$F_h=U_0/(L_0N)$ is the Froude number, 
$N=\sqrt {-{\alpha }g/{\rho }_0}$ is the Brunt–Väisälä frequency, 
$g$ is the acceleration of gravity, 
$\nu$ is the kinematic viscosity, and 
$U_0$ and 
$L_0$ are a characteristic velocity and a length scale. In the following the values are scaled by 
$U_0$ and 
$L_0$ unless explicitly stated otherwise.
2.2. Numerical procedure
The equations (2.4)–(2.6) were solved numerically by the Fourier spectral method (Peyret Reference Peyret2010) assuming periodic boundary conditions in all three directions. This is valid for the 2-D Taylor–Green vortices, while some care must be taken for the Stuart vortices, which are periodic in 
$x$ but not in 
$y$. In the following we place the Stuart vortices at 
$y=nL_y$ and those with the opposite-signed vorticity at 
$y=(n+1/2)L_y$, where 
$n$ is an integer, so that the base flow is periodic in 
$y$. The spatial period 
$L_y$ should be large enough to minimize the effects of the periodic boundary condition, which will be checked in § 4.1. The time marching was performed by the fourth-order Runge–Kutta method.
Since the base flow is 2-D, the time evolution of disturbances is separable in the vertical direction. Thus, we set
\begin{equation} \boldsymbol{u}' = \exp({{\textrm{{i}}}k_z z}) \sum_{k_x=-K_1}^{K_1}\sum_{k_y=-K_2}^{K_2} \tilde{\boldsymbol{u}}_{k_x,k_y} \exp({\textrm{i} (k_x x+k_y y)})\end{equation}
with similar expressions for 
$p'$ and 
$\rho '$.
For the 2-D Taylor–Green vortices, the base flow consists of only four Fourier modes: 
$\exp ({2{\rm \pi} {\textrm {{i}}}(\pm x/L_x \pm y/L_y)})$. Thus the products of the base flow and the disturbances in (2.5) and (2.6) were taken directly in the Fourier space. For the Stuart vortices, on the other hand, the products were calculated in the physical space to reduce computational time since all Fourier modes are non-zero in general.
2.3. Parameters of Krylov method
The growth rate and frequency were obtained by the Krylov subspace method (Edwards et al. Reference Edwards, Tuckerman, Friesner and Sorensen1994; Julien, Ortiz & Chomaz Reference Julien, Ortiz and Chomaz2004; Donnadieu et al. Reference Donnadieu, Ortiz, Chomaz and Billant2009). Starting from randomized initial conditions, (2.4)–(2.6) were integrated for a certain long time. Intermediate states 
$\{(\pmb {u}'(T_0), \rho '(T_0)), (\pmb {u}'(T_0+{\rm \Delta} T), \rho '(T_0+{\rm \Delta} T)), \ldots , (\pmb {u}'(T_0+(N-1){\rm \Delta} T), \rho '(T_0+(N-1){\rm \Delta} T))\}$ were used as generators of the Krylov subspace. Then the eigenvalues and the eigenmodes were obtained in the 
$N$-dimensional Krylov subspace.
In this method, the error of an eigenvalue 
$\lambda$ of a linear operator 
${\mathcal {L}}$ can be evaluated by
\begin{equation} \delta = \frac{\| {\mathcal{L}}\boldsymbol{v}-\lambda\boldsymbol{v} \|}{\| \boldsymbol{v} \|}, \end{equation}
where 
$\boldsymbol {v}$ is the corresponding approximate eigenvector. The error 
$\delta$ depends on the initial time of the data 
$T_0$, the interval between the data 
${\rm \Delta} T$ and the dimension of the Krylov subspace 
$N$. In order to obtain eigenvalues accurately, several Krylov subspaces were generated from different sets of parameters and the eigenvalue with the smallest error for each eigenmode was chosen. The actual values of the parameters were chosen after trial and error. For the 2-D Taylor–Green vortices, the number of data 
$N$ was fixed to 
$10$, the start time of the data was 
$T_0=195$ or 
$T_0=245$, and the interval between the data 
${\rm \Delta} T$ was fixed to 
$5$. For the Stuart vortices, 
$N$ was either 
$10$ or 
$20$ and 
$T_0=92, 112, \ldots , 192$, while 
${\rm \Delta} T$ was fixed to 
$2$. Typically the error of the eigenvalue is 
$\delta =O(10^{-10})$ for the largest eigenvalue for a fixed wavenumber 
$k_z$, while it increases for subdominant eigenmodes. In the following, we discarded the eigenmodes with 
$\delta \geqslant 10^{-3}$.
3. Stability of stratified 2-D Taylor–Green vortices
3.1. Base flow and parameter values
In this section we study the stability of the stratified 2-D Taylor–Green vortices (figure 1). The vorticity is
\begin{equation} \omega(x,y) = \omega_0 \sin \frac{x}{L_x} \sin \frac{y}{L_y}, \end{equation}
where 
$\omega _0$ is the maximum vorticity and 
$L_x$ and 
$L_y$ are spatial periods in the 
$x$ and 
$y$ directions, respectively. The stream function is 
$\varPsi = B \sin ({x}/{L_x}) \sin ({y}/{L_y})$, where 
$\omega _0 = (L_x^{-2}+L_y^{-2})B$. The hyperbolic points are 
$(x,y)=(mL_x/2, nL_y/2)$, where 
$m$ and 
$n$ are integers, while 
$(x,y)=((m/2+1/4)L_x, (n/2+1/4)L_y)$ are elliptic stagnation points. The strain rate at the elliptic stagnation points normalized by 
$\omega _0$ is
\begin{equation} \varepsilon_e=\frac{|{L_x}^{2}-{L_y}^{2}|}{2({L_x}^{2}+{L_y}^{2})}. \end{equation}
We choose 
$L_0=(L_xL_y)^{1/2}$ and 
$U_0=\omega _0L_0/(2{\rm \pi} )$ as the characteristic length and velocity, respectively, so that the time scale is 
$L_0/U_0=2{\rm \pi} /\omega _0$.
The Froude number is set to 
$F_h^{-1}=0, 5$ and 
$10$. The Reynolds number is fixed to 
$\textit {Re}=10^{5}$, except that 
$\textit {Re}=10^{4}$ is also considered for 
$F_h^{-1}=5$ to see the effect of viscosity on the growth rate. The strain rate at the elliptic stagnation points is set to 
$\varepsilon _e=0$ and 
$\varepsilon _e=0.2$. The range of the wavenumber 
$k_z$ is 
$0 < k_z < 60$. The number of Fourier modes is 
$500\times 500$.
Since the velocity field of the base flow is mirror-symmetric with respect to 
$x=mL_x/2$ and 
$y=nL_y/2$, the unstable modes can be classified into four different symmetry types, namely 
$S_xS_y$, 
$A_xS_y$, 
$S_xA_y$ and 
$A_xA_y$, where the mode is of type 
$S_x$ when it is symmetric with respect to 
$x=0$, i.e.
\begin{equation} u'(-x,y,z)=-u'(x,y,z), \quad v'(-x,y,z)=v'(x,y,z), \quad w'(-x,y,z)=w'(x,y,z), \end{equation}
while it is of type 
$A_x$ when it is antisymmetric with respect to 
$x=0$, i.e.
\begin{equation} u'(-x,y,z)=u'(x,y,z), \quad v'(-x,y,z)=-v'(x,y,z), \quad w'(-x,y,z)=-w'(x,y,z). \end{equation}
The types 
$S_y$ and 
$A_y$ are defined in the same way by exchanging 
$\{x, u'\}$ and 
$\{y, v'\}$. It is pointed out that these symmetry types are preserved under the linearized equations (2.4)–(2.6).
3.2. Types of instability and their mode structures
The unstable modes are classified into four types: (i) pure hyperbolic instability (PH), (ii) strato-hyperbolic instability (SH), (iii) mixed hyperbolic instability (MH) and (iv) elliptic instability (E). Examples of each instability type are shown below.
Mode structures of the pure hyperbolic instability are shown in figure 3. This mode exists for both the non-stratified and stratified cases, although the symmetry type is restricted to 
$A_xA_y$. For 
$F_h^{-1}=0$ and 
$k_zL_0=1.3$ strong vorticity is concentrated near the cell boundaries. The region of strong vorticity is wider for 
$F_h^{-1}=5$ and 
$k_zL_0=25.1$; the branch of the pure hyperbolic instability mode is further discussed in § 3.3.
Figure 3. Structures of pure hyperbolic instability modes. 2-D Taylor–Green vortices with 
$\varepsilon _e=0$. (a,b) Iso-surface of the magnitude of vorticity coloured by 
$\omega^{\prime}_z$; (c,d) contours of 
$\omega^{\prime}_z$ on the 
$xy$ plane; (e,f) contours of 
$\omega '_y$ on 
$y/L_0=0.25$ near the cell boundaries. The contours of 
$\omega '_i$ (
$i=y, z$) are drawn for 
$|\omega '_i|/|\omega '_i|_{{max}} = 0.1, 0.3, 0.5, 0.7$ and 
$0.9$; the red and blue lines correspond to positive and negative values, respectively. (a,c,e) 
$F_h^{-1}=0, k_zL_0=1.3$, 
$\sigma =0.445$, 
$\omega =1.860$, 
$A_xA_y$. (b,d,f) 
$F_h^{-1}=5, k_zL_0=25.1$, 
$\sigma =0.378$, 
$\omega =1.351$, 
$A_xA_y$. The vertical domain size of the iso-surface and the magnitude of vorticity on the iso-surface are (a) 
$L_z=(2{\rm \pi} /k_z)L_0, |\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.5$ and (b) 
$4L_z=(8{\rm \pi} /k_z)L_0, |\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.4$. The colour bar of 
$\omega '_z$ shown in (a,c,e) is common to all iso-surface figures. The contours of 
$\omega '_z$ are drawn on 
$z/L_z=0.29$ in (a,c,e) and 
$z/L_z=0$ in (b,d,f).
Mode structures of the strato-hyperbolic instability are shown in figure 4. This mode exists only for the stratified case; all symmetry types are found as confirmed later. Figure 4(a) shows a stationary mode of 
$S_xS_y$ symmetry with 
$F_h^{-1}=5$ and 
$k_zL_0=28.9$. Tubular structures of strong vorticity are observed near the cell boundaries, although the 
$z$ component of vorticity 
$\omega '_z$ vanishes at the boundaries. Cellular structures symmetric in the 
$z$ direction are observed on the 
$xz$ plane. Figure 4(b) shows an oscillatory mode of 
$A_xS_y$ symmetry with 
$F_h^{-1}=5$ and 
$k=41.5$. The vorticity distribution on the 
$xy$ plane is similar to figure 4(a), while the sizes of the tubular structures are smaller. The 
$z$ component of vorticity 
$\omega '_z$ vanishes at 
$y=nL_y/2$ because of 
$S_y$ symmetry; however, it does not vanish at 
$x=nL_x/2$. Vorticity distribution on the 
$xz$ plane is asymmetric in the 
$z$ direction; the complex conjugate counterpart of this mode has the reverse structure.
Figure 4. Structures of strato-hyperbolic instability modes. 2-D Taylor–Green vortices with 
$\varepsilon _e=0$, same as in figure 3: (a,c,e) 
$F_h^{-1}=5, k_zL_0=28.9$, 
$\sigma =0.608$, 
$\omega =0$, 
$S_xS_y$; (b,d,f) 
$F_h^{-1}=5, k_zL_0=41.5$, 
$\sigma =0.396$, 
$\omega =1.607$, 
$A_xS_y$. The vertical domain size of the iso-surface is 
$4L_z=(8{\rm \pi} /k_z)L_0$. The magnitude of vorticity on the iso-surface is 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.25$ in (a) and 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.2$ in (b). The contours of 
$\omega '_z$ are drawn on 
$z/L_z=0$.
Mode structures of the mixed hyperbolic instability are shown in figure 5. This mode appears mostly for the stratified case, although there exist a few modes for the non-stratified case; symmetry types other than 
$S_xS_y$ are found. The vorticity distribution is a mixture of the pure hyperbolic instability and the strato-hyperbolic instability. The mode with 
$F_h^{-1}=5$ and 
$k_zL_0=10.1$ shown in figure 5(a) is similar to the pure hyperbolic instability mode shown in figure 3(b) along 
$x=L_x/2$, while it is similar to the strato-hyperbolic instability along 
$y=L_y/2$. The mode with 
$F_h^{-1}=5$ and 
$k_zL_0=31.4$ shown in figure 5(b) has strong vorticity not only in regions near the cell boundaries but also at the cell boundaries. The difference between the strato-hyperbolic and the mixed hyperbolic instabilities will be clarified later in this section (figure 7) and in §§ 3.3 and 3.6.
Figure 5. Structures of mixed hyperbolic instability modes. 2-D Taylor–Green vortices with 
$\varepsilon _e=0$, same as in figure 3: (a,c,e) 
$F_h^{-1}=5, k_zL_0=10.1$, 
$\sigma =0.687$, 
$\omega =0$, 
$A_xS_y$; (b,d,f) 
$F_h^{-1}=5, k_zL_0=31.4$, 
$\sigma =0.458$, 
$\omega =2.628$, 
$A_xA_y$. The vertical domain size of the iso-surface and the magnitude of vorticity on the iso-surface are (a) 
$2L_z=(4{\rm \pi} /k_z)L_0, |\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.2$ and (b) 
$4L_z=(8{\rm \pi} /k_z)L_0, |\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.25$. The contours of 
$\omega '_z$ are drawn on 
$z/L_z=0$.
Mode structures of the elliptic instability are shown in figure 6. All the elliptic instability modes are found for 
$\varepsilon _e=0.2$. This mode is stabilized by strong stratification; no elliptic instability modes are found for 
$F_h^{-1}=5$ or 
$F_h^{-1}=10$. All symmetry types are found for 
$\varepsilon _e=0.2$ as confirmed later. The mode with 
$F_h^{-1}=0$ and 
$k_zL_0=6.1$ (figure 6a) is a typical one found in previous studies (Waleffe Reference Waleffe1990; Sipp & Jacquin Reference Sipp and Jacquin1998). Strong vorticity is concentrated near the elliptic stagnation points at higher wavenumbers (figure 6b).
Figure 6. Structures of elliptic instability modes. 2-D Taylor–Green vortices with 
$\varepsilon _e=0.2$, same as in figure 3 except that the contours of 
$\omega '_y$ near the cell boundaries are not included: (a,c) 
$F_h^{-1}=0, k_zL_0=6.1$, 
$\sigma =0.618$, 
$\omega =0$, 
$S_xA_y$; (b,d) 
$F_h^{-1}=0, k_zL_0=23.4$, 
$\sigma =0.514$, 
$\omega =2.794$, 
$S_xS_y$. The vertical domain size of the iso-surface and the magnitude of vorticity on the iso-surface are (a) 
$2L_z=(4{\rm \pi} /k_z)L_0, |\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.25$ and (b) 
$4L_z=(8{\rm \pi} /k_z)L_0, |\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.25$. The contours of 
$\omega '_z$ are drawn on 
$z/L_z=0$.
The classification of the hyperbolic instability modes into the pure hyperbolic, mixed hyperbolic and strato-hyperbolic instabilities is shown in figure 7. We define two quantities which characterize the structures of the unstable modes as
\begin{equation} s_1 = \frac{\int_{|\varPsi|\leqslant 0.1B} {\omega'}_z^{2} \,\textrm{d}x\,\textrm{d}y\,\textrm{d}z}{\int {\omega'}_z^{2} \,\textrm{d}x\,\textrm{d}y\,\textrm{d}z}, \quad s_2 = \frac{\int \left(\boldsymbol{\nabla}_h \boldsymbol{\cdot} \boldsymbol{u}'_h\right)^{2} \textrm{d}x\,\textrm{d}y\,\textrm{d}z}{\int {\omega'}_z^{2} \,\textrm{d}x\, \textrm{d}y\,\textrm{d}z}. \end{equation}
Here, 
$s_1$ measures the ratio of 
${\omega '}_z^{2}$ in the region near the cell boundaries 
$|\varPsi |/B \leqslant c_b$, where 
$B$ is the maximum of the stream function of the base flow, while 
$s_2$ is the ratio of the horizontal divergence, defined by 
$\boldsymbol {\nabla }_h \boldsymbol {\cdot } \boldsymbol {u}'_h={\partial u'}/{\partial x}+{\partial v'}/{\partial y}$, to the vertical component of vorticity; the constant 
$c_b$ was set to 
$0.1$ for 
$F_h^{-1}=5$. Figure 7 shows that the strato-hyperbolic instability modes are clustered in 
$s_1 \lesssim 0.03$ as 
$\omega '_z$ is small near the cell boundaries. In contrast 
$s_1$ is larger than 
$0.2$ for the pure hyperbolic instability modes and for most of the mixed hyperbolic instability modes. The horizontal divergence is small as 
$s_2 < 0.06$ for the pure hyperbolic instability modes. The mixed hyperbolic instability modes are separated from the other two instabilities, although some of them are close to the strato-hyperbolic instability modes and the pure hyperbolic instability modes, suggesting that the corresponding instability has larger amplitude than the other in its mode structure. The difference between the three hyperbolic instabilities will be elucidated in § 3.3.
Figure 7. Classification of the type of instability of unstable modes. 2-D Taylor–Green vortices with 
$\varepsilon _e=0$, 
$F_h^{-1}=5$. Blue: pure hyperbolic instability. Red: strato-hyperbolic instability. Green: mixed hyperbolic instability. See the text for the definitions of 
$s_1$ and 
$s_2$.
3.3. Branches of unstable modes and the growth rate: 
$\varepsilon _e=0$
Figure 8 shows the growth rate 
$\sigma$ and frequency 
$\omega$, which are the real and imaginary parts of the eigenvalue 
$\lambda =\sigma +{\textrm {{i}}} \omega$, plotted against the vertical wavenumber 
$k_z$ for 
$\varepsilon _e=0$. For the non-stratified case (figure 8a), there is only one pure hyperbolic instability branch, which is oscillatory and takes maximum at 
$k_zL_0=0$, and one branch which is stationary and broad, extending from 
$k_zL_0=7.54$ to 
$k_zL_0 \geqslant 60$; this branch is categorized as the mixed hyperbolic instability, as explained below. The branch of the pure hyperbolic instability survives without significant change of the growth rate. It does not change very much for 
$F_h^{-1}=1$ and 
$F_h^{-1}=2$ (figure 8b,c). However, it extends to high wavenumber for 
$F_h^{-1}=5$ and 
$F_h^{-1}=10$ (figure 8d,e). The maximum growth rate shifts to 
$k_zL_0=12.6$ for 
$F_h^{-1}=10$.
Figure 8. Growth rate (a–e) and frequency (f–j). 2-D Taylor–Green vortices, 
$\textit {Re}=10^{5}, \varepsilon _e=0$; 
$F_h^{-1}$ equals (a,f) 
$0$, (b,g) 
$1$, (c,h) 
$2$, (d,i) 
$5$ and (e,j) 
$10$. Symbols: 
$\square$ (blue) denotes pure hyperbolic instability; 
$\circ$ (red) denotes strato-hyperbolic instability; 
$\triangle$ (green) denotes mixed hyperbolic instability. Solid symbols and open symbols stand for stationary (
$\omega =0$) and oscillatory (
$\omega \neq 0$) modes, respectively. Only the oscillatory modes with positive frequencies are shown.
The two branches found for 
$F_h^{-1}=0$ are tracked by changing the strength of stratification. Figure 9(a) shows the pure hyperbolic instability branch for 
$F_h^{-1}=2, 4, 5$; also included in this figure is a mixed hyperbolic instability for 
$F_h^{-1}=4$ which merges with the pure hyperbolic instability at 
$F_h^{-1}=5$. As a result the pure hyperbolic instability branch extends to high wavenumber for 
$F_h^{-1}=5$ and 
$F_h^{-1}=10$, as we observed in figures 8(d) and 8(e). Figure 9(b) tracks the other branch of 
$F_h^{-1}=0$ in the range 
$0 \leqslant F_h^{-1} \leqslant 5$. This branch is tentatively called a ‘spiral’ mode since the vorticity distribution exhibits spiral structures (figure 10a), but it is later categorized as the mixed hyperbolic instability for the reason described below. The growth rate of the ‘spiral’ mode increases monotonically in 
$0 \leqslant F_h^{-1} \leqslant 5$ except at low wavenumbers 
$k_0L_z \lesssim 6$. For the non-stratified case, the mode possess spiral structures as shown in figures 10(a) and 10(b). The mode at 
$k_zL_0=15.1$ has maximum vorticity component 
$\omega '_z$ near the elliptic stagnation points, where the structures resemble those of the elliptic instability. The mode at 
$k_zL_0=28.9$ exhibits similar structures, although 
$\omega '_z$ is large near the cell boundary 
$x/L_x=0.5$ and is small near the elliptic stagnation points. This change is also observed at 
$k_zL_0=15.1$ as stratification becomes strong (figure 10c,d); for 
$F_h^{-1}=2$ spiral structures are weak and the mode shows the features of the pure hyperbolic instability. Since this mode is connected smoothly to the strongest mixed hyperbolic instability at 
$F_h^{-1}=5$, we categorize this branch for 
$0 \leqslant F_h^{-1} < 5$ as the mixed hyperbolic instability.
Figure 9. (a) Pure hyperbolic instability mode and (b) ‘spiral’ mode (mixed hyperbolic instability mode). 2-D Taylor–Green vortices, 
$\varepsilon _e=0.$ Panel (a) shows the pure hyperbolic instability mode for 
$F_h^{-1}=2, 4, 5$ and a mixed hyperbolic instability mode for 
$F_h^{-1}=4$ which merges with the pure hyperbolic instability mode. The symbols are the same as in figure 8. In panel (b), the branch other than the pure hyperbolic instability at 
$F_h^{-1}=0$ is tracked for 
$F_h^{-1}=0, 0.5, 1, 1.5, 2, 3$ and 
$5$. The growth rate for 
$20 \lesssim k_zL_0 \lesssim 40$ increases with 
$F_h^{-1}$.
Figure 10. Structures of ‘spiral’ modes (mixed hyperbolic instability mode). 2-D Taylor–Green vortices with 
$\varepsilon _e=0$, 
$A_xS_y$: (a) 
$F_h^{-1}=0, k_zL_0=15.1$, 
$\sigma =0.327$, 
$\omega =0$; (b) 
$F_h^{-1}=0, k_zL_0=28.9$, 
$\sigma =0.355$, 
$\omega =0$; (c) 
$F_h^{-1}=1, k_zL_0=15.1$, 
$\sigma =0.422$, 
$\omega =0$; (d) 
$F_h^{-1}=2, k_zL_0=15.1$, 
$\sigma =0.571$, 
$\omega =0$.
One more mixed hyperbolic instability appears at 
$F_h^{-1}=1$, while little change is observed for the pure hyperbolic instability branch. At 
$F_h^{-1}=2$ another mixed hyperbolic instability appears; the strato-hyperbolic instability modes with small growth rates also appear. The growth rate of the strato-hyperbolic instability mode increases and the number of the branches also increases at 
$F_h^{-1}=5$ and 
$F_h^{-1}=10$ (figure 8d,e). The lowest wavenumbers of the branches of the strato-hyperbolic instability are 
$k_zL_0 \approx 10 \sim 20$ for 
$F_h^{-1}=5$, while they shift to 
$k_zL_0 \approx 22 \sim 40$ for 
$F_h^{-1}=10$. This shift is in accordance with the self-similarity of strongly stratified flows (Billant & Chomaz Reference Billant and Chomaz2001), which states that the vertical scale is 
$U_0/N=F_hL_0$. In fact, the resonance condition for the strato-hyperbolic instability obtained by local stability analysis (Suzuki et al. Reference Suzuki, Hirota and Hattori2018) is
\begin{equation} \frac{1}{F_h} \int_{mT/2}^{(m+1)T/2} \sin \theta = 2n{\rm \pi}, \end{equation}
where 
$T$ is the period of the fluid particle motion, 
$\theta$ is the angle between the 
$z$ axis and the wavevector, and 
$m$ and 
$n$ are integers; according to (3.6), when 
$F_h^{-1}$ increases, 
$\theta$ decreases, which implies that 
$k_z$ increases for a fixed horizontal wavenumber. The growth rate increases with 
$k_z$, although viscous damping decreases the growth rate slightly for 
$k_zL_0 \gtrsim 50$ for 
$F_h^{-1}=5$. The growth rate is expected to converge in the inviscid limit. The largest value 
$\sigma \sim 0.63$ is reasonably close to the value 
$0.76$ obtained by local stability analysis (Suzuki et al. Reference Suzuki, Hirota and Hattori2018). Since the growth rates of the other modes decrease for large 
$k_z$, the strato-hyperbolic instability is dominant when the wavenumber is bounded from below as 
$k_z \geqslant k_{z0}$ with large 
$k_{z0}$. The frequencies of the oscillatory modes decrease with 
$k_z$.
The modes of the mixed hyperbolic instability exist in the entire range of 
$k_z$ for 
$F_h^{-1}=5$ and 
$F_h^{-1}=10$. The first branch of the mixed hyperbolic instability gives the largest growth rate among all modes: 
$\sigma =0.69$ at 
$k_zL_0=10.1$ for 
$F_h^{-1}=5$ and 
$\sigma =0.77$ at 
$k_zL_0=12.6$ for 
$F_h^{-1}=10$. All branches of the mixed hyperbolic instability take maximum at moderate wavenumbers (
$k_zL_0\approx 9 \sim 44$) and decrease for large wavenumbers. These features are intermediate between the pure hyperbolic instability and the strato-hyperbolic instability. The frequencies of the oscillatory modes behave similarly to those of the strato-hyperbolic instability. For a fixed wavenumber 
$k_z$, the higher frequencies are close to integer multiples of the lowest frequency of the oscillatory modes.
The branches shown in figure 8 involve different symmetry types. Figure 11 shows branches decomposed by symmetry types 
$A_xA_y$, 
$A_xS_y$ and 
$S_xS_y$ for 
$F_h^{-1}=5$; note that the branches of the symmetry type 
$S_xA_y$ coincide with those of 
$A_xS_y$ for 
$\varepsilon _e=0$, thanks to mirror reflection symmetry with respect to 
$y=x$. This shows that the pure hyperbolic instability exists only for 
$A_xA_y$, the strato-hyperbolic instability exists for all symmetry types and the mixed hyperbolic instability exists for 
$A_xA_y$ and 
$A_xS_y$. The mode of the largest growth rate belongs to 
$A_xS_y$. Although the number of branches of the strato-hyperbolic instability depends on the symmetry type, both the growth rate and frequency are similar across the three symmetry types at high wavenumbers. This is because the modes of the strato-hyperbolic instability do not have large amplitude at the boundaries, so that they are less affected by the symmetry conditions at the cell boundaries.
Figure 11. Growth rate and frequency decomposed by symmetry type. 2-D Taylor–Green vortices; 
$\textit {Re}=10^{5}, \varepsilon _e=0, F_h^{-1}=5$: (a,d) 
$A_xA_y$, (b,e) 
$A_xS_y$ and (c,f) 
$S_xS_y$. The symbols are the same as in figure 8.
3.4. Effects of viscosity and density diffusion
Figure 12 shows the growth rate and frequency plotted against 
$k_z$ for 
$F_h^{-1}=5$ and 
$\textit {Re}=10^{4}$. Compared with figure 8(b) for 
$F_h^{-1}=5$ and 
$\textit {Re}=10^{5}$, the growth rate decreases with 
$k_z$ at high wavenumbers, while the curves of the frequency are similar. This implies that essentially the same branches are captured with reduced growth rate due to viscous damping which is proportional to 
$k^{2}$.
Figure 12. Growth rate (a) and frequency (b). 2-D Taylor–Green vortices; 
$\textit {Re}=10^{4}, \varepsilon _e=0, F_h^{-1}=5$. The symbols are the same as in figure 8.
In order to check the effects of density diffusion on the growth rate and frequency, the diffusion term 
$1/(Sc\textit {Re}) \nabla ^{2} \rho '$ was added to (2.6) with 
$F_h^{-1}=5$, 
$\textit {Re}=10^{5}$ and the Schmidt number 
$Sc=1$. The results are shown in figure 13, where the growth rate and frequency are plotted against 
$k_z$. It is similar to figure 8(d) with 
$Sc=\infty$; all branches observed in figure 8(d) are identified with similar growth rates and frequencies in figure 13. The maximum growth rate of the most unstable mixed hyperbolic instability in figure 13 is only 
$0.07\,\%$ smaller than that of figure 8(d), while the growth rate of the most unstable strato-hyperbolic instability in figure 13 is 
$5.0\,\%$ smaller than that of figure 8(d) at 
$k_zL_0=60$. This confirms that the effects of density diffusion are small at high Reynolds numbers and Schmidt number of 
$O(1)$ or larger.
Figure 13. Growth rate (a) and frequency (b). 2-D Taylor–Green vortices; 
$Sc=1, \textit {Re}=10^{5}, \varepsilon _e=0, F_h^{-1}=5$. The symbols are the same as in figure 8.
3.5. Correspondence between the unstable modes and local solutions
It is of interest to check the correspondence between the unstable modes of the strato-hyperbolic instability and the local stability analysis. This has been shown in Suzuki et al. (Reference Suzuki, Hirota and Hattori2018) for one mode by comparing the distributions on the streamline where the amplitude of the mode is nearly the largest with the corresponding local solution. Here, we check it to clarify the structures of not only the dominant modes but also the subdominant modes. The local stability analysis is performed in the same way as in Suzuki et al. (Reference Suzuki, Hirota and Hattori2018); see appendix A for a brief summary of the local stability analysis.
Figure 14 compares the distributions of horizontal divergence 
$\boldsymbol {\nabla }_h \boldsymbol {\cdot } \boldsymbol {u}'_h$, vertical component of vorticity 
$\omega '_z$ and density 
$\rho '$ of the disturbance with those of the corresponding variables 
$p$, 
$q$ and 
$s$ of the solution obtained by the local stability analysis. Here, 
$p \propto \boldsymbol {k}_\perp \boldsymbol {\cdot } \boldsymbol {a}_\perp$, 
$q \propto \boldsymbol {k}_\perp \times \boldsymbol {a}_\perp$ and 
$s \propto r$, where 
$\boldsymbol {k}_\perp$ is the horizontal local wavevector, 
$\boldsymbol {a}_\perp$ is the horizontal amplitude wavevector and 
$r$ is the amplitude of density disturbance. The streamlines of nearly the largest amplitude are chosen as in Suzuki et al. (Reference Suzuki, Hirota and Hattori2018); they are parametrized by one parameter 
$\beta$ which specifies the initial condition in the local stability analysis as
\begin{equation} \boldsymbol{x}(0)=\left(\frac{L_x}{2}(1-\beta), \frac{L_y}{2}, 0\right)^\textrm{T}, \end{equation}
with 
$\beta =0$ and 
$\beta =1$ corresponding to the elliptic stagnation points and the cell boundaries. We focus on stationary modes (
$\omega =0$) since it is difficult to make direct comparisons for oscillating modes. The four modes shown in figure 14 are chosen from different branches: (a) the dominant mode of 
$S_xS_y$ type; (c) the subdominant stationary mode of 
$S_xS_y$ type; (e) the stationary mode of 
$A_xS_y$ type; (g) the stationary mode of 
$A_xA_y$ type. It is confirmed that the time evolution of the corresponding variables in the modal stability analysis is in good agreement with that in the local stability analysis; the phase of the oscillation and the rate of change along the streamline are nearly the same, although there are small differences in the ratio of the amplitudes.
Figure 14. Comparison between modal and local stability analysis, for the unstable eigenmode obtained by modal stability analysis for the stratified 2-D Taylor–Green vortices (
$F_h^{-1}=5$); 
$\varepsilon _e=0, \textit {Re}=10^{5}$: (a,b) 
$S_xS_y$, 
$k_zL_0=41.5$, 
$\sigma =0.625$, 
$\omega =0$, 
$\beta =0.95$; (c,d) 
$S_xS_y$, 
$k_zL_0=60.3$, 
$\sigma =0.413$, 
$\omega =0$, 
$\beta =0.75$; (e,f) 
$A_xS_y$, 
$k_zL_0=44.0$, 
$\sigma =0.471$, 
$\omega =0$, 
$\beta =0.725$; (g,h) 
$A_xA_y$, 
$k_zL_0=52.8$, 
$\sigma =0.556$, 
$\omega =0$, 
$\beta =0.9$. Panels (a,c,e,g) show the values of 
$\boldsymbol {\nabla }_h \boldsymbol {\cdot } \boldsymbol {u}'_h$, 
$\omega '_z$ and 
$\rho '$, and panels (b,d,f,h) those of 
$p$, 
$q$ and 
$s$ on the streamline where 
$\omega '_z$ is maximum, plotted as a function of time. (a,c,e,g) Modal stability analysis; (b,d,f,h) local stability analysis.
It should be pointed out that all local solutions are obtained for the first resonance or the most unstable band. In other words, at least the four branches of strato-hyperbolic instability chosen here are essentially of the same nature. The modes corresponding to higher resonances are not found in the present analysis. This is possibly because unstable regions of higher resonance obtained by local analysis are limited near the cell boundaries (Suzuki et al. Reference Suzuki, Hirota and Hattori2018). However, these modes would exist and can be detected if we can obtain many more modes by analysis with higher resolution. It is emphasized that agreement is obtained both for the streamlines near the boundary (
$\beta =0.9$ and 
$\beta =0.95$, figure 14a,g) and for those inside the cell (
$\beta =0.725$ and 
$\beta =0.75$, figure 14c,e).
3.6. Branches of unstable modes and the growth rate: $\varepsilon _e=0.2$
Figure 15 shows the growth rate and frequency plotted against the vertical wavenumber 
$k_z$ for 
$\varepsilon _e=0.2$. As in the case 
$\varepsilon _e=0$, there is one branch of the pure hyperbolic instability for all values of 
$F_h^{-1}$; the features are similar to the case 
$\varepsilon _e=0$ with the wavenumber range extending to high wavenumber for the stratified case. For the non-stratified case (figure 15a), several branches of the elliptic instability appear since the strain rate at the elliptic stagnation points is non-zero. Some of them form narrow peaks of parametric resonance at low wavenumbers, while they are broad at high wavenumbers. The frequency of the oscillatory modes does not depend very much on the wavenumber. The elliptic instability is stabilized and does not appear for 
$F_h^{-1}=5$ or 
$F_h^{-1}=10$.
Figure 15. Growth rate (a–c) and frequency (d–f). 2-D Taylor–Green vortices with 
$\textit {Re}=10^{5}$, 
$\varepsilon _e=0.2$; 
$F_h^{-1}=0$ in panels (a,d), 
$F_h^{-1}=5$ in panels (b,e), 
$F_h^{-1}=10$ in panels (c,f). The symbols are the same as in figure 8, with 
$\diamond$ (violet) added for the elliptic instability.
The strato-hyperbolic and mixed hyperbolic instabilities appear for the stratified cases as in the case of 
$\varepsilon _e=0$. Their major features are also similar to the case 
$\varepsilon _e=0$, although the wavenumbers are higher than those for 
$\varepsilon _e=0$; for 
$F_h^{-1}=5$ the maximum growth rate of the mixed hyperbolic instability is at 
$k_zL_0=12.7$, which is about 
$1.4$ times the value 
$k_zL_0=8.8$ for 
$\varepsilon _e=0$, while the modes of the strato-hyperbolic instability appear at 
$k_zL_0 \approx 12$, which is larger than the value 
$k_zL_0 \approx 10$ for 
$\varepsilon _e=0$. One of the reasons for this difference is scaling. Although we take the geometric mean of the sides of the cell as the length scale 
$L_0$, the actual scale of the instability is controlled by the shorter side of the cell, which is 
$0.81$ times the length scale 
$L_0$ for 
$\varepsilon _e=0.2$. It is also pointed out that the maximum growth rate of the mixed hyperbolic instability is 
$\sigma =0.85$ at 
$k_zL_0=24.4$ for 
$F_h^{-1}=10$, which is larger than that for 
$\varepsilon _e=0$.
Figures 16 and 17 show the growth rate decomposed by the symmetry type for 
$F_h^{-1}=0$ and 
$F_h^{-1}=5$, respectively. Since there is no mirror symmetry with respect to 
$y=x$ for 
$\varepsilon _e=0.2$, there are differences between 
$A_xS_y$ and 
$S_xA_y$. For 
$F_h^{-1}=0$, the pure hyperbolic instability exists only for 
$A_xA_y$. The narrow peak which gives the maximum growth rate belongs to 
$S_xA_y$. Besides these, however, the branches are quite similar across the four symmetry types, particularly at high wavenumbers. Some modes are missing at high wavenumbers as numerical errors are above the threshold.
Figure 16. Growth rate decomposed by symmetry type. 2-D Taylor–Green vortices; 
$\textit {Re}=10^{5}$, 
$\varepsilon _e=0.2$, 
$F_h^{-1}=0$: (a) 
$A_xA_y$, (b) 
$A_xS_y$, (c) 
$S_xA_y$ and (d) 
$S_xS_y$. The symbols are the same as in figure 15.
Figure 17. Growth rate decomposed by symmetry type. 2-D Taylor–Green vortices; 
$\textit {Re}=10^{5}$, 
$\varepsilon _e=0.2$, 
$F_h^{-1}=5$: (a) 
$A_xA_y$, (b) 
$A_xS_y$, (c) 
$S_xA_y$ and (d) 
$S_xS_y$. The symbols are the same as in figure 8.
For 
$F_h^{-1}=5$, we observe that the stationary mixed hyperbolic instability, which is stronger than the oscillatory mixed hyperbolic instability, appears only for 
$A_xS_y$ and 
$S_xA_y$, and the growth rate of 
$S_xA_y$ is larger than that of 
$A_xS_y$. This is because for the modes of 
$S_xA_y$ the pure hyperbolic instability works on the longer sides of the cells. The strato-hyperbolic instability does not depend on the symmetry type very much, although the number of detected branches is large for 
$S_xS_y$.
4. Stability of stratified Stuart vortices
4.1. Base flow and parameter values
In this section we study the stability of the stratified Stuart vortices (figure 2). The stream function of the Stuart vortices is
\begin{equation} \psi = \log{\left(C\cosh{y}+\sqrt{C^{2}-1}\cos{x}\right)}, \end{equation}
where 
$C \geqslant 1$ is a constant. The vorticity is
\begin{equation} \omega(x,y) = -\frac{1}{\left(C\cosh{y}+\sqrt{C^{2}-1}\cos{x}\right)^{2}}. \end{equation}
The elliptic stagnation points are located at 
$(x, y)=((2n+1){\rm \pi} , 0)$, while the hyperbolic stagnation points are at 
$(x, y)=(2n{\rm \pi} , 0)$. The streamlines are closed inside the heteroclinic streamlines connecting the hyperbolic stagnation points shown by the dot-dashed lines; the heteroclinic streamlines are called separatrices in the present paper, since the mode of the fluid particle motion changes across them. The velocity tends to 
${\boldsymbol {u}}_b \to \pm \boldsymbol {e}_x$ when 
$y \to \pm \infty$, so that the Stuart vortices can be regarded as a model of rolled-up vortices in a mixing layer. We choose the half of the velocity jump as the velocity scale: 
$U_0=1$. The length scale is set to the spatial period in the 
$x$ direction: 
$L_0=2{\rm \pi}$.
The constant 
$C$ is set to 
$C=1.2$ and 
$C=2$. Higher resolution is required for the Stuart vortices than for the 2-D Taylor–Green vortices. Thus, the Reynolds number is fixed at 
$10^{4}$, which is lower than the value 
$10^{5}$ used for most of the cases of 2-D Taylor–Green vortices. The Froude number 
$F_h^{-1}$ is 
$0$, 
$5$ or 
$7.5$. The number of Fourier modes is 
$1024\times 4096$. The values of the Froude number are limited since numerical resolution was not sufficient for 
$F_h^{-1}>7.5$.
As described in § 2.2, the spatial period 
$L_y$ should be large enough to minimize the effects of the periodic boundary condition and the opposite-signed vortices at 
$y=L_y/2$. Table 1 shows how the growth rates and frequencies of two chosen modes depend on 
$L_y$. The differences in the growth rate and the frequency are less than 
$0.15\,\%$ when 
$L_y/L_0 \geqslant 3$. Thus we set 
$L_y/L_0=4$ in the following.
Table 1. Dependence of growth rate (real part) and frequency (imaginary part) on the size of the domain. Stuart vortices with 
$C=1.2$; 
$\textit {Re}=10^{4}$, 
$F_h^{-1}=5$, 
$k_zL_0=88.0$. Mode 1 is the strongest stationary mode; mode 2 is the strongest oscillatory mode.
4.2. Types of instability and their mode structures
For the Stuart vortices, the modes of the strato-hyperbolic instability were not identified. The unstable modes are classified into three types: (i) pure hyperbolic instability (PH), (ii) mixed-type instability (M) and (iii) elliptic instability (E). Note that the mixed-type instability is the mixture of the pure hyperbolic instability and the elliptic instability, which is different from the mixed hyperbolic instability observed for the 2-D Taylor–Green vortices. Examples of each instability type are shown below.
Now we explain how the types of modes are identified, since classification is not so obvious as in the case of the 2-D Taylor–Green vortices. Since the mode structures can change continuously along a branch, as found in Potylitsin & Peltier (Reference Potylitsin and Peltier1999), we should introduce a criterion for classification of the modes. Figure 18 shows the ratio of enstrophy for a branch along which the type of instability changes from the mixed-type instability to the pure hyperbolic instability. This branch is dominant at all wavenumbers for 
$C=1.2, F_h^{-1}=0$ (figure 22a). We define the core region, which includes the elliptic stagnation point and strong vorticity, by 
$(\psi -\psi _e)/(\psi _h-\psi _e) \leqslant 0.5$, where 
$\psi _e$ and 
$\psi _h$ are the values of the stream function at the elliptic and hyperbolic stagnation points, respectively; the separatrix region is defined as 
$0.75\leqslant (\psi -\psi _e)/(\psi _h-\psi _e) < 1.25$, so that it includes the hyperbolic stagnation points and the separatrix. Enstrophy of the disturbance is then integrated in each region to give the ratio divided by the total enstrophy of the disturbance; enstrophy ratios of the core region and the separatrix region are denoted by 
$\chi _c$ and 
$\chi _s$, respectively. In figure 18 the enstrophy ratio of the core region is about 
$0.5$ at low wavenumbers, which implies that the mode has large vorticity due to the elliptic instability; it decreases to 
$0.1 \sim 0.4$ at 
$25 \leqslant k_zL_0 \leqslant 35$, showing that vorticity near the separatrix increases due to the pure hyperbolic instability; and it almost vanishes at high wavenumbers as the mode is dominated by the pure hyperbolic instability. We define the elliptic, mixed and pure hyperbolic instabilities by the enstrophy ratios of the core region and the separatrix region: the mode is named (i) the mixed-type instability when 
$\chi _c \geqslant 0.1$ and 
$\chi _s \geqslant 0.2$; (ii) the elliptic instability when 
$\chi _c \geqslant 0.1$ and 
$\chi _s < 0.2$; and (iii) the pure hyperbolic instability when 
$\chi _c < 0.1$ and 
$\chi _s \geqslant 0.2$. There was no mode which could not be classified into any of the three instabilities above. Although these values are chosen just for convenience, they lead to a reasonable classification, as confirmed below.
Figure 18. (a) Ratio of enstrophy along a branch. Stuart vortices with 
$C=1.2$; 
$F_h^{-1}=0$. Magenta: core region. Blue: separatrix region. The dashed line separates the type of instability along this branch: (M) indicates mixed-type instability, (PH) pure hyperbolic instability. (b) Core and separatrix regions. Three red lines show the streamlines where 
$(\psi -\psi _e)/(\psi _h-\psi _e)=0.5$, which is the boundary of the core region (denoted by 
$c$), and 
$(\psi -\psi _e)/(\psi _h-\psi _e)=0.75$ and 
$(\psi -\psi _e)/(\psi _h-\psi _e)=1.25$, which bound the separatrix region (denoted by 
$s$). The separatrix is shown by the dashed line.
Figure 19 shows examples of the pure hyperbolic instability modes. They look similar to those found in Potylitsin & Peltier (Reference Potylitsin and Peltier1998) and Potylitsin & Peltier (Reference Potylitsin and Peltier1999), although the resolution is better in the present results. Thin regions of strong vorticity are observed near the separatrix. They are wider for the stratified case than for the non-stratified case. This is an effect of stratification; the ratio of the vertical scale to the horizontal scale of the hyperbolic instability is decreased owing to stratification, so that the horizontal scale becomes large for the same vertical wavenumber 
$k_zL_0=88.0$. Both modes have symmetry under rotation by 
${\rm \pi}$ about the elliptic stagnation point: 
$\omega '_z(-(x-L_x/2), y, z)=-\omega '_z(x-L_x/2, y, z)$. The structures are rather simple along a streamline, as the phase change 
$\omega '_z$ is 
$2{\rm \pi}$ in one period of fluid particle motion.
Figure 19. Structures of pure hyperbolic instability modes. Stuart vortices with 
$C=1.2$. (a,b) Iso-surface of the magnitude of vorticity coloured by 
$\omega '_z$; (c,d) contours of 
$\omega '_z$ on the 
$xy$ plane; (e,f) contours of 
$\omega '_x$ on 
$x/L_0=0.5$ near the separatrix, whose position is denoted by 
$y_{{sp}}$. The vertical domain sizes are 
$4L_z$ and 
$2L_z$ for the iso-surface and the contours of 
$\omega '_x$, respectively. The thick black dashed lines in the contours of 
$\omega '_z$ show the separatrix of the base flow. (a,c,e) 
$F_h^{-1}=0$, 
$k_zL_0=88.0$, 
$\sigma =1.896$, 
$\omega =0$. (b,d,f) 
$F_h^{-1}=5$, 
$k_zL_0=88.0$, 
$\sigma =1.645$, 
$\omega =0$. The magnitude of vorticity on the iso-surface is 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.25$ in (a) and 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.3$ in (b). The contours of 
$\omega '_z$ are drawn on 
$z/L_z=0$.
Figure 20 shows examples of the mixed-type instability modes. Strong vorticity is observed in the core region as well as near the separatrix. The structures in the core region are due to the elliptic instability as shown in figure 21, while those along the separatrix are due to the pure hyperbolic instability as shown in figure 19. Thus these modes are classified as the mixed-type instability.
Figure 20. Structures of mixed-type instability modes. Stuart vortices with 
$C=1.2$. (a,b) Iso-surface of the magnitude of vorticity coloured by 
$\omega'_{z}$; (c,d) contours of 
$\omega'_{z}$ on the 
$xy$ plane. The vertical domain size for the iso-surface is 
$3L_z$. The thick black dashed lines in the contours of 
$\omega'_{z}$ show the separatrix of the base flow: (a,c) 
$F_h^{-1}=0$, 
$k_zL_0=30.2$, 
$\sigma =1.846$, 
${\omega =0}$; (b,d) 
$F_h^{-1}=5$, 
$k_zL_0=30.2$, 
$\sigma =1.379$, 
$\omega =0$. The magnitude of vorticity on the iso-surface is 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.2$ in (a) and 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.25$ in (b). The contours of 
$\omega '_z$ are drawn on 
$z/L_z=0$.
Figure 21. Structures of elliptic instability modes. Stuart vortices with 
$C=1.2$, same as figure 20 except that the vertical domain size for the iso-surface is 
$2L_z$: (a,c) 
$F_h^{-1}=0$, 
$k_zL_0=15.1$, 
$\sigma =1.525$, 
$\omega =0$; (b,d) 
$F_h^{-1}=5$, 
$k_zL_0=40.2$, 
$\sigma =1.053$, 
$\omega =0$. The magnitude of vorticity on the iso-surface is 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.5$ in (a) and 
$|\boldsymbol {\omega }'|/|\boldsymbol {\omega }'|_{{max}}=0.2$ in (b). The contours of 
$\omega '_z$ are drawn on 
$z/L_z=0$.
Figure 21 shows examples of the elliptic instability modes. Vorticity distributions in the core region are similar to those of the elliptic instability modes of the 2-D Taylor–Green vortices (figure 6). Vorticity is weak outside the core region. The iso-surface of the magnitude of vorticity in figure 21(a) shows typical structures of bending waves.
4.3. Branches of unstable modes and the growth rate: $C=1.2$
Figure 22 shows the growth rate and frequency plotted against the vertical wavenumber 
$k_z$ for 
$C=1.2$. As explained in the previous subsection, the type of mode can change along a branch, in contrast to the case of the 2-D Taylor–Green vortices. The elliptic instability exists mostly at low wavenumbers. It forms narrow bands due to parametric resonance, while a wide band is also observed for 
$F_h^{-1}=5$ and 
$F_h^{-1}=7.5$. It is worth noting that the elliptic instability is not completely stabilized for the stratified case, although the growth rate is reduced; it would be stabilized for stronger stratification. There is an oscillatory mode of the elliptic instability which does not change the type for all values of 
$F_h^{-1}$.
Figure 22. Growth rate (a–c) and frequency (d–f). Stuart vortices with 
$C=1.2$; 
$\textit {Re}=10^{4}$: (a,d) 
$F_h^{-1}=0$; (b,e) 
$F_h^{-1}=5$; (c,f) 
$F_h^{-1}=7.5$. Symbols: 
$\square$ (blue) denotes pure hyperbolic instability; 
$\triangle$ (green) denotes mixed-type instability; 
$\diamond$ (violet) denotes elliptic instability. Solid symbols and open symbols stand for stationary (
$\omega =0$) and oscillatory (
$\omega \neq 0$) modes, respectively. Only the oscillatory modes with positive frequencies are shown.
The dominant branch changes type from the mixed-type instability to the pure hyperbolic instability for 
$F_h^{-1}=0$ and 
$F_h^{-1}=5$, while the elliptic instability is missing for 
$F_h^{-1}=7.5$. At high wavenumbers the pure hyperbolic instability is dominant; no mode of the elliptic or mixed-type instabilities is observed at 
$k_zL_0 \geqslant 100$. The growth rate decreases at high wavenumbers due to viscous effects. This decrease would not be observed for higher Reynolds numbers, although the verification of this requires much higher numerical costs and is not addressed in the present paper. In the inviscid limit, the growth rate is expected to converge to a finite value depending on the branch. The growth rate of the dominant branch decreases slightly owing to the effects of stratification. The number of the obtained branches also decreases with stratification as the branches shift to higher wavenumbers at which viscous damping is strong; the decrease is partly due to limited numerical resolution. The frequencies of the pure hyperbolic instability modes do not depend on stratification significantly except that they are scattered for 
$F_h^{-1}=7.5$.
4.4. Branches of unstable modes and the growth rate: $C=2$
Figure 23 shows the growth rate and frequency plotted against the vertical wavenumber 
$k_z$ for 
$C=2$. The most important difference from the 
$C=1.2$ case is that the elliptic instability is stronger, so that the narrow peaks are more pronounced, and the growth rate does not decrease very much when 
$F_h^{-1}=5$; as a result the most unstable mode is the elliptic instability at 
$k_zL_0=45.2$ for 
$F_h^{-1}=5$. There are also branches of the elliptic instability occupying the whole wavenumber range considered in the present paper. The reason for this difference is that the Froude number based on the vortex core is larger for 
$C=2$ than for 
$C=1.2$. In fact, we may define 
$\tilde {F}_h = \omega _{{max}}/N$, where 
$\omega _{{max}}$ is the maximum of the magnitude of vorticity, as the Froude number based on the vortex core; then 
$\tilde {F}_h/F_h=21.8$ for 
$C=1.2$ and 
$\tilde {F}_h/F_h=87.5$ for 
$C=2$, implying that effect of stratification with same value of 
$F_h$ is weaker in the vortex core for 
$C=2$ than for 
$C=1.2$.
Figure 23. Growth rate (a,c) and frequency (b,d). Stuart vortices with 
$C=2$; 
$\textit {Re}=10^{4}$: (a,c) 
$F_h^{-1}=0$; (b,d) 
$F_h^{-1}=5$. The symbols are the same as in figure 22.
4.5. Differences between 2-D Taylor–Green vortices and Stuart vortices
Now we discuss the differences between the 2-D Taylor–Green vortices and the Stuart vortices in more detail. The strato-hyperbolic instability was found for the 2-D Taylor–Green vortices, while it was not found for the Stuart vortices; or, equivalently, only the dominant instability near the cell boundaries or the separatrix was observed, except the single branch of the pure hyperbolic instability which survives in the stratified cases. For the 2-D Taylor–Green vortices, although the strato-hyperbolic instability forms multiple unstable bands according to the local stability analysis (Suzuki et al. Reference Suzuki, Hirota and Hattori2018), all the observed modes of the strato-hyperbolic instability correspond to the first and dominant unstable band (figure 14). For the Stuart vortices, the strato-hyperbolic instability is weaker than the pure hyperbolic instability which survives for the stratified case; all modes of the hyperbolic instability are those of the pure hyperbolic instability. It is not easy to capture the modes corresponding to higher resonances, since both the unstable region and the growth rate are smaller than for the modes corresponding to the first resonance according to the local stability analysis; these modes would be much more localized near the cell boundaries and would consist of high-wavenumber Fourier modes. We expect to observe these modes if we increase the spatial resolution and the Reynolds number, although the numerical costs are high.
The elliptic instability is stabilized for 
$F_h^{-1}=5$ and 
$F_h^{-1}=10$ for the 2-D Taylor–Green vortices, while it is not completely stabilized for the Stuart vortices within the range of Froude numbers studied in the present work. However, stronger stratification would stabilize the elliptic instability of the Stuart vortices as shown by local stability analysis: when 
$C=1.2$, it exists for 
$F_h^{-1}=5$ but does not for 
$F_h^{-1}=10$ (Suzuki et al. Reference Suzuki, Hirota and Hattori2018). In fact, the growth rates for 
$F_h^{-1}=7.5$ are smaller than those for 
$F_h^{-1}=0$ in figure 22. Stronger stratification is required for 
$C=2$ to stabilize the elliptic instability completely, since the strain rate at the elliptic stagnation points is larger than 
$C=1.2$.
5. Concluding remarks
We have studied the linear stability of a periodic array of vortices in stratified fluid by modal stability analysis. Two base flows were considered: the 2-D Taylor–Green vortices and the Stuart vortices. In the case of the 2-D Taylor–Green vortices, four types of instability were identified: the pure hyperbolic instability, the strato-hyperbolic instability, the mixed hyperbolic instability and the elliptic instability. There is only one branch of the pure hyperbolic instability. Although it is the most unstable in the non-stratified case, it is surpassed by the strato-hyperbolic instability and the mixed hyperbolic instability in the stratified case. The strato-hyperbolic instability and the mixed hyperbolic instability are caused by stratification. The mixed hyperbolic instability is a mixture of the pure hyperbolic instability, which has large amplitude along the cell boundaries, and the strato-hyperbolic instability, which has large amplitude near the boundaries but inside the cells. The mixed hyperbolic instability gives the maximum growth rate when 
$F_h^{-1}=5$ and 
$F_h^{-1}=10$, while the strato-hyperbolic instability is dominant at large wavenumbers (e.g. 
$k_zL_0 \gtrsim 30$ for 
$F_h^{-1}=5$ and 
$\varepsilon _e=0$). It should be emphasized that the growth rate of the strato-hyperbolic instability tends to a constant along each branch at the large-wavenumber limit, while that of the mixed hyperbolic instability decreases with the wavenumber. Many branches of the elliptic instability appear when 
$\varepsilon _e \neq 0$; however, they are stabilized when 
$F_h^{-1}=5$ and 
$F_h^{-1}=10$.
In the case of the Stuart vortices, the unstable modes were classified into three types: the pure hyperbolic instability, the elliptic instability and the mixed-type instability, which is a mixture of the pure hyperbolic instability and the elliptic instability. In contrast to the 2-D Taylor–Green vortices, the type of instability for the Stuart vortices can change along a branch; it is elliptic or mixed-type at low wavenumbers, while it becomes pure hyperbolic at large wavenumbers. The pure hyperbolic instability modes have large amplitude along the separatrix of the base flow, while the elliptic instability modes have large amplitude in the core region; the mixed-type instability modes have both features simultaneously. The growth rate of the elliptic instability decreases with stratification for 
$C=1.2$, although it is not completely stabilized. For 
$C=2$, on the other hand, the elliptic instability is not affected very much by stratification at 
$F_h^{-1}=5$; stronger stratification is expected to stabilize the elliptic instability, although this has not been confirmed due to limited numerical resolution.
It should be pointed out that the strato-hyperbolic instability has not been found for the Stuart vortices. As discussed in the previous section, however, it would be possible to capture modes of the strato-hyperbolic instability by increasing spatial resolution. This also applies to the modes of the strato-hyperbolic instability corresponding to higher resonance in the case of 2-D Taylor–Green vortices. Selected modes of the strato-hyperbolic instability have been shown to correspond to solutions of local stability analysis which are due to the first resonance. It is also pointed out that the strato-hyperbolic instability in the 2-D Taylor–Green vortices and the pure hyperbolic instability in the Stuart vortices are not stabilized by strong stratification; Potylitsin & Peltier (Reference Potylitsin and Peltier1998) claimed that the three-dimensional instability would appear to be entirely suppressed by strong density stratification, but this must be corrected. It is possible that stratification makes the vertical wavenumber high so that the growth rate is reduced by viscous damping; however, the instability is not stabilized in the inviscid limit.
The present results imply that both the pure hyperbolic instability and the strato-hyperbolic instability are important in stably stratified flows of geophysical or planetary scale. Since the vertical scale is normally much smaller than the horizontal scale in these flows, the vertical wavenumber 
$k_z$ is large. Then the dominant instability can be the pure hyperbolic instability or the strato-hyperbolic instability when there is an array of vortices. The present results on linear stability are indispensable for understanding the dynamics of the stratified vortices. When these instability modes grow into nonlinear stages, the core of the vortices can survive as turbulence is localized near the cell boundaries or the separatrix where these modes are concentrated. In the present paper, we have assumed that the Brunt–Väisälä frequency 
$N$ is constant, which is not always the case in actual stratified flows of geophysical or planetary scale. There can arise interesting phenomena in non-uniformly stratified flows. Even if 
$N^{2}$ varies in the vertical direction, however, both the pure hyperbolic instability and the strato-hyperbolic instability are expected to occur since they are short-wave instabilities for which WKB renormalization would work. They can appear as disturbances localized in the vertical direction, while global modes are possible for small variations of 
$N^{2}$, since the instability branches exist in wide ranges of the vertical wavenumber.
Some topics for future work are listed below. Unstable modes due to higher resonance, i.e. strato-hyperbolic instability modes of second or higher bands in the 2-D Taylor–Green vortices and strato-hyperbolic instability modes in the Stuart vortices, should be searched for with increased numerical resolution. The effects of stratification on other flows possessing hyperbolic stagnation points, which include vortex pairs and wake vortices such as the von Kármán vortex street, are of interest. The case of non-uniform stratification is of practical importance. The effects of rotation should be also investigated since they are often important in geophysical and planetary problems. Finally, the role of the strato-hyperbolic and the pure hyperbolic instabilities in nonlinear dynamics of arrays of vortices is of great interest. If unstable modes develop into turbulence near the cell boundaries or the separatrix, strong dissipation due to turbulence may result in concentration of vorticity (Hattori Reference Hattori2016, Reference Hattori2018), which may be one of the mechanisms preventing the breakdown of long-lived vortices.
Acknowledgements
This work was supported by JSPS KAKENHI 17K05561. Numerical calculations were performed on the UV1000, UV2000 and PRIMERGY at the Institute of Fluid Science, Tohoku University.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Summary of local stability analysis
Here we briefly summarize the method of the local stability analysis used in § 3.5; see Suzuki et al. (Reference Suzuki, Hirota and Hattori2018) for the details.
In the local stability analysis the disturbance is assumed to be in the form of wave packet
\begin{equation} \boldsymbol{u}'=\left(\hat{\boldsymbol{u}}_0+\delta\hat{\boldsymbol{u}}_1+\cdots \right)\exp\left(\frac{\textrm{i}}{\delta}\varPhi\right) \end{equation}
with similar expressions for 
$p'$ and 
$\rho '$, where 
$\delta$ is a small parameter and 
$\varPhi$ is eikonal satisfying 
$\partial \varPhi /\partial t+ \boldsymbol {u}_b \boldsymbol {\cdot } \boldsymbol {\nabla }\varPhi =0$. Substituting these expressions into (2.4)–(2.6) yields a set of ordinary differential equations at the leading order:
\begin{gather} \boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{k}=0, \end{gather}
\begin{gather}\frac{\textrm{d}\boldsymbol{x}}{\textrm{d}t}=\boldsymbol{u}_b, \end{gather}
\begin{gather}\frac{\textrm{d}\boldsymbol{k}}{\textrm{d}t}=-{\boldsymbol{\mathsf{L}}}^{T}\boldsymbol{k}, \end{gather}
\begin{gather}\frac{\textrm{d}\boldsymbol{a}}{\textrm{d}t}= \left(2\frac{\boldsymbol{k}\boldsymbol{k}^{T}}{|\boldsymbol{k}|^{2}}- {\boldsymbol{\mathsf{I}}}\right){\boldsymbol{\mathsf{L}}}\boldsymbol{a} +\left(\frac{\boldsymbol{k}\boldsymbol{k}^{T}}{|\boldsymbol{k}|^{2}}-{\boldsymbol{\mathsf{I}}}\right)r\boldsymbol{e}_3, \end{gather}
\begin{gather}\frac{\textrm{d}r}{\textrm{d}t}= \frac{1}{F_h^{2}}a_3, \end{gather}
where 
$\boldsymbol {k}=\boldsymbol {\nabla }\varPhi$, 
$\boldsymbol {a}=\hat {\boldsymbol {u}}_0$, 
$r=\hat {\rho }_0$, 
${\mathsf{L}}_{ij}=\partial u_{b,i}/\partial x_j$ and 
${\boldsymbol{\mathsf{I}}}$ is the identity matrix (Friedlander & Vishik Reference Friedlander and Vishik1991).
We consider periodic orbits of fluid particles. We also assume that the wavevector 
$\boldsymbol {k}$ is time-periodic, which is a necessary condition for exponential instability on the periodic orbits. It is known that 
$\boldsymbol {k}$ is time-periodic if it is perpendicular to the streamline initially:
\begin{equation} \boldsymbol{k}(0)\boldsymbol{\cdot}\boldsymbol{u}_b(\boldsymbol{x}(0))=0. \end{equation}
Then the time evolution of amplitude is described by a Floquet matrix 
${\boldsymbol{\mathsf{F}}}$, since the matrices which appear in (A 5) are also time-periodic:
\begin{equation} \{ \boldsymbol{a}, r\}(t+T)={\boldsymbol{\mathsf{F}}}(T)\{ \boldsymbol{a}, r\}(t), \end{equation}
where 
$T$ is the period of 
$\boldsymbol {k}$, which coincides with that of the particle motion 
$\boldsymbol {x}$. The growth rate is obtained as
\begin{equation} \sigma_i=\frac{\log{|\mu_i|}}{T}, \end{equation}
where 
$\mu _i$ is an eigenvalue of 
${\boldsymbol{\mathsf{F}}}(T)$.
It is useful to introduce
\begin{equation} p=\frac{|\boldsymbol{k}|}{|\boldsymbol{k}_\perp|}\boldsymbol{k}_\perp \boldsymbol{\cdot} \boldsymbol{a}_\perp, \quad q=\frac{|\boldsymbol{k}|}{|\boldsymbol{k}_\perp|}(\boldsymbol{k}_\perp \times \boldsymbol{a}_\perp) \boldsymbol{\cdot} \boldsymbol{e}_3, \quad s=\frac{|\boldsymbol{k}|}{|\boldsymbol{k}_\perp|}r \end{equation}
as in Bayly, Holm & Lifschitz (Reference Bayly, Holm and Lifschitz1996), where 
$\boldsymbol {k}_\perp =(k_1, k_2)^\textrm {T}$ and 
$\boldsymbol {a}_\perp =(a_1, a_2)^\textrm {T}$. Then (A 5) and (A 6) are reduced to
\begin{equation} \frac{\textrm{d}}{\textrm{d}t}\left( \begin{array}{c}p\\q\\s\end{array}\right) = \left( \begin{array}{ccc} \dfrac{\textrm{d}}{\textrm{d}t}\log{\dfrac{|\boldsymbol{k}_\perp|}{|\boldsymbol{k}|}} & \dfrac{2k^{2}_3{\boldsymbol{\mathsf{H}}}\boldsymbol{k}_\perp\boldsymbol{\cdot} \boldsymbol{k}_\perp}{|\boldsymbol{k}|^{2}|\boldsymbol{k}_\perp|^{2}} & \dfrac{|\boldsymbol{k}_\perp|^{2}}{|\boldsymbol{k}|^{2}}k_3\\ -W & -\dfrac{\textrm{d}}{\textrm{d}t}\log{\dfrac{|\boldsymbol{k}_\perp|}{|\boldsymbol{k}|}} & 0\\ -\dfrac{1}{F_h^{2} k_3} & 0 & -\dfrac{\textrm{d}}{\textrm{d}t}\log{\dfrac{|\boldsymbol{k}_\perp|}{|\boldsymbol{k}|}} \end{array}\right) \left(\begin{array}{c} p\\q\\s\end{array}\right), \end{equation}
where 
$W={\mathsf{L}}_{12}-{\mathsf{L}}_{21}$ is the vorticity of the base flow, 
${\boldsymbol{\mathsf{L}}}_\perp =\left (\begin {smallmatrix}{\mathsf{L}}_{11} & {\mathsf{L}}_{12}\\ {\mathsf{L}}_{21} & {\mathsf{L}}_{22}\end {smallmatrix}\right )$, 
${\boldsymbol{\mathsf{H}}}={\boldsymbol{\mathsf{L}}}_\perp \left (\begin {smallmatrix}0 & 1\\-1 & 0\end {smallmatrix}\right )$ and 
${\boldsymbol{\mathsf{L}}}=\left (\begin {smallmatrix}{\boldsymbol{\mathsf{L}}}_\perp & \boldsymbol {0}\\ \boldsymbol {0}^{T} & 0\end {smallmatrix}\right )$. It is pointed out that the number of variables is reduced from four (
$a_1$, 
$a_2$, 
$a_3$, 
$r$) to three (
$p$, 
$q$, 
$s$) by virtue of the incompressibility condition (A 2). Since the horizontal derivative 
$\boldsymbol {\nabla }_h$ is replaced by 
${\textrm {{i}}} \boldsymbol {k}_\perp$ for the disturbance in the form of (A 1), the horizontal divergence 
$\boldsymbol {\nabla }_h \boldsymbol {\cdot }\boldsymbol {u}'_h$ becomes 
${\textrm {{i}}} \boldsymbol {k}_\perp \boldsymbol {\cdot } \boldsymbol {a}_\perp \propto p$ and the vertical component of vorticity 
$\omega '_z$ becomes 
${\textrm {{i}}} (\boldsymbol {k}_\perp \times \boldsymbol {a}_\perp ) \boldsymbol {\cdot } \boldsymbol {e}_3 \propto q$, while 
$\rho ' \propto s$ by the definition (A 10a–c).
