2.1. Governing equations

We consider the linear stability of the Stuart vortices in stably stratified and rotating fluids. The Stuart vortices are a one-dimensional array of periodic vortices with vorticity of the same sign (figure 1). The streamfunction and the vorticity are

(2.1)$$\begin{gather} \psi = \log{\left(C\cosh{y}+\sqrt{C^2-1}\cos{x}\right)}, \end{gather}$$

(2.2a,b)$$\begin{gather}\boldsymbol{\omega}_b = \omega_b(x,y)\,\boldsymbol{e}_z, \quad \omega_b(x,y) ={-}\frac{1}{\left(C\cosh{y}+\sqrt{C^2-1}\cos{x}\right)^2}, \end{gather}$$

where $C \geqslant 1$ is a constant. In the previous works (Leblanc & Cambon Reference Leblanc and Cambon1998; Godeferd et al. Reference Godeferd, Cambon and Leblanc2001), a different expression, $\psi = \log {(\cosh {y}+\rho _s\cos {x})}$, is used instead of the above streamfunction; the two are related by $\rho _s=\sqrt {1-C^{-2}}$, with a constant that does not affect the flow field. Two adjacent vortices are connected by a hyperbolic point. The vorticity is parallel to the vertical direction, which coincides with the direction of gravity force. The effects of density stratification are taken into account by the Boussinesq approximation. Then the Stuart vortices are steady in the absence of diffusion (Hattori & Hirota Reference Hattori and Hirota2023).

Figure 1.

Stuart vortices: (a) streamlines and (b) contours of vorticity distribution for $C=1.2$; (c) comparison of vorticity distribution on $y=0$ between $C=1.2$ (solid line), $1.6$ (dotted line) and $2$ (dashed line). In (a), the red, green and blue lines show the closed streamlines, the separatrix and the open streamlines, respectively.

Viscosity is taken into account in general, while diffusion of density is neglected since its effects are negligible (Hattori et al. Reference Hattori, Suzuki, Hirota and Khandelwal2021). The base flow is assumed steady because the growth of instabilities is much faster than the time evolution of the base flow due to viscous diffusion in the high-Reynolds-number flows considered in this paper. The velocity, pressure and density fields are decomposed as

(2.3)$$\begin{gather} \boldsymbol{u}=\boldsymbol{u}_b+\boldsymbol{u}', \end{gather}$$

(2.4)$$\begin{gather}p=p_b+p', \end{gather}$$

(2.5)$$\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'_x, u'_y, u'_z, p', \rho ')$ are the base flow and the disturbance, the direction of the gravity force is taken as $-\boldsymbol {e}_z$, and the base density is assumed to be $\rho _b=\rho _0+\alpha z$, with $\alpha =\partial \rho _b/\partial z<0$ being a constant. The magnitude of the disturbance is infinitesimally small. Then the governing equations in non-dimensionalized form are

(2.6)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'=0, \end{gather}$$

(2.7)$$\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}' - \frac{1}{Ro}\,\boldsymbol{e}_z \times \boldsymbol{u}' ={-}\boldsymbol{\nabla}{p'}-\rho'\boldsymbol{e}_z+\frac{1}{Re}\,\nabla^2 \boldsymbol{u}', \end{gather}$$

(2.8)$$\begin{gather}\frac{\partial{\rho'}}{\partial{t}}+(\boldsymbol{u}_b\boldsymbol{\cdot}\boldsymbol{\nabla}){\rho'}-\frac{1}{F_h^2}\,u'_z =0, \end{gather}$$

where $Ro =U_0/(2 \varOmega _0 L_0)$ is the Rossby number, $Re =U_0L_0/\nu$ is the Reynolds number, $F_h=U_0/(L_0N)$ is the Froude number based on the horizontal scale, $N=\sqrt {-{\alpha }g/{\rho }_0}$ is the Brunt–Väisälä frequency, $g$ is the acceleration of gravity, $\nu$ is the kinematic viscosity, $\varOmega _0$ is the angular velocity, and $U_0$ and $L_0$ are a characteristic velocity and a length scale, respectively. As in our previous works (Suzuki et al. Reference Suzuki, Hirota and Hattori2018; Hattori et al. Reference Hattori, Suzuki, Hirota and Khandelwal2021), we set $L_0=2{\rm \pi}$, which is the spatial period in the $x$ direction, while $U_0$ is set to the maximum velocity $1$. In the following, the values are scaled by $U_0$ and $L_0$ unless stated explicitly. It is pointed out that the minus sign of the Coriolis force in (2.7) is introduced because the vorticity distribution of the Stuart vortices in (2.2a,b) is negative; $Ro>0$ and $Ro<0$ correspond to the cyclonic and anticyclonic cases, respectively.

In the local stability analysis, the disturbance is assumed to be in the form of a wave packet:

(2.9)$$\begin{gather} \boldsymbol{u}'=\left(\hat{\boldsymbol{u}}_0+\delta\hat{\boldsymbol{u}}_1+\cdots \right)\exp \left(\frac{{\rm i}}{\delta}\,\varPhi\right), \end{gather}$$

(2.10)$$\begin{gather}p'=\left(\hat{p}_0+\delta\hat{p}_1+\cdots \right)\exp\left(\frac{{\rm i}}{\delta}\,\varPhi\right), \end{gather}$$

(2.11)$$\begin{gather}\rho'=\left(\hat{\rho}_0+\delta\hat{\rho}_1+\cdots \right)\exp\left(\frac{{\rm i}}{\delta}\,\varPhi\right), \end{gather}$$

where $\delta$ is a small parameter, and $\varPhi$ is eikonal, satisfying ${\rm D}\varPhi /{\rm D}t=0$ (where ${\rm D}/{\rm D}t=\partial /\partial t+ \boldsymbol {u}_b \boldsymbol {\cdot } \boldsymbol {\nabla }$). Viscosity is neglected in the local stability analysis. Substituting the above expressions into (2.6)–(2.8) yields a set of ordinary differential equations at the leading order:

(2.12)$$\begin{gather} {\frac{{\rm{d}} {{\boldsymbol{X}}}}{{\rm{d}} {t}}} = {\boldsymbol{U}}({\boldsymbol{X}}), \end{gather}$$

(2.13)$$\begin{gather}{\frac{{\rm{d}} {{\boldsymbol{k}}}}{{\rm{d}} {t}}} ={-}{\boldsymbol{\mathsf{L}}}^{\rm T} {\boldsymbol{k}}, \end{gather}$$

(2.14)$$\begin{gather}{\frac{{\rm{d}} {{\boldsymbol{a}}}}{{\rm{d}} {t}}} = \left( 2{\hat{\boldsymbol{k}}}{\hat{\boldsymbol{k}}}^{\rm T} - {\boldsymbol{\mathsf{I}}} \right) {\boldsymbol{\mathsf{L}}}{\boldsymbol{a}} +\left({\hat{\boldsymbol{k}}}{\hat{\boldsymbol{k}}}^{\rm T}-{\boldsymbol{\mathsf{I}}}\right)r\boldsymbol{e}_z - \frac{1}{Ro} \left( {\hat{\boldsymbol{k}}}{\hat{\boldsymbol{k}}}^{\rm T} - {\boldsymbol{\mathsf{I}}} \right) {\boldsymbol{e}_z} \times {\boldsymbol{a}}, \end{gather}$$

(2.15)$$\begin{gather}\frac{{\rm d}r}{{\rm d}t}= \frac{1}{F_h^2}\,a_z, \end{gather}$$

where ${\boldsymbol{\mathsf{L}}}_{ij} = \partial U_i/\partial x_j$ and $\hat {\boldsymbol {k}}=\boldsymbol {k}/|{\boldsymbol {k}}|$ (Friedlander & Vishik Reference Friedlander and Vishik1991; Lifschitz & Hameiri Reference Lifschitz and Hameiri1991; Leblanc Reference Leblanc1997). Here, ${\boldsymbol {X}}$ is the position of the fluid particle, and $\boldsymbol {k}=\boldsymbol {\nabla }\varPhi$ is the local wavevector, while ${\boldsymbol {a}}=\hat {\boldsymbol {u}}_0$ and $r=\hat {\rho }_0$ are the amplitudes of the disturbance corresponding to velocity and density, respectively. The incompressibility condition (2.6) leads to ${\boldsymbol {a}} \boldsymbol {\cdot } {\boldsymbol {k}} = 0$, which is satisfied for $t>0$ if it holds at $t=0$. The base flow is unstable if the amplitude $\{{\boldsymbol {a}}, r\}$ grows without bound.

The number of variables in (2.14) and (2.15) can be reduced from four to three by introducing

(2.16a–c)\begin{equation} p = \frac{k}{k_\perp}\,{\boldsymbol{k}}_\perp \boldsymbol{\cdot}{\boldsymbol{a}}_\perp{=}{-}\frac{kk_z}{k_\perp}\,a_z, \quad q = \left(\frac{k}{k_\perp}\,{\boldsymbol{k}}_\perp{\times} {\boldsymbol{a}}_\perp\right)\boldsymbol{\cdot}{\boldsymbol{e}}_z, \quad s = \frac{k}{k_\perp}\,r \end{equation}

as in Bayly, Holm & Lifschitz (Reference Bayly, Holm and Lifschitz1996), where $\boldsymbol {k}_\perp =(k_x, k_y)^{\rm T}$, $\boldsymbol {a}_\perp =(a_x, a_y)^{\rm T}$ and

(2.17a,b) \begin{equation} {\boldsymbol{\mathsf{L}}} = \left(\begin{array}{@{}cc@{}} {\boldsymbol{\mathsf{L}}}_\perp & 0 \\ 0 & 0 \end{array} \right), \quad {\boldsymbol{\mathsf{H}}} = {\boldsymbol{\mathsf{L}}}_\perp \left( \begin{array}{@{}cc@{}} 0 & 1 \\ -1 & 0 \end{array} \right). \end{equation}

Then the equations (2.14) and (2.15) reduce to

(2.18) \begin{equation} \dfrac{{\rm d}}{{\rm d}t}\left( \begin{array}{@{}c@{}}p\\q\\s\end{array}\right) = \left(\begin{array}{@{}ccc@{}} \dfrac{{\rm d}}{{\rm d}t}\log{\dfrac{|\boldsymbol{k}_\perp|}{|\boldsymbol{k}|}} & \dfrac{2k^2_z\boldsymbol{\mathsf{H}}\boldsymbol{k}_\perp \boldsymbol{\cdot}\boldsymbol{k}_\perp}{|\boldsymbol{k}|^2\,|\boldsymbol{k}_\perp|^2}+\dfrac{k_z^2}{Ro\,k^2} & \dfrac{|\boldsymbol{k}_\perp|^2}{|\boldsymbol{k}|^2}\,k_z\\ -\omega_z-{Ro^{{-}1}} & -\dfrac{{\rm d}}{{\rm d}t}\log{\dfrac{|\boldsymbol{k}_\perp|}{|\boldsymbol{k}|}} & 0\\ -\dfrac{1}{F_h^2 k_z} & 0 & -\dfrac{{\rm d}}{{\rm d}t}\log{\dfrac{|\boldsymbol{k}_\perp|}{|\boldsymbol{k}|}} \end{array}\right)\left( \begin{array}{@{}c@{}} p\\q\\s\end{array}\right). \end{equation}

2.2. Instability condition and estimate of growth rates

Here, we summarize the condition and the estimate of the growth rate for each instability obtained in the framework of local stability analysis; see Hattori & Hirota (Reference Hattori and Hirota2023) for the details. It is pointed out that our aim here is to give compact and useful expressions for the instability condition and the growth rate for each instability, which are not always rigorous but allow us to interpret the results in § 4 without difficulties. In the following, $C_{PH}$, $C_{SH}$, $C_{RH}$, $C_{C}$ and $C_{E}$ are $O(1)$ coefficients that depend on the parameters in general; the actual dependence will be checked numerically in § 4.1 (figure 8).

1. (i) Pure-hyperbolic instability. The pure-hyperbolic (PH) instability is due to stretching near the hyperbolic stagnation points. It occurs when

   (2.19)\begin{equation} |{Ro^{{-}1}}| < \varepsilon_h. \end{equation}
   The growth rate is estimated as
   (2.20)\begin{equation} \sigma = C_{PH} \left(\varepsilon_h^2 -Ro^{{-}2}\right)^{1/2}. \end{equation}

2. (ii) Strato-hyperbolic instability. The strato-hyperbolic (SH) instability is a variant of the pure-hyperbolic instability under stratification effects. It occurs when the exponential growth near the hyperbolic stagnation points is connected with phase shift due to the inertia-gravity waves in favour of exponential growth. It occurs when

   (2.21a,b)\begin{equation} |{Ro^{{-}1}}| < \varepsilon_h \quad\mbox{and} \quad {F_h^{{-}1}} \gtrsim \omega_{max}/2. \end{equation}
   The growth rate is estimated as
   (2.22)\begin{equation} \sigma = C_{SH} \left(\varepsilon_h^2 -Ro^{{-}2}\right)^{1/2}. \end{equation}

3. (iii) Rotational-hyperbolic instability. The rotational-hyperbolic (RH) instability is also a variant of the pure-hyperbolic instability. It occurs when

   (2.23)\begin{equation} |{Ro^{{-}1}}| \gtrsim \varepsilon_h. \end{equation}
   The growth rate is estimated as
   (2.24)\begin{equation} \sigma = C_{RH} \varepsilon_h. \end{equation}

4. (iv) Centrifugal instability. The essential mechanism of the centrifugal (C) instability was given by Rayleigh (Reference Rayleigh1917). It occurs when positive energy is released by interchanging fluid elements under conservation of angular momentum. Assuming a monotonically decreasing vorticity distribution, the instability condition turns out to be

   (2.25)\begin{equation} {-}\omega_{max} < {Ro^{{-}1}} < 0. \end{equation}
   The growth rate is estimated as
   (2.26)\begin{equation} \sigma = C_{C} \omega_{{max}}. \end{equation}

5. (v) Elliptic instability. The elliptic (E) instability is caused by resonance due to strain at elliptic stagnation points. It occurs when

   (2.27a,b)\begin{equation} {F_h^{{-}1}}<\tfrac{1}{2}\omega_{max} \quad\mbox{and} \quad {Ro^{{-}1}} <{-}\tfrac{3}{2}\omega_{max} \mbox{or} \ {Ro^{{-}1}} >{-}\tfrac{1}{2}\omega_{max}, \end{equation}
   or
   (2.28a,b)\begin{equation} {F_h^{{-}1}}>\tfrac{1}{2}\omega_{max} \quad \mbox{and} \quad {-\tfrac{3}{2}}\omega_{max}<{Ro^{{-}1}} \lesssim 0. \end{equation}
   The growth rate is estimated as
   (2.29)\begin{equation} \sigma = C_{E} \varepsilon_e. \end{equation}
   The coefficient $C_{E}$ can be expressed as a function of $F_h$ and $Ro$ at the elliptic stagnation point (Leblanc Reference Leblanc2000).

The values of the strain rates $\varepsilon _h$ and $\varepsilon _e$ at the hyperbolic and elliptic stagnation points, respectively, and the maximum magnitude of vorticity $\omega _{max}$, which appear in the expressions above, are listed in table 1, which includes those for the 2-D Taylor–Green vortices considered in Hattori & Hirota (Reference Hattori and Hirota2023) for comparison purposes. It is pointed out that $\omega _{max}$ and $\varepsilon _e$ of the Stuart vortices are larger than values for the 2-D Taylor–Green vortices, while $\varepsilon _h$ is not much different between the two base flows.

Table 1.

Strain rates at hyperbolic and elliptic stagnation points, and maximum vorticity of the Stuart vortices considered in the present paper. The values for the 2-D Taylor–Green vortices studied in Hattori & Hirota (Reference Hattori and Hirota2023) (HH2023) are also included for comparison.

3.1. Local stability analysis

The numerical method for local stability analysis is same as that in Hattori & Hirota (Reference Hattori and Hirota2023). Equations (2.12)–(2.15) were integrated by the fourth-order Runge–Kutta method. We consider periodic orbits of fluid particles throughout this paper. 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 (Lifschitz & Hameiri Reference Lifschitz and Hameiri1993; Hattori & Fukumoto Reference Hattori and Fukumoto2003):

(3.1)\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 that appear in (2.14) are also time-periodic:

(3.2)\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}$ that coincides with that of the particle motion $\boldsymbol {X}$ (Lifschitz & Hameiri Reference Lifschitz and Hameiri1993; Hattori & Fukumoto Reference Hattori and Fukumoto2003). Our task is to calculate the eigenvalues $\{ \mu _i \}$ of $\boldsymbol{\mathsf{F}}(T)$, which determines the growth rate as

(3.3)\begin{equation} \sigma_i=\frac{\log{|\mu_i|}}{T}. \end{equation}

The initial conditions should be specified to have particular solutions for a given set of the Rossby number $Ro$ and the Froude number $F_h$. One parameter, which is denoted by $\beta$ in the following subsections, is required for $\boldsymbol {X}(0)$ to identify a streamline in a 2-D flow. We set

(3.4)\begin{equation} \boldsymbol{X}(0)=\left({\rm \pi}, \beta y_c, 0\right)^{\rm T},\end{equation}

where $y_c=\cosh ^{-1}(1+2\sqrt {1-C^{-2}})$ is the maximum of $y$ on the separatrix. The elliptic stagnation point corresponds to $\beta =0$, while $\beta =1$ corresponds to the separatrix. The open streamlines corresponding to $\beta >1$ are also considered.

Another parameter is required for $\boldsymbol {k}(0)$ to specify the direction of the wavevector that satisfies (3.1); we take the angle between $\boldsymbol {e}_z$ and $\boldsymbol {k}(0)$, which is denoted by $\theta _0$. It should be pointed out that the magnitude of $\boldsymbol {k}(0)$ is arbitrary since the right-hand side of (2.14) depends only on the direction of $\boldsymbol {k}$ and is independent of the magnitude after taking the short-wave limit. For the amplitudes $\boldsymbol {a}(0)$ and $r(0)$, three independent initial conditions satisfying the incompressibility condition $\boldsymbol {a}(0)\boldsymbol {\cdot } \boldsymbol {k}(0)=0$ are considered; the results do not depend on the choice of the initial conditions since the space spanned by the three initial conditions is common. As a result, we obtain the largest growth rate $\sigma$ as a function of $\beta$, $\theta _0$, $Ro$ and $F_h$, namely $\sigma =\sigma (\beta, \theta _0, Ro, F_h)$.

3.2. Modal stability analysis

In the modal stability analysis, (2.6)–(2.8) were solved numerically by the Fourier spectral method (Peyret Reference Peyret2010), assuming periodic boundary conditions in all three directions as in Hattori et al. (Reference Hattori, Suzuki, Hirota and Khandelwal2021) and Hattori & Hirota (Reference Hattori and Hirota2023). Since the Stuart vortices are periodic in $x$ but not in $y$, the vortices are placed at $y=nL_y$, and those with the opposite-signed vorticity are placed at $y=(n+1/2)L_y$, for $n=0, \pm 1, \pm 2, \ldots\,$, to make the base flow periodic in $y$. The spatial period $L_y$ is fixed at $L_y=4L_0=8{\rm \pi}$, which is large enough to make the effects of the periodic boundary condition in the $y$ direction negligible (Hattori et al. Reference Hattori, Suzuki, Hirota and Khandelwal2021). 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

(3.5)\begin{equation} \boldsymbol{u}' = \exp({{\rm{i}} [\mu (x/L_x) +k_z z]}) \sum_{k_x={-}K_x}^{K_x}\sum_{k_y={-}K_y}^{K_y} \tilde{\boldsymbol{u}}_{k_x,k_y} \exp({{\rm i} [k_x (x/L_x)+k_y (\kern0.7pt y/L_y)]}), \end{equation}

with similar expressions for $p'$ and $\rho '$. In the above equation, we have included the Floquet exponent ${\rm {i}}\mu$ to consider Floquet modes in general. The number of the Fourier modes is $1024 \times 4096$, as in Hattori et al. (Reference Hattori, Suzuki, Hirota and Khandelwal2021).

The growth rate and frequency were obtained by the method of Krylov subspaces (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; Hattori et al. Reference Hattori, Suzuki, Hirota and Khandelwal2021; Hattori & Hirota Reference Hattori and Hirota2023). Starting from randomized initial conditions, (2.6)–(2.8) were integrated for a certain long time. Intermediate states $\{(\boldsymbol {u}'(T_0), \rho '(T_0))$, $(\boldsymbol {u}'(T_0\,{+}\,\Delta T), \rho '(T_0\,{+}\,\Delta T)), \ldots, (\boldsymbol {u}'(T_0+(N_K\,{-}\,1)\,\Delta T), \rho '(T_0+(N_K-1)\,\Delta T))\}$ were used as generators of the Krylov subspace. Then the eigenvalues and the eigenmodes were obtained in the $N_K$-dimensional Krylov subspace. In this method, the error of an eigenvalue $\lambda$ of a linear operator ${\boldsymbol{\mathsf{L}}}$ can be evaluated by

(3.6)\begin{equation} \epsilon = \frac{\| {\boldsymbol{\mathsf{L}}}\boldsymbol{v}-\lambda\boldsymbol{v} \|}{\| \boldsymbol{v} \|}, \end{equation}

where $\boldsymbol {v}$ is the corresponding approximate eigenvector. The error $\epsilon$ depends on the initial time of the data $T_0$, the interval between the data $\Delta T$, and the dimension of the Krylov subspace $N_K$. 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; $N_K$ was $5$, $10$ or $20$, and $T_0=92, 112, \ldots, 192$, while $\Delta T$ was fixed at $2$. Typically, the error of the eigenvalue is $\epsilon =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 $\epsilon \geqslant 10^{-3}$.

It turned out that the growth rates of the cyclonic case are sometimes difficult to obtain because they are smaller than those of the anticyclonic case. Therefore, we applied a filter that damps the anticyclonic modes and keeps the cyclonic modes unchanged, to obtain the growth rates of the cyclonic case.

3.3. Realizability as a mode

As we will see later, in § 4, the instabilities that are found by local stability analysis are not always found in modal stability analysis at finite Reynolds numbers since high-wavenumber modes are damped by viscous damping. In this case, the corresponding region of the instability in the $(\beta, \theta _0)$ plane is often thin so that it is difficult to construct an unstable mode. We use the realizability introduced in Hattori & Hirota (Reference Hattori and Hirota2023),

(3.7)\begin{equation} \mathcal{R} = \int_{S} \sigma \sin\theta_0 \,{\rm d}\beta \,{\rm d}\theta_0, \end{equation}

where $S$ is the region of an instability in the $(\beta, \theta _0)$ plane, to quantify the realizability as a mode of each instability. The idea behind the above definition is that the eigenmode corresponding to a wider unstable region in the $(\beta, \theta _0)$ plane is less affected by viscous damping because the corresponding mode also has a large spatial width and thereby the radial wavenumber of the mode is small.