2.1. Governing equations

We use the incompressible Navier–Stokes equations under the Boussinesq approximation

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

(2.2)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \left( \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \right) \boldsymbol{u} ={-}\boldsymbol{\nabla} \left( \frac{p} {\rho_0} \right) + b \boldsymbol{e_z} - 2\boldsymbol{\varOmega_b} \times \boldsymbol{u} + \nu {\nabla}^{2}\boldsymbol{u}, \end{gather}$$

(2.3)$$\begin{gather} \frac{\partial b}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} b + N^{2} u_z = \kappa \nabla^{2}b, \end{gather}$$

where $\boldsymbol {u}$ is the velocity field, $p$ is the pressure and $b = -g \rho /\rho _0$ is the buoyancy, $g$ is the gravity, $\rho$ is the density perturbation and $\rho _0$ is a constant reference density. Here $\boldsymbol {e_z}$ is the vertical unit vector and $\boldsymbol {\varOmega _b}$ is the background rotation vector. It is assumed to have not only a vertical component but also a horizontal component along the $y$ direction: $2\boldsymbol {\varOmega _b} = \tilde {f} \boldsymbol {e_y} + f \boldsymbol {e_z}$, where $f = 2\varOmega _b \sin {(\phi )}$ and $\tilde {f} = 2 \varOmega _b \cos {(\phi )}$ with $\phi$ the latitude or, equivalently, the angle between the background rotation vector and the unit vector in the $y$ direction, $\boldsymbol {e}_y$ (figure 1). Here $\nu$ and $\kappa$ are the viscosity of the fluid and diffusivity of the stratifying agent, respectively. The total density field $\rho _t$ reads $\rho _t(\boldsymbol {x},t)=\rho _0+\bar {\rho }(z)+\rho (x,t)$, where $\bar {\rho }$ is the mean density profile along the $z$ axis. The Brunt–Väisälä frequency

(2.4)\begin{equation} N = \sqrt{-\frac{g}{\rho_0}\frac{{\rm d}\bar{\rho}}{{\rm d}z}} \end{equation}

is assumed to be constant.

Figure 1. Sketch of the initial vortex in a stratified fluid and in the presence of a background rotation $\varOmega _b$ inclined with an angle $\phi$.

2.2. Initial conditions

A single vertical vortex with a Lamb–Oseen profile is taken as initial conditions. Its vorticity field reads

(2.5)\begin{equation} \boldsymbol{\omega} (\boldsymbol x,t=0) = \zeta \boldsymbol{e_z} = \frac{\varGamma}{{\rm \pi} a_0^{2}}\exp({-r^{2}/a_0^{2}}) \boldsymbol{e_z}, \end{equation}

where $\varGamma$ is the circulation, $a_0$ is the radius and $r$ is the radial coordinate. The geometry of the flow is sketched in figure 1.

2.3. Non-dimensionalization

Equations (2.1)–(2.3) are non-dimensionalized by using $2{\rm \pi} a_0^{2}/\varGamma$ and $a_0$ as time and length units,

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

(2.7)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \left( \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \right) \boldsymbol{u} ={-}\boldsymbol{\nabla}p + b \boldsymbol{e_z} - 2\left ( \frac{1}{Ro} \boldsymbol{e_z} + \frac{1}{\widetilde{Ro}} \boldsymbol{e_y} \right) \times \boldsymbol{u} + \frac{1}{Re} {\nabla}^{2} \boldsymbol{u}, \end{gather}$$

(2.8)$$\begin{gather}\frac{\partial b}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} b + \frac{1}{F_h^{2}} u_z = \frac{1}{ReSc} \nabla^{2}b, \end{gather}$$

where the same notation has been kept for the non-dimensional variables. Note that $\rho _0$ has been eliminated by redefining the pressure $p$. The Reynolds, Froude, Rossby and Schmidt numbers are defined as

(2.9a–e)\begin{equation} Re = \frac{\varGamma}{2{\rm \pi}\nu}, \quad F_h = \frac{\varGamma}{2{\rm \pi} a_0^{2} N}, \quad Ro = \frac{\varGamma}{{\rm \pi} a_0^{2}\,f}, \quad \widetilde{Ro} = \frac{\varGamma}{{\rm \pi} a_0^{2}\,\tilde{f}}, \quad Sc=\frac{\nu}{\kappa}. \end{equation}

Note that two Rossby numbers $Ro$ and $\widetilde {Ro}$ are defined based on the two components of the rotation vector. The Schmidt number is always set to unity. In the following, all results will be reported in non-dimensional form.

2.4. Numerical method

A pseudo-spectral method with periodic boundary conditions is used to integrate (2.1)–(2.3) in space (Deloncle, Billant & Chomaz Reference Deloncle, Billant and Chomaz2008). Time integration is performed with a fourth-order Runge–Kutta scheme. Most of the aliasing is removed by truncating the top one third of the modes along each direction. The viscous and diffusive terms are integrated exactly. The horizontal sizes of the computational domain have been set to $l_x=l_y=18$. As shown by Bonnici (Reference Bonnici2018) and Billant & Bonnici (Reference Billant and Bonnici2020), this is sufficiently large to minimize the effects of the periodic boundary conditions and to give results independent of the box sizes. In particular, the periodic boundary condition imposes that the net circulation over the domain is zero implying that the initial vorticity is not exactly (2.5) but $\zeta ^{\prime } = \zeta - \varGamma /(l_xl_y)$. However, with $l_x=l_y=18$, the artificial background vorticity $\varGamma /(l_xl_y)$ is weak and represents only $1\,\%$ of the maximum vorticity of the vortex. The horizontal resolution has been varied from $n_x=n_y=512$ for $Re=2000$ up to $n_x=n_y=1024$ for $Re=10\,000$. As mentioned in the introduction, preliminary simulations were fully three dimensional with a resolution $n_z$ and a vertical size $l_z$ similar to the horizontal ones. However, the flow is always observed to remain independent of the vertical coordinate. It can be seen indeed from (2.6)–(2.8) that if $\partial / \partial z = 0$ at $t=0$ then the flow will remain independent of the vertical coordinate for all time. Therefore, only two-dimensional simulations but with three components of velocity will be presented in the following. Several tests using different horizontal resolutions have been performed in order to verify the accuracy of the computations. For $Re=2000$, the velocity has been found to differ by less than $0.1\,\%$ when the resolution is increased from $512 \times 512$ to $1024 \times 1024$. Similarly, for $Re=10\,000$, the relative variation of the velocity is less than $0.5\,\%$ when the resolution is increased from $1024 \times 1024$ to $1536 \times 1536$.

3.1. Illustrative example of the vortex dynamic

To get an overview of the effect of the complete Coriolis force, we start by presenting the vortex evolution for the sample set of parameters $Re=2000$, $F_h=10$, $Ro=23.1$, $\phi =60^{\circ }$ ($\widetilde {Ro}=40$). Figures 2 and 3 show the evolution of the vertical velocity and vertical vorticity at six different times (a movie is available in the supplementary material available online at https://doi.org/10.1017/jfm.2022.812). Initially, the vortex is completely axisymmetric (figure 3a) and the vertical velocity is zero (figure 2a). As time goes on, a vertical velocity field with an azimuthal wavenumber $m=1$ develops (figure 2b). This structure tends to intensify and becomes more and more concentrated at a particular radius (figure 2c). Concomitantly, a ring of negative vertical vorticity appears and grows at the same radius (figure 3b,c). Later, the vertical velocity structure and the ring of anomalous vertical vorticity are no longer perfectly circular (figure 3d–f). Two negative vortices appear on the vertical vorticity ring and revolve around the vortex centre (figure 3e, f). Simultaneously, the shape of the vertical velocity structure is deformed similarly (figure 2e, f). As already mentioned, preliminary three-dimensional simulations with various vertical sizes $l_z$ of the computational domain and resolutions $n_z$ have shown that the velocity and vorticity fields remain always completely independent of the vertical coordinate, i.e. the same evolution is observed in any horizontal cross-section. It is also important to stress that this phenomenon occurs only in the presence of the complete Coriolis force. Indeed, if $\widetilde {Ro} = \infty$, the vertical velocity remains identically zero while the vertical vorticity simply decays by viscous diffusion.

Figure 2. Vertical velocity at different times: $(a)$ $t=1$ , $(b)$ $t=30$ , $(c)$ $t=60$, $(d)$ $t=75$, $(e)$ $t=80$ and $(f)$ $t=90$ for $Re=2000$, $F_h=10$, $Ro=23.1$, $\phi =60^{\circ }$ ($\widetilde {Ro}=40$).

Figure 3. Vertical vorticity at different times: $(a)$ $t=1$ , $(b)$ $t=30$ , $(c)$ $t=60$, $(d)$ $t=75$, $(e)$ $t=80$ and $(f)$ $t=90$ for $Re=2000$, $F_h=10$, $Ro=23.1$, $\phi =60^{\circ }$ ($\widetilde {Ro}=40$).

From figure 3 we can distinguish two phases in the evolution of the vortex. First a circular ring of anomalous vertical vorticity develops and then this ring becomes non-axisymmetric. A more in-depth analysis of these two phases will be discussed later. Let us first examine the effects of the control parameters on this phenomenon.

3.2. Effects of the stratification

When the Froude number is decreased from $F_h=10$ to $F_h=2$, the same evolution of the vertical velocity (figure 4a–d) – see the movie in the supplementary material – and the vertical vorticity (figure 4e–h) is observed but at a smaller radius. However, the anomaly of the vertical vorticity does not become negative this time and the maximum vertical velocity is also lower than for $F_h=10$. Figure 5(a) shows the evolution of the vertical velocity maximum $u_{zm}$. In fact, three phases can be distinguished. First, $u_{zm}$ increases linearly with time with some small oscillations superimposed that will be later attributed to inertia-gravity waves. The propagation of waves can be also seen in the vertical velocity field at the beginning of the movies. Then, $u_{zm}$ saturates and tends to slightly decrease. When $t\gtrsim 60$, large oscillations arise when the vortex becomes fully non-axisymmetric.

Figure 4. Vertical velocity (a–d) and vertical vorticity (e–h) at different times: (a,e) $t=10$, (b, f) $t=40$, (c,g) $t=65$ and (d,h) $t=80$ for $Re=2000$, $F_h=2$, $Ro=23.1$, $\phi =60^{\circ }$ ($\widetilde {Ro}=40$).

Figure 5. Maximum vertical velocity as a function of time for $F_h=2$ and (a) $Re=2000$, $Ro=23.1$, $\phi =60^{\circ }$ ($\widetilde {Ro}=40$), (b) $Re=2000$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$) and (c) $Re=10\,000$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$).

If the Froude number is below unity, such evolution is no longer observed. In § 4 it will be shown that this phenomenon is due to the presence of a critical layer where the angular velocity of the vortex is equal to the non-dimensional Brunt–Väisälä frequency, i.e. $\varOmega = 1/F_h$, where $\varOmega$ is the non-dimensional angular velocity of the vortex corresponding to the vorticity field (2.5).

3.3. Effects of the Rossby numbers

When the traditional Rossby number $Ro$, which is based on the vertical component of the background rotation, is varied, while keeping the other numbers fixed ($Re$, $F_h$, $\widetilde {Ro}$), the vortex evolution remains strictly identical. This is consistent with the fact that the dynamics is independent of the vertical coordinate so that the Coriolis force associated with $Ro$ can be eliminated from (2.7) by redefining the pressure.

In contrast, varying the non-traditional Rossby number $\widetilde {Ro}$, which is based on the horizontal component of the background rotation, has important effects on the evolution of the vortex. Since the Rossby number $Ro$ has no effect, the effect of $\widetilde {Ro}$ has been studied by varying the latitude $\phi$ while keeping the background rotation rate $\varOmega _b$ constant. Hence, both $\widetilde {Ro}$ and $Ro$ varies. Figure 6 shows the evolution of the vertical velocity and vertical vorticity when the latitude is increased from $\phi = 60^{\circ }$ to $\phi =80^{\circ }$ ($\widetilde {Ro}$ is increased from $\widetilde {Ro}=40$ to $\widetilde {Ro}=115.2$) while keeping the other parameters fixed (a movie is available in the supplementary material). The initial evolution (figure 6a,b,e, f) is similar to that in figure 4 but, later (figure 6c,d,g,h), the fields keep the same shape, i.e. no significant asymmetric deformations can be seen in contrast to figure 4(d,h). As seen in figure 5(b), only two phases are then present in the evolution of the maximum vertical velocity. First, a phase where $u_{zm}$ grows linearly with weak oscillations and, second, a phase where $u_{zm}$ remains approximately constant. In addition, the maximum vertical velocity is lower (figure 6b,c) and the anomalous vorticity ring is weaker (figure 6g) than in figure 4. At late time $t=200$ (figure 6d,h), we can see that the vertical velocity has decreased by viscous diffusion while the vertical velocity field has moved towards the centre of the vortex. This is consistent with the critical layer's interpretation since, as the angular velocity decays, the radius where $\varOmega = 1/F_h$ decreases.

Figure 6. Vertical velocity (a–d) and vertical vorticity (e–h) at different times: (a,e) $t=10$, (b, f) $t=50$, (c,g) $t=150$ and (d,h) $t=200$ for $Re=2000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$).

3.4. Effects of the Reynolds number

Figure 7 shows the evolution of the vertical velocity and vertical vorticity when the Reynolds number is increased from $Re=2000$ to $Re=10\,000$ while keeping the other parameters as in figure 6. A movie is also available in the supplementary material. The maximum vertical velocity is almost doubled and the vertical velocity field is much thinner and focused near a given radius (figure 7b) than in figure 6. Furthermore, the ring of anomalous vertical vorticity (figure 7f) is more intense. Later, asymmetric deformations of this ring and of the vertical velocity field are clearly visible (figure 7c,d,g,h) in contrast to figure 6. Oscillations are then visible in the evolution of $u_{zm}$ (figure 5c).

Figure 7. Vertical velocity (a–d) and vertical vorticity (e–h) at different times: (a,e) $t=10$, (b, f) $t=100$, (c,g) $t=150$ and (d,h) $t=200$ for $Re=10\,000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$).

3.5. Combined effects of $\widetilde {Ro}$ and $Re$

In the previous sections we have seen that the vertical vorticity becomes fully non-axisymmetric in a second stage if $\phi$ is sufficiently lower than $90^{\circ }$, i.e. if $\widetilde {Ro}$ is not too large (figure 4) or if the Reynolds number is large enough (figure 7), otherwise the vertical vorticity field remains quasi-axisymmetric (figure 6). Figure 8 summarizes several other simulations for various $\widetilde {Ro}$ and Reynolds numbers $Re$ keeping the Froude number equal to $F_h=2$. Yellow and blue symbols indicate simulations where the vertical vorticity remains quasi-axisymmetric or becomes non-axisymmetric, respectively. We can see that the critical Rossby number $\widetilde {Ro}_c$ above which the vortex remains quasi-axisymmetric increases with the Reynolds number from $\widetilde {Ro}_c \simeq 100$ for $Re=2000$ to $\widetilde {Ro}_c \simeq 400$ for $Re=10\,000$.

Figure 8. Map of the simulations in the parameter space $(Re,\widetilde {Ro})$ for $F_h=2$. The yellow and blue circles represent simulations where the vertical vorticity remains quasi-axisymmetric or not, respectively. The solid and dashed lines represent the criterion (6.14) for different values of (a,c): $(\infty,0)$ and $(\infty,0.4)$, respectively. The number near some points indicate the figure numbers where the corresponding simulation is shown.

4.1. Unsteady inviscid analysis

Accordingly, we first consider (4.1a)–(4.1e) in the inviscid limit $Re = \infty$, but keep the time derivatives. At leading order, the equations reduce to

(4.7)$$\begin{gather} -\frac{u_{\theta 0}^{2}}{r} = \frac{\partial p_0}{\partial r} + \frac{2u_{\theta 0}}{Ro}, \end{gather}$$

(4.8)$$\begin{gather}\frac{\partial u_{\theta 0}}{\partial t} = 0, \end{gather}$$

so that $u_{\theta 0} = r \varOmega$ as before. At first order, (4.1d)–(4.1e) become

(4.9a)$$\begin{gather} \frac{\partial u_{z1}}{\partial t} + \varOmega \frac{\partial u_{z1}}{\partial \theta} = b_{1} - r\varOmega \sin{(\theta)}, \end{gather}$$

(4.9b)$$\begin{gather}\frac{\partial b_{1}}{\partial t} + \varOmega \frac{\partial b_{1}}{\partial \theta} ={-} \frac{u_{z1}}{F_h^{2}}. \end{gather}$$

The only difference with (4.6) is the presence of the time derivatives. By imposing ${u_{z1} = b_1 = 0}$ at $t = 0$, the solutions can be found in the form

(4.10a)$$\begin{gather} u_{z1} = {u_z}_{p} \mathrm{e}^{{\rm i}\theta} + {u^{*}_z}_{p} \mathrm{e}^{-{\rm i}\theta}, \end{gather}$$

(4.10b)$$\begin{gather}{b}_{1}= {b}_{p} \mathrm{e}^{{\rm i}\theta} + {b^{*}}_{p} \mathrm{e}^{-{\rm i}\theta}, \end{gather}$$

where the star denotes the complex conjugate and

(4.11a)$$\begin{gather} {u_z}_{p} = \frac{r \varOmega }{4} \left[- \frac{1}{\alpha} \left( 1 - \mathrm{e}^{{\rm i} \alpha t} \right) + \frac{1}{\beta} \left( 1 - \mathrm{e}^{-{\rm i} \beta t} \right) \right], \end{gather}$$

(4.11b)$$\begin{gather}{b}_{p} ={-} {\rm i} \frac{r \varOmega }{4 F_h} \left[ \frac{1}{\alpha} \left( 1 - \mathrm{e}^{{\rm i} \alpha t} \right)+ \frac{1}{\beta} \left( 1 - \mathrm{e}^{-{\rm i} \beta t} \right) \right], \end{gather}$$

with

(4.12a,b)\begin{equation} \alpha = \frac{1-F_h\varOmega}{F_h}, \quad \beta =\frac{1+F_h\varOmega}{F_h}. \end{equation}

Compared with the steady solution (4.6), the additional terms present in (4.11) correspond to waves generated at $t=0$ to satisfy the initial conditions. These waves oscillate at the frequencies $1/F_h - \varOmega$ and $-1/F_h - \varOmega$, i.e. the non-dimensional Brunt–Väisälä frequency with an additional Doppler shift coming from the azimuthal motion of the vortex. Hence, they correspond to inertia-gravity waves with a zero vertical wavenumber. In contrast to (4.6), we see now that the vertical velocity and buoyancy (4.10a)–(4.10b) are no longer singular at the radius $r_c$ where $\varOmega (r_c) = 1/F_h$. Indeed, we have $(1-e^{i\alpha t})/\alpha \simeq -i t$ when $\alpha \to 0$ so that (4.11a)–(4.11b) remain finite at $r=r_c$.

Following Wang & Balmforth (Reference Wang and Balmforth2020), the behaviour of these solutions in the vicinity of the critical radius can be studied more precisely by introducing the variable $\eta = r - r_c$. When $\eta \ll 1$, the vertical velocity approximates to

(4.13)\begin{align} u_{z1} &= \left[\! \left(\!\frac{r_c \varOmega_c}{2\eta \varOmega^{\prime}_c} + \frac{\varOmega_c}{2\varOmega^{\prime}_c} + \frac{r_c}{2} - \frac{r_c \varOmega^{\prime\prime}_c \varOmega_c }{{4\varOmega^{\prime}_c}^{2}} \!\right) \left(\!1\!-\cos{(\eta \varOmega^{\prime}_c t)}\!\right) + \frac{r_c}{4} \left(\!1\!-\cos{\left(\!\frac{2}{F_h}t\!\right)}\!\right) \!\right] \cos{(\theta)} \nonumber\\ &\quad - \left[ \left(\frac{r_c \varOmega_c}{2\eta \varOmega^{\prime}_c} + \frac{\varOmega_c}{2\varOmega^{\prime}_c} + \frac{r_c}{2} - \frac{r_c \varOmega^{\prime\prime}_c \varOmega_c }{4{\varOmega^{\prime}_c}^{2}} \right) \sin{(\eta \varOmega^{\prime}_c t)} + \frac{r_c}{4} \sin{\left(\frac{2}{F_h}t\right)} \right] \sin{(\theta)} \!+\! \textit{O} (\eta), \end{align}

where the subscript $c$ indicates a value taken at $r_c$. The terms involving $\eta \varOmega ^{\prime }_ct$ in the sin and cos functions have not been expanded in (4.13) since this quantity can be large when $t$ is large even for small $\eta$. The approximation (4.13) is therefore uniformly valid whatever $t$. Taking into account only the leading order, (4.13) can be simplified and rewritten as

(4.14)\begin{equation} u_{z1} = \frac{r_c \varOmega_c t}{2} \left[ \left(\frac{1-\cos{U}}{U}\right) \cos{(\theta)} - \frac{\sin{U}}{U}\sin{(\theta)} \right] + {O} (1), \end{equation}

where $U = \eta \varOmega ^{\prime }_c t$. This expression shows that the radial profile of the vertical velocity in the vicinity of $r_c$ depends only on the self-similar variable $U$. This implies that the vertical velocity will be more and more concentrated around $r_c$ as time increases. In addition, (4.14) shows that the amplitude of $u_{z1}$ will increase linearly with time for a fixed value of $U$. In practice, the approximation (4.14) will be accurate only when $t$ is large since the neglected terms are of order unity. Since $\varOmega _c = 1/F_h$, this will occur increasingly later when $F_h$ increases. We emphasize that the inertia-gravity wave oscillating at frequency $1/F_h+\varOmega _c=2/F_h$ is neglected in (4.14) unlike in (4.13).

4.2. Unsteady viscous analysis

Here, viscous and diffusive effects are taken into account in addition to the time evolution. Following Boulanger et al. (Reference Boulanger, Meunier and Le Dizes2007) and Wang & Balmforth (Reference Wang and Balmforth2020, Reference Wang and Balmforth2021), we assume that the Reynolds number is large $Re \gg 1$ and consider only the vicinity of the critical radius by introducing a rescaled radius such that $\tilde {r} = Re^{1/3}(r-r_c)$. We also assume that the evolution occurs over the slow time $T=Re^{-1/3}t$. Since $Re \gg 1$, the leading-order solution of (4.1a)–(4.1e) is still $u_{\theta 0} = r \varOmega$. At order $\varepsilon$, (4.1d)–(4.1e) become

(4.15a) \begin{align} & \frac{1}{Re^{1/3}}\frac{\partial u_{z1}}{\partial T}+ \left( \varOmega_c + \frac{\tilde{r} \varOmega_c^{\prime}}{Re^{1/3}} + {O} \left(\frac{1}{Re^{2/3}}\right) \right) \frac{\partial u_{z1}}{\partial \theta}\nonumber\\ &\quad = b_1 \!- \left( r_c \varOmega_c \!+ \frac{\tilde{r}(\varOmega_c+r_c\varOmega_c^{\prime})}{Re^{1/3}} + {O} \left(\frac{1}{Re^{2/3}}\right) \right) \sin{(\theta)} \!+ \frac{1}{Re^{1/3}} \frac{\partial^{2} u_{z1}}{\partial\tilde{r}^{2}} + {O} \left(\frac{1}{Re^{2/3}}\right) , \end{align}

(4.15b)\begin{align}& \frac{1}{Re^{1/3}} \frac{\partial b_{1}}{\partial T} \!+ \left(\!\varOmega_c + \frac{\tilde{r} \varOmega_c^{\prime}}{Re^{1/3}} + {O} \left(\!\frac{1}{Re^{2/3}}\!\right)\!\right) \frac{\partial b_1}{\partial\theta}={-} \frac{u_{z1}}{F_h^{2}} \!+ \frac{1}{Sc} \left[\!\frac{1}{Re^{1/3}} \frac{\partial^{2}{b}_1}{\partial\tilde{r}^{2}} + {O} \left(\!\frac{1}{Re^{2/3}}\!\right) \!\!\right] . \end{align}

These equations can be solved by expanding $u_{z1}$ and $b_1$ as

(4.16a)$$\begin{gather} {u_{z1}} = Re^{1/3} \left[ \tilde{u}_{z1}(\tilde{r},T) \mathrm{e}^{{\rm i}\theta} +{\rm c.c.} \right]+ \tilde{u}_{z2}(\tilde{r},T) \mathrm{e}^{{\rm i}\theta} +{\rm c.c.} + \dots, \end{gather}$$

(4.16b)$$\begin{gather}b_1 = Re^{1/3} \left[ \tilde{b}_1(\tilde{r},T) \mathrm{e}^{{\rm i}\theta} + {\rm c.c.} \right] + \tilde{b}_2(\tilde{r},T) \mathrm{e}^{{\rm i}\theta} +{\rm c.c.} + \dots\ . \end{gather}$$

By substituting (4.16a)–(4.16b) in (4.15a)–(4.15b), we get at order ${O}(Re^{1/3})$,

(4.17a)$$\begin{gather} {\rm i} \varOmega_c \tilde{u}_{z1} = \tilde{b}_1 , \end{gather}$$

(4.17b)$$\begin{gather}{\rm i} \varOmega_c \tilde{b}_1 ={-} \frac{ \tilde{u}_{z1}}{F_h^{2}}, \end{gather}$$

which both yield

(4.18)\begin{equation} \tilde{b}_1 = {\rm i} \varOmega_c \tilde{u}_{z1}. \end{equation}

The equations at order ${O} (1)$ are

(4.19a)$$\begin{gather} \frac{\partial\tilde{u}_{z1}}{\partial T} + {\rm i}\varOmega_c \tilde{u}_{z2} + {\rm i} \tilde{r}\varOmega_c^{\prime} \tilde{u}_{z1}= \tilde{b}_2 - \frac{r_c\varOmega_c}{2i} + \frac{{\rm d}^{2}\tilde{u}_{z1}}{{\rm d}\tilde{r}^{2}} , \end{gather}$$

(4.19b)$$\begin{gather}\frac{\partial\tilde{b}_{1}}{\partial T}+ {\rm i}\varOmega_c \tilde{b}_2 + {\rm i} \tilde{r}\varOmega_c^{\prime}\tilde{b}_1 ={-}\frac{1}{F_h^{2}} \tilde{u}_{z2} + \frac{1}{Sc} \frac{{\rm d}^{2}\tilde{b}_1}{{\rm d}\tilde{r}^{2}}. \end{gather}$$

By using (4.18), the solvability condition to find $\tilde {b}_2$ and $\tilde {u}_{z2}$ from (4.19a)–(4.19b) requires $\tilde {u}_{z1}$ to satisfy

(4.20)\begin{equation} \frac{\partial\tilde{u}_{z1}}{\partial T}+ {\rm i} \tilde{r}\varOmega_c^{\prime} \tilde{u}_{z1} = \frac{{\rm i}}{4} r_c\varOmega_c + \frac{1}{2} \left (1+\frac{1}{Sc}\right) \frac{{\rm d}^{2}\tilde{u}_{z1}}{{\rm d}\tilde{r}^{2}}. \end{equation}

As shown by Wang & Balmforth (Reference Wang and Balmforth2020, Reference Wang and Balmforth2021), the solution is

(4.21)\begin{equation} \tilde{u}_{z1} = {\rm i} A \frac{1}{\rm \pi} \int_0^{\left|{\varOmega_c^{\prime}}\right|T/ \gamma} \exp \left ( - \frac{z^{3}}{3} + {\rm i} \gamma \tilde{r} z \right) \mathrm{d}z, \end{equation}

where

(4.22a,b)\begin{equation} A=\frac{{\rm \pi} r_c \varOmega_c}{2 \left|{2\varOmega_c^{\prime}}\right|^{2/3}\left(1+\dfrac{1}{Sc}\right)^{1/3}}, \quad \gamma = \frac{ \left|{2\varOmega_c^{\prime}}\right|^{1/3}}{\left(1+\dfrac{1}{Sc}\right)^{1/3}}. \end{equation}

When $\left |{\varOmega _c^{\prime }}\right |T/ \gamma \gg 1$, the upper bound in the integral can be replaced by infinity so that (4.21) recovers the steady solution (Boulanger et al. Reference Boulanger, Meunier and Le Dizes2007)

(4.23)\begin{equation} \tilde{u}_{z1} = {\rm i} A \,{Hi}( {\rm i} \gamma \tilde{r}), \end{equation}

where ${Hi}$ is the Scorer's function (Abramowitz & Stegun Reference Abramowitz and Stegun1972). We will see that the time interval where both the unsteadiness and viscous effects are important is short so that (4.23) turns out to be reached quickly after the inviscid regime.

4.3. Effect on the vertical vorticity

We now turn to the study of the effect of the vertical velocity on the vertical vorticity. Since the flow is invariant along the vertical, the equation for the vertical vorticity (4.2) reduces to

(4.24)\begin{align} \frac{\partial\zeta}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\zeta = \frac{\varepsilon}{r}\left[ \frac{\partial }{\partial r} \left (r {u_z} \sin{(\theta)}\right) + \frac{\partial}{\partial\theta} \left({u_z} \cos{(\theta)} \right) \right] + \frac{1}{Re} \varDelta_h \zeta. \end{align}

This equation shows that the vertical velocity generated at order ${O}(\varepsilon )$ will in turn force a vertical vorticity field at order ${O}(\varepsilon ^{2})$. In order to compute this second-order horizontal flow, the vertical vorticity and streamfunction can be expanded as

(4.25)$$\begin{gather} \zeta = \zeta_0 + \varepsilon^{2} \zeta_2 + \dots, \end{gather}$$

(4.26)$$\begin{gather}\psi = \psi_0 + \varepsilon^{2} \psi_2 + \dots, \end{gather}$$

where $(\zeta _0,\psi _0)$ are the non-dimensional vertical vorticity and streamfunction of the base flow (2.5) and $\zeta _2 = \varDelta _h \psi _2$. The second-order vertical vorticity is governed by

(4.27)\begin{gather} \frac{\partial \zeta_{2}}{\partial t} + \varOmega \frac{\partial \zeta_2}{\partial \theta} - \frac{1}{r}\frac{\partial\psi_2}{\partial\theta} \frac{\partial\zeta_0}{\partial r} = \frac{1}{r}\left[ \frac{\partial}{\partial r} \left (r u_{z1} \sin{(\theta)}\right) + \frac{\partial}{\partial\theta} \left(u_{z1} \cos{(\theta)} \right) \right]+ \frac{1}{Re} \varDelta_h \zeta_2 . \end{gather}

Since the first-order vertical velocity $u_{z1}$ has an azimuthal wavenumber $m=1$, $\zeta _2$ and $\psi _2$ can be sought in the form

(4.28)$$\begin{gather} \zeta_2 = \zeta_{20} (r,t) + \left[ \zeta_{22}(r,t) {\rm e}^{2{\rm i}\theta} + {\rm c.c.} \right], \end{gather}$$

(4.29)$$\begin{gather}\psi_2 = \psi_{20} (r,t) + \left[ \psi_{22}(r,t) {\rm e}^{2{\rm i}\theta} + {\rm c.c.} \right]. \end{gather}$$

In the following, we will determine only the axisymmetric part $\zeta _{20}$ since we will show that the component $\zeta _{22}$ grows slower than $\zeta _{20}$.

4.3.1. Inviscid evolution of $\zeta _{20}$

When $u_{z1}$ follows the unsteady inviscid expression (4.10a), the axisymmetric vertical vorticity $\zeta _{20}$ is given by

(4.30)\begin{equation} \frac{\partial\zeta_{20}}{\partial t} ={-}\frac{{\rm i}}{2r} \frac{\partial}{\partial r} \left( r(u_{zp}^{*}-u_{zp}) \right). \end{equation}

By assuming that $\zeta _{20}=0$ at $t=0$, the solution reads

(4.31)\begin{align} {\zeta}_{20} &= \frac{ \zeta_0 }{4}\left( \frac{\cos{(\alpha t)}-1}{\alpha ^{2}} + \frac{\cos{(\beta t)}-1}{\beta^{2}} \right)\nonumber\\ &\quad + \frac{r\varOmega \varOmega^{\prime}}{4} \left[ \frac{\alpha t \sin{(\alpha t)} + 2(\cos{(\alpha t)}-1) }{\alpha^{3}} - \frac{ \beta t \sin{(\beta t)} + 2(\cos{(\beta t)}-1) }{\beta ^{3}} \right]. \end{align}

Considering the vicinity of the critical radius $\eta = r - r_c$, this solution approximates at leading order to

(4.32)\begin{equation} {\zeta}_{20} = \frac{r_c\varOmega_c}{{4\varOmega^{\prime}_c}^{2}\eta^{3}} \left[ 2 \left(1-\cos{(\eta\varOmega^{\prime}_c t)}\right) - \eta\varOmega_c^{\prime}t \sin{( \eta\varOmega_c^{\prime}t)} \right] + {O}\left(\frac{1}{\eta^{2}}\right), \end{equation}

where $\eta t$ has been considered finite as in (4.13) so that the approximation remains valid even for long time. In terms of the similarity variable $U$, (4.32) can be rewritten as

(4.33)\begin{equation} {\zeta}_{20} = \frac{r_c \varOmega_c \varOmega_c^{\prime} t^{3}}{4} \left[ \frac{2(1-\cos{U})-U\sin{U}}{U^{3}} \right] + {O} \left(\frac{1}{\eta^{2}}\right). \end{equation}

This expression shows that $\zeta _{20}$ remains finite and even vanishes when $U \to 0$, but $\zeta _{20}$ is more and more concentrated in the vicinity of the critical layer as time increases. Furthermore, its amplitude increases like $t^{3}$ at leading order.

Using the same approach for $\zeta _{22}$, it can be shown that its amplitude grows like $t^{2}$ instead of $t^{3}$. Indeed, the forcing term of $\zeta _{22}$ is proportional to $t^{2}$ like in the case of $\zeta _{20}$, but the left-hand side of (4.27) is dominated by the term $2 {\rm i} \varOmega _c \zeta _{22}$ instead of $\partial \zeta _{20} / \partial t$ for the axisymmetric component. Thus, $\zeta _{22}$ is proportional to $t^{2}$, at least initially. This explains why the vertical vorticity is observed in the DNS to remain quasi-axisymmetric during the first two phases.

4.3.2. Viscous evolution of $\zeta _{20}$

When the vertical velocity is given by the unsteady viscous solution (4.16a) and (4.21), it is possible to also obtain its effect on the second-order axisymmetric vertical vorticity $\zeta _{20}$. Expressing first (4.27) in terms of $\tilde {r}$ and $T$ gives at leading order when $Re \gg 1$,

(4.34)\begin{equation} \frac{1}{Re^{1/3}} \frac{\partial\zeta_{20}}{\partial T} = \frac{-{\rm i}}{2}Re^{2/3} \frac{\partial}{\partial\tilde{r}} \left(\tilde{u}_{z1}^{*}-\tilde{u}_{z1} \right) + {O} (Re^{1/3}) + \frac{1}{Re^{1/3}} \frac{\partial^{2}\zeta_{20}}{\partial\tilde{r}^{2}} + \dots\ . \end{equation}

As shown by Wang & Balmforth (Reference Wang and Balmforth2021), the exact solution of (4.34) can be found by means of a Fourier transform in $\tilde {r}$ using (4.21). This gives

(4.35)\begin{align} \zeta_{20} ={-}{\rm i} \frac{Re A}{2\gamma{\rm \pi}} \int_0^{\left|{\varOmega_c^{\prime}}\right|T} \exp \left( \frac{-q^{3}}{3\gamma^{3}} + {\rm i} q \tilde{r} \right) \left( \frac{1-\exp\left( q^{3}/\left|{\varOmega_c^{\prime}}\right| - q^{2}T \right)}{q} \right ) \mathrm{d}q + {\rm c.c.}.\end{align}

In the appendix an approximation of (4.35) valid for large time $T \gg 1$ is found in the form

(4.36)\begin{equation} \zeta_{20} = \zeta_{20}^{(1)} + \zeta_{20}^{(2)}, \end{equation}

where

(4.37a)$$\begin{gather} \zeta_{20}^{(1)} (\tilde{r}) = Re \frac{A}{2} \int_{0}^{\tilde{r}} \left( {{Hi}}^{*}( {\rm i} \gamma u) + {Hi}( {\rm i} \gamma u) \right) \mathrm{d}u, \end{gather}$$

(4.37b)$$\begin{gather}\zeta_{20}^{(2)} (\tilde{r}, t) ={-}Re \frac{A}{2\gamma} \mathrm{erf}\left( \frac{\tilde{r}}{\sqrt{4tRe^{{-}1/3}}} \right). \end{gather}$$

This approximation shows that $\zeta _{20}$ saturates and tends towards (4.37a) for $t Re^{-1/3} \rightarrow \infty$ for $\tilde {r} = {O}(1)$. However, when $|{\tilde {r}}| \rightarrow \infty$ and $tRe^{-1/3}$ is finite, $\zeta _{20}$ vanishes. The approximation (4.36) will be compared with numerical solutions of (4.34) as well as DNS in § 5.2.

To summarize, the axisymmetric vertical vorticity correction at order ${O}(\varepsilon ^{2})$, $\zeta _{20}$, follows two distinct regimes. First, it evolves purely inviscidly and grows like $t^{3}$ according to (4.31). Then, it tends to saturate and follows the approximation (4.36) for large times. This shows that $\zeta _{20}$ saturates towards the steady solution (4.37a) for finite $\tilde {r}$.

4.4. Nonlinear analysis of the critical layer

The viscous linear analysis of the critical layer in §§ 4.2 and 4.3 has shown that the vertical velocity scales as $u_z = {O}(\varepsilon Re^{1/3})$. This creates a vorticity correction of the order $\delta \zeta = {O}(\varepsilon ^{2}Re)$ (see (4.25), (4.28), (4.36)–(4.37)). Since $r = r_c + \tilde {r} / Re^{1/3}$ in the critical layer, the corresponding angular velocity correction $\delta \varOmega$ is given at leading order by $\delta \zeta \simeq r_c Re^{1/3} \partial \delta {\varOmega } / \partial {\tilde {r}}$ so that $\delta \varOmega = \textit {O}(\varepsilon ^{2} Re^{2/3})$. In turn, we see that this angular velocity correction would have the same order as the other terms of order ${O}(1/Re^{1/3})$ in (4.15) if $\varepsilon ^{2}Re^{2/3} = {O}(1/Re^{1/3})$, i.e. if

(4.38)\begin{equation} Re = \frac{\widetilde{Re}}{\varepsilon^{2}}, \end{equation}

where $\widetilde {Re}$ is of order unity. For this distinguished limit, there is therefore a nonlinear feedback of the horizontal flow on the evolutions of the vertical velocity and buoyancy, as considered by Wang & Balmforth (Reference Wang and Balmforth2020, Reference Wang and Balmforth2021). Using the scaling (4.38), the typical order of magnitudes of the different variables can be expressed in terms of $\varepsilon$ only: $u_z={O} ( \varepsilon ^{1/3} )$, $b= {O} ( \varepsilon ^{1/3} )$, $\delta \zeta = {O} ( 1 )$, $\delta \varOmega = {O} ( \varepsilon ^{2/3} )$, $u_r = {O} ( \varepsilon ^{4/3} )$ and the radius and slow time are $\tilde {r} = \varepsilon ^{-2/3}(r-r_c)$ and $T = \varepsilon ^{2/3}t$, respectively. We also assume that $\partial /\partial z=0$. Accordingly, we expand the variables as

(4.39a)$$\begin{gather} u_z = \varepsilon^{1/3} \left[ \tilde{u}_{z1} \mathrm{e}^{{\rm i}\theta} + {\rm c.c.} \right] +\varepsilon \left[ \tilde{u}_{z2} \mathrm{e}^{{\rm i}\theta} + {\rm c.c.} \right] + \dots, \end{gather}$$

(4.39b)$$\begin{gather}b = \varepsilon^{1/3} \left[ \tilde{b}_{1} \mathrm{e}^{{\rm i}\theta} + {\rm c.c.} \right] +\varepsilon \left[ \tilde{b}_{2} \mathrm{e}^{{\rm i}\theta} + {\rm c.c.} \right] + \dots, \end{gather}$$

(4.39c)$$\begin{gather}\varOmega = \varOmega_0 + \varepsilon^{2/3} \varOmega_1 +\dots, \end{gather}$$

(4.39d)$$\begin{gather}u_r = \varepsilon^{4/3} u_{r1} + \dots, \end{gather}$$

(4.39e)$$\begin{gather}\zeta = \zeta_0 + \zeta_1 + \varepsilon^{2/3} \zeta_2 + \dots, \end{gather}$$

where $\varOmega _0$ and $\zeta _0$ are the non-dimensional angular velocity and vorticity corresponding to (2.5). In the vicinity of $r_c$, they can be expanded as

(4.40a)$$\begin{gather} \varOmega_0 =\varOmega_c + \tilde{r} \varOmega_c^{\prime} \varepsilon^{2/3} + \dots, \end{gather}$$

(4.40b)$$\begin{gather}\zeta_0 =\zeta_c + \tilde{r} \zeta_c^{\prime} \varepsilon^{2/3} + \dots\ . \end{gather}$$

Note that there are only components of the form $\exp (\pm {\rm i}\theta )$ in (4.39a)–(4.39b) because the flow is invariant along the vertical and because third harmonics arise only at higher order. Then, (4.1d)–(4.1e) become at leading order

(4.41a)\begin{align} & \varepsilon \frac{\partial\tilde{u}_{z1}}{\partial T} + \varepsilon u_{r1} \frac{\partial \tilde{u}_{z1}}{\partial \tilde{r}} + {\rm i} \left( \varOmega_c + \varOmega_c^{\prime}\tilde{r}\varepsilon^{2/3} \right) \varepsilon^{1/3} \tilde{u}_{z1} + {\rm i} \varepsilon \varOmega_1 \tilde{u}_{z1} + {\rm i} \varepsilon \varOmega_c \tilde{u}_{z2}\nonumber\\ &\quad = \varepsilon^{1/3} \tilde{b}_{1} + \varepsilon \tilde{b}_{2} - \varepsilon \frac{r_c\varOmega_c}{2{\rm i}} + \frac{\varepsilon}{\widetilde{Re}} \frac{\partial^{2}\tilde{u}_{z1}}{\partial \tilde{r}^{2}} + {O} \left( \varepsilon^{5/3} \right), \end{align}

(4.41b)\begin{align} & \varepsilon \frac{\partial \tilde{b}_{1}}{\partial T} + \varepsilon u_{r1} \frac{\partial\tilde{b}_{1}}{\partial\tilde{r}} + {\rm i} \left( \varOmega_c + \varOmega_c^{\prime}\tilde{r}\varepsilon^{2/3} \right) \varepsilon^{1/3} \tilde{b}_{1} + {\rm i} \varepsilon \varOmega_1 \tilde{b}_{1} + {\rm i} \varepsilon \varOmega_c \tilde{b}_{2}\nonumber\\ &\quad ={-} \frac{\varepsilon^{1/3}}{F_h^{2}}\tilde{u}_{z1} - \frac{\varepsilon}{F_h^{2}}\tilde{u}_{z2} + \frac{\varepsilon}{\widetilde{Re}Sc} \frac{\partial^{2}\tilde{b}_{1}}{\partial\tilde{r}^{2}} + {O} \left( \varepsilon^{5/3} \right). \end{align}

Similarly, (4.24) becomes

(4.42)\begin{align} & \varepsilon^{2/3} \frac{\partial\zeta_{1}}{\partial T} + \left( \varOmega_c + \varOmega_c^{\prime}\tilde{r}\varepsilon^{2/3} \right)\frac{\partial\zeta_1}{\partial\theta} + \varepsilon^{2/3} \varOmega_1 \frac{\partial\zeta_1}{\partial\theta} + \varepsilon^{2/3} u_{r1} \frac{\partial\zeta_1}{\partial\tilde{r}} + \varepsilon^{2/3} \varOmega_c \frac{\partial\zeta_2}{\partial\theta} \nonumber\\ & \quad = \varepsilon^{2/3} \frac{\partial}{\partial\tilde{r}} \left( \left( \tilde{u}_{z1} \mathrm{e}^{{\rm i}\theta} + \tilde{u}_{z1}^{*} \mathrm{e}^{-{\rm i}\theta} \right) \sin{(\theta)} \right) + \frac{\varepsilon^{2/3}}{\widetilde{Re}} \frac{\partial^{2}\zeta_1}{\partial\tilde{r}^{2}} + {O} \left( \varepsilon^{4/3} \right). \end{align}

At leading order, (4.41)–(4.42) become

(4.43a)$$\begin{gather} \varepsilon^{1/3} {\rm i} \varOmega_c \tilde{u}_{z1} = \varepsilon^{1/3} \tilde{b}_{1}, \end{gather}$$

(4.43b)$$\begin{gather}\varepsilon^{1/3} {\rm i} \varOmega_c \tilde{b}_{1} ={-}\frac{\varepsilon^{1/3}}{F_h^{2}} \tilde{u}_{z1}, \end{gather}$$

(4.43c)$$\begin{gather}\varOmega_c \frac{\partial\zeta_1}{\partial\theta} = 0. \end{gather}$$

The first two equations are identical to (4.17) and the third one implies that $\zeta _1 \equiv \zeta _1(\tilde {r},T)$ and $u_{r1}=0$. At the next order, we have

(4.44a)$$\begin{gather} \frac{\partial \tilde{u}_{z1}}{\partial T} + {\rm i} \varOmega_c^{\prime}\tilde{r} \tilde{u}_{z1} + {\rm i} \varOmega_1 \tilde{u}_{z1} + {\rm i} \varOmega_c \tilde{u}_{z2} = \tilde{b}_{2} - \frac{r_c\varOmega_c}{2i} + \frac{1}{\widetilde{Re}}\,\frac{\partial^{2}\tilde{u}_{z1}}{\partial\tilde{r}^{2}}, \end{gather}$$

(4.44b)$$\begin{gather}\frac{\partial \tilde{b}_{1}}{\partial T} + {\rm i} \varOmega_c^{\prime}\tilde{r} \tilde{b}_{1} + {\rm i} \varOmega_1 \tilde{b}_{1} + {\rm i} \varOmega_c \tilde{b}_{2} ={-} \frac{\tilde{u}_{z2}}{F_h^{2}} + \frac{1}{\widetilde{Re}Sc}\,\frac{\partial^{2}\tilde{b}_{1}}{\partial \tilde{r}^{2}}, \end{gather}$$

(4.44c)$$\begin{gather}\frac{\partial \zeta_{1}}{\partial T} + \varOmega_c \frac{\partial \zeta_2}{\partial \theta}={-}\frac{{\rm i}}{2} \frac{\partial }{\partial \tilde{r}} \left( \tilde{u}_{z1}^{*} - \tilde{u}_{z1} + \tilde{u}_{z1} \mathrm{e}^{2{\rm i}\theta} - \tilde{u}_{z1}^{*} \mathrm{e}^{{-}2{\rm i}\theta} \right) + \frac{1}{\widetilde{Re}}\, \frac{\partial ^{2}\zeta_1}{\partial \tilde{r}^{2}} . \end{gather}$$

Equations (4.44a) and (4.44b) are identical to (4.19) except for the presence of the terms involving $\varOmega _1$. They can be combined to give

(4.45)\begin{equation} \frac{\partial \tilde{u}_{z1}}{\partial T} + {\rm i} \varOmega_c^{\prime}\tilde{r} \tilde{u}_{z1} + {\rm i} \varOmega_1 \tilde{u}_{z1} = \frac{i}{4}r_c \varOmega_c + \frac{1}{2\widetilde{Re}} \left(1 + \frac{1}{Sc} \right) \frac{\partial ^{2}\tilde{u}_{z1}}{\partial \tilde{r}^{2}}, \end{equation}

whereas (4.44c) splits into

(4.46a)$$\begin{gather} \frac{\partial \zeta_{1}}{\partial T} ={-}\frac{{\rm i}}{2} \frac{\partial }{\partial \tilde{r}} \left( \tilde{u}_{z1}^{*} - \tilde{u}_{z1} \right) + \frac{1}{\widetilde{Re}}\, \frac{\partial ^{2}\zeta_1}{\partial \tilde{r}^{2}}, \end{gather}$$

(4.46b)$$\begin{gather}\varOmega_c \frac{\partial \zeta_2}{\partial \theta} ={-}\frac{{\rm i}}{2} \frac{\partial }{\partial \tilde{r}} \left( \tilde{u}_{z1} \mathrm{e}^{2{\rm i}\theta} - \tilde{u}_{z1}^{*} \mathrm{e}^{{-}2{\rm i}\theta} \right). \end{gather}$$

Since $\zeta _1 = r_c \partial \varOmega _1 / \partial \tilde {r}$, (4.45)–(4.46a) form a closed system of equations. An equation for $\varOmega _1$ can be obtained by integrating (4.46a) as follows:

(4.47)\begin{equation} \frac{\partial \varOmega_{1}}{\partial T} ={-} \frac{{\rm i}}{2r_c} \left( \tilde{u}_{z1}^{*} - \tilde{u}_{z1} \right) + \frac{1}{\widetilde{Re}} \,\frac{\partial^{2}\varOmega_1}{\partial \tilde{r}^{2}}. \end{equation}

Equation (4.46b) shows that the azimuthal wavenumbers $m=\pm 2$ are generated in the vorticity only at higher order explaining again why the vorticity remains quasi-axisymmetric in the DNS during the first two phases. For this reason, components of the form $\exp (\pm 3 {\rm i}\theta )$ arise in the vertical velocity and density (4.39a)–(4.39b) only at the higher order $\varepsilon ^{5/3}$.

5.1. Vertical velocity

Figure 9 shows the evolution of the maximum vertical velocity $u_{zm}(\theta,t)$ (solid line) for $\theta =0$ (figure 9a) and $\theta ={\rm \pi} /2$ (figure 9b) for the set of parameters $Re=10\,000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$). Two different angles are considered since the theoretical vertical velocity has the form $u_z = u_{zc} (r,t) \cos {(\theta )} + u_{zs} (r,t) \sin {(\theta )}$. Hence, the plots for $\theta =0$ and $\theta ={\rm \pi} /2$ allow us to check the predictions for $u_{zc}$ and $u_{zs}$, respectively.

Figure 9. Comparison between the maximum vertical velocity in the DNS (black solid line), predicted by the unsteady inviscid solution (4.10a) (red dashed line), by the unsteady viscous solution (4.16a), (4.21) (yellow dashed line), by the steady viscous solution (4.16a), (4.23) (green dashed line) and by the nonlinear equations (4.45), (4.47) (blue dashed line) for $(a)$ $\theta =0$ and $(b)$ $\theta ={\rm \pi} /2$ for $Re=10\,000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ and for $(c)$ $\theta =0$ and $(d)$ $\theta ={\rm \pi} /2$ for $Re=2000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$).

We can see that the inviscid theoretical solution (4.10a) (red dashed line in figure 9a,b) predicts very well the initial linear increase of the maximum vertical velocity in the DNS for both angles. The unsteady viscous solution (4.16a), (4.21) (yellow dashed line) increases also linearly initially and is in good agreement with the DNS except that it lacks the small oscillations. They are indeed due to inertia-gravity waves oscillating at frequency $1/F_h + \varOmega = 2/F_h$ near $r_c$ (see the last term of (4.11a)) and of the two lines of (4.13) which are neglected in (4.14) and § 4.2.

When the growth of $u_{zm}$ is no longer linear, the time-dependent viscous solution (4.21) remains in very good agreement with the DNS and describes perfectly the transition towards the steady viscous solution (4.23) (green dashed line). The latter solution is close to the levels of saturation of $u_{zm}(\theta,t)$ for $\theta =0$ and $\theta = {\rm \pi}/2$ in the DNS, although the agreement is not as good as for the initial regime. However, the nonlinear equations (4.45), (4.47) (blue dashed lines) are in better agreement with the DNS indicating that nonlinear effects are also active in the saturation. In essence, none of the theoretical solutions can exhibit oscillations associated with the late non-axisymmetric evolution observed in the DNS.

Figure 9(c,d) shows a similar comparison when the Reynolds number is reduced to ${Re=2000}$, keeping the other parameters fixed. In this case, $u_{zm}(\theta =0,t)$ and $u_{zm}(\theta ={\rm \pi} /2,t)$ do not oscillate at late time in the DNS (black solid line). The agreement with (4.10a) and (4.21) or (4.23) are then excellent in the initial and saturation regimes, respectively. In addition, we see that the predictions of the nonlinear equations (4.45), (4.47) remain very close to the linear ones, i.e. (4.21), showing that nonlinear effects are weak in this case.

Figure 9 shows that the transition from the unsteady inviscid solution (4.10a) to the steady viscous one, (4.23), occurs in a short time range. Therefore, the unsteady viscous solution (4.21) can be well approximated by (4.10a) for $t \leqslant \mathcal {T}$ and (4.23) for $t \geqslant \mathcal {T}$, where

(5.1)\begin{equation} \mathcal{T} = 2 {\rm \pi}\frac{Re^{1/3}{Hi}(0)}{\left|{2 \varOmega_c^{\prime}}\right|^{2/3}(1+1/Sc)^{1/3}} \end{equation}

is the time when the overall maximum given by (4.10a) and (4.23) becomes equal. This time is independent of $\widetilde {Ro}$ and depends only on the Reynolds number, the Froude number via $\varOmega _c^{\prime }$ and the Schmidt number.

Figure 10 displays a detailed comparison of the radial profile of $u_z$ for $\theta =0$ and $\theta ={\rm \pi} /2$ predicted by (4.10a) (red dashed line) and observed in the DNS (black solid line) for different instants in the inviscid regime, i.e. $t \leqslant \mathcal {T}$ for the same parameters as figure 9(a,b). We see that the agreement is excellent even when $t \simeq \mathcal {T}$ (figure 10c, f). The approximation (4.21) and the solution of the nonlinear equations (4.45), (4.47) for the vertical velocity are also represented by yellow and blue dashed lines, respectively, in figure 10(b,c,e, f). The approximation (4.14) is almost identical to (4.21) in this time range and not represented. As expected, the agreement between the DNS and (4.21) or (4.45), (4.47) is very good near $r_c$ but deteriorates away from $r_c$. We can note that the blue and yellow dashed lines are superposed everywhere except close to the critical radius $r_c$ for $t=40$ (figure 10c, f). In this region the blue dashed lines are in very good agreement with the DNS indicating that nonlinear effects are important there. However, away from $r_c$, it is (4.10a) (red dashed line) which better agrees with the DNS. For $t=5$ (figure 10a,d), the approximations (4.21) and (4.45), (4.47) are not accurate and not shown. This is because the profile of $u_z$ is not yet sufficiently localized around $r_c$ at this early time and, therefore, it cannot be well described by a local approximation near $r_c$. Indeed, we can see in figure 10(a,d) that the profiles of $u_z$ are quite different from those in figure 10(b,c,e, f).

Figure 10. Comparison between the vertical velocity at $\theta =0$ (a–c) and $\theta ={\rm \pi} /2$ (d–f) in the DNS (black solid line), predicted by the unsteady inviscid solution (4.10a) (red dashed line), by the unsteady viscous solution (4.21) (yellow dashed line) and by the nonlinear equations (4.45), (4.47) (blue dashed line) at (a,d) $t=5$, (b,e) $t=25$ and (c, f) $t=40$ for $Re=10\,000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$). The circle symbols represent the location of the critical radius in the unsteady inviscid solution (4.10a).

Figure 11 shows again the vertical velocity profiles observed in the DNS for the same parameters but, for $t \geqslant \mathcal {T}$, this time they are compared with the unsteady and steady viscous solutions (4.21) (yellow dashed lines) and (4.23) (green dashed lines) as well as the predictions of the nonlinear equations (4.45), (4.47) (blue dashed lines). At $t=40$ (figure 11a,d), the steady viscous solution (4.23) (green dashed lines) is already in good agreement with the DNS since $t=40$ is close to the time $\mathcal {T}$ where the transition from (4.10a) to (4.23) occurs (figure 9a,b). Nevertheless, it departs slightly from the DNS near $r_c$ unlike the unsteady viscous solution (4.21) and nonlinear predictions from (4.45), (4.47). At longer times, $t=60$ (figure 11b,e) and $t=80$ (figure 11c, f), (4.21) and (4.23) become almost identical and remain in satisfactory agreement with the vertical velocity profiles of the DNS. Nevertheless, we can see a shift between the numerical and theoretical profiles. In contrast, the solution of the nonlinear equations (4.45), (4.47) remains in very good agreement with the DNS and does not exhibit such a shift. This indicates that the shift is due to the viscous and nonlinear variations of the angular velocity. For example, this makes the critical radius where $\varOmega (r_c) = 1/F_h$ to move towards the vortex centre, as seen in figure 11. This phenomenon is absent from the linear equations (4.21), (4.23) since they do not take into account any variation of the angular velocity.

Figure 11. Comparison between the vertical velocity at $\theta =0$ (a–c) and $\theta ={\rm \pi} /2$ (d–f) in the DNS (black solid line), predicted by the unsteady viscous solution (4.16a), (4.21) (yellow dashed line), by the steady viscous solution (4.16a), (4.23) (green dashed line) and by the nonlinear equations (4.45), (4.47) (blue dashed line) at (a,d) $t=40$, (b,e) $t=60$, (c, f) $t=80$ for $Re=10\,000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$). The circle symbols represent the location of the critical radius in the unsteady viscous solution (4.16a).

5.2. Vertical vorticity

The asymptotic analyses have shown that the axisymmetric component of the vertical velocity is given by $\zeta = \zeta _0(r) + \varepsilon ^{2} \zeta _{20}(r,t) + \dots$, where $\zeta _{20}$ follows (4.31) and (4.35) in the inviscid and viscous regimes, i.e. $t \leqslant \mathcal {T}$ and $t\geqslant \mathcal {T}$, respectively. A global view of these two regimes and the associated approximations can be gained by plotting $\partial \zeta _{20} / \partial r$ at $r=r_c$ (figure 12). The black solid line shows the evolution of $\partial \zeta _{20} / \partial r (r_c,t)$ computed numerically from (4.35). Here $\partial \zeta _{20} / \partial r (r_c,t)$ increases initially like $t^{4}$ in agreement with the approximation (4.33) (red dashed line). Subsequently, for $t \gg \mathcal {T}$, $\partial \zeta _{20} / \partial r (r_c,t)$ increases more slowly and saturates towards $\partial \zeta _{20}^{(1)} / \partial r (r_c)$ (black dashed line) for $t \rightarrow \infty$. The approximation (4.36) (blue dashed line) is in good agreement with the solution (4.35) (black line) in this regime. It will allow us to derive a theoretical criteria for the onset of non-axisymmetry in § 6.4.

Figure 12. Comparison between the evolution of $\partial \zeta _{20}/\partial r(r_c,t)$ from the theoretical expressions: the unsteady viscous solution (4.35) (black solid line), the unsteady inviscid solution (4.31) (red dashed line), the viscous solutions (4.36) (blue dashed line) and (4.37a) (black dashed line) for $Re=10\,000$, $F_h=2$.

Figure 13 compares the radial profile of the theoretical vorticity $\zeta = \zeta _0(r) + \varepsilon ^{2} \zeta _{20}(r,t) + \dots$, with $\zeta _{20}$ given by the unsteady inviscid solution (4.31) (red dashed line) to the DNS when $t \leqslant \mathcal {T}$. The agreement is excellent and predicts very well the deformation of the vertical vorticity profile near $r_c \simeq 1.2$. A similar comparison is displayed in figure 14 for three different times such that $t \gtrsim \mathcal {T}$ and $\zeta _{20}$ is given by the unsteady viscous solution (4.35). The agreement continues to be very good even for $t=85$ (figure 14c), confirming the validity of the solution (4.35). A slight shift can however be noticed near $r_c$ and near the vortex axis. In contrast, the predictions of the nonlinear equations (4.45), (4.47) (blue dashed lines) are in very good agreement with the DNS near $r_c$. This shows again that nonlinear effects are not negligible in the critical layer.

Figure 13. Comparison between the vertical vorticity at $\theta ={\rm \pi} /2$ in the DNS (black solid line) and the asymptotic expressions $\zeta = \zeta _0 + \varepsilon ^{2}\zeta _{20}$ where $\zeta _{20}$ follows the unsteady inviscid solution (4.31) (red dashed line) and $\zeta = \zeta _0 + \zeta _1$ with $\zeta _1$ given by the nonlinear equations (4.45), (4.47) (blue dashed line) at $(a)$ $t=25$, $(b)$ $t=35$ and $(c)$ $t=40$ for $Re=10\,000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$). The circle symbols represent the location of the critical radius in the unsteady inviscid solution (4.31).

Figure 14. Comparison between the vertical vorticity at $\theta ={\rm \pi} /2$ in the DNS (black solid line) and the asymptotic expressions $\zeta = \zeta _0 + \varepsilon ^{2}\zeta _{20}$ where $\zeta _{20}$ follows the unsteady viscous solution (4.35) (red dashed line) and $\zeta = \zeta _0 + \zeta _1$ with $\zeta _1$ given by the nonlinear equations (4.45), (4.47) (blue dashed line) at $(a)$ $t=50$, $(b)$ $t=65$ and $(c)$ $t=85$ for $Re=10\,000$, $F_h=2$, $Ro=20.3$, $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$). The circle symbols represent the location of the critical radius in the unsteady inviscid solution (4.35).

Finally, the prediction for late time $t \gg \mathcal {T}$ based on the steady viscous solution (4.36) has been compared with a DNS for a larger non-traditional Rossby number $\widetilde {Ro} = 500$, the other parameters being identical. The vortex remains then quasi-axisymmetric in the DNS (figure 8) allowing us to see the late evolution of the vorticity without the non-axisymmetric perturbations. Figure 15 shows that the deformation of the vorticity profile near $r_c$ is weak but well predicted by (4.36), although there is a shift for $t \geqslant 300$. In this case, the predictions from the nonlinear equations (4.45), (4.47) (not shown) are identical to those from (4.36) indicating that the nonlinear effects are weak. However, if we take also into account the global viscous decay of the leading-order vorticity $\zeta _0$ which is significant for these large times, then the agreement with the DNS is perfect (yellow dashed line).

Figure 15. Comparison between the vertical vorticity at $\theta ={\rm \pi} /2$ in the DNS (black solid line) and the asymptotic expressions $\zeta = \zeta _0 + \varepsilon ^{2}\zeta _{20}$ where $\zeta _{20}$ follows the viscous solution (4.36) (red dashed line) and $\zeta = \zeta _0 + \zeta _1$ with $\zeta _1$ given by the nonlinear equations (4.45), (4.47) and the viscous decay of $\zeta _0$ also taken into account (yellow dashed line) at $(a)$ $t=200$, $(b)$ $t=300$ and $(c)$ $t=400$ for $Re=10\,000$, $F_h=2$, $Ro=20.0$, $\phi =87.7^{\circ }$ ($\widetilde {Ro}=500$). The circle symbols represent the location of the critical radius in the viscous solution (4.36).

Besides, a feature of high interest is that the vertical vorticity profile exhibits two extrema near $r_c$, ${\rm d}\zeta / {\rm d}r = 0$, as soon as $t \geqslant 40$ (figures 13(c) and 14). According to Rayleigh's inflectional criterion (Rayleigh Reference Rayleigh1880), this is a necessary condition for the shear instability. However, only the first extremum, where $\zeta$ has a local minimum, satisfies the stiffer (Fjørtoft Reference Fjørtoft1950) instability condition $[\varOmega (r)-\varOmega (r_I)] \partial \zeta / \partial r < 0$, where $r_I$ is the extremum. This gives us hindsight on the possible origin of the late non-axisymmetric evolution of the vertical vorticity. The next section will investigate whether or not this hypothesis is correct.

6.1. Azimuthal decomposition of $\zeta$ and $u_z$

In order to better understand the onset of non-axisymmetric vertical vorticity, we have first decomposed $\zeta$ thanks to an azimuthal Fourier transform

(6.1)\begin{equation} \hat{\zeta}_m (r,t) = \int_0^{2{\rm \pi}} \zeta (r,\theta,t) {\rm e}^{-{\rm i}m\theta} \,\mathrm{d}\theta, \end{equation}

where $\zeta (x,y,t)$ has been first interpolated on a grid of cylindrical coordinates $(r,\theta )$. The same transform has been applied to $u_z$ in order to obtain $\hat {u}_{z_m}$. The mean power in each azimuthal wavenumber is then defined as

(6.2)\begin{equation} E_{\zeta}(m,t)=\frac{ \displaystyle\int_0^{l_x/2}\hat{\zeta}_m^{2}(r,t) r\,\mathrm {d}r}{\displaystyle\int_0^{l_x/2} r\,\mathrm {d}r} = \frac{8}{l_x^{2}} \int_0^{l_x/2}\hat{\zeta}_m^{2}(r,t) r\,\mathrm {d}r,\end{equation}

and

(6.3)\begin{equation} E_{u_z}(m,t)= \frac{8}{l_x^{2}} \int_0^{l_x/2} \hat{u}_{z_m}^{2}(r,t) r\,\mathrm {d}r. \end{equation}

Figure 16(a) displays the evolution of the logarithm of the power of the first three azimuthal wavenumbers of the vertical vorticity. These are only even, i.e. $m=0$, $m=2$, $m=4$. Similarly, figure 16(b) shows the logarithm of the power of the first three azimuthal wavenumbers of $u_z$. They are odd in this case: $m=1$, $m=3$, $m=5$. It can be seen that $E_{\zeta }(0,t)$ and $E_{u_z}(1,t)$ (black solid lines) remain approximately constant except that $E_{u_z}(1,t)$ sustains large oscillations after $t \simeq 120$.

Figure 16. Evolution of the power $(a)$ $E_{\zeta }(m,t)$ for the azimuthal wavenumbers $m=0$ (black solid line), $m=2$ (red dashed line) and $m=4$ (green dashed line) and the power $(b)$ $E_{u_z}(m,t)$ for the azimuthal wavenumbers $m=1$ (black solid line), $m=3$ (red dashed line) and $m=5$ (green dashed line) for $Re=10\,000$, $F_h=2$, $Ro=20.3$ and $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$).

It is also worth pointing out that $E_{\zeta }(2,t)$ starts to grow at the beginning of the simulation (figure 16a) since $\zeta _{22}$ increases like $t^{2}$ due to the forcing by the vertical velocity as detailed in § 4.3.1. However, when $t \leqslant 80$, its power remains negligible compared with that of the axisymmetric mode, $E_\zeta (0,t)$. In contrast, after $t \simeq 80$, $E_{\zeta }(2,t)$ grows exponentially before saturating at $t \geqslant 120$. There is therefore a clear transition towards an exponential growth, a feature consistent with the instability hypothesis. Nevertheless, we can see that the azimuthal mode $m=3$ of $u_z$ (figure 16b) grows also exponentially at the same time. The higher modes, $m=4$ of $\zeta$ and $m=5$ of $u_z$, start to increase also exponentially but somewhat later. Thus, it is unclear if the growth of $E_{u_z}(3,t)$ is a consequence of the growth of $E_{\zeta }(2,t)$, if it is the opposite or if the exponential growth is due to a coupling between $E_{u_z}(3,t)$ and $E_{\zeta }(2,t)$.

6.2. Truncated model

To answer the latter question, we have derived a truncated model taking into account only the first azimuthal wavenumbers of each quantity. More precisely, the different variables have been written as

(6.4a)$$\begin{gather} u_r = \hat{u}_{r_{2c}} (r,t) \cos{(2\theta)} + \hat{u}_{r_{2s}} (r,t) \sin{(2\theta)}, \end{gather}$$

(6.4b)$$\begin{gather}u_\theta = \hat{u}_{\theta_0} (r,t) + \hat{u}_{\theta_{2c}} (r,t) \cos{(2\theta)} + \hat{u}_{\theta_{2s}} (r,t) \sin{(2\theta)}, \end{gather}$$

(6.4c)$$\begin{gather}p (r,\theta) = \hat{p}_0 (r,t) + \hat{p}_{2c} (r,t) \cos{(2\theta)} + \hat{p}_{2s} (r,t) \sin{(2\theta)}, \end{gather}$$

(6.4d)$$\begin{gather}\zeta = \hat{\zeta}_0 (r,t) + \hat{\zeta}_{2c} (r,t) \cos{(2\theta)} + \hat{\zeta}_{2s} (r,t) \sin{(2\theta)}, \end{gather}$$

(6.4e)$$\begin{gather}u_z = \hat{u}_{z_{1c}} (r,t) \cos{(\theta)} + \hat{u}_{z_{1s}} (r,t) \sin{(\theta)}, \end{gather}$$

(6.4f)$$\begin{gather}b = \hat{b}_{1c} (r,t) \cos{(\theta)} + \hat{b}_{1s} (r,t) \sin{(\theta)}. \end{gather}$$

These decompositions have been introduced in (4.1a)–(4.1e) and the following governing equations have been obtained for the vertical velocity and buoyancy by truncating all the higher modes:

(6.5a)\begin{align} \frac{\partial \hat{u}_{z_{1c}}}{\partial t} &={-}\frac{1}{2} \left ( \hat{u}_{r_{2c}} \frac{\partial \hat{u}_{z_{1c}}}{\partial r} + \hat{u}_{r_{2s}} \frac{\partial \hat{u}_{z_{1s}}}{\partial r} \right) - \frac{1}{r} \left (\hat{u}_{\theta_0} \hat{u}_{z_{1s}} + \frac{1}{2} ( \hat{u}_{\theta_{2c}} \hat{u}_{z_{1s}} - \hat{u}_{\theta_{2s}} \hat{u}_{z_{1c}} ) \right)\nonumber\\ &\quad + \hat{b}_{1c} + \frac{1}{Re} \nabla ^{2} \hat{u}_{z_{1c}}+ \frac{1}{\widetilde{Ro}}( \hat{u}_{r_{2c}} - \hat{u}_{\theta_{2s}} ), \end{align}

(6.5b)\begin{align} \frac{\partial \hat{u}_{z_{1s}}}{\partial t} &= \frac{1}{2} \left ( \hat{u}_{r_{2c}} \frac{\partial \hat{u}_{z_{1s}}}{\partial r} - \hat{u}_{r_{2s}} \frac{\partial \hat{u}_{z_{1c}}}{\partial r} \right) + \frac{1}{r} \left ( \hat{u}_{\theta_0} \hat{u}_{z_{1c}} - \frac{1}{2} ( \hat{u}_{\theta_{2c}} \hat{u}_{z_{1c}} + \hat{u}_{\theta_{2s}} \hat{u}_{z_{1s}} ) \right) \nonumber\\ &\quad + \hat{b}_{1s} + \frac{1}{Re} \nabla ^{2} \hat{u}_{z_{1s}} + \frac{1}{\widetilde{Ro}} ( \hat{u}_{r_{2s}} + \hat{u}_{\theta_{2c}} -2\hat{u}_{\theta_0} ) , \end{align}

(6.5c)\begin{align} \frac{\partial \hat{b}_{1c}}{\partial t} &={-}\frac{1}{2} \left ( \hat{u}_{r_{2c}} \frac{\partial \hat{b}_{1c}}{\partial r} + \hat{u}_{r_{2s}} \frac{\partial \hat{b}_{1s}}{\partial r} \right) - \frac{1}{r} \left( \hat{u}_{\theta_0} \hat{b}_{{1s}} + \frac{1}{2} ( \hat{u}_{\theta_{2c}} \hat{b}_{1s} - \hat{u}_{\theta_{2s}} \hat{b}_{1c} ) \right )\nonumber\\ &\quad - \frac{\hat{u}_{z_{1c}}}{F_h^{2}} + \frac{1}{ReSc} \nabla ^{2} \hat{b}_{{1c}}, \end{align}

(6.5d)\begin{align} \frac{\partial \hat{b}_{{1s}}}{\partial t} &= \frac{1}{2} \left ( \hat{u}_{r_{2c}} \frac{\partial \hat{b}_{1s}}{\partial r} - \hat{u}_{r_{2s}} \frac{\partial \hat{b}_{1c}}{\partial r} \right) + \frac{1}{r} \left ( \hat{u}_{\theta_0} \hat{b}_{{1c}} - \frac{1}{2} ( \hat{u}_{\theta_{2c}} \hat{b}_{1c} + \hat{u}_{\theta_{2s}} \hat{b}_{1s} ) \right) \nonumber\\ &\quad - \frac{\hat{u}_{z_{1s}}}{F_h^{2}} + \frac{1}{ReSc} \nabla ^{2} \hat{b}_{{1s}}. \end{align}

Applying the same truncation approach for the vertical vorticity gives

(6.6a)\begin{align} \frac{\partial \hat{\zeta}_0}{\partial t} &={-} \frac{1}{2} \left( \hat{u}_{r_{2c}} \frac{\partial \hat{\zeta}_{2c}}{\partial r} + \hat{u}_{r_{2s}} \frac{\partial \hat{\zeta}_{2s}}{\partial r} \right)- \frac{1}{r} \left( \hat{u}_{\theta_{2c}} \hat{\zeta}_{2s} - \hat{u}_{\theta_{2s}} \hat{\zeta}_{2c} \right) + \frac{1}{Re}\nabla ^{2} \hat{\zeta}_0 \nonumber\\ &\quad + \frac{1}{ \widetilde{Ro}} \left( \frac{\partial \hat{u}_{z_{1s}}}{\partial r} + \frac{\hat{u}_{z_{1s}}}{r} \right), \end{align}

(6.6b)$$\begin{align} \frac{\partial \hat{\zeta}_{2c}}{\partial t} &={-} \hat{u}_{r_{2c}} \frac{\partial \hat{\zeta}_0}{\partial r} - \frac{2 \hat{u}_{\theta_0}}{r} \hat{\zeta}_{2s} + \frac{1}{ Re} \nabla ^{2} \hat{\zeta}_{2c} - \frac{1}{\widetilde{Ro}} \left( \frac{\partial \hat{u}_{z_{1s}}}{\partial r} -\frac{\hat{u}_{z_{1s}}}{r} \right), \end{align}$$

(6.6c)$$\begin{align}\frac{\partial \hat{\zeta}_{2s}}{\partial t} &={-} \hat{u}_{r_{2s}} \frac{\partial\hat{\zeta}_0}{\partial r} + \frac{2 \hat{u}_{\theta_0}}{r} \hat{\zeta}_{2c} + \frac{1}{Re} \nabla ^{2} \hat \zeta_{2s} + \frac{1}{\widetilde{Ro}}\left( \frac{\partial \hat{u}_{z_{1c}}}{\partial r} -\frac{\hat{u}_{z_{1c}}}{r} \right), \end{align}$$

where

(6.7a)$$\begin{gather} \hat{\zeta}_0 = \frac{1}{r} \,\frac{\partial r \hat{u}_{\theta_0} }{\partial r}, \end{gather}$$

(6.7b)$$\begin{gather}\hat{\zeta}_{2c} = \frac{1}{r} \,\frac{\partial r \hat{u}_{\theta_{2c}} }{\partial r} - \frac{2}{r} \hat{u}_{r_{2s}}, \end{gather}$$

(6.7c)$$\begin{gather}\hat{\zeta}_{2s} = \frac{1}{r} \,\frac{\partial r \hat{u}_{\theta_{2s}} }{\partial r} + \frac{2}{r} \hat{u}_{r_{2c}}. \end{gather}$$

The divergence equation also implies that

(6.8a)$$\begin{gather} \frac{1}{r} \,\frac{\partial r \hat{u}_{r_{2c}} }{\partial r} + \frac{2}{r} \hat{u}_{\theta_{2s}} = 0, \end{gather}$$

(6.8b)$$\begin{gather}\frac{1}{r} \,\frac{\partial r \hat{u}_{r_{2s}} }{\partial r} - \frac{2}{r} \hat{u}_{\theta_{2c}} = 0. \end{gather}$$

Such a truncated model can be seen as a heuristic extension of the asymptotic analyses. Indeed, it takes into account both time dependence and diffusive effects in the evolution of the vertical velocity and buoyancy (6.5). Moreover, the modifications of the axisymmetric flow field ($\hat {u}_{\theta _0}$) and the generated $m=2$ mode ($\hat {u}_{r_{2c}}$, $\hat {u}_{r_{2s}}$, $\hat {u}_{\theta _{2c}}$, $\hat {u}_{\theta _{2s}}$) are also taken into account. The latter is governed by (6.6b)–(6.6c) and it appears also in the evolution of the axisymmetric flow field (6.6a). However, as in the asymptotic analyses, only the first azimuthal wavenumber of $u_z$ and $b$ are considered, i.e. the modes $m=3$, $m=5,\ldots$ are neglected. Similarly, the higher modes $m=4$, $m=6,\ldots$ are neglected in the evolution of the horizontal flow field (6.6).

Figure 17(a,b) compares the power $E_{\zeta } (2,t)$ and $E_{u_{z}} (1,t)$, respectively, obtained in the DNS and in the truncated model. There are some differences for $t \geqslant 100$ but, qualitatively, the same type of evolution as in the DNS is obtained with the truncated model. This is remarkable since the truncated model crudely neglects many azimuthal modes and, in particular, the azimuthal mode $m=3$ in the vertical velocity and buoyancy fields. Hence, this proves that the growth of $E_{u_z} (3,t)$ in figure 16(b) does not play a key role in the onset of non-axisymmetry in the vertical vorticity.

Figure 17. Evolution of the power $(a)$ $E_{\zeta }(2,t)$ and $(b)$ $E_{u_z}(1,t)$ in the DNS (black solid line) and truncated model (red dashed line) for $Re=10\,000$, $F_h=2$, $Ro=20.3$ and $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$).

We can take a step further in the understanding of this phenomenon by freezing the axisymmetric velocity field. In other words, the time evolution of $(\hat {u}_{\theta _0},\hat {\zeta }_0)$ is suppressed after a given time $t_f$. We also set the forcing terms due to the non-traditional Coriolis force to be zero in (6.6b) and (6.6c). Hence, the latter equations become

(6.9a)$$\begin{gather} \frac{\partial \hat{\zeta}_{2c}}{\partial t} ={-} \hat{u}_{r_{2c}} \frac{\partial \hat{\zeta}_0 (r,t_f)}{\partial r} - \frac{2 \hat{u}_{\theta_0}(r,t_f)}{r} \hat{\zeta}_{2s} + \frac{1}{Re} \nabla ^{2} \hat{\zeta}_{2c}, \end{gather}$$

(6.9b)$$\begin{gather}\frac{\partial \hat{\zeta}_{2s}}{\partial t} ={-} \hat{u}_{r_{2s}} \frac{\partial \hat{\zeta}_0 (r,t_f) }{\partial r} + \frac{2 \hat{u}_{\theta_0}(r,t_f)}{r} \hat{\zeta}_{2c} + \frac{1}{Re}\nabla ^{2} \hat{\zeta}_{2s}. \end{gather}$$

These equations describe simply the linear evolution of perturbations with azimuthal wavenumber $m=2$ on a steady axisymmetric vortex with azimuthal velocity $\hat {u}_{\theta _0}(r,t_f)$. The perturbations $(\hat {u}_{r_{2c}},\hat {u}_{\theta _{2c}})$ and $(\hat {u}_{r_{2s}},\hat {u}_{\theta _{2s}})$ are initialized by a white noise whose amplitude is adjusted so as to have a power of the same order as $E_{\zeta }(2,t_f)$. Figure 18 shows the evolution of the power $E_{\zeta }(2,t)$ for different freezing times $t_f$ ($t_f=40$, yellow dashed line; $t_f=50$, blue dashed line; $t_f=65$, red dashed line; $t_f=85$, green dashed line) compared with the evolution of $E_{\zeta }(2,t)$ in the truncated model (black solid line). Strikingly, we see that $E_{\zeta }(2,t)$ grows also exponentially regardless of the value of $t_f$ investigated. Furthermore, the growth rate, i.e. the slope, increases with $t_f$. Most interestingly, the growth rate for $t_f=85$ is close to that observed in the truncated model (black solid line). This demonstrates that the onset of non-axisymmetry in the vertical vorticity is due to an instability of the vortex profile. When the anomaly of vertical vorticity is sufficient to have an extremum, a shear instability with an azimuthal wavenumber $m=2$ develops. Subsequently, this triggers the growth of higher azimuthal modes through coupling with the non-traditional Coriolis force.

Figure 18. Evolution of the power $E_{\zeta }(2,t)$ in the truncated model (black solid line) or using (6.9) for different freezing times: $t_f=85$ (green dashed line), $t_f=65$ (red dashed line), $t_f=50$ (blue dashed line) and $t_f=40$ (yellow dashed line) for $Re=10\,000$, $F_h=2$, $Ro=20.3$ and $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$).

6.3. Equivalent vortex with piecewise uniform vorticity

A simple model of the instability can be obtained by considering the inviscid limit and by using a vortex with piecewise uniform vorticity with four concentric regions, as considered by Carton & Legras (Reference Carton and Legras1994) and Kossin, Schubert & Montgomery (Reference Kossin, Schubert and Montgomery2000). As shown by two examples in figure 19, the vorticity profile in the DNS can be crudely approximated by four levels of constant vorticity, i.e.

(6.10)\begin{equation} \zeta = \left\{ \begin{array}{ll} \zeta_1 = 2, & 0 < r < r_1, \\ \zeta_2 = \zeta_c - \delta_v/2, & r_1 < r < r_2, \\ \zeta_3 = \zeta_c + \delta_v/2, & r_2 < r < r_3, \\ \zeta_4 = 0, & r_3 < r, \end{array} \right. \end{equation}

where $r_1=r_c - \delta _h$, $r_2=r_c$ and $r_3= r_c + \delta _h$, with $\delta _v$ and $\delta _h$ the amplitude and size of the vorticity anomaly in the vicinity of the critical radius $r_c$. More explicitly, $\delta _v$ is the difference between the local maximum and minimum of the vorticity and $\delta _h$ is the distance between these two extrema. The corresponding angular velocity of the vortex is continuous and given by

(6.11)\begin{gather} \varOmega (r) = \tfrac{1}{2}\left\{ \begin{array}{ll} \zeta_1, & 0 < r < r_1, \\ \zeta_2 - (\zeta_2 -\zeta_1)(r_1/r)^{2}, & r_1 < r < r_2, \\ \zeta_3 - (\zeta_2 -\zeta_1)(r_1/r)^{2} - (\zeta_3 -\zeta_2)(r_2/r)^{2}, & r_2 < r < r_3, \\ - (\zeta_2 -\zeta_1)(r_1/r)^{2} - (\zeta_3 -\zeta_2)(r_2/r)^{2} + \zeta_3 (r_3/r)^{2}, & r_3 < r. \end{array} \right. \end{gather}

For a given Froude number $F_h$, the position of the critical radius $r_c$ and the value of $\zeta _c$ are fixed. Hence, the problem has only two control parameters: $\delta _v$ and $\delta _h$. The stability of such a vortex with respect to perturbations of the form $\psi \exp ({{\rm i}m\theta +\sigma t})$ is governed by the eigenvalue problem (Carton & Legras Reference Carton and Legras1994; Kossin et al. Reference Kossin, Schubert and Montgomery2000)

(6.12)\begin{align} &\begin{pmatrix} m \varOmega(r_1) + \frac{1}{2} (\zeta_2 - \zeta_1) & \frac{1}{2} (\zeta_2 - \zeta_1)(r_1/r_2)^{m} & \frac{1}{2} (\zeta_2 - \zeta_1)(r_1/r_3)^{m}\\ \frac{1}{2} (\zeta_3 - \zeta_2)(r_1/r_2)^{m} & m \varOmega(r_2) + \frac{1}{2} (\zeta_3 - \zeta_2) & \frac{1}{2} (\zeta_3 - \zeta_2)(r_2/r_3)^{m} \\ -\frac{1}{2} \zeta_3 (r_1/r_3)^{m} & -\frac{1}{2} \zeta_3 (r_2/r_3)^{m} & m \varOmega(r_3) - \frac{1}{2} \zeta_3 \end{pmatrix} \begin{pmatrix} \psi_1 \\ \psi_2 \\ \psi_3 \end{pmatrix} \nonumber\\ &\quad = \sigma \begin{pmatrix} \psi_1 \\ \psi_2 \\ \psi_3 \end{pmatrix} . \end{align}

Figure 19. Examples of the piecewise uniform vorticity (red line) fitting the continuous vertical vorticity profiles (black line) at $(a)$ $t=40$, $(b)$ $t=50$ for $Re=10\,000$, $F_h=2$, $Ro=20.3$ and $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$). The circle symbols represent the location of the critical radius.

Figure 20(a) shows the growth rate contours for the $m=2$ perturbations for $F_h=2$ as a function of $(\delta _v,\delta _h)$. We see that the growth rate is positive only when $\delta _v$ and $\delta _h$ are sufficiently away from zero in the ranges investigated. The symbols in figure 20(a) indicate the parameters $(\delta _v,\delta _h)$ estimated by fitting (6.10) to the vorticity field $\hat {\zeta }_0(r,t_f)$ at different times $t_f$ for $F_h=2$, for two different latitudes $\phi =80^{\circ }$ (red circles) and $\phi =75^{\circ }$ (black squares). For example, for $\phi =80^{\circ }$, the time $t_f$ varies from $t_f=45$ (leftmost point) to $t_f=85$ (rightmost point). The size of the vorticity anomaly $\delta _h$ does not vary very much and is around $\delta _h \simeq 0.2$ for both latitudes. In contrast, the amplitude of the anomaly $\delta _v$ increases with $t_f$ as expected. Figure 20(b) displays the growth rate as a function of $\delta _v$ (circle and square symbols). The dashed lines show the corresponding growth rate computed from (6.9), i.e. by considering the continuous vorticity profile $\hat {\zeta }_0(r,t_f)$ and with the perturbations initialized by white noise. This shows that the piecewise vortex model is able to predict quite well the growth rate of the instability observed in the truncated model, which is itself in good agreement with the DNS.

Figure 20. (a) Growth rate contours of the piecewise vortex model as a function of $\delta _h$ and $\delta _v$ for $F_h=2$. The contour interval is 0.03. The bold line indicates the growth rate $\sigma = 0$. The symbols correspond to the values of $\delta _h$ and $\delta _v$ estimated at different freezing times for $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$) (red circles) and $\phi =75^{\circ }$ ($\widetilde {Ro}=77.27$) (black squares) for $Re=10\,000$, $F_h=2$. (b) Growth rates as a function of $\delta _v$ obtained by the truncated model at different freezing times (symbols) and given by (6.12) for $m=2$ (dashed lines) for $\phi =75^{\circ }$ ($\widetilde {Ro}=77.27$) (black) and $\phi =80^{\circ }$ ($\widetilde {Ro}=115.2$) (red), $Re=10\,000$ and $F_h=2$.

Using (6.12), we have also computed the growth rate of higher azimuthal wavenumbers $m>2$. The results show that the most unstable wavenumber is not $m=2$ but is between $m=3$ and $m=5$ for the parameters indicated by the symbols in figure 20(a). Three reasons might explain the actual dominance of $m=2$. First, the velocity jumps in the piecewise vortex model could favour larger wavenumbers compared with a continuous vorticity profile. Second, we have seen from (4.27)–(4.28) that the non-traditional Coriolis force generates not only an axisymmetric vorticity at order $1/\widetilde {Ro}^{2}$ but also a vorticity field with an azimuthal wavenumber $m=2$. Figure 17(a) shows that the latter is weak before the onset of the instability. However, it is not zero and, therefore, this small amplitude could favour its dominance over more unstable higher wavenumbers whose initial amplitudes are much lower (see $m=4$ in figure 17a). Third, the vortex profile is continuously evolving with time while, in the stability problems (6.9) or (6.12), we have frozen this evolution. Hence, the $m=2$ wavenumber could be selected first when the vortex becomes slightly unstable. This early selection would then ensure its subsequent dominance even if it is no longer the most unstable wavenumber. Such an effect has been evidenced by Wang & Balmforth (Reference Wang and Balmforth2021) in their study of the evolution of the wavenumber selection as the critical layer becomes finer.

6.4. Theoretical criterion

Even if the Rayleigh–Fjørtoft criterion is only a necessary condition for the shear instability in inviscid fluids, we can try to use it to establish a theoretical criterion for the onset of the shear instability in the DNS for finite Reynolds numbers. Since the radial derivative of the vorticity is maximum at $r=r_c$ and is negative away from the critical radius, a necessary condition ensuring that there exists extrema, ${\rm d} \zeta / {\rm d} r = 0$, reads

(6.13)\begin{equation} \frac{{\rm d}\zeta (r_c,t_s)}{{\rm d}r} \geqslant c, \end{equation}

where ${c=0}$ is the minimum requirement for the existence of an inflection point, but we have explored also the consequences of larger values of $c$. Besides, the time $t_s$ will be set as $t_s=a\mathcal {T}$, where $a$ is a constant larger than unity. Indeed, the onset of the shear instability always occurs after the time $\mathcal {T}$. Therefore, the vorticity $\zeta$ will be taken as the asymptotic axisymmetric vorticity $\zeta = \zeta _0 + \varepsilon ^{2} \zeta _{20}$, where $\zeta _{20}$ is given by (4.36). As seen in figure 12, $\partial \zeta _{20}/\partial r(r_c,t)$ is indeed well predicted by (4.36) for large times. Then, (6.13) becomes

(6.14)\begin{equation} \frac{Re^{2/3}}{\widetilde{Ro}} \geqslant \left|{2 \varOmega_c^{\prime}}\right|^{1/3} \left (1+\frac{1}{Sc}\right)^{1/6} \sqrt{ \frac{\displaystyle { c - 3 \varOmega_c^{\prime} - r_c \varOmega_c^{\prime\prime}}}{\displaystyle { 2 {\rm \pi}r_c \varOmega_c \left( {Hi}(0) - \left(\frac{1+1/Sc}{8a{\rm \pi}^{2}{Hi}(0)}\right)^{1/2} \right)}}}. \end{equation}

Remarkably, the right-hand side depends only on the Froude number through $r_c$, the Schmidt number $S_c$ and the constants $a$ and $c$. The criterion (6.14) when ${c=0}$ and $a=\infty$ is represented by a solid line in figure 8. It delimits quite well the quasi-axisymmetric/non-axisymmetric domains observed in the DNS, except for the lowest Reynolds number investigated $Re=2000$. Such a difference for moderate Reynolds and Rossby numbers is not surprising since the asymptotics have been derived for a high Reynolds number and large Rossby number $\widetilde {Ro}$. Furthermore, viscous effects might damp the shear instability growth when the Reynolds number is moderate.

As seen from the piecewise vortex model, the shear instability for $m=2$ does not appear when $\delta _h \simeq 0.2$ as soon as $\delta _v >0$. It can be roughly estimated that the instability arises only when $\delta _v$ such that $\delta_v/\delta_h \gtrsim 0.4$ (figure 20a). Therefore, we can estimate $\partial \zeta /\partial r =$ $c \simeq \delta _v / \delta _h = 0.4$. The criterion (6.14) with this value of $c$ and $a = \infty$ is represented by a dashed line in figure 8. The agreement with the DNS is as good as the criterion (6.14) with ${c=0}$. The actual threshold is likely to be in between these two curves, i.e. in the hatched region (figure 8). Finally, we stress that the criterion (6.14) applies only to the shear instability due to an inflection point and not to other types of instability that may exist in viscous shear flows.