2.1 Equations of motion

The Boussinesq approximation of the equations of motion for 3-D rotating, stratified, inviscid flows in the Cartesian coordinates $(x,y,z)$ , as observed in a frame rotating with angular velocity $(f/2)\hat{\boldsymbol{z}}$ , is (Vallis Reference Vallis2006)

(2.1a,b ) $$\begin{eqnarray}\displaystyle \frac{\text{D}\boldsymbol{v}}{\text{D}t}=-\frac{\unicode[STIX]{x1D735}p}{\unicode[STIX]{x1D70C}_{o}}+\boldsymbol{v}\times f\hat{\boldsymbol{z}}+b\hat{\boldsymbol{z}},\quad \frac{\text{D}b}{\text{D}t}=-\bar{N}^{2}v_{z},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & & \displaystyle\end{eqnarray}$$

where the operator $\text{D}/\text{D}t\equiv \unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the material derivative, $t$ denotes time, $\boldsymbol{v}=(v_{x},v_{y},v_{z})$ is the 3-D velocity vector, $f$ is the Coriolis frequency (constant in our study) and $g$ is the acceleration of gravity. The total pressure and the total density of the fluid are $p_{tot}\equiv \bar{p}(z)+p(x,y,z,t)$ and $\unicode[STIX]{x1D70C}_{tot}\equiv \bar{\unicode[STIX]{x1D70C}}(z)+\unicode[STIX]{x1D70C}(x,y,z,t)$ , where $\bar{\unicode[STIX]{x1D70C}}(z=0)=\unicode[STIX]{x1D70C}_{o}$ . We define the buoyancy as $b(x,y,z,t)\equiv -g\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{o}$ . Quantities with a bar are properties of the equilibrium background flow (i.e. far from the vortex where $\boldsymbol{v}\rightarrow 0$ , $b\rightarrow 0$ , $\unicode[STIX]{x1D70C}\rightarrow 0$ and $p\rightarrow 0$ ). The background pressure $\bar{p}$ and density $\bar{\unicode[STIX]{x1D70C}}$ are in hydrostatic balance $\text{d}\bar{p}/\text{d}z=-\bar{\unicode[STIX]{x1D70C}}g$ . The background Brunt–Väisälä frequency $\bar{N}\equiv \sqrt{-(g/\unicode[STIX]{x1D70C}_{o})(\text{d}\bar{\unicode[STIX]{x1D70C}}/\text{d}z)}$ is assumed to be constant, so that $\bar{\unicode[STIX]{x1D70C}}(z)=\unicode[STIX]{x1D70C}_{o}(1-\bar{N}^{2}z/g)$ .

In the above equations, we have ignored viscosity in the momentum equations and diffusion in the density equation, which are reasonable approximations for atmospheric and oceanic flows. Furthermore, we have dropped the planetary centrifugal term from the momentum equations, assuming that the rotational Froude number $f^{2}d/g$ is small (Barcilon & Pedlosky Reference Barcilon and Pedlosky1967), where $d$ is the distance between the centre of the vortex and the planetary rotation axis.

2.2 Numerical method

A pseudo-spectral initial-value solver is developed to solve (2.1) in a triply periodic domain with 256 or 512 Fourier modes in each direction. In numerical simulations of strongly rotating stratified flows, resolving the fast inertia–gravity waves can substantially limit the size of the time step $\unicode[STIX]{x0394}t$ and thus increase the computational cost. Here we use the semi-analytic method developed by Barranco & Marcus (Reference Barranco and Marcus2006) for rotating stratified flows, which enables us to accurately and efficiently deal with large $f\unicode[STIX]{x0394}t$ and $\bar{N}\unicode[STIX]{x0394}t$ .

A vortex in the middle of a periodic domain interacts with its periodic images. To minimize this interaction and its potential impact on the stability of the vortex (and to simulate having an unbounded flow) the computational domain size is chosen to be large compared to the vortex size: the domain size in the $x$ and $y$ directions, i.e. the values of $L_{x}$ and $L_{y}$ are 7.5 (or more often 15) times larger than the initial vortex diameter $(2L)$ , and, similarly, the domain size in the $z$ direction $L_{z}$ is 7.5 (or more often 15) times larger than the initial vortex height $(2H)$ . There are two reasons for sometimes making the domain size very large. First, we wanted to ensure that the periodic boundary conditions had no perceptible effects on the flow dynamics; secondly, in the follow-up paper to this one (see our discussion § 7) we show that unstable vortices often fragment with pieces of the initial vortices becoming widely separated, so that the calculations require a large domain. To help simulate an unbounded flow, we also added a cylindrical sponge layer near the boundaries of the computational domain (see appendix A). The sponge layer, implemented as Rayleigh drag and Newtonian cooling in (2.1), damps $\boldsymbol{v}$ and $\unicode[STIX]{x1D70C}$ outside a cylindrical surface of diameter $24L$ and height $24H$ (for the large domain calculations) or $12L$ and height $12H$ (for the small domain calculations) around the centre of the domain. Another advantage of adding the sponge layer is that it damps the reflection of the outgoing inertia–gravity waves, and occasional detached filaments back to the domain at the periodic boundaries. One more advantage of the axisymmetric sponge layer is that we find that it prevents the (non-axisymmetric) periodic boundary conditions in $x$ and $y$ from adding any significant non-axisymmetric perturbations to the initial vortices. The latter is important when computing the stability of the vortices. One way of determining if the domain size is too small is to compute the ratio of the magnitude of each component of the velocity and density of a numerically computed eigenmode at a damped location just inside the sponge layer to the maximum value of that component over the entire domain. With the domain sizes presented here, that ratio is always of order $10^{-4}$ (or smaller), but the ratio increases to values with orders as large as $10^{-2}$ when the computational domain is reduced to $(10L)\times (10L)\times (10H)$ and a sponge layer with diameter of $8L$ and height  $8H$ .

Hyperviscosities and hyperdiffusivities are added to our otherwise inviscid and non-diffusive calculations to stabilize the code. See Barranco & Marcus (Reference Barranco and Marcus2006) for more details.

2.3 Initial equilibria: Gaussian vortices

In this study we focus on 3-D axisymmetric baroclinic vortices that are initially in horizontal cyclo-geostrophic balance and vertical hydrostatic balance, and hence they are in gradient–wind balance (Vallis Reference Vallis2006). The initial vortex is centred at $r=0$ and $z=0$ , where $r$ denotes the radial coordinate. A widely used model for oceanic and laboratory vortices is that of an axisymmetric vortex with a Gaussian pressure distribution (e.g. McWilliams Reference McWilliams1985)

(2.2) $$\begin{eqnarray}\displaystyle p=p_{o}\unicode[STIX]{x1D712}(r,z), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D712}(r,z)\equiv \exp [-(r/L)^{2}-(z/H)^{2}]$ . Using (2.2) and the definitions presented in § 2.1, an exact, steady, axisymmetric equilibrium solution to the Boussinesq equations in (2.1) is the vortex

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle v_{\unicode[STIX]{x1D719}}(r,z)={\displaystyle \frac{fr}{2}}\left(-1+\sqrt{1-(8p_{o}\unicode[STIX]{x1D712}(r,z))/(\unicode[STIX]{x1D70C}_{o}f^{2}L^{2})}\right),\quad v_{r}=v_{z}=0, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle b(r,z)=-{\displaystyle \frac{2p_{o}z}{\unicode[STIX]{x1D70C}_{o}H^{2}}}\unicode[STIX]{x1D712}(r,z), & \displaystyle\end{eqnarray}$$

where the cylindrical coordinate is used for convenience ( $v_{\unicode[STIX]{x1D719}}$ is the azimuthal velocity). For any vortex, whether or not it is Gaussian, we shall define a quantity written with a subscript ‘ $c$ ’ to mean that the quantity is to be evaluated at the vortex centre, so $N_{c}$ is the Brunt–Väisälä frequency at the centre of a vortex, or $N_{c}^{2}\equiv \bar{N}^{2}+(\unicode[STIX]{x2202}b/\unicode[STIX]{x2202}z)_{c}$ . For the Gaussian vortex described by (2.3)–(2.4),

(2.5) $$\begin{eqnarray}\displaystyle N_{c}^{2}=\bar{N}^{2}-2p_{o}/(\unicode[STIX]{x1D70C}_{o}H^{2}). & & \displaystyle\end{eqnarray}$$

As discussed in the next section, for some values of $p_{o}$ , $N_{c}^{2}<0$ , which means that the density distribution is locally unstable at the vortex centre with heavy fluid over light fluid (i.e. statically unstable). It is convenient to define the Rossby number $Ro$ , which by definition has $Ro>0$ for a cyclone and $Ro<0$ for an anticyclone, in terms of the maximum (or minimum) value of a vortex’s vertical vorticity $\unicode[STIX]{x1D714}_{E}$ , such that $Ro\equiv \unicode[STIX]{x1D714}_{E}/2f$ . For the Gaussian vortices described above $\unicode[STIX]{x1D714}_{E}=\unicode[STIX]{x1D714}_{c}$ , and

(2.6) $$\begin{eqnarray}Ro=\unicode[STIX]{x1D714}(r=0,z=0)/2f=-1/2+\sqrt{1/4-2p_{o}/(\unicode[STIX]{x1D70C}_{o}f^{2}L^{2})}.\end{eqnarray}$$

Note that the Gaussian vortex has an aspect ratio of

(2.7) $$\begin{eqnarray}\left(\frac{H}{L}\right)^{2}=\frac{-Ro(1+Ro)f^{2}}{\bar{N}^{2}[1-(N_{c}/\bar{N})^{2}]},\end{eqnarray}$$

in accord with the universal scaling law of Aubert et al. (Reference Aubert, Le Bars, Le Gal and Marcus2012) and Hassanzadeh et al. (Reference Hassanzadeh, Marcus and Le Gal2012), which is valid for all vortices that are in cyclo-geostrophic and hydrostatic balance. This can be seen by simply replacing $2p_{o}/\unicode[STIX]{x1D70C}_{o}$ in (2.6) with $H^{2}(\bar{N}^{2}-N_{c}^{2})$ using (2.5), and then solving for $H/L$ .

The three independent dimensional parameters in the governing equations (2.1) are $f$ , $\bar{N}$ and $\unicode[STIX]{x1D70C}_{o}$ . The sizes of the computational domain $L_{x}\times L_{y}\times L_{z}$ have no effect (on the dimensional analysis), due to the fact that the cylindrical sponge layer is far from the vortices and that the net circulations of the flow are zero, which makes the velocity due to the vortices fall off exponentially fast and be effectively zero at the sponge layer. (See the definition of shielded below and in the appendices.) The equilibrium Gaussian vortices in (2.2)–(2.4) introduce three additional dimensional parameters $H$ , $L$ and $p_{o}$ . Thus, there are three independent, dimensionless parameters that describe the dynamics of Gaussian vortices. The choice of these parameters is not unique, but in this paper we choose $Ro$ , $f/\bar{N}$ , and

(2.8) $$\begin{eqnarray}Bu\equiv \left(\frac{\bar{N}}{f}\frac{H}{L}\right)^{2}=(L_{r}/L)^{2},\end{eqnarray}$$

where the latter is the Burger number and $L_{r}\equiv H\bar{N}/f$ is the deformation radius. It should be noted that whether the vortices studied here are big or small depends on the inverse of their Burger number, which is the square of the vortex radius over $L_{r}$ . Big vortices have small $Bu$ and vice versa.

Figure 1.

Vertical vorticity $\unicode[STIX]{x1D714}(r,z)$ and potential vorticity $Q(r,z)$ for Gaussian cyclones defined by (2.2)–(2.4). Panel (a) shows $\unicode[STIX]{x1D714}$ at $z=0$ as a function of $r$ as a solid curve (blue, in colour) for $Ro=0.2$ , $Bu=0.1$ and $f/\bar{N}=0.1$ . The thin and thick dashed vertical lines show the boundaries of the core and the shield (see appendix B). The solid vertical lines at large radii show, with increasing thickness from left to right, where the boundary damping function $f_{bd}$ (see appendix A) reaches values of 0.01, 0.5 and 0.99, respectively. Panel (b) shows $\unicode[STIX]{x1D714}(r,z)$ in units of $2f$ in the $r$ – $z$ plane of the vortex in (a). Dashed lines indicate the boundaries of the core and shield. (c–e) Show the potential vorticity $Q(r,z)$ for a Gaussian vortex with $Ro=0.2$ and $f/\bar{N}=0.1$ , for $Bu=0.1$ , $Bu=1.0$ and $Bu=2.0$ , respectively. For larger values of $Bu$ , the distributions of $Q$ and $\unicode[STIX]{x1D714}$ are similar.

The Gaussian vortices defined in the above model are shielded. Here we define a shielded flow as one in which the circulation computed with the $z$ -component of the vorticity over the entire $(x,y)$ -plane for any fixed value of $z$ is zero. In addition, the circulation computed with the $x$ -component of the vorticity over the $(y,z)$ -plane for any fixed value of $x$ is zero; and the circulation computed with the $y$ -component of the vorticity over the $(x,z)$ -plane for any fixed value of $y$ is zero. (Note that figure 1(b) does not violate our definition of shielded because the figure shows the vertical component of the vorticity in an $x$ – $z$ plane, not an $x$ – $y$ plane.) Our governing equations and boundary conditions show that if the initial flow is shielded, then the flow is shielded for all time. In practical terms, a shielded isolated vortex is one in which the central core of the vortex is surrounded, or partially surrounded, by a region (shield) of opposite vorticity and that the circulation quickly vanishes outside the shield. For an arbitrary (i.e. not necessarily Gaussian) cyclonic vortex, the core of a cyclone is a contiguous region at and near the vortex centre where the vertical component of its vorticity $\unicode[STIX]{x1D714}$ is greater than or equal to zero. The shield is a region around the core (usually looking like a shell or annular ring) located not too far from the core, where $\unicode[STIX]{x1D714}<0$ . The precise definitions that we use for core and shield are in appendix B. The core and shield of an example Gaussian vortex are illustrated in figure 1(a,b). The definitions of the core and shield of an anticyclone are analogous to those of the cyclone. For Gaussian vortices and many other types of shielded cyclones, outside the shield the amplitude of the vorticity decays exponentially with the radial distance $r$ (or $r^{p}$ with $p\geqslant 2$ ) from the vortex centre. In our calculations, the circulation due to the vertical component of the vorticity $\int \unicode[STIX]{x1D714}(x,y,z)\,\text{d}x\,\text{d}y$ (where the integral is over the entire $x$ – $y$ computational domain) at each value of $z$ must remain zero due to the periodic boundary conditions.

Commonly, in the studies of oceanic and atmospheric vortices, potential vorticity (PV) is used to describe the vortices, instead of vertical vorticity, due to its conservation property (Hoskins et al. Reference Hoskins, Mcintyre and Robertson1985; Morel & McWilliams Reference Morel and McWilliams1997). Ertel’s PV in figure 1 is defined (Ertel Reference Ertel1942) as

(2.9) $$\begin{eqnarray}Q\equiv [\unicode[STIX]{x1D74E}+f\hat{\boldsymbol{z}}]\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D735}b+\bar{N}^{2}\hat{\boldsymbol{z}}}{f\bar{N}^{2}}\right)-1,\end{eqnarray}$$

where $\unicode[STIX]{x1D74E}\equiv \unicode[STIX]{x1D735}\times \boldsymbol{v}$ is the vorticity vector as observed in the rotating frame. To provide a better sense about the PV structure of the vortices studied here, $Q(r,z)$ for a Gaussian vortex with $Ro=0.2$ and three values of $Bu=0.1,1$ and $2$ are depicted in figure 1(c–e), showing that the PV structure can significantly change with $Bu$ (see Morel & McWilliams (Reference Morel and McWilliams1997) for a discussion of potential vorticity of Gaussian vortices). Our purpose for showing the PV of Gaussian vortices is to allow the reader the ability to make comparisons of the vortex model with what is used in some other stability studies such as Tsang & Dritschel (Reference Tsang and Dritschel2015) who model the initial vortex with uniform patches of PV.

Finally it should be noted that there is a restriction on the equilibrium of anticyclones in the Gaussian model (2.2)–(2.4); there is no equilibrium for anticyclones for $Ro<-0.5$ . This is because (2.3) and (2.6) show that $v_{\unicode[STIX]{x1D719}}$ does not have a real solution for $Ro<-0.5$ , as noted, for example, by McWilliams (Reference McWilliams1985) and Olson (Reference Olson1991).

2.4 Eigenmodes

The symmetries of the governing equations in (2.1) linearized around the equilibrium vortex (2.2)–(2.4) are presented in dimensionless form in § 5 in (5.10)–(5.14). These equations and their boundary conditions show that the eigenfunctions are either symmetric or antisymmetric with respect to the $z=0$ horizontal plane and have an $m$ -fold azimuthal symmetry about the $z$ -axis. We use the labels S $m$ or A $m$ for each eigenmode, to identify it as symmetric (or antisymmetric) with respect to the $z=0$ horizontal plane and with $m$ -fold symmetry.

The complex eigenvalues $\unicode[STIX]{x1D706}$ and eigenfunctions are of the form

(2.10) $$\begin{eqnarray}\text{e}^{\unicode[STIX]{x1D706}t}\boldsymbol{g}_{eig}(r,z)\text{e}^{\text{i}m\unicode[STIX]{x1D719}}=\text{e}^{\unicode[STIX]{x1D70E}t}\boldsymbol{g}_{eig}(r,z)\text{e}^{\text{i}m(\unicode[STIX]{x1D719}-ct)},\end{eqnarray}$$

where the eigenvector has three velocity components, a density component and a pressure component:

(2.11) $$\begin{eqnarray}\boldsymbol{g}_{eig}(r,z)\equiv [v_{r,eig},v_{\unicode[STIX]{x1D719},eig},v_{z,eig},\unicode[STIX]{x1D70C}_{eig},p_{eig}].\end{eqnarray}$$

The three velocity components are with respect to cylindrical coordinates, where $m$ is the integer azimuthal wavenumber, $\unicode[STIX]{x1D70E}$ is a real growth (or decay) rate and $c$ is a real azimuthal phase speed. By taking the complex conjugate of the linearized equation, we can show that if $\unicode[STIX]{x1D706}$ is an eigenvalue with eigenfunction given by (2.11), then $\unicode[STIX]{x1D706}^{\dagger }$ is also an eigenvalue with eigenfunction $\boldsymbol{g}_{eig}^{\dagger }(r,z)\equiv [v_{r,eig}^{\dagger },v_{\unicode[STIX]{x1D719},eig}^{\dagger },v_{z,eig}^{\dagger },\unicode[STIX]{x1D70C}_{eig}^{\dagger },p_{eig}^{\dagger }]$ , with $m$ replaced by $-m$ , $c$ unchanged and where the superscript $\dagger$ denotes complex conjugate. Or in other words, the eigenvalues $\unicode[STIX]{x1D706}$ when plotted in the complex plane are symmetric with respect to the real axis. Because the equations are non-dissipative, replacing $t$ with $-t$ in the linearized equations shows that if $\unicode[STIX]{x1D706}\equiv \unicode[STIX]{x1D70E}-\text{i}mc$ is an eigenvalue with eigenfunction given by (2.11), then $\unicode[STIX]{x1D706}^{\prime }\equiv -\unicode[STIX]{x1D70E}-\text{i}m^{\prime }c^{\prime }$ is an eigenvalue that corresponds to $\boldsymbol{g}_{eig}^{\prime }(r,z)\equiv [v_{r,eig},-v_{\unicode[STIX]{x1D719},eig},v_{z,eig},-\unicode[STIX]{x1D70C}_{eig},-p_{eig}]$ , with $m^{\prime }=-m$ and $c^{\prime }=c$ . Or in other words, the eigenvalues $\unicode[STIX]{x1D706}$ when plotted in the complex plane are symmetric with respect to the imaginary axis and for each eigenfunction with a positive growth rate there is one with a negative growth rate and vice versa. The flow can never be linearly stable with all of its eigenmodes having decay rates. The flow can either be unstable or be neutrally stable with all of its eigenmodes on the imaginary axis with $\unicode[STIX]{x1D70E}=0$ . For the Gaussian vortices, the two symmetries of the linearized equations combine and therefore the eigenvalues appear as quartets of the form $\pm a\pm \text{i}b$ , with all four possible combinations of the signs, and where $a$ and $b$ are real functions of $m$ and of the parameters of the unperturbed vortex $Ro$ , $Bu$ and $f/\bar{N}$ . For Hamiltonian systems (Ozorio de Almeida Reference Ozorio de Almeida1988), it can be shown that the quartet of eigenvalues is of a more specialized form:

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\pm \sqrt{A}\pm \text{i}B,\end{eqnarray}$$

with all four possible combinations of the signs, and where $A$ and $B$ are real functions of the control parameters of the system. For many non-dissipative flows, e.g. unidirectional shears flows with vortex sheets and/or vortex layers made up of piecewise-constant vorticity (Drazin & Reid Reference Drazin and Reid2004), it can be shown that the quartets of the eigenvalues are of the form of (2.12). Consider a system with eigenvalue quartets such as those in (2.12). When $A>0$ , the eigenvalues in the quartets are symmetric about the real and imaginary axes and each quartet has two unstable and two stable eigenmodes. If a control parameter changes such that $A$ decreases, then eigenvalues symmetrically approach the imaginary axis and collide when $A=0$ . For that parameter value, there are two pairs of degenerate, neutrally stable eigenmodes with all four eigenvalues on the imaginary axis. If the control parameter is further changed such that $A$ continues to decrease and becomes negative, then the eigenvalues are no longer degenerate, but they remain on the imaginary axis and all four eigenmodes remain neutrally stable, regardless of how negative $A$ becomes. Although we cannot prove that the eigenvalue quartets of the linear eigenmodes of the Gaussian vortex have the form of (2.12), all of our numerical simulations are consistent with (2.12). (See § 3.)

Note that although we are studying the stability of axisymmetric vortices, we solve (2.1) in the Cartesian coordinates rather than in the cylindrical coordinates. A numerical solver in the Cartesian coordinates avoids the difficulties of handling the singularity at the origin ( $r=0$ ), which requires using special polynomial basis functions (Matsushima & Marcus Reference Matsushima and Marcus1995). However, our main reason for using Cartesian coordinates is that future studies can include background shear flows, so that the stability of vortices in planetary atmospheres and protoplanetary disks can be examined, as discussed in the Introduction. To minimize the effect of the square computational domain, we have used a circular sponge layer as described in § 2.2. In order to find the eigenmodes with various classes of azimuthal (and vertical) symmetry in the Cartesian coordinates, we use our initial-value solver as an eigenvector/eigenvalue solver and additionally use a spatial symmetrizer (see appendix C for details). Using the spatial symmetrizer, the eigenmodes can be restricted to be symmetric or antisymmetric in the vertical direction, while in the azimuthal direction we can enforce one of the following classes of symmetry: $m$ odd; $m$ even not divisible by $4$ ; and $m$ even and divisible by $4$ . We use these specific symmetry groups to apply the azimuthal symmetry directly in the Cartesian coordinates, which greatly speeds up the convergence of the calculations, and also avoids introducing additional errors due to transformation between Cartesian and cylindrical coordinates (see appendix C).

4.1 Spatial symmetries and growth rates of the eigenmodes

The parameter map of stability in the $Ro{-}Bu$ space is shown in figure 8(a). Gaussian anticyclones do not exist with $Ro<-0.5$ (see § 2.3). The region to the lower left of the thick dashed black curve corresponds to equilibrium Gaussian vortices for which $N_{c}^{2}<0$ (or $Bu<-Ro(1+Ro)$ according to (2.7)). These vortices are not unphysical, but near their cores they have heavy fluid above light fluid (i.e. $\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}z>0$ at the vortex centre).

As shown in figure 8(a), the most unstable eigenmodes (i.e. those with the largest growth rates) of the vortices generally have either S2 or A1 symmetries. A few points in the figure correspond to vortices for which the fastest-growing eigenmode is A2, A3 or A4. We found that no vortex had a fastest-growing eigenmode with a symmetry different from those just listed. To our surprise, only four out of the 130 vortices that we examined were neutrally stable. All the neutrally stable vortices were cyclones with $0.02\lesssim Ro\lesssim 0.05$ and $0.8\lesssim Bu\lesssim 1$ . The neutrally stable eigenmodes are denoted in figure 8(a) as solid circles in the region circumscribed by a small rectangle. The rectangle is to ‘guide the eye’ and is used to denote the approximate boundary of the region of neutral stability. Computing the actual boundary between the regions where vortices are all neutrally stable and where they are unstable would be expensive and rather pointless given how small the neutrally stable region is. Anticyclones have linear growth rates that are slow and would not destroy a vortex in less than 50 vortex turnaround times if $0.5\lesssim Bu\lesssim 1.3$ . For nearly geostrophic cyclones with $|Ro|<0.05$ , linear growth rates are slow and would not destroy a vortex in less than 50 vortex turnaround times if $0.7\lesssim Bu\lesssim 1.2$ . As $Ro$ increases, the growth rates of large diameter cyclones (i.e. with $Bu\lesssim 1.05$ or $L\gtrsim 0.98L_{r}$ ) becomes faster.

Figure 8.

Parameter map of the stability of Gaussian vortices in the $Ro{-}Bu$ space. No equilibrium Gaussian vortices exist with $Ro<-0.5$ . The thick black dashed line in the lower left corner indicates the locus over which $N_{c}^{2}=0$ , i.e. $Bu=-Ro(1+Ro)$ , with $N_{c}^{2}<0$ for vortices with smaller $Ro$ or $Bu$ . Panel (a): the thick solid (blue, in colour) and thick dotted (red, in colour) lines indicate the iso-surface where $\unicode[STIX]{x1D70E}$ of the fastest-growing eigenmode is $0.02\unicode[STIX]{x1D70F}^{-1}$ . The region bounded by this iso-surface, the thick black dashed curve (but see the caveat in the text describing figure 10) and the bottom of the figure has $\unicode[STIX]{x1D70E}<0.02\unicode[STIX]{x1D70F}^{-1}$ (the iso-contour is to guide the eye and is approximated by interpolating among the growth rates calculated at the locations of the discrete symbols). The symbols denote the spatial symmetry of the fastest-growing eigenmode, with diamonds (blue, in colour) as S2, solid triangles (red, in colour) as A1, squares (green, in colour) as A2, hollow triangles as A3 and hollow circles as A4. Black solid circles correspond to vortices for which the most unstable eigenmodes have growth rates slower than $0.02\unicode[STIX]{x1D70F}^{-1}$ . Panel (b): four iso-contours (approximated as in (a)) of growth rate $\unicode[STIX]{x1D70E}$ of the fastest-growing eigenmode. Each contour consists of one solid curve (blue, in colour) and one or two dotted curves (red, in colour). The fastest-growing eigenmodes along the dotted curves (red, in colour) have A1 symmetry and along the solid curves (blue, in colour) have S2 symmetry. The small rectangular box near $Bu=1$ is to guide the eye and shows the approximate, very small, region where all of the eigenmodes of the cyclones are neutrally stable. The $\unicode[STIX]{x1D70E}$ and the symmetries of the most unstable eigenmodes with $\unicode[STIX]{x1D70E}>0.02\unicode[STIX]{x1D70F}^{-1}$ for vortices with $N_{c}^{2}>0$ in (a) are given in appendix D.

Considering the smallness of the region of neutral stability, clearly, linear stability cannot be used to explain the differences between the numbers of observed cyclones and anticyclones in the oceans or in planetary atmospheres. On the other hand, ocean vortices can survive for more than 50 of their own turnaround times, $\unicode[STIX]{x1D70F}$ . So, one plausible explanation of the cyclonic/anticyclonic asymmetry in the frequency of observation of mesoscale oceanic eddies and of planetary vortices might depend on the differences of the growth rates of the linear instabilities, rather than just the fact that some vortices are not linearly unstable and others are. For example, if there are physical processes (such as turbulence, interactions with other vortices or currents or boundaries) that are likely to destroy a vortex after $50\unicode[STIX]{x1D70F}$ , which is more than ${\sim}1/2$ year for ocean Meddies (McWilliams Reference McWilliams1985; Armi et al. Reference Armi, Hebert, Oakey, Price, Richardson, Rossby and Ruddick1989; Hebert, Oakey & Ruddick Reference Hebert, Oakey and Ruddick1990; Pingree & Le Cann Reference Pingree and Le Cann1993; D’Asaro, Walker & Baker Reference D’Asaro, Walker and Baker1994; Prater & Sanford Reference Prater and Sanford1994; Paillet et al. Reference Paillet, Le Cann, Carton, Morel and Serpette2002), then a vortex need not be neutrally stable to be observed, it needs only have growth rates less than ${\sim}1/50\unicode[STIX]{x1D70F}^{-1}$ . So, it is plausible that the asymmetry between the numbers of observed cyclones and anticyclones depends upon the relative amount of area in $Ro{-}Bu$ parameter space for which the fastest-growing eigenmodes grow slower than ${\sim}1/50\unicode[STIX]{x1D70F}^{-1}$ , or some other critical growth rate. For Gaussian vortices, the region in $Ro{-}Bu$ parameter space where the growth rate of the fastest-growing eigenmode is less than $1/50\unicode[STIX]{x1D70F}^{-1}$ (i.e. the ‘slow growth region’ for linear instability) is the region bounded above by the solid (blue, in colour) and dotted (red, in colour) curves in figure 8(a) and to the lower left by the thick dashed curve. Along the solid curve (blue, in colour), the fastest-growing eigenmode has S2 symmetry, whereas along the dotted curve (red, in colour) it is A1. The solid and dotted curves are drawn to guide the eye, and the vortices corresponding to the black solid circles have $\unicode[STIX]{x1D70E}<1/50\unicode[STIX]{x1D70F}^{-1}$ . In general, for large $Bu$ , the fastest-growing eigenmodes have A1 symmetry, while for smaller $Bu$ , they have S2 symmetry. However, for cyclones with $Bu\lesssim 0.4$ , some of the fastest-growing eigenmodes also have A1 symmetry, or even A2, A3 or A4 symmetry, and the growth rates are often faster than $1\unicode[STIX]{x1D70F}^{-1}$ . There are two regions in the $Ro{-}Bu$ parameter space where the fastest-growing eigenmodes of the cyclones have A1 symmetry. In the region with higher $Bu$ , the growth rate of the fastest-growing modes is smaller than that in the lower $Bu$ region, and, as discussed previously and elaborated on in § 6, the radial structures of the fastest-growing A1 eigenmodes in the large and small $Bu$ regions differ as well.

Of course, our choice of $50\unicode[STIX]{x1D70F}$ to define the ‘slow growth region’ for linear instability is arbitrary, so figure 8(b) shows how the slow growth region changes when we change our choice from $50\unicode[STIX]{x1D70F}$ to $20\unicode[STIX]{x1D70F}$ , $10\unicode[STIX]{x1D70F}$ or $6.67\unicode[STIX]{x1D70F}$ . That is, the two sets of (solid/broken) curves are iso-surfaces in $Ro{-}Bu$ parameter space where $\unicode[STIX]{x1D70E}$ is 0.02, 0.05, 0.1 and 0.15 in units of $\unicode[STIX]{x1D70F}^{-1}$ . For the iso-surface for the growth rate of $0.15\unicode[STIX]{x1D70F}^{-1}$ in figure 8(b), the fastest-growing eigenmode has S2 symmetry for $0.5\lesssim Bu\lesssim 1.8$ , otherwise the fastest-growing eigenmode has A1 symmetry. Note that the iso-surfaces for the growth rates of $0.10\unicode[STIX]{x1D70F}^{-1}$ and $0.15\unicode[STIX]{x1D70F}^{-1}$ are very close to each other for $Bu\lesssim 0.3$ . Most of the fastest-growing eigenmodes in the $Ro{-}Bu$ parameter space shown in figure 8 have $\unicode[STIX]{x1D70E}<0.2\unicode[STIX]{x1D70F}^{-1}$ . However, cyclones in the upper left corner of figure 8 can have $\unicode[STIX]{x1D70E}$ of order one $\unicode[STIX]{x1D70F}^{-1}$ . As shown below, anticyclones to the lower left of the thick dashed curve in the lower left side of figure 8 (with $N_{c}^{2}<0$ ) can have much larger  $\unicode[STIX]{x1D70E}$ .

The growth rates of the three fastest-growing eigenmodes for $Ro=0.05$ as functions of $Bu$ are plotted in figure 9(a) (combining figures 2, 4 and 6) showing that the fastest-growing eigenmode is A1 for $Bu\lesssim 0.2$ ; is S2 for $0.2\lesssim Bu\lesssim 0.8$ ; and is A1 for $1\lesssim Bu\lesssim 2.1$ . However, for $0.8\lesssim Bu\lesssim 1$ , the eigenmodes are all neutrally stable. This region of neutral stability is consistent with the neutrally stable region shown in figure 8. The change in the spatial symmetry from A1 to S2 back to A1 of the fastest-growing eigenmode as $Bu$ increases was discussed in § 3 and it was shown that (i) the family of eigenmodes continues to exist even after the eigenmodes become neutrally stable and (ii) the A1 modes at small and large $Bu$ belong to two different families of eigenmodes. Similar changes in the symmetries of the most unstable mode are observed at $Ro=0.2$ (figure 9 b); however, at $Ro=0.2$ there is not a region where the vortex is neutrally stable to all eigenmodes. Similar to $Ro=0.05$ , the two families of A1 eigenmodes shown in figure 9(b) with triangles and with filled circles are distinct families with different radial structures. How these results, particularly at the small $Ro$ of 0.05, compare with those obtained from analysing the QG equations is discussed in § 4.2.

Figure 9.

Growth rates $\unicode[STIX]{x1D70E}$ (in units of $\unicode[STIX]{x1D70F}^{-1}$ ) of the three fastest-growing modes as functions of $Bu$ for fixed $Ro$ . $f/\bar{N}=0.1$ . Triangles, filled circles and diamonds, respectively, indicate the fastest-growing eigenmodes at low $Bu$ (which have A1 symmetry), the fastest-growing eigenmodes at high $Bu$ (which also have A1 symmetry) and the fastest-growing eigenmodes for intermediate $Bu$ (which have S2 symmetry). The three lines connecting the three sets of symbols are to guide the eye to show the three families of eigenmodes. (a) $Ro=0.05$ ; in this case as $Bu$ increases, the fastest-growing mode changes from A1 to S2; then all modes are linearly neutrally stable; then the fastest-growing mode is A1. (b) $Ro=0.2$ ; the fastest-growing mode changes from A1 to S2 and again to A1 as $Bu$ increases.

The growth rates for the region with statically unstable vortex cores, i.e. with $N_{c}^{2}<0$ , are shown in figure 10. Eigenmodes for this region have A4 symmetry and the growth rates can be as large as ${\sim}100\unicode[STIX]{x1D70F}^{-1}$ . The $\unicode[STIX]{x1D70E}$ as a function of $Bu$ (for fixed $Ro$ ), and as a function of $Ro$ (for fixed $Bu$ ) for vortices with $N_{c}^{2}<0$ are shown in figures 11(a) and 11(b), respectively. In each of the eight panels, the value of the horizontal coordinate axis on the right-hand side of the panel corresponds to a vortex with $N_{c}^{2}=0$ (i.e. a point on the thick dashed curve in figure 8 or in the broken curve in figure 10). The figure shows that $\unicode[STIX]{x1D70E}$ increases rapidly as a function of distance from the $N_{c}^{2}=0$ boundary. Due to this rapid growth in $\unicode[STIX]{x1D70E}$ , for all practical purposes we can consider the thick dashed line at $N_{c}^{2}=0$ to be the left boundary of the region in figure 8(b) in $Ro{-}Bu$ for which $\unicode[STIX]{x1D70E}<0.02\unicode[STIX]{x1D70F}^{-1}$ , and also the boundary for the region $0.02\unicode[STIX]{x1D70F}^{-1}\leqslant \unicode[STIX]{x1D70E}<0.05\unicode[STIX]{x1D70F}^{-1}$ , and for the region $0.05\unicode[STIX]{x1D70F}^{-1}\leqslant \unicode[STIX]{x1D70E}<0.10\unicode[STIX]{x1D70F}^{-1}$ .

Figure 10.

Blow up of the lower left corner of figure 8(a), showing details of the eigenvalues in the region where the Gaussian vortices have $N_{c}^{2}<0$ . The axes of the figure, line styles and symbols have the same meaning as they do in figure 8(a). In the lower left region, below the broken line, numbers rather than symbols are used to indicate where in parameter space we have carried out linear stability calculations. The numbers are the values $\unicode[STIX]{x1D70E}$ (in units of $\unicode[STIX]{x1D70F}^{-1}$ ) of the fastest-growing eigenmode (which in all cases has an A4 symmetry).

Figure 11.

Growth rates (in units of $\unicode[STIX]{x1D70F}^{-1}$ ) for the most unstable eigenmode of vortices with $N_{c}^{2}<0$ for fixed $Ro$ and $Bu$ . $f/\bar{N}=0.1$ . For all vortices examined in this region, the fastest-growing eigenmode has A4 symmetry; (a) $\unicode[STIX]{x1D70E}$ as a function of $Bu$ for $Ro=-0.4$ , $-0.3$ , $-0.2$ and $-0.1$ ; (b) $\unicode[STIX]{x1D70E}$ as a function of $Ro$ for $Bu=0.225$ , 0.19, 0.14 and 0.075. In each panel, the value of the horizontal coordinate axis on the right side of the panel corresponds to a vortex with $N_{c}^{2}=0$ . The dotted lines are to guide the eye.

4.2 Comparison with previous studies

As discussed in § 1, this paper extends the analyses of previous studies by using the full 3-D non-hydrostatic Boussinesq equations and by employing the 3-D Gaussian vortex model, which has continuous velocity and density (and PV) fields and is initially in exact equilibrium. The latter is necessary for a rigorous linear stability analysis. A comparison of our results with those of many previous studies is not straightforward because various different vortex models and flow models have been used. Below we compare our parameter map of stability with the results of the most relevant study in the QG limit (Nguyen et al. Reference Nguyen, Hua, Schopp and Carton2012) and with the results of several relevant studies using multilayer models. We also discuss the results of Yim et al. (Reference Yim, Billant and Ménesguen2016), who used the full Boussinesq equations but studied a different family of vortices.

In the limit of vanishing $Ro$ , the most relevant study to ours is that of Nguyen et al. (Reference Nguyen, Hua, Schopp and Carton2012), who numerically calculated the unstable modes of a Gaussian vortex using the QG equations. They found that the fastest-growing mode changes from S2 to A1 around $Bu=1$ , which along with the general dependence of the growth rate of the fastest-growing mode on $Bu$ in their figure 1(a) agrees overall with the results of current study (see figure 9 a which is for $Ro=0.05$ ). However, they also found that for $Bu$ as small as 0.05, modes with higher $m$ dominate. In our results, for anticyclones, as $Bu$ decreases, the most unstable mode changes from S2 to A4 once the vortex becomes statically unstable (this instability is not considered in the QG framework used by Nguyen et al. (Reference Nguyen, Hua, Schopp and Carton2012)). For cyclones, as $Bu$ decreases, the most unstable mode changes from S2 to A1 for small $Ro$ and to A2 or A3 for moderate $Ro$ (see figure 8).

There are a number of studies which have used the shallow-water equations with the Gaussian vortex model and are relevant to current work. Consistent with the results of our analysis, these studies find that anticyclones become more stable as the absolute value of the Rossby number increases, whereas for cyclones the growth rates decrease with decreasing the Rossby number (Stegner & Dritschel Reference Stegner and Dritschel2000; Baey & Carton Reference Baey and Carton2002; Benilov & Flanagan Reference Benilov and Flanagan2008). (In this section note that our results are only for vortices with stably stratified interiors.)

How the growth rates in these studies vary with the Burger number, however, shows a strong dependence on the vertical structure of the vortex and the background flow. Stegner & Dritschel (Reference Stegner and Dritschel2000) studied the stability of isolated Gaussian vortices using a $1{-}1/2$ layer model and found that for vortices with small Rossby numbers, the growth rate decreases with decreasing the Burger number. This is consistent with our results only for $Bu\gtrsim 1$ . Benilov & Flanagan (Reference Benilov and Flanagan2008) used a two-layer model to examine the stability of the ‘compensated’ (i.e. $\boldsymbol{v}=\mathbf{0}$ in the bottom layer) Gaussian vortices, and also Gaussian vortices with uniform PV in the lower layer. They found that compensated vortices are neutrally stable for intermediate Burger numbers, while vortices with uniform PV in the lower layer are neutrally stable for Burger numbers smaller than a critical value of order 1. Baey & Carton (Reference Baey and Carton2002) studied two-layer Gaussian vortices and found, in contrast to the previous results and those of ours, that the growth rate decreases with Burger number for both cyclones and anticyclones and the eigenmodes are stable for Burger numbers larger than a critical value. It is apparent that identifying a unique stability behaviour with Burger number in these studies is difficult and the behaviour is highly dependent on the vertical structure of the flow/vortex. An example of such dependence is given by Sutyrin (Reference Sutyrin2015), who examined two- and three-layer compensated shallow-water vortices and showed that the addition of a third middle layer with uniform PV weakens the coupling between the upper and lower layers and enhances the stability of vortices. Considering these results, comparing the Burger number dependence of the stability behaviour of 3-D vortices in continuously stratified Boussinesq flows and vortices in shallow water and layer models is not particularly useful.

Only a few studies have used the full Boussinesq equations, and even those have focused on very different vortex models such as barotropic Taylor columns (Smyth & McWilliams Reference Smyth and McWilliams1998), evolving (out-of-equilibrium) 3-D vortices interacting with large-scale internal waves (Brunner-Suzuki et al. Reference Brunner-Suzuki, Sundermeyer and Lelong2012), out-of-equilibrium, ellipsoidal 3-D vortices with discontinuous PV profiles (Tsang & Dritschel Reference Tsang and Dritschel2015) and 3-D equilibrium vortices with Gaussian angular velocity (Yim et al. Reference Yim, Billant and Ménesguen2016). Here we focus on the latter, because the main difference between our analysis and that of Yim et al. (Reference Yim, Billant and Ménesguen2016) is in the vortex model: Gaussian pressure anomaly in the current study versus their Gaussian angular velocity (also note that the flow in their study is not inviscid). Such comparison provides some understanding of how the stability properties depend on the vortex profile.

Yim et al. (Reference Yim, Billant and Ménesguen2016) conducted a linear stability analysis of 3-D equilibrium vortices with Gaussian angular velocity in unbounded, rotating, stratified flows for a wide range of Rossby numbers, $|Ro|\leqslant 20$ . Here we only focus on their results for $|Ro|\leqslant 0.5$ and inviscid and non-diffusive flows, which are relevant to the present study. Consistent with our results, for $Bu\gtrsim 1$ , they found A1 as the most unstable mode for both cyclones and anticyclones (their figure 39d and f), which they attributed to the instability mechanism of Gent & McWilliams (Reference Gent and McWilliams1986) (this is also consistent with the results of Smyth & McWilliams (Reference Smyth and McWilliams1998) for Taylor columns). For $Bu\lesssim 1$ , Yim et al. (Reference Yim, Billant and Ménesguen2016) found anticyclones neutrally stable for $0.5\lesssim Bu\lesssim 1$ (while we found them weakly unstable), and they found S2 as the most unstable mode for anticyclones between the statically unstable region and $Bu\sim 0.4{-}0.5$ (depending on $Ro$ ), which is consistent with our results. For cyclones with $Bu\lesssim 1$ , Yim et al. (Reference Yim, Billant and Ménesguen2016) found a neutrally stable region between $0.5\lesssim Bu\lesssim 1$ (variable with $Ro$ ), which is much larger than (and encompasses) the neutrally stable region we found; they also found that as $Bu$ decreases from one, modes with $m=2$ become the most unstable ones before modes with $m=1$ also becoming unstable at lower $Bu$ , which is overall consistent with our results. At $Bu$ as low as 0.3, the family of vortices studied by Yim et al. (Reference Yim, Billant and Ménesguen2016) can have statically unstable cyclones, while cyclones in the family of vortices we studied are always statically stable. The comparison of the results of the current study and those of Yim et al. (Reference Yim, Billant and Ménesguen2016), as summarized above, suggests that for these two vortex families, while the linear stability properties are not sensitive to the vortex profile for $Bu\gtrsim 1$ , the stability properties strongly depend on the vortex profile for $Bu\lesssim 1$ . Whether this behaviour is generic or not requires further studies with other vortex families.