4.1. Flow geometry

Vortices are flows combining rotation and shear, and the majority of them possess a localized distribution of vorticity. The Lamb–Oseen vortex is perhaps the simpler solution of the Navier–Stokes equation containing these ingredients, thus is adopted in the rest of this paper. It expresses

(4.1a,b)\begin{equation} U_{{b},r}(r,t) = 0,\quad U_{{b},\theta}(r,t)=\frac{1-\exp({-r^2/(1+4t/ )})}{r}. \end{equation}

The associated vorticity field, $W_z$, has a Gaussian radial profile

(4.2)\begin{equation} W_z(r,t) = \frac{U_{{b},\theta}}{r} + \frac{\partial U_{{b},\theta}}{\partial _r} = \frac{2}{1+4t/{\textit {Re}}}\exp({-r^2/(1+4t/{\textit {Re}})}), \end{equation}

that diffuses with time for finite ${\textit {Re}}$ values. For the linear velocity perturbation, we select the azimuthal wavenumber $m=2$, since it is the wavenumber associated with the subcritical bifurcation towards the tripolar state described in the introduction (see for instance Rossi et al. Reference Rossi, Lingevitch and Bernoff1997). Only the response to an initial condition with $m=2$ will be considered in this paper.

4.2. Numerical methods

In practice, the Poisson equation (2.4) for the pressure is never solved directly, the latter being included in the state space instead:

(4.3)\begin{equation} \boldsymbol{\mathsf{M}}\partial_t\begin{bmatrix} \hat{\boldsymbol{u}}\\ \hat{p} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\mathsf{A}}_{m}(t) & -\hat{\boldsymbol{\nabla}}_{m} \\ \hat{\boldsymbol{\nabla}}_{m}\boldsymbol{\cdot} & 0 \end{bmatrix}\begin{bmatrix} \hat{\boldsymbol{u}}\\ \hat{p} \end{bmatrix} \doteq \boldsymbol{\mathsf{B}}_{m}(t)\begin{bmatrix} \hat{\boldsymbol{u}}\\ \hat{p} \end{bmatrix},\quad \text{with}\ \boldsymbol{\mathsf{M}} \doteq \begin{bmatrix} \boldsymbol{I} & 0\\ 0 & 0 \end{bmatrix}. \end{equation}

Systems (4.3) and (2.2) lead to the same solution for the velocity field, and (2.2) was selected for the analytical development only because it does not contain the singular mass matrix, which lightens the writing. To produce the linear and weakly nonlinear results, (4.3) is discretized on Matlab by means of the pseudospectral Chebyshev collocation method. Specifically, the variables are collocated at the $N$ Gauss–Lobatto nodes $s=\cos (\,j{\rm \pi} /(N-1))$, with $j=0,1,\ldots,N-1$, which includes the endpoints $s=-1$ and $s=1$. This grid is mapped on the physical domain $0\leq r \leq R_{max}$ using an algebraic mapping with domain truncation $r=L(1+s)/(s_{max}-s)$, where $L$ is a mapping parameter to cluster the points close to the origin, and is set equal to $3$. The parameter $s_{max}$ is defined as $s_{max} \doteq (2L+R_{max})/R_{max}$ (see Canuto et al. Reference Canuto, Quarteroni, Hussaini and Zang2007; Viola, Arratia & Gallaire Reference Viola, Arratia and Gallaire2016). In the present work, $R_{max}=50$ and $N=300$ proved to be sufficient for the convergence of all results. The lines corresponding to $r=0$ and $r=R_{max}$ in the discretized version of (4.3) are replaced by those enforcing the boundary conditions there; this also avoids the problem of the singularity in $r=0$.

The action of the propagator $\boldsymbol {\varPsi }(t,0)$ is computed in practice by marching (4.3) in time using a Crank–Nicholson scheme. Specifically, if the pressure and the velocity field are known and equal to $[\hat {\boldsymbol {u}}^{(n)},\hat {p}^{(n)}]^{\rm T}$ at a time $t_n$, their values at $t_{n+1}=t_n+\Delta t$ are obtained upon inverting the system

(4.4)\begin{equation} \left(\frac{\boldsymbol{\mathsf{M}}}{\Delta t} - \frac{1}{2}\boldsymbol{\mathsf{B}}_{m}(t_{n+1})\right)\begin{bmatrix} \hat{\boldsymbol{u}}^{(n+1)}\\ \hat{p}^{(n+1)} \end{bmatrix} = \left(\frac{\boldsymbol{\mathsf{M}}}{\Delta t} + \frac{1}{2}\boldsymbol{\mathsf{B}}_{m}(t_{n}) \right)\begin{bmatrix} \hat{\boldsymbol{u}}^{(n)}\\ \hat{p}^{(n)} \end{bmatrix}, \end{equation}

which is done by means of the command backslash on Matlab. The solutions to the other linear time-dependent systems at higher orders are also approximated upon discretizing in time with a Crank–Nicholson scheme. The adjoint system to (4.3) reads

(4.5)\begin{equation} -\boldsymbol{\mathsf{M}}\partial_t\begin{bmatrix} \hat{\boldsymbol{u}}^{{{\dagger}}}\\ \hat{p}^{{{\dagger}}} \end{bmatrix} = \boldsymbol{\mathsf{B}}_{m}(t)^{{{\dagger}}} \begin{bmatrix} \hat{\boldsymbol{u}}^{{{\dagger}}}\\ \hat{p}^{{{\dagger}}} \end{bmatrix}, \end{equation}

where we used $\boldsymbol{\mathsf{M}}^{{\dagger} } = \boldsymbol{\mathsf{M}}$, and where the expression of $\boldsymbol{\mathsf{B}}_{m}(t)^{{\dagger} }$ and the boundary conditions for the adjoint field are derived in Appendix E. The action of the adjoint propagator $\boldsymbol {\varPsi }(t,0)^{{\dagger} }$, necessary for the linear optimization, is determined by time marching (4.5) backward. That is, from known adjoint velocity and pressure fields at a time $t_{n+1}$, their evolution at $t_{n}=t_{n+1}-\Delta t$ are obtained by solving

(4.6)\begin{equation} \left(\frac{\boldsymbol{\mathsf{M}}}{\Delta t} - \frac{1}{2}\boldsymbol{\mathsf{B}}^{{{\dagger}}}_{m}(t_{n}) \right)\begin{bmatrix} \hat{\boldsymbol{u}}^{{{\dagger}},(n)}\\ \hat{p}^{{{\dagger}},(n)} \end{bmatrix} = \left(\frac{\boldsymbol{\mathsf{M}}}{\Delta t} + \frac{1}{2}\boldsymbol{\mathsf{B}}^{{{\dagger}}}_{m}(t_{n+1}) \right)\begin{bmatrix} \hat{\boldsymbol{u}}^{{{\dagger}},(n+1)}\\ \hat{p}^{{{\dagger}},(n+1)} \end{bmatrix}. \end{equation}

A time step of $\Delta t = 0.02$ was selected and was found sufficient to guarantee convergence of the results. The fully nonlinear direct numerical simulations (DNS) were performed using the spectral element solver Nek5000. A two-dimensional square grid $(x,y)\in [-R_{max},R_{max}]\times [-R_{max},R_{max}]$ was used, with a particularly high density of elements near the vortex core region where the flow gradients are intense. Convergence of the results by refining the mesh and extending the size of the computational domain has been checked.

4.3. Linear results

The results of the linear transient growth analysis for $m=2$ and varying ${\textit {Re}} \in [1.25,2.5,5,10]\times 10^{3}$ are presented below. The optimal gain $G_o(t_o)$ defined in (2.13), also called ‘envelope’, is shown as a function of $t_o$ in figure 1(a) and for different ${\textit {Re}}$. For all considered ${\textit {Re}}$ values, $G_o(t_o)$ reaches a clear maximum at relatively small $t_o$, before decreasing and plateauing by further increasing $t_o$. The values $G_o(t_{o,{max}})$ of these maxima, and of the corresponding temporal horizons $t_{o,{max}}$, seem to increase monotonically with the ${\textit {Re}}$.

Figure 1. Linear transient growth in the two-dimensional, time-dependent Lamb–Oseen vortex flow for varying ${\textit {Re}}$ $\in [1.25,2.5,5,10]\times 10^{3}$, larger ${\textit {Re}}$ corresponding to lighter colours. The perturbation has wavenumber $m=2$. (a) Optimal gain as defined in (2.13) as a function of the temporal horizon $t_o$. The star marker corresponds to its maximum value over $t_o$ and for a given ${\textit {Re}}$, that is, to $G_o(t_{o,{max}})=1/\epsilon _{o,{max}}$; (b) $\epsilon _{o,{max}}$ multiplied by ${\textit {Re}}^{1/3}$ and plotted as a function of ${\textit {Re}}$; (c) $t_{o,{max}}$ multiplied by ${\textit {Re}}^{-1/3}$ and shown as a function of ${\textit {Re}}$.

The type of dependence is made explicit in figure 1(b,c), where we propose in (b) to plot $\epsilon _{o,{max}}$ multiplied by ${\textit {Re}}^{1/3}$ as a function of ${\textit {Re}}$, demonstrating that $G_o(t_{o,{max}}) \propto {\textit {Re}}^{1/3}$. In (c), $t_{o,{max}}$ multiplied by ${\textit {Re}}^{-1/3}$ is also shown as a function of ${\textit {Re}}$: while the data align slightly less well on a constant line, the agreement remains correct and stating $t_{o,{max}} \propto {\textit {Re}}^{1/3}$ is a fair approximation. All these results are in excellent agreement with those presented in chapter 3 of Antkowiak (Reference Antkowiak2005). In particular, the power-law dependencies of both $G_o(t_{o,{max}})$ and $t_{o,{max}}$ have already been observed and interpreted physically, and discussed next.

In figure 2, by fixing ${\textit {Re}}=5000$, the optimal gain $G_o(t_o)$ over $t_o$ is reproduced from figure 1, as well as the gain $G(t)$ associated with the linear trajectory optimized for $t_o=t_{o,{max}}=35$ specifically. We check that $G(0)=1$, that $G(35)=G_o(35)$ and that $G$ is below $G_o$ everywhere else. At the times corresponding to the black dots, the axial vorticity structure, expressed as

(4.7)\begin{equation} \hat{{\omega}} = (1/r+\partial_r)\hat{u}_{\theta}-(\text{i} m/r)\hat{u}_r, \end{equation}

is shown in figure 3. For panel (a), we observe that the optimal initial condition consists of an arrangement of vorticity sheets, or ‘spirals’, orientated in counter-shear with respect to the base flow. As time increases, the sheets unfold under the effect of advection by the base flow, in analogy with the Orr mechanism. Antkowiak & Brancher (Reference Antkowiak and Brancher2004) and Antkowiak (Reference Antkowiak2005) also argue that a Kelvin wave (corresponding to a core mode) is transitorily excited in the heart of the vortex, by induction of the vorticity spirals. By ‘induction’ is meant that a radial velocity perturbation is induced in the heart of the vortex as the spirals unfold, and advects the base vorticity which is large here, thus acting as a source for the vorticity perturbation. The core mode is particularly visible in the heart of the vortex in panel (d), corresponding to $t=45$. However, it quickly disappears for later times, as the vorticity spirals roll up in the other direction, this time in line with the base advection (see panel e for $t=60$). Doing so, the vorticity perturbation decays with time due to a shear-diffusion averaging mechanism studied by Rhines & Young (Reference Rhines and Young1983) for a passive scalar and in Lundgren (Reference Lundgren1982) for vorticity. In particular, Lundgren (Reference Lundgren1982) concludes that, as long as the vorticity perturbations are rapidly radially varying, the shear-diffusion homogenization mechanism is qualitatively identical to the passive scalar case. This was confirmed numerically a few years later in Bernoff & Lingevitch (Reference Bernoff and Lingevitch1994), where is further demonstrated that the decay of the perturbation vorticity spirals acts on a ${\textit {Re}}^{1/3}$ time scale. From here, Antkowiak (Reference Antkowiak2005) argues that, in the limit of vanishing viscosity, the evolution equation for the perturbation vorticity $\hat {{\omega }}$ becomes time reversible; therefore, the unfolding (Orr) amplification mechanism can be seen as the ‘mirror’ of the shear diffusion, a damping one. Indeed, the curve for $G(t)$ is fairly symmetric around $t=t_o$ in figure 2, and would be even more for larger ${\textit {Re}}$ values. In this sense, the Orr mechanism has an ‘anti-mixing’ effect. This explains why the temporal horizon leading to the largest amplification is also in ${\textit {Re}}^{1/3}$, because it is the ‘natural’ amplification time scale of the vortex. It should be noted that this conception of the Orr and shear-diffusion mechanism as ‘mirroring’ each other was already present in Farrell (Reference Farrell1987) in the context of plane shear flows. As a side comment, Antkowiak (Reference Antkowiak2005) also derives the ${\textit {Re}}^{1/3}$ dependence of $t_{o,{max}}$ by simply balancing the unfolding and the diffusion (in the radial direction) time scales.

Figure 2. The continuous line is the reproduction of the envelope shown in figure 1(a) for ${\textit {Re}}=5000$. The dash-dotted line is the gain $G(t)$ associated with the linear trajectory optimized for $t_o=t_{o,{max}}=35$, defined in (2.16). By definition, both curves collapse at $t=35$. The black dots correspond to $t=0,15,30,\ldots,105$, for which the corresponding vorticity structures are shown in figure 3.

Figure 3. Temporal evolution of the vorticity structure $\hat {{\omega }}(r,t)$ of the optimal linear response, shown at the specific times corresponding to by the black dots in figure 2. Each panel shows only $[x,y]\in [-4,4]\times [-4,4]$. The plus sign denotes the origin, the dotted circle is the unit circle, and the dashed circle highlights the radius $r_q$ as defined in (4.9).

We insist that, apart from a tenuous transient excitation of a core mode by vortex induction, the amplification mechanism of the perturbation is essentially the Orr mechanism. Farrell & Ioannou (Reference Farrell and Ioannou1993) have shown that this mechanism was associated with the scaling law $G_o(t_o)\propto t_o$ (note that the gain defined in Farrell & Ioannou (Reference Farrell and Ioannou1993) is the square of the gain presently defined). From here, the scaling law $G_o(t_{o,{max}}) \propto {\textit {Re}}^{1/3}$ follows immediately.

The shear-diffusion (thus the Orr) mechanism is in essence non-viscous. The only role played by viscosity is to select the radial length scale of the initial vorticity structure (which would tend to be infinitely thin in the absence of viscosity, and the optimal gain and associated temporal horizon follow). Indeed, vorticity sheets that are too thin will be smoothed out by viscosity. If both the optimal gain and the decay time are $\propto {\textit {Re}}^{1/3}$, then the decay rate is clearly independent of ${\textit {Re}}$. In fact, the decay of the vorticity moments associated with the structures observable in figure 3 after $t=35$, can also be interpreted in terms of the Landau damping phenomenon, which is purely inviscid (Schecter et al. Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000), as developed in the introduction.

More precisely, the exponential decay rate of the vorticity moments can be obtained from a Landau pole of the vortex base profile, as first demonstrated in Briggs, Daugherty & Levy (Reference Briggs, Daugherty and Levy1970). For this, one first needs to define the $m$th multipole moment of the vorticity perturbation as

(4.8)\begin{equation} Q^{(m)}(t)=\int_{0}^{\infty} r^{m+1}\hat{{\omega}}(r,t)\,\mathrm{d} r = [r^{m+1}(\hat{u}_{\theta}-\text{i}\hat{u}_r)]_{r \rightarrow \infty}. \end{equation}

Then, to quote Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000), ‘A Landau pole is a complex frequency $\omega _q-\text {i} \gamma$ at which the Laplace transform of $Q^{(m)}(t)$ is singular [$\ldots$]. A Landau pole contributes a term to $Q^{(m)}(t)$ of the form $\exp ({-\gamma t - \text {i} \omega _q t})$’. A Landau pole, like an eigenvalue, is a property of some specific linear operator acting on a perturbation and, in this sense, depends on the base profile and on the wavenumber of the perturbation, but not on the specific radial shape of the latter. A Landau pole is generally not unique, and its calculation amounts to solving an eigenvalue problem along a radial contour that was deformed in the complex plane (see the appendix in Schecter et al. Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000); however, one can reasonably expect one of these potentially multiple Landau poles to dominates over the others.

In considering the response of an inviscid Gaussian vortex to a generic external impulse with $m=2$, Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000) indeed find that the quadrupole moment $Q^{(2)}(t)$ of the vorticity perturbation ‘decays exponentially, and the decay rate is given by a Landau pole’. In addition, is also mentioned that ‘the vorticity perturbation [$\ldots$] is poorly described by a single damped wave [$\ldots$] Rather, the vorticity perturbation is rapidly dominated by spiral filament’. This is also the structure observed from $t=60$ in figure 3, which closely resembles the figure 9 in Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000). In other terms, the response excites a very large number of structurally different eigenmodes, yet it does not invalidate the relevance of the dominant Landau pole.

By writing the quadrupole in terms of amplitude and phase, $Q^{(2)}(t)=|Q^{(2)}(t)| \exp ({{\rm i} \phi (t)})$, we display in figure 4(a) the temporal evolution of the phase $\phi (t)$ normalized by $2{\rm \pi}$ and corresponding to the response shown in figure 3. The amplitude $|Q^{(2)}(t)|$ is also shown in figure 4(b). The best fit of the form of a pure Landau damping $Q^{(2)}(t)=\exp ({-\text {i} \omega _q t-\gamma t})$ is also proposed. Clearly, the phase changes at a constant rate in figure 4(a) for $15 \leq t \leq 65$, and $|Q^{(2)}(t)|$ decays exponentially in figure 4(b) for $40 \leq t \leq 65$. Therefore, we confirm that the quadrupole moment $Q^{(2)}(t)$, associated with the linear response in figure 3, is well approximated by a Landau pole in its decaying part, despite the fact that we consider a viscous flow over a time-dependent base vortex, which is not the case in Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000). The Reynolds number ${\textit {Re}} = 5000$, however, seems sufficiently large for the conclusions drawn in Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000) to also hold here. In the following, reference will be made to the Landau pole in order to qualitatively interpret the nonlinear effects, even if the fact that ${\textit {Re}}$ is finite possibly makes it less rigorous. The precise value of $\gamma$ in figure 4(b) is unimportant, for the comparison with any analytical prediction of the Landau pole is outside of the scope of this paper. Note that the exponential decay is followed by an algebraic one for large times, which is also observed in Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000). The fitted value $\omega _q=0.42$ in figure 4(a) is, however, important, for it gives the location $r_q$ of the radius where the angular velocity $\omega _q/m$ associated with the Landau pole is equal to that of the base flow

(4.9)\begin{equation} m \varOmega(r_q(t),t) = \omega_q. \end{equation}

The radius solving (4.9) at each time is highlighted by a dashed circle in figure 3. It takes the value $r_q \approx 2.16$, very weakly dependent on time.

Figure 4. Temporal evolution of the vorticity moment $Q^{(2)}(t)$ corresponding to the linear response shown in figures 2 and 3; the black dots again denote the specific times where the vorticity field is shown in figure 3. (a) Phase $\phi (t)$ divided by $2{\rm \pi}$. (b) Amplitude $|Q^{(2)}(t)|$. The black dashed lines correspond to a pure Landau damping $Q^{(2)}(t) = \exp ({-\text {i} \omega _q t-\gamma t})$, where we use fitted values for $\omega _q$ and $\gamma$.

It can be shown that the decay rate of the Landau pole is linked to the radial derivative of the base vorticity at $r_q$ specifically. In fact, as demonstrated in Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000), Turner & Gilbert (Reference Turner and Gilbert2007) and Turner et al. (Reference Turner, Gilbert and Bassom2008), the exponential damping mechanism can be removed from any vortex, in particular a Gaussian one, by cancelling such derivative at $r_q$.

The poor modal description of a perturbation evolution over a Gaussian vortex, mentioned in Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000) and confirmed in figure 3, is the reason for which non-modal analytical tools are deployed in the present study, in particular for studying nonlinear effects. We address these latter now.

4.4. Fully nonlinear results

We introduce the fully nonlinear results obtained from DNS in the present section. Let us define $\boldsymbol {u}_p$ as being the perturbation between the fully nonlinear solution $\boldsymbol {U}$ obtained from a DNS, and the reference Lamb–Oseen solution $\boldsymbol {U}_{{b}}$ in (4.1a,b); in other terms, $\boldsymbol {u}_p \doteq \boldsymbol {U}-\boldsymbol {U}_{{b}}$. We recall that this perturbation is initialized along the wavenumber $m=2$ with the linear optimal initial condition with an amplitude $U_0$, such that

(4.10)\begin{equation} \boldsymbol{u}_p|_{t=0} = U_0 \hat{\boldsymbol{u}}_o \,{\rm e}^{\text{i} m \theta} + \text{c.c.}. \end{equation}

Owing to nonlinear effects, $\boldsymbol {u}_p$ generally does not oscillate purely along $m$ for strictly positive times. For this reason we need to further extract $\hat {\boldsymbol {u}}^{(1)}_p$, the component of $\boldsymbol {u}_p$ oscillating at the fundamental wavenumber $m$, hence the superscript ‘$(1)$’. For this, a Fourier series is naturally used as

(4.11)\begin{align} \left.\begin{gathered} \hat{\boldsymbol{u}}^{(1)}_p(r) \doteq \frac{1}{2{\rm \pi}}\int_{0}^{2{\rm \pi}} \boldsymbol{u}_p(r,\theta) \cos(m\theta)\,\mathrm{d} \theta - \frac{\text{i}}{2{\rm \pi}}\int_{0}^{2{\rm \pi}} \boldsymbol{u}_p(r,\theta) \sin(m\theta)\,\mathrm{d} \theta,\quad \text{such that} \\ \boldsymbol{u}^{(1)}_p(r,\theta) \doteq \hat{\boldsymbol{u}}^{(1)}_p(r)\,{\rm e}^{\text{i} m\theta} + \text{c.c.}. \end{gathered}\right\} \end{align}

From here, the fully nonlinear gain (associated with the fundamental pair) is immediately defined as

(4.12)\begin{equation} G_{{DNS}} \doteq \frac{\|\hat{\boldsymbol{u}}^{(1)}_p\|}{U_0}. \end{equation}

By fixing ${\textit {Re}}=5000$ and $t_o=35$, we report in figure 5(a) the fully nonlinear evolution of the perturbation amplitude $\|\hat {\boldsymbol {u}}^{(1)}_p\|$, for different amplitudes of the initial condition $U_0$. The associated gain, that is, the same data divided by the corresponding $U_0$, is shown in figure 5(b).

Figure 5. (a) Amplitude of the fully nonlinear fundamental perturbation $\hat {\boldsymbol {u}}_p^{(1)}$ as defined in (4.11), initialized with the linear optimal condition for $({\textit {Re}},m,t_o)=(5000,2,35)$. Larger $U_0$ values correspond to lighter colours. (b) The same data are divided by their corresponding $U_0$, yielding the nonlinear gains as defined in (4.12). For the largest considered $U_0 = 2.82 \times 10^{-2}$, black dots are located at $t=0,30,60,\ldots,210$, for which the corresponding vorticity structures are shown in figure 6.

Figure 6. Temporal evolution of the vorticity structure of the fully nonlinear perturbation for $U_0 = 2.82\times 10^{-2}$, shown at the specific times highlighted by the black dots in figure 5. Each panel shows only $[x,y]\in [-4,4]\times [-4,4]$, the plus sign denotes the origin, and the dashed line is the unit circle.

Figure 5(a) illustrates well a bypass transition, also called ‘bootstrapping effect’ in Trefethen et al. (Reference Trefethen, Trefethen, Reddy and Driscoll1993), where linear transient growth and nonlinear mechanisms interact to bring about a transition to a new state, distinct from the reference one. Indeed, as we increase the amplitude of the initial condition (linearly in log scale), the transient growth of the perturbation, initiated by the linear Orr mechanism as we have seen, becomes sufficient to trigger nonlinear terms; for the two largest considered $U_0$, this prevents the flow to re-axisymmetrize to the Lamb–Oseen vortex. Interestingly, for the third and fourth largest considered $U_0$, the flow seems to temporarily approach a new state but nevertheless relaxes towards the reference one.

For the largest considered $U_0=2.82 \times 10^{-2}$ and the set of times highlighted by the black dots, we report in figure 6 the structure of the vorticity perturbation. From $t=60$, the flow has reached a new nonlinear quasi-equilibrium state, that diffuses very slowly due to viscous effects, and rotates counterclockwise with a period of approximately $45$ units of times. The corresponding total vorticity field, which is obtained by adding to figure 6 the Gaussian reference vorticity (4.2), takes a tripolar shape. This is shown in figure 7 for $t\geq 120$, where the heart of the vortex takes an elliptical shape that is surrounded by two satellite vortices of low and negative vorticities (corresponding to the two blues spots in figure 6 from $t=60$). As developed in the introduction, this tripolar state was already largely reported in a variety of contexts. In particular, its spatial structure compares well with the one reported in figure 2 of Antkowiak & Brancher (Reference Antkowiak and Brancher2007) for ${\textit {Re}}=1000$, or with the experimental visualization in Kloosterziel & van Heijst (Reference Kloosterziel and van Heijst1991).

Figure 7. Same as in figure 6 but the reference vorticity field (4.2) has been added, so as to visualize the total, fully nonlinear, vorticity field; some isolines are also shown to better expose the elliptical deformation of the vortex core.

In terms of the gain in figure 5(b), we notice that all the curves corresponding to different $U_0$ collapse for small times, confirming the linearity of the initial growth mechanism. However, they significantly depart from each other at larger times. If nonlinearities seem to saturate the gain for times around $t_o=35$, they clearly increase the latter up to three orders of magnitude for large times; as said, that is because the flow bifurcated to another state, thereby the perturbation that is measured around the reference vortex remains large.

The weakly nonlinear model is now employed to assess the gain evolution with $U_0$ shown in figure 5(b).

5.1. Transient change from nonlinear saturation to nonlinear amplification

The real part of the weakly nonlinear coefficient, defined in (3.29), is shown in figure 8(a) for ${\textit {Re}}=5000$ and $t_o=t_{o,{max}}=35$. It is further split into the sum of a mean flow distortion contribution, $\mu _r^{(0)}$, arising from $\boldsymbol {u}_2^{(0)}$ only, and second harmonic contribution, $\mu _r^{(2)}$, arising from $\hat {\boldsymbol {u}}_2^{(2)}$ only. In this manner, $\mu _r = \mu _r^{(0)} + \mu _r^{(2)}$.

Figure 8. (a) Real part of the weakly nonlinear coefficient $\mu$ defined in (3.29) (black continuous line). It is further split as the sum of a contribution from of a mean flow distortion (red dash-dotted line), plus a contribution from the second harmonic (blue dotted line). In the grey zone, the weakly nonlinear gain $G_{{w}}$ defined in (3.30) increases monotonically with $U_0$, since $\mu _r > 0$. In the white zone, it decreases monotonically. (b) The coefficient $\mu _r$ is compared with its reconstruction from DNS data (magenta dashed line) according to (5.2).

Remarkably, after having significantly decreased until reaching a minimum at $t=59$, the coefficient $\mu _r$ increases again and changes sign at $t=87$. It then reaches its maximum at $t=107$, and decreases again until plateauing around zero for $t \geq 150$. A change of sign in $\mu _r$ brings diversity to the behaviours reported in Ducimetière et al. (Reference Ducimetière, Boujo and Gallaire2022), concerning transient growths in the streamwise-invariant Poiseuille flow, and in the flow past a backward-facing step. There, only negative coefficients, corresponding to saturating nonlinearities, were observed. In the present work, however, nonlinearities appear to reinforce the gain for some times. When it is negative, the behaviour of $\mu _r$ seems largely dominated by the contribution from the mean flow distortion $\mu _r^{(0)}$. When the former is positive, however, $\mu _r^{(0)}$ is reduced and the second harmonics appear to contribute as much to the sum.

Upon taking the derivative of the square the weakly nonlinear gain in (3.34) with respect to $U_0^2$, we obtain

(5.1)\begin{equation} \frac{\partial (G_{{w}}^2)}{\partial (U_0^2)} = \frac{2\mu_r}{\epsilon_o^2-2U_0^2\mu_r}G_{{w}}^2. \end{equation}

Therefore, the coefficient $\mu _r$ relates directly to the rate of change of $G_{{w}}^2$ with respect to $U_0^2$ in the limit where $U_0$ tends towards zero as

(5.2)\begin{equation} \mu_r = \frac{\epsilon_o^2}{2} \left(\lim_{U_0 \rightarrow 0} \frac{1}{G_{{w}}^2 }\frac{\partial (G_{{w}}^2)}{\partial (U_0^2)}\right). \end{equation}

By using (5.2), the coefficient $\mu _r$ can be reconstructed directly from DNS data. For this, the derivative in (5.2) is estimated by using a first-order finite difference approximation between DNS data for the gain squared, corresponding to the two lowest considered $U_0= 2.26\times 10^{-4}$ and $U_0= 2.82 \times 10^{-4}$. The result is shown as the magenta dotted line in figure 8(b), and compared with the actual $\mu _r$. The good agreement between both curves for $t\leq 130$ a posteriori validates for these times our weakly nonlinear expansion, a least in the limit of small $U_0$. However, both curves depart significantly after $t=130$, presumably due to the violation of the condition (3.11) ensuring that the asymptotic expansion is well posed. Indeed, the response to the initial condition optimized for $t=t_o=35$ was shown to rapidly decay in amplitude for larger times in figure 2; therefore the norm of the perturbation operator, $\|\boldsymbol{\mathsf{P}}(t)\|=1/\|\hat {\boldsymbol {l}}(t)\|$, increases accordingly.

5.2. Comparison of fully and weakly nonlinear gains

In figure 9 the weakly nonlinear gain, associated with the coefficient in figure 8, is compared with the fully nonlinear one. For the moment, only short times where the gains are large are shown. The later evolution, being associated with orders of magnitude smaller gains, will be considered in figure 10 in log scale. On the time interval chosen for figure 9, the coefficient $\mu _r$ is negative, thereby the model predicts the gain to decrease monotonically with $U_0$. For values of $U_0$ small enough so that the linear gain (black curves in figure 9) is only slightly modulated, the agreement with DNS data is excellent. By increasing $U_0$, the agreement degrades only slowly and remains very good for these relatively small times, for instance over $0 \leq t \leq 40$ in the panel corresponding to ${\textit {Re}}=10^4$. This corresponds to times when the nonlinear structure remains symptomatic of the linear one, in the sense that the flow has not yet reached the tripolar state (temporarily or not, as shown in figure 5). Indeed, the amplitude equation (3.33) is independent of space, which condemns the response to be structurally close to the linear one. For later times and large $U_0$, as soon as the fully nonlinear gain begins to oscillate and to bear a non-monotonic behaviour, the agreement with our weakly nonlinear model rapidly degrades. The reason is precisely that, the response whose nonlinear evolution is approached by our model, is not selected by the flow anymore, which has reached another state, temporarily or not, and that is structurally completely different; about this new state, the amplitude equation has no information.

Figure 9. Weakly and fully nonlinear gains as defined in (3.30) and (4.12), respectively. Larger $U_0 \in [0.03,0.45,1.12,2.82]\times 10^{-2}$ correspond to lighter colours (direction of increasing $U_0$ is also indicated by the arrow). The continuous line stands for DNS data whereas the dash-dotted is for the weakly nonlinear model. Temporal horizons are $t_o = t_{o,{max}} = [25,30,35,40]$ (times at which a black star is shown) for ${\textit {Re}}=[1.25,2.5,5,10]\times 10^3$, respectively.

Figure 10. Weakly and fully nonlinear gains for ${\textit {Re}}=5000$ and increasing amplitude of the initial condition $U_0\in [4.5,7.1,11.2,17.8]\times 10^{-3}$. Each panel corresponds to a different $U_0$. The grey box denotes the time interval (5.3) over which the weakly nonlinear gain is undefined. The latter is singular at the boundaries of this interval.

Nevertheless, while the proposed amplitude equation fails to predict the gain and the structure of the flow when at the tripolar state, it does predict a bifurcation threshold. This is illustrated in figure 10, where the same data as in figure 9 for ${\textit {Re}}=5000$ are shown, although the time interval is extended until $t=150$ so as to include the time interval where the real part of the weakly nonlinear coefficient becomes positive, indicating ‘anti-saturation’. The weakly nonlinear gain is undefined on the time interval where

(5.3)\begin{equation} \mu_r(t) \geq \frac{1}{2}\left(\frac{\epsilon_o}{U_0} \right)^2 > 0, \end{equation}

which widens by increasing $U_0$. The gain is singular (tends to infinity) at the times corresponding to the boundaries of this interval. The latter is highlighted by the grey zones in figure 10.

5.3. Bifurcation thresholds

By looking at the panels corresponding to $U_0=1.12\times 10^{-2}$ and $U_0=1.78\times 10^{-2}$, we observe that, over the time interval where the equation has no solution, i.e. the grey zone, the fully nonlinear simulation seems to have reached the tripolar state (although for $U_0=1.12\times 10^{-2}$ it is only temporary as it eventually relaxes towards the reference state). In this sense, a loss of solution in the amplitude equation could indicate that the DNS has left the state around which the weakly nonlinear expansion was constructed. The minimum $U_0$ for which the weakly nonlinear gain becomes singular may then be considered as an approximation of the actual bifurcation threshold. Such minimal $U_0$, named $U_0^{s}$ is what follows, reads

(5.4)\begin{equation} U_0^{s} = \frac{\epsilon_o}{\sqrt{2 \mu_r(t_s)}},\quad \text{where}\ \mu_r(t_s)=\max_{t}\mu_r(t), \end{equation}

and is defined if and only if $\mu _r(t_s)>0$.

A bifurcation threshold in the DNS solutions must also be defined and compared directly with (5.4). To the knowledge of the authors, there is no clear and universal subcritical bifurcation criterion, and the choice made is always arguably arbitrary. Nevertheless, we will opt for the criterion proposed in Antkowiak (Reference Antkowiak2005), based on the observation that the tripolar state is characterized by a deformed, non-axisymmetric, heart. Therefore a characteristic aspect ratio $\lambda$ of the heart is established as

(5.5)\begin{equation} \lambda^2 \doteq \frac{J+R}{J-R},\quad \text{with}\ J=(J_{20}+J_{02}),\quad R=\sqrt{(J_{20}-J_{02})^2+4J_{11}^2}, \end{equation}

and $J_{mn}\doteq \int _{\varOmega }x^m y^n \hat {{\omega }}(x,y)\,\mathrm {d}\kern 0.06em x\,\mathrm {d} y$ a vorticity moment. A characteristic eccentricity is further defined as $e\doteq \sqrt {1-1/\lambda ^2}$. From here, Antkowiak (Reference Antkowiak2005) computes the typical half-life time (that is where the criterion is rather arbitrary) of the heart deformation as

(5.6)\begin{equation} \tau,\quad \text{such that}\ \int_{0}^{\tau}e(t)\,\mathrm{d} t = \frac{1}{2}\int_{0}^{t_f} e(t)\,\mathrm{d} t, \end{equation}

where $t_f=500$ is the final time of the fully nonlinear simulations.

We report the half-life time as a function of $U_0$ and for different ${\textit {Re}}$ values in figure 11(a). For each ${\textit {Re}}$, we also highlight an inflection point in $\tau$ at a certain $U_0$, and we declare the latter as being the subcritical bifurcation threshold. It is reported directly as a function of ${\textit {Re}}$ in figure 11(b), and compared with the weakly nonlinear prediction $U_0^s$, in both linear and log scales. Both approaches clearly highlight a decreasing power-law dependence of the subcritical bifurcation threshold with the ${\textit {Re}}$. For the fully nonlinear data, the fitted exponent $-0.88$ found in the present study, agrees relatively well with the one reported in Antkowiak & Brancher (Reference Antkowiak and Brancher2007), which is $-0.8$. However, since little information is given regarding how the threshold amplitude was computed in Antkowiak & Brancher (Reference Antkowiak and Brancher2007), the agreement with our results is perhaps fortunate. In any case, the negative power-law dependence with ${\textit {Re}}$ implies that the threshold amplitude above which the flow goes into the basin of attraction of the tripolar state, and in the direction given by the linear optimal, vanishes for increasing ${\textit {Re}}$. This may explain the formation and persistence of elliptical vortex eyewalls in some tropical cyclones (Kuo et al. Reference Kuo, Williams and Chen1999; Reasor et al. Reference Reasor, Montgomery, Marks and Gamache2000). We insist that we only consider perturbations in the direction of the linear optimal, whereas nonlinear optimal perturbations can also be found by relying on the techniques presented in Pringle & Kerswell (Reference Pringle and Kerswell2010). The latter would be associated with a threshold amplitude $U_0$ for a subcritical bifurcation even smaller than the one shown in figure 11. The fitted exponent for the weakly nonlinear model is found as being $-0.66$, thereby the threshold amplitudes between both models differ significantly by decreasing ${\textit {Re}}$ in figure 11(b).

Figure 11. (a) Typical half-life time of the heart deformation defined in (5.6), as a function of the amplitude of the initial condition $U_0$. Larger ${\textit {Re}}=[1.25,2.5,5,10]\times 10^3$ correspond to lighter colours. A star highlights an inflection point, for which the corresponding $U_0$ is declared as being the threshold amplitude for the subcritical bifurcation. Such thresholds $U_0$ are reported in (b) as a function of ${\textit {Re}}$ (also with a star symbol). The prediction from the weakly nonlinear model, $U_0^s$ defined in (5.4) is also shown. The thin continuous line is a power law fitted on the DNS data for the first three considered ${\textit {Re}}$, with ${\propto }Re^{-0.88}$, whereas the thin dashed line is fitted on the weakly nonlinear data with ${\propto }Re^{-0.66}$. The inset shows the same in log–log scale.

This discrepancy between exponents may possibly be explained as follows. For the DNS, our bifurcation criterion aims at computing the threshold above which the flow has definitely bifurcated, whereas the loss of solution in the amplitude equation refers to a specific time interval. The importance of this conceptual difference is well illustrated in figure 10(c) for $U_0=1.12\times 10^{-2}$. Here, the amplitude equation predicts that no solution exists over some time interval; thus, according to the criterion (5.4), the flow has bifurcated. However, if the DNS data seem indeed to have reached the tripolar state over this time interval, it relaxes to the reference state at longer times; thus, the criterion based on $\tau$ concludes that the flow has not yet bifurcated.

For this reason, we propose another criterion, rather artificial, for which both ‘bifurcations’ refer to the same specific time $t_s$, the first time for which the amplitude equation predicts a loss of solution (see (5.4)). The threshold $U_0$ in the DNS is decreed as being the one corresponding to an inflection point in $G_{{DNS}}(t=t_s)$, as shown in figure 12(a), despite a quite coarse resolution. The agreement with $U_0^s$ improves significantly in figure 12(b) (as compared with figure 11b). Note that, for both models, the power-law dependence of the threshold amplitude cannot be only explained by the power-law dependence of the linear optimal transient gain, for the latter was shown to be in $-1/3$.

Figure 12. Same as in figure 11, although the inflection point is sought in $G_{{DNS}}(t=t_s)$, where $t_s$, defined in (5.4), is the first time for which the amplitude equation predicts a loss of solution. In (a), the inset shows the same data in linear scale. In (b), the thin continuous line is a power law fitted on the DNS data with ${\propto }Re^{-0.69}$, whereas the thin dashed line ${\propto }Re^{-0.66}$ is similar to figure 11.

5.4. Physical interpretation of the nonlinear ‘anti-saturation’: mean flow distortion and inversion of the vorticity gradient

We now propose to interpret physically the behaviour of the coefficient associated with the (azimuthal) mean flow distortion, $\mu _r^{(0)}$ shown in figure 8, and attempt to discuss the reason why it changes sign and leads to the subcritical behaviour discussed above. In the following discussion, we will focus on the case ${\textit {Re}}=5000$. First, the energy of the mean flow distortion $\boldsymbol {u}_2^{(0)}$ is shown as the red curve in figure 13(a).

Figure 13. (a) For ${\textit {Re}}=5000$, the evolution of the energy $\|\ast \|^2$ of the linear response $\hat {\boldsymbol {l}}$ (black dashed line), of the second harmonic $\hat {\boldsymbol {u}}_2^{(2)}$ (blue dashed-dotted line) and, mainly, of the mean flow distortion $\boldsymbol {u}_2^{(0)}$ (red continuous line). The red dots correspond to the specific times $t=30,40,\ldots,100$, for which the vorticity structure of $\boldsymbol {u}_2^{(0)}$ is reported in figure 14. Two horizontal dotted lines are drawn at $t=35$ and $t=85$. (b) Slope of the vorticity at the radius $r_q$, i.e. $\partial _r \omega _2^{(0)}|_{r=r_q}$, as a function of time.

Figure 14. Temporal evolution of $\omega _2^{(0)}$, the vorticity structure of $\boldsymbol {u}_{2}^{(0)}$, the mean flow distortion induced by the Reynolds stress of the linear response shown at the specific times corresponding to the red dots in figure 13(a). Each panel shows only $[x,y]\in [-4,4]\times [-4,4]$. The plus sign denotes the origin, the dotted circle is the unit circle, and the dashed circle highlights the radius $r_q$ solving (4.9) with $\omega _q=0.42$.

For comparison, the energies of the linear response and of the second harmonic are also shown. As written in (3.17), the mean flow distortion is forced by

(5.7)\begin{equation} \boldsymbol{f}_2^{(0)}=-\boldsymbol{\mathsf{C}}_{m} \big[\hat{\boldsymbol{l}}^* , \hat{\boldsymbol{l}}\big] + \text{c.c.}, \end{equation}

which is the Reynolds stress of the linear response. Under its action, the energy of $\boldsymbol {u}_2^{(0)}$ increases until reaching a maximum around $t=t_o=35$, before decaying until $t=65$. From here, the energy rebounds very slightly, then decays extremely slowly for $t\geq 85$. It is striking to notice in figure 13(a) that $\boldsymbol {u}_2^{(0)}$ can persist for extremely long times, even when the linear response, whose nonlinear interactions force $\boldsymbol {u}_2^{(0)}$, has vanished. This suggests that $\boldsymbol {u}_2^{(0)}$ can persist even in the absence of sustained forcing. This is also illustrated in figure 14, where we show $\omega _2^{(0)}$, the vorticity structure associated with $\boldsymbol {u}_2^{(0)}$, for the different times highlighted by the red dots in figure 13(a). If we observe a transient regime until the panel for $t=70$, the vorticity field does not seem to evolve over time afterward, except by slow diffusion.

It will simplify the rest of the analysis to notice that $\boldsymbol {u}_2^{(0)}$ only possesses an azimuthal component, i.e. $\boldsymbol {u}_2^{(0)}= u_{2,\theta }^{(0)} \boldsymbol {e}_{\theta }$. That is because $\boldsymbol {u}_2^{(0)}$ is associated with the wavenumber $m=0$, for which the equation for the radial perturbation is $\partial _r(r u_{2,r}^{(0)}) = 0$ from the continuity. This leads to $u_{2,r}^{(0)}=0$ and a forced diffusion equation for $u_{2,\theta }^{(0)}$, reading

(5.8)\begin{equation} \partial_t u_{2,\theta}^{(0)} = {\textit {Re}}^{-1}(\varDelta_{0}-1/r^2)u_{2,\theta}^{(0)} + f_{2,\theta}^{(0)}. \end{equation}

The panels in the first row in figure 15 show the evolution of $f_{2,\theta }^{(0)}$ (left) and $u_{2,\theta }^{(0)}$ (right) over $0 \leq t \leq 35$. The Reynolds stress forcing $f_{2,\theta }^{(0)}$ grows in amplitude (since $\hat {\boldsymbol {l}}$ does) while roughly conserving its shape, which $u_{2,\theta }^{(0)}$ seems to imitate closely, although logically with some delay. This direct structural link between the forcing and the response is certainly due to the absence of an advection term in (5.8).

Figure 15. Temporal evolution of the forcing $f_{2,\theta }^{(0)}$ (panels on the left) and of the velocity response $u_{2,\theta }^{(0)}$ (panels on the right). Panels on the top and the bottom correspond to two different and successive time intervals. The top line is for $t=0,5,\ldots,35$ (on the left of the first vertical line in figure 13a) and the second for $t=35,40,\ldots,85$ (between the two vertical lines in figure 13a). Larger times correspond to lighter colours.

The panels in the second row in figure 15 show the evolution over the times $35 \leq t \leq 85$. Notice that the forcing $f_{2,\theta }^{(0)}$ has shapes at $t=30$ and $t=40$ that are similar but have opposite signs. As developed in § 4.3, this is because the spiral structures in $\hat {\boldsymbol {l}}$, that unroll in one direction until $t=35$, roll up symmetrically in the other direction afterward. Under the action of this forcing that is now adverse, the velocity $u_{2,\theta }^{(0)}$ tends to invert in figure 15(c), which implies its momentary flattening, therefore $u_{2,\theta }^{(0)}$ rapidly loses energy. This energy would have grown again (which it does shortly after $t=65$ in figure 13a) if the forcing could maintain its intensity, but the latter decays rapidly with the one of $\hat {\boldsymbol {l}}$.

After $t = 85$ the velocity $u_{2,\theta }^{(0)}$ evolves freely, for $f_{2,\theta }^{(0)}$ is negligible. The persistence of $u_{2,\theta }^{(0)}$ is then easily understood once it is realized that the diffusion operator in (5.8) possesses a continuum of orthogonal Bessel eigenfunctions (the domain is infinite). They are associated with eigenvalues spanning the strictly negative part of the imaginary axis. Accordingly, after discretization, we found an extremely dense packing of eigenvalues all along the negative part of the imaginary axis, increasingly dense with the number of discretization points. The velocity structure left by the forcing that has vanished at $t=85$, projects over a significant number of these Bessel eigenfunctions, including some with very low damping rates of the order of $10^{-4}$. It is therefore not surprising to see this structure subject to very slow diffusion in figures 13(a) and 14.

After having proposed some elements to understand the evolution of $u_{2,\theta }^{(0)}$, let us study now how it feeds back on the response. We show in figure 13(b) the slope of $\omega _2^{(0)}$ at the specific radius $r_q$ where the angular velocity $\omega _q/m$ associated with the Landau pole is equal to that of the base flow. We recall this specific radius to be found by solving (4.9) with the fit value $\omega _q = 0.42$, resulting in $r_q\approx 2.16$, very weakly time dependent.

The slope is negative until $t = 60$, where it changes sign. As a consequence, from $t = 60$ onward, the addition of $\omega _2^{(0)}$ to the reference vorticity has a positive contribution to the total vorticity slope at $r_q$. In parallel, we recall that $\mu _r^{(0)}$ can be directly interpreted as the sensibility at a time $t$ of the transient gain to the addition of $u_{2,\theta }^{(0)}$ to the reference flow. For this, $\mu _r^{(0)}$ computes the integrated effect of this addition between $0$ and $t$ (see formula (C10)). Therefore, the decreasing tendency of $\mu _r^{(0)}$ until $t=60$ in figure 8, results from the negativity of the slope of $\omega _2^{(0)}$ at $r_q$, which enhances the Landau damping rate of the response, whose integrated effect is to reduce the gain. Such an effect of the mean vorticity slope at $r_q$ on the Landau damping was reported in Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000), Turner & Gilbert (Reference Turner and Gilbert2007) and Turner et al. (Reference Turner, Gilbert and Bassom2008). On the contrary, from $t=60$ onward the coefficient $\mu _r^{(0)}$ increases, for the presence of $\omega _2^{(0)}$ now reduces the Landau damping rate; the integrated effect of such mitigation of the Landau damping rate over time leads to a coefficient that becomes positive, traducing a weakly nonlinear gain larger than the linear one. In other words, we interpret the subcritical behaviour of our amplitude equation as being partially due to the fact that the mean flow distortion reminiscence, shown in figure 14 for large time, tends to erase the vorticity slope at the specific radius $r_q$, therefore letting the perturbation persist. This conclusion is comparable to the one drawn in Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young2001). In addition, the present amplitude equation leads to the conclusion that, when $\mu _r$ is at its maximum, the effect of the second harmonic is just as important as the one of the mean flow, although it is harder to interpret.

We insist that the conclusions drawn below relate only to predictions of the amplitude equation, thereby possibly explaining DNS behaviour only when both approaches agree. However, precisely because they do not at large $U_0$ and after increasingly short times in figures 9 and 10, the departure between weakly and fully nonlinear responses there remains to be interpreted. This can still be done by using the amplitude equation, although in an indirect manner.

In figure 16, we compare the linear stability to $m=2$ perturbations of the (azimuthal) mean flow extracted from the DNS, and the one predicted by the weakly nonlinear approach. This is done for $U_0=2.82 \times 10^{-2}$, the largest considered initial condition amplitude in figure 9. The weakly nonlinear mean flow is simply obtained by evaluating the weakly nonlinear expansion as $U_{{b},\theta } + a^2 u_{2,\theta }^{(0)}$, where $a$ solves the amplitude equation. The linear stability is assessed by replacing the reference flow by the mean one in (2.2), then assuming a temporal dependence as $\hat {\boldsymbol {u}}(r,t) \doteq \breve {\boldsymbol {u}}(r)\ \textrm {e}^{\sigma t}$, and solving the subsequent eigenvalue problem. The azimuthal mean is linearly unstable if and only if there exists an eigenvalue $\sigma$ such that ${Re}(\sigma ) = \sigma _r > 0$, and the corresponding eigenmode grows exponentially to the extent that the growth rate is much smaller than the rate of change of the base flow. The mean vorticity profiles are also shown in figure 17.

Figure 16. (a) Maximum growth rate of the eigenvalues of the Navier–Stokes operator at ${\textit {Re}}=5000$ and linearized around the (azimuthal) mean flow, in order to describe $m=2$ perturbations. Mean flows corresponding to $U_0=2.82\times 10^{-2}$ are obtained either from DNS data (continuous line linking red dots) or by evaluating the weakly nonlinear expansion (dash-dotted line linking blue diamonds). Positive growth rate values imply linear instability. (b) Shows the eigenmode corresponding to the most unstable eigenvalue of the operator linearized around the DNS mean flow, at $t=20$. The radius of the dotted circle is equal to one, the radius of the dashed circle is equal to $r_q$, and the radii of the continuous-line circle highlight the extrema of the mean vorticity profile. The colour scale is arbitrary.

Figure 17. Mean (axisymmetric) vorticity profiles corresponding to $U_0=2.82\times 10^{-2}$, extracted from DNS data (red continuous) and reconstituted from the weakly nonlinear expansion as $W_z + a^2\omega _2^{(0)}$, where $a$ solves the amplitude equation (blue dashed-dotted line). The largest dot (respectively diamond) is located at the radius for which $-m\varOmega =\sigma _i$ where $\sigma _i$ is the imaginary part of the most unstable eigenvalue of the DNS (respectively weakly nonlinear) profile. This also corresponds to the critical radius of the most unstable mode. The respectively smaller dot or diamond (if exists) stands for the second most unstable mode.

For $10 \leq t \leq 20$ in figure 16, both DNS and weakly nonlinear mean flows appear unstable, with an excellent agreement between the growth rates, and between the corresponding vorticity profiles in figure 17. For $t=15$ and $t=20$ in figure 17, these profiles consist of a central region of large vorticity, surrounded by a ring of locally enhanced but relatively smaller vorticity. A local minimum is to be noticed at $r=1.6$ between these two regions, as well as a second one a little further at $r=2.5$. Note the qualitative resemblance with the profiles considered in Kossin et al. (Reference Kossin, Schubert and Montgomery2000) (see their figures 1 and 12). The structure of the most unstable eigenmode supported by the DNS mean vorticity profile at $t=20$ is shown at the right of figure 16. It bears several characteristics of a shear instability. Specifically, its phase velocity is very close to the angular velocity of the mean flow at the radius of the vorticity extremum (see the largest dot in figure 17). This radius, highlighted by the first black continuous circle in figure 16, is also the zero amplitude isoline of the eigenmode structure. Tightly around, the eigenmode reaches its largest amplitude under the form of four external lobes, rotated of ${\approx }{\rm \pi} /2$ in the anticlockwise direction with respect to four internal ones. This structure is similar to the one shown in figure 2 of Carton & Legras (Reference Carton and Legras1994), who interpret it under the Rossby-wave interaction paradigm. In words, the mean vorticity profile can be thought of as supporting Rossby (vorticity) waves at the ‘edges’ (sharp slope) at one side and the other of the local minima. These two vorticity waves are precisely the internal and external series of lobes of the eigenmode structures. The velocity perturbation induced by a vorticity wave has a phase lag of ${\rm \pi} /2$ with respect to the latter. Therefore, since the lobe series are also rotated of ${\approx }{\rm \pi} /2$ with respect to each other, the velocity perturbation induced by one is in phase with the vorticity perturbation of the other, thus reinforcing its growth. In addition, Carton & Legras (Reference Carton and Legras1994) mentions that both waves interact in such a way as to travel with the same phase speed, thus the reinforcement of the one by the other is preserved in time, eventually leading to instability.

From $20 \leq t$ in figure 16, growth rates determined with both approaches depart significantly from each other, and so do the mean vorticities profiles in figure 17. This is because the unstable mode grows and evolves in the DNS and feeds back on the mean flow, whereas in our weakly nonlinear approach, the structure of the mean flow distortion is fully determined by the Reynolds stress of the linear response.

Overall, it is plausible to think of the departure of the fully nonlinear solution from the weakly nonlinear one, observable in figure 9 for large $U_0$, as resulting from a shear instability. This instability reaches its maximum growth rate around $t=20$, and requires some time for the unstable mode to gain in amplitude and for its nonlinear evolution, taking a tripolar shape, to dominate the energy of the response from $t=40$ and for $U_0 = 2.82 \times 10^{-2}$ in figure 9. This goes in the sense of Carton & Legras (Reference Carton and Legras1994) and Kossin et al. (Reference Kossin, Schubert and Montgomery2000), who have also shown shear instability modes to saturate into a tripolar state. The weakly nonlinear expansion does not predict an amplitude for the unstable mode, for the latter only declares as a ‘secondary’ mode, on the top of the optimal response. Therefore the predictions from the amplitude equation, for the optimal response, depart from DNS results after the time needed for the shear-driven unstable mode to become dominant.