A weakly nonlinear amplitude equation approach to the bypass transition in the two-dimensional Lamb–Oseen vortex

Abstract We analytically derive an amplitude equation for the weakly nonlinear evolution of the linearly most amplified response of a non-normal dynamical system. The development generalizes the method proposed in Ducimetière et al. (J. Fluid Mech., vol. 947, 2022, A43), in that the base flow now arbitrarily depends on time, and the operator exponential formalism for the evolution of the perturbation is not used. Applied to the two-dimensional Lamb–Oseen vortex, the amplitude equation successfully predicts the nonlinearities to weaken or reinforce the transient gain in the weakly nonlinear regime. In particular, the minimum amplitude of the linear optimal initial perturbation required for the amplitude equation to lose a solution, interpreted as the flow experiencing a bypass (subcritical) transition, is found to decay as a power law with the Reynolds number. Although with a different exponent, this is recovered in direct numerical simulations, showing a transition towards a tripolar state. The simplicity of the amplitude equation and the link made with the sensitivity formula permits a physical interpretation of nonlinear effects, in light of existing work on Landau damping and on shear instabilities. The amplitude equation also quantifies the respective contributions of the second harmonic and the spatial mean flow distortion in the nonlinear modification of the gain.


Introduction
The two-dimensional axisymmetric Lamb-Oseen (Gaussian) vortex flow, the vorticity of which is a strictly decreasing radial function, is linearly stable: the eigenvalues of L, the linearized Navier-Stokes operator around this flow, all possess a positive or null damping rate.In fact, even in the absence of viscosity where all damping rates are null, linear perturbations experience an inviscid exponential decay; this phenomenon, called 'Landau damping', was observed and interpreted for instance in Schecter et al. (2000).In particular, it was analytically derived that the Landau damping rate can be related to the vorticity gradient, at the specific radius where the angular velocity associated with the dominant Landau pole is equal to that of the base flow.
The Landau damping can be interpreted in mathematical terms.In an inviscid framework, Schecter et al. (2000) consider at first a 'top-hat' base vortex, for which the vorticity decreases slowly then rapidly drops to zero when the radial coordinate r is greater than a specific value r v .This vortex supports a continuum of modes whose critical layers are located at some r ≤ r v , where they are singular since the base vorticity gradient is non-zero there; it also supports a discrete ('Kelvin', according to the terminology of, for instance, Balmforth, Llewellyn Smith & Young 2001) mode whose critical layer is located at r c ≥ r v where the base vorticity gradient is exactly zero: for this reason, the discrete mode is smooth, regular and can be interpreted as a classical eigenmode.Secondly, Schecter et al. (2000) consider a base vortex that is equivalent to the first, with the addition of a low-vorticity 'skirt' that extends radially to a new, larger, r v ≥ r c .This introduces a non-zero base vorticity gradient at r c which ruins the regularity of the previously discrete mode.However, Schecter et al. (2000) argue that symptoms of the original discrete mode remain; some of the continuum modes closely resemble the latter in terms of structure and frequency and combine to form what Schecter et al. (2000) call a 'quasi-mode'.Therefore, if the discrete mode is used as an initial condition, it will excite the continuum following a Lorentzian distribution that is peaked around the discrete mode frequency (see figure 3 in Schecter et al. 2000).As time evolves, the continuum modes disperse, and their superposition behaves like an exponentially damped version of the original discrete mode (hence the appellation 'quasi-mode').Schecter et al. (2000) also consider the response of a Gaussian vortex, such as the one considered here, to a generic external impulse.Although the response does not behave as a single damped wave but projects well on a very large number of structurally different modes, the Landau damping is still found to be relevant and to dominate the initial decay of the perturbation.
Nevertheless, Rossi, Lingevitch & Bernoff (1997) evidenced that the Gaussian vortex flow, despite its linear stability, could relax to a new non-axisymmetric, called 'tripolar', state when subject to an arbitrary perturbation of sufficiently large amplitude.Such phenomenology is symptomatic of a subcritical bifurcation.This tripolar structure is well described for instance in Nolan & Farrell (1999) as a vortex for 'which the low vorticity of the moat pools into two satellites of an elliptically deformed central vortex, with the whole structure rotating cyclonically'.It has been observed in the laboratory experiments of Denoix, Sommeria & Thess (1994) as well as in Van Heijst & Kloosterziel (1989), Kloosterziel & van Heijst (1991) and Van Heijst, Kloosterziel & Williams (1991), and described in this last article as being 'a very stable structure, even persisting in a highly sheared fluid environment'.
This corroborates its importance in geophysical contexts, where tropical cyclones sometimes show rotating elliptical eyes, as was reported for instance by Kuo, Williams & Chen (1999) for Typhoon Herb, which occurred in Taiwan in 1996.A radar located on the Wu-Feng mountain could measure the horizontal distribution of maximum reflectivity for the Typhoon, from which we observe an elliptical eye (see their figures 1 and 2).The eye was rotating cyclonically with a period of 144 min.Concerned with Hurricane Olivia, which was part of the 1994 Pacific hurricane season, the work of Reasor et al. (2000) also documented an elliptical eye (see their figure 16).The ratio of minor to major axis was approximately 0.7, and the period of (cyclonic) rotation was found to be around 50 min.
A substantial body of theoretical work has therefore been devoted to the apparition and persistence of the tripolar state.Some of them reflect on the problem in terms of the Landau damping, for instance Schecter et al. (2000), Le Dizès (2000), Balmforth et al. (2001), Turner & Gilbert (2007) and Turner, Gilbert & Bassom (2008).Common to all these latter works is the following idea: if the perturbation is large enough to nonlinearly feedback on the mean vorticity (averaged in the azimuthal direction), in such a way that the mean vorticity gradient vanishes near the radius where the angular velocity associated with the dominant Landau pole equates to that of the base flow, then the Landau damping is deactivated.Indeed, it was shown in Schecter et al. (2000), Turner & Gilbert (2007) and Turner et al. (2008) that such cancelling of the Landau damping goes with the appearance of an undamped Kelvin mode.The effects of flattening the mean vorticity distribution have been thoroughly studied in Schecter et al. (2000) and Turner et al. (2008), although it was introduced rather artificially in order to a posteriori mimic nonlinear effects.It has, however, been rigorously quantified using a matched asymptotic expansion in Balmforth et al. (2001), the small parameter being directly linked to the amount of vorticity near the critical radius r c of the neutral Kelvin mode of a compact vortex, whose vorticity there would be zero otherwise.The developments result in an amplitude equation for the weakly nonlinear quasi-mode, that can predict a secondary instability for a sufficiently strong disturbance amplitude.
Under certain conditions, the vorticity tripole has also been shown to be the nonlinear fate of a shear instability (also sometimes called 'barotropic' instability, in opposition to 'baroclinic' instability, this latter requiring density stratification).This was illustrated clearly for instance in Carton, Flierl & Polvani (1989), Carton & Legras (1994), Carnevale & Kloosterziel (1994) and Kossin, Schubert & Montgomery (2000), and many other works.If the mean vorticity profile (averaged along the azimuthal direction) presents a local extremum at some radius, a necessary condition for shear instability is satisfied according to a generalization of the Rayleigh theorem of inflectional point by Billant & Gallaire (2005).This is, for instance, the case of the family of shielded monopoles, where the vortex core of positive vorticity is surrounded by a ring of negative vorticity.By increasing the intensity of the shear, the shielded vortex becomes unstable with a maximum growth rate for perturbations of azimuthal wavenumber m = 2; by increasing the shear further, m = 3 becomes even more unstable (see figure 7 in Carnevale & Kloosterziel 1994).Inspired by velocity measurements of Hurricane Gilbert that occurred in 1988, Kossin et al. (2000) considered the vorticity of a piecewise-constant vorticity profile.The latter is constituted of four distinct regions of vorticity: an inner region of very high vorticity, a moat region of relatively low vorticity, an annular ring of positive vorticity and an irrotational far field.Kossin et al. (2000) then show that, by narrowing the moat, a shear instability appears with a maximum growth rate at a wavenumber m = 2 (see their figure A1).This shear instability can be conceptualized as resulting from a wave interaction across the moat, between two Rossby waves riding respectively the inner edge of the annular ring, and the outer edge of the central vortex.Note that what we designate here as a 'Rossby' wave, according to the terminology of Kossin et al. (2000), is similar to the 'Kelvin' mode discussed until now.Nonlinear simulations have then shown this instability to saturate into a tripolar state.
Furthermore, an important ingredient to the subcritical transition towards the tripolar state was found to be the non-normality of the linear operator L.An operator is non-normal if it does not commute with its adjoint, the expression of the latter being relative to the choice of an inner product.The eigenvectors of a non-normal operator do not form an orthogonal basis under this inner product, thus a negative growth rate for all eigenvalues is not a guarantee for the associated norm to decay monotonically.In fact, some small-amplitude perturbations may experience a large transient amplification.In particular, the linear transient growth analyses conducted for instance in Antkowiak & Brancher (2004), Antkowiak (2005), Pradeep & Hussain (2006), Heaton & Peake (2007) and Mao & Sherwin (2012), have revealed the Lamb-Oseen (and Batchelor) vortex flow to support such strong transient energy growth.Therefore, perturbations to the base flow field, such as free-stream turbulence or acoustic disturbances, can be amplified strongly enough through non-normal, linear, mechanisms to lead to a regime where nonlinearities come into play; the flow may then escape from its linearly stable solution.This conjunction of non-normality then nonlinearity corresponds to the 'by-pass' scenario proposed in Trefethen et al. (1993) and contextualized to the Lamb-Oseen vortex flow in Antkowiak (2005).Recently, the transient growth analysis of the Lamb-Oseen vortex has been numerically extended in the fully nonlinear regime via a Lagrangian optimization in Navrose et al. (2018).
The present work aims at reconciling the simplicity of a weakly nonlinear model (such as in Sipp 2000;Balmforth et al. 2001), in the sense that it is easier to solve and interpret than the original equation, with non-normality.Specifically, the objective is to construct an amplitude equation that is not restricted to the description of close-to-neutral modes or quasi-modes, but that extends to responses associated with optimal transient growth.
The amplitude equation analytically derived in this paper will not restrict the shape of the base flow in order for the latter to support a close-to-neutral mode or a quasi-mode, and will tolerate arbitrary temporal dependence of this base flow; this is precisely because these complexities are already incorporated in the optimal transient growth analysis of which we study the weakly nonlinear continuation.This makes possible the weakly nonlinear analysis of optimal responses on vortices with more realistic vorticity profiles from field measurements.For instance, this could be applied to the profile reported in figure 1 of Kossin et al. (2000), from flight-level radar measurement of Hurricane Gilbert.
The derivation of a weakly nonlinear reduced-order model will make it possible to distinguish the regimes where weak nonlinearities reinforce the transient gain, from the regimes where they cause it to decrease.It will also provide a rough criterion for the minimum amplitude of the initial condition required to trigger a bypass transition away from the axisymmetric state.It will also help us quantify the importance of the distortion of the flow averaged in the azimuthal direction, called 'mean flow', with respect to the importance of the second harmonic in nonlinear effects at stake.
After a brief derivation on the linear formulation in § 2, the method advanced to derive the amplitude equation is outlined in § 3; specifically, we vary the amplitude of a given initial condition and predict, at low numerical cost, the gain of the response at a selected time t = t o .The (general) method is then illustrated with the two-dimensional Lamb-Oseen vortex flow with azimuthal wavenumber m = 2 in § 5, exhibiting large gain and subcritical bifurcation.

Linear formulation
In the following subsection, the formulation of the linear transient growth problem is briefly recalled.In the rest of the paper, we shall consider a purely two-dimensional flow, invariant and with zero velocity in the axial direction.This two-dimensionality implicitly assumes the flow not to be subject to any three-dimensional instabilities, or any kind of spontaneous axial variations.This assumption is often made in considering vortex 976 A10-4 flows, and will not be discussed further in the present paper.Note that the restriction to a two-dimensional flow is not intrinsic to the analytical development proposed in the following, but is simply made to ease the computations thus concentrating our efforts on the analysis of the subcritical transition toward the tripolar state.
Let U b (r, t) = [0, U b,θ ] T (r, t) denote a reference vortex flow, satisfying the Navier-Stokes equations exactly (without body force) and supporting a small-amplitude perturbation field of the form r, t).The invariance of the base flow along the azimuthal (θ) coordinate justifies the Fourier mode expansion of the perturbation in this direction; m ∈ Z denotes the wavenumber in the azimuthal direction.Linearizing the Navier-Stokes equations around U b (r, t) leads to an equation for the temporal evolution of the velocity field û(r, t) where ∇m The letter Ω denotes the angular velocity of the base flow, where W z is its vorticity along the z-axis.The operator imΩ is the advection by the base flow, and Δ m is the Laplacian associated with viscous diffusion: it is therefore systematically multiplied by 1/Re, with Re the Reynolds number.The presence of the pressure term in (2.2) ensures the velocity field û to be divergence-free for all times.Indeed, the pressure is linearly linked to the velocity field according to the Poisson equation obtained by taking the divergence of the momentum equations, then enforcing the continuity (∂ r + 1/r)û r + imû θ /r = 0.The perturbation fields are subject to the following boundary conditions over r ∈ [0; +∞[, valid for all times: (2.5a,b) as well as lim r→∞ û = 0.As shown in the appendix of Kerswell & Davey (1996), by imposing the parity conditions (2.5a,b), 'the correct axial behaviour automatically follows without need to explicitly impose the regularity conditions'.
Only the temporal dependence of the operators will be made explicit in the rest of the paper; for instance, L m (t), whose temporal dependence is inherited from the base flow, is actually L m (t) = L m (r, t; Re).Precisely due to the fact that it depends on time, the operator exponential formalism cannot be used to solve (2.2).Instead, and given the value of the perturbation field at a time t i , its temporal evolution at a further time t according to (2.2) formally reads û(t) = Ψ (t, t i ) û(t i ).The operator Ψ (t, t i ) is traditionally named the propagator, for its action directly maps û(t i ) onto û(t) (Farrell & Ioannou 1996).If t i = 0 in particular, û(t) = Ψ (t, 0) û(0). (2.6) The propagator responds to the semigroup properties and 2), and enforcing the equation to be satisfied for all û(t i ) leads to an evolution equation for the propagator (2.8) By evaluating (2.7) at t = t i and s = t, it follows that (2.9) Eventually, taking the temporal derivative of the first equation in (2.9) results in an evolution equation for the inverse propagator (2.10) Note that in presence of a forcing term f (t) at the right-hand side of the system (2.2), its solution (2.6) generalizes into (2.11)This formula will turn out to be useful in a moment.
The transient growth analysis amounts to constructing an orthonormal basis for the initial flow field, with the particular feature that elements of this basis are ranked according to how much they are amplified at a given temporal horizon t o > 0. The first element of this basis is the initial condition that is most amplified at t = t o , the second is the most amplified under the constraint that it must be orthogonal to the first, etc. Due to the non-normality of L m (t), the first few elements often lead to amplifications that are much larger than all the others.In the case where the initial condition projects equivalently on each element of the basis, the flow response at t = t o is dominated by those stemming from the few first elements of the basis only.In other terms, the key idea behind a transient growth (and/or a resolvent) analysis is that the response of a system to external perturbations is intrinsic to the operator, thus the precise form of these external perturbations matters little.
The first step is to define according to which measure the amplification is sought for.In the following, we choose the norm induced by the Hermitian inner product (2.12) the superscript 'H' denoting the Hermitian transpose.Note that û| û = û 2 is directly the kinetic energy of the perturbation.The largest linear amplification (or 'gain') between an initial flow structure and its evolution at t = t o > 0 (subscript o for 'optimal') reads (2.13) The optimal gain does not depend on the time itself, but only on the temporal horizon t o .In the following, G o alone implies G o (t o ).The propagator formalism (2.6) is used to re-write the maximization problem (2.13) as an eigenvalue problem where the operator Ψ (t o , 0) † denotes the adjoint of Ψ (t o , 0) under the inner product (2.12).From (2.14), it is clear that G 2 o is also the largest (necessary real) eigenvalue of the self-adjoint operator Ψ (t o , 0) † Ψ (t o , 0), and the associated eigenvector, named ûo in what follows, is the optimal initial condition.The latter is normalized such that ûo | ûo = 1.Smaller eigenvalues and corresponding eigenmodes constitute sub-optimal gains and initial conditions, and the family of eigenvectors is orthonormal.Let lo be the unit-norm response to the optimal initial condition ûo at t = t o , such that lo .= Ψ (t o , 0) ûo /G o and lo | lo = 1.From this definition, and using that o ûo and that the inverse of the adjoint is the adjoint of the inverse, two relations follow where we recall that o is the inverse of the optimal gain.Note that these two last relations also indicate the optimal initial condition and its associated response to be respectively the right and left singular vectors of Ψ (t o , 0), associated with its largest singular value.
Since the operator L m (t) is assumed in the rest of the study to be strongly non-normal, the structures ûo and lo generally project poorly each of its direct and adjoint eigenmodes (except in the limit t o → ∞).Of all the temporal horizons t o , the one leading to the largest optimal gain will be highlighted with the subscript 'max' such that max . (2.17)

Weakly nonlinear formulation
In the linear paradigm, the gain is independent of û(0) , for the latter is assumed to be arbitrarily small.In reality, û(0) may be sufficiently large for the nonlinear corrections to the response not to be negligible anymore, thus for the transient gain to depart from its linear prediction.Building on our previous work (Ducimetière, Boujo & Gallaire 2022), we propose thereafter a method for capturing weakly nonlinear effects on the transient gain over a time-varying base flow.
In the rest of the present study, we subject the two-dimensional, unforced, Navier-Stokes equations governing the flow field U(t) to a small-amplitude initial perturbation around the reference vortex flow U b (t).The initial perturbation has an azimuthal wavenumber m and its radial structure ûh , with ûh = 1, is for now arbitrary.The strong non-normality assumption justifies further assuming the transient gain to be large, that is to say, o 1, such that o constitutes a natural choice of small parameter with which to scale the amplitude of the initial perturbation.Specifically, where the amplitude of the initial condition, U 0 .= α √ o 3 , can vary through the real prefactor α = O(1).We intend at capturing the variation of the transient gain by increasing U 0 .
For this purpose, the total flow field U is developed as a multiple-scale asymptotic expansion in terms of powers of √ o around the reference vortex flow U b (t) The slow time scale T .= o t was introduced, aiming at capturing the slow modulation of the linear trajectory as nonlinearities progressively set in.The reason for which the expansion and the scaling of the initial condition are made in terms of powers of √ o shall be specified later.Injecting (3.2) in the Navier-Stokes equations, then expanding each u j as a Fourier series in space according to where the nonlinear advection operator (3.5) has been defined.The fundamental field, corresponding to n = 1, has been isolated at each order in (3.4), and the harmonics have been incorporated in sj with  1) 1) We recall that the application of the operator Ψ (0, t) is equivalent to integrating the system backward from t to 0, and that it maps the optimal response lo on a field of very small amplitude o 1 (by assumption) in (2.15a,b).The main idea behind our method is to take advantage of this by perturbing the operator Ψ (0, t) for all t ≥ 0 according to and where the Heaviside distribution H(t) is defined as H(0) = 0 and H(t > 0) = 1.The operator l(t)| * applied to any field ĝ simply results in l(t)|ĝ .The operator Φ(0, t) satisfies Φ(0, 0) = I and Φ(0, t) l(t) = 0 for t > 0, and therefore is singular for all strictly positive times; the linear trajectory l(t) constitutes its non-trivial (time-evolving) kernel.We further show in Appendix A that the non-trivial kernel b(t) of the adjoint operator , where the term in parenthesis is an orthogonal projection operator: it is self-adjoint, linear, idempotent and its application projects into the subspace that is complementary to l(t).Therefore (3.9), by stating that Ψ (0, t) ≈ Φ(0, t), implies that applying Ψ (0, t), or applying Φ(0, t) that first removes the component on l(t) then applies Ψ (0, t), both lead to very similar initial states.That is precisely because Ψ (0, t) maps l(t) on o ûo of very small amplitude o 1.In principle, expansion (3.9) requires Ψ (0, t) o P(t) = o / l(t) = 1/ Ψ (t, 0) ûo ; by using that the norm of the inverse of an operator is the inverse of its minimum singular value In other terms, the operator perturbation a priori holds for the times t such that the response to the initial condition that is the least amplified at t, is much smaller than the response to ûo that is the most amplified at t o .This is certainly verified for times around t = t o , which is appropriate, for the response l(t) is expected to be nonlinearly modified at first at these times.By then perturbing the operator Ψ (0, t) in (3.8) according to (3.9), we are left with 1) 1) 1) ū( 1) Note that the perturbation operator o P(t), modifying Ψ (0, t) in Φ(0, t) at leading order 3 ) specifically.This was made purposely and explains a posteriori why the asymptotic expansion was outlined in terms of integer powers of √ o , instead of being for instance in terms of integer powers of o .More precisely, we anticipated the forcing term C [u 1 , u 2 ] + C [u 2 , u 1 ] appearing at third order, to have a non-zero component on the fundamental wavenumber m, due to the bi-linearity of the operator C [ * , * ]; therefore, making the term in P(t) appear at third order as well, was a way of putting all forcing terms oscillating at m at the same order.
Note also that only the propagator associated with the wavenumber m was perturbed, although harmonics may equally lead to significant transient gains.This selective perturbation is justified a priori by the fact these harmonics are not externally excited like the fundamental pair is through the initial condition, but are only generated nonlinearly at higher orders.If, however, the harmonic responses are a posteriori found to endanger the asymptotic hierarchy, they can always be included in the kernel of their associated propagators in the manner of (3.9).
Terms in (3.12) are collected at each order in √ o , yielding a cascade of linear problems.At order This leads to ū(n) 1 = 0 for all times and n / = 1.In addition, multiplying (3.13) by Ψ (0, t), then integrating in time and using Φ(0, 0) = I, gives Φ(0, t) ū( 1) where A(T) ∈ C is a slowly varying scalar amplitude verifying ∂ t A = 0. Eventually, the general solution at √ o is written In the linear regime, A would be constant over time; therefore, by continuity, we expect its variation to be small in the weakly nonlinear regime.This is what is implied in the fact that A depends on the slow time T, At order o , the solution is for t ≥ 0, where and The advection operator in the Fourier space C m x, ŷ is as (3.5), excepted that ∇ is replaced by ∇m ; it computes the transport, by the field x, of the field ŷ oscillating at m.In principle, the Heaviside distribution H(t) multiplies the forcing terms in (3.17).However, it can be ignored here without loss of generality, for the forcing terms appear inside an integral between 0 and t in the formal expression of the solution, where the presence of a Heaviside distribution is unimportant.In (3.16), the homogeneous solution A 2 (T)H(t) l(t), solving the system Φ(0, t) ū(1) 2 = 0 for t ≥ 0, was ignored.In this manner, the component of the overall solution on l(t) is fully embedded in A, without needing to be corrected.
At order √ o 3 in (3.12) are gathered two equations: the first yields the Fourier component of the solution oscillating at m subject to ū(1) 3 | t=0 = α ûh and where we defined depending only on the fast time scale t.The second equation governs the Fourier component oscillating at 3m Integration over the time t, detailed in Appendix B, leads to where The operator Φ(0, t) being singular for strictly positive time, (3.22) admits a non-diverging particular solution if and only if its right-hand side is orthogonal to b(t) for all t > 0 (Fredholm alternative).This condition results in for t > 0. The adjoint field b(t) = Ψ (t, 0) † l(t) tends towards b(0) = l(0) = o ûo when t tends towards 0. Therefore, in this limit, (3.24) reduces to where we also used that the limits of both the left-hand side and the nonlinear term in (3.24) are 0 when t goes to 0; the former because of the presence of the factor t, and the latter because ũ(0) = 0. and where we defined the nonlinear coefficient From here, the weakly nonlinear transient gain corresponding to the wavenumber m for t > 0 is simply The linear regime corresponds to the limit U 0 / o → 0, or simply U 0 → 0, since if the amplitude of the initial condition is much smaller than the linear gain, then the amplitude of its response is much smaller than one.In this limit, the nonlinear terms in (3.28) becomes negligible, and the solution tends towards which is the properly normalized projection on l(t) of the response to the initial condition U 0 ûh , as expected at a linear level since a(t) is the amplitude along l(t).
Note that in (3.9) any other perturbation operator of the form P = H ûo x| * / x| l , with x any trajectory, would also have led to a singular operator with l as a non-trivial kernel, but with b = Ψ (t, 0) † x as a non-trivial kernel of its adjoint.Choosing x = l implying b = Ψ (t, 0) † l, as was done so far, is the only choice leading to a coherent result in (3.31), thereby guaranteeing the uniqueness of P in (3.9); as a side comment, it should be noted that among all x, x = l also leads to the P of the smallest possible norm.
For the gain at t = t o , (3.31) corresponds to lim The result (3.32) also is as expected from the linear prediction.In particular, we recover lim when the optimal initial condition is applied ( ûh = ûo ).It also predicts the gain to be null when ûh is orthogonal to ûo , indicating the linear response at t = t o to be orthogonal to lo without taking into account the gains associated with sub-optimal forcings.Therefore, these latter should be O(1/ √ o ) to be mathematically consistent.For the rest of the paper, we set ûh = ûo .In physical terms, this choice, although very specific, is found to be the most natural in the absence of information regarding the structure of the actual initial condition ûh .Note, however, that ûh might project very poorly on ûo , and in the latter case reducing the dynamics to the response to ûo would give inaccurate results.
Further expressing a(t) in terms of an amplitude |a(t)| ∈ R + and a phase φ(t the subscripts 'r' and 'i' denoting the real and imaginary parts, respectively.The equation for |a|, in particular, bears the following analytical solution: For any time t, the gain decreases strictly monotonically with U 0 if μ r (t) < 0, thus is maximum in the linear regime.Moreover, the gain stays constant with U 0 if μ r (t) = 0, and increases strictly monotonically if μ r (t) > 0. In the latter case, by increasing U 0 the gain eventually becomes singular (tending towards +∞ at time t) for In the following, we call (3.34) the weakly nonlinear non-normal transient (WNNt) model.In Appendix C we show that at first order in the gain variation, (3.34) partly reduces to the sensitivity of the gain G(t) to the axisymmetric base flow modification +a 2 l u (0) is the amplitude of the response in the linear regime.In this sense, our model states that small gain modifications by increasing U 0 , are partly due to the fact that the perturbation now evolves over a base flow that is distorted through the action of the Reynolds stress forcing −C m a l l * , a l l + c.c. of this very same perturbation (nonlinear feedback).In addition, (3.34) also includes the effects of the second harmonic oscillating at 2m.
In Appendix D, a second method that is more immediate, although perhaps less rigorous, is proposed to derive the amplitude equation (3.24).This method does not rely on the singularization of the propagator, and was named the 'pseudo-forcing method'.

Linear and fully nonlinear transient growth in the diffusing, two-dimensional
Lamb-Oseen vortex Equation (3.34) is employed thereafter for the analysis of the transient gain experienced by a two-dimensional Gaussian vortex in a weakly nonlinear regime.
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 The associated vorticity field, W z , has a Gaussian radial profile that diffuses with time for finite 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. 1997).Only the response to an initial condition with m = 2 will be considered in this paper.

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: 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π/(N − 1)), with j = 0, 1, . . ., N − 1, which includes the endpoints s = −1 and s = 1.This grid is mapped on the physical domain 0 ≤ r ≤ 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 Canuto et al. 2007;Viola, Arratia & Gallaire 2016).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 Ψ (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 [ û(n) , p(n) ] T at a time t n , their values at t n+1 = t n + t are obtained upon inverting the system 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 where we used M † = M, and where the expression of B m (t) † and the boundary conditions for the adjoint field are derived in Appendix E. The action of the adjoint propagator Ψ (t, 0) † , 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 − t are obtained by solving A time step of 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 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.

Linear results
The results of the linear transient growth analysis for m = 2 and varying Re ∈ [1.25, 2.5, 5, 10] × 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 The type of dependence is made explicit in figure 1(b,c), where we propose in (b) to plot o,max multiplied by Re 1/3 as a function of Re, demonstrating that G o (t o,max ) ∝ Re 1/3 .In (c), t o,max multiplied by Re −1/3 is also shown as a function of Re: while the data align slightly less well on a constant line, the agreement remains correct and stating t o,max ∝ Re 1/3 is a fair approximation.All these results are in excellent agreement with those presented in chapter 3 of Antkowiak (2005).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 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 (35) and that G is below G o everywhere else.At the times corresponding to the black dots, the axial vorticity structure, expressed as 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 (2004) and Antkowiak (2005) 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 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).
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 (1983) for a passive scalar and in Lundgren (1982) for vorticity.In particular, Lundgren (1982) 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 (1994), where is further demonstrated that the decay of the perturbation vorticity spirals acts on a Re 1/3 time scale.From here, Antkowiak (2005) argues that, in the limit of vanishing viscosity, the evolution equation for the perturbation vorticity ω 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 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 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 (1987) in the context of plane shear flows.As a side comment, Antkowiak (2005) also derives the Re 1/3 dependence of t o,max by simply balancing the unfolding and the diffusion (in the radial direction) time scales.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 (1993) have shown that this mechanism was associated with the scaling law G o (t o ) ∝ t o (note that the gain defined in Farrell & Ioannou (1993) is the square of the gain presently defined).From here, the scaling law G o (t o,max ) ∝ 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 ∝ Re 1/3 , then the decay rate is clearly independent of 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. 2000), 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 (1970).For this, one first needs to define the mth multipole moment of the vorticity perturbation as Then, to quote Schecter et al. (2000), 'A Landau pole is a complex frequency ω q − iγ at which the Laplace transform of of the form exp(−γ t − iω 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. 2000); 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. (2000) 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 [. ..] is poorly described by a single damped wave [. ..]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. (2000).In other terms, the response excites a very , where we use fitted values for ω q and γ .
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(iφ(t)), we display in figure 4(a) the temporal evolution of the phase φ(t) normalized by 2π 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(−iω q t − γ t) is also proposed.Clearly, the phase changes at a constant rate in figure 4(a) for 15 ≤ t ≤ 65, and |Q (2) (t)| decays exponentially in figure 4(b) for 40 ≤ t ≤ 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. (2000).The Reynolds number Re = 5000, however, seems sufficiently large for the conclusions drawn in Schecter et al. (2000) 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 Re is finite possibly makes it less rigorous.The precise value of γ 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. (2000).The fitted value ω q = 0.42 in figure 4(a) is, however, important, for it gives the location r q of the radius where the angular velocity ω q /m associated with the Landau pole is equal to that of the base flow mΩ(r q (t), t) = ω q . (4.9) The radius solving (4.9) at each time is highlighted by a dashed circle in figure 3. It takes the value r q ≈ 2.16, very weakly dependent on time.
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. (2000), Turner & Gilbert (2007) and Turner et al. (2008), 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. (2000) and confirmed in figure 3, is the reason for which 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 × 10 −2 , black dots are located at t = 0, 30, 60, . . ., 210, for which the corresponding vorticity structures are shown in figure 6. non-modal analytical tools are deployed in the present study, in particular for studying nonlinear effects.We address these latter now.

Fully nonlinear results
We introduce the fully nonlinear results obtained from DNS in the present section.Let us define u p as being the perturbation between the fully nonlinear solution U obtained from a DNS, and the reference Lamb-Oseen solution U b in (4.1a,b); in other terms, u p .= U − 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 u p | t=0 = U 0 ûo e imθ + c.c.. (4.10) Owing to nonlinear effects, u p generally does not oscillate purely along m for strictly positive times.For this reason we need to further extract û(1) p , the component of u p oscillating at the fundamental wavenumber m, hence the superscript '(1)'.For this, a Fourier series is naturally used as From here, the fully nonlinear gain (associated with the fundamental pair) is immediately defined as (4.12) By fixing Re = 5000 and t o = 35, we report in figure 5(a) the fully nonlinear evolution of the perturbation amplitude û(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) illustrates well a bypass transition, also called 'bootstrapping effect' in Trefethen et al. (1993), 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 × 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 ≥ 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 (2007) for Re = 1000, or with the experimental visualization in Kloosterziel & van Heijst (1991).
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).
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.

Weakly nonlinear transient growth in the diffusing, two-dimensional Lamb-Oseen vortex
In the weakly nonlinear paradigm, the perturbation field u p is approached by p is approximated by ), therefore the fully nonlinear gain G DNS defined in (4.12) is expected to reduce to the weakly nonlinear one G w defined in (3.30), since o 1.These two quantities are compared in this section.Note that we checked that for all considered Re values the first sub-optimal transient gain was O(1/ √ o ), in contrast to 1/ o for the optimal one, which mathematically justifies focusing on the response to ûo only.

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 (2) r .Remarkably, after having significantly decreased until reaching a minimum at t = 59, the coefficient μ 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 ≥ 150.A change of sign in μ r brings diversity to the behaviours reported in Ducimetière et al. (2022), 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 μ r seems largely dominated by the contribution from the mean flow distortion μ 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 μ r > 0. In the white zone, it decreases monotonically.(b) The coefficient μ r is compared with its reconstruction from DNS data (magenta dashed line) according to (5.2). positive, however, μ (0)  r 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 2 0 , we obtain (5.1) Therefore, the coefficient μ r relates directly to the rate of change of G 2 w with respect to U 2 0 in the limit where U 0 tends towards zero as .
(5.2) By using (5.2), the coefficient μ 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 × 10 −4 and U 0 = 2.82 × 10 −4 .The result is shown as the magenta dotted line in figure 8(b), and compared with the actual μ r .The good agreement between both curves for t ≤ 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, P(t) = 1/ l(t) , increases accordingly.

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 μ 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 ≤ t ≤ 40 in the panel corresponding to 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.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 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 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.

Bifurcation thresholds
By looking at the panels corresponding to U 0 = 1.12 × 10 −2 and U 0 = 1.78 × 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 × 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 s 0 is what 976 A10-25 , 10] × 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 Re (also with a star symbol).The prediction from the weakly nonlinear model, U s 0 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 Re, with ∝Re −0.88 , whereas the thin dashed line is fitted on the weakly nonlinear data with ∝Re −0.66 .The inset shows the same in log-log scale.
follows, reads and is defined if and only if μ 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 (2005), based on the observation that the tripolar state is characterized by a deformed, non-axisymmetric, heart.Therefore a characteristic aspect ratio λ of the heart is established as and J mn .
= Ω x m y n ω(x, y) dx dy a vorticity moment.A characteristic eccentricity is further defined as e .= 1 − 1/λ 2 .From here, Antkowiak (2005) computes the typical half-life time (that is where the criterion is rather arbitrary) of the heart deformation as 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 Re values in figure 11(a).For each Re, we also highlight an inflection point in τ at a certain U 0 , and we declare the latter as being the subcritical bifurcation threshold.It is reported directly as a function of Re in figure 11(b), and compared with the weakly nonlinear prediction U s 0 , in both linear and log scales.Both approaches clearly highlight a decreasing power-law dependence of the subcritical bifurcation threshold with the 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 (2007), which is −0.8.However, since little information is given  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 ∝Re −0.69 , whereas the thin dashed line ∝Re −0.66 is similar to figure 11.
regarding how the threshold amplitude was computed in Antkowiak & Brancher (2007), the agreement with our results is perhaps fortunate.In any case, the negative power-law dependence with 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 Re.This may explain the formation and persistence of elliptical vortex eyewalls in some tropical cyclones (Kuo et al. 1999;Reasor et al. 2000).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 (2010).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 Re in figure 11 (b).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 × 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 τ 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 s 0 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.

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, μ (0)  r 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 Re = 5000.First, the energy of the mean flow distortion u (0)  2 is shown as the red curve in figure 13(a).For comparison, the energies of the linear response and of the second harmonic are also shown.As written in (3.17 2 can persist for extremely long times, even when the linear response, whose nonlinear interactions force u (0) 2 , has vanished.This suggests that u (0) 2 can persist even in the absence of sustained forcing.This is also illustrated in figure 14, where we show ω (0) 2 , the vorticity structure associated with u (0) 2 , 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 u (5.8)  −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 ω q = 0.42.
The panels in the first row in figure 15 show the evolution of f 2,θ grows in amplitude (since l does) while roughly conserving its shape, which u (0)  2,θ 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).
The panels in the second row in figure 15 show the evolution over the times 35 ≤ t ≤ 85.Notice that the forcing f (0) 2,θ 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 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 (0) 2,θ tends to invert in figure 15(c), which implies its momentary flattening, therefore u (0) 2,θ 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 l.
After t = 85 the velocity u (0) 2,θ evolves freely, for f 2,θ is negligible.The persistence of u (0)  2,θ 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,θ (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, . . ., 35 (on the left of the first vertical line in figure 13a) and the second for t = 35, 40, . . ., 85 (between the two vertical lines in figure 13a).Larger times correspond to lighter colours.specific radius r q where the angular velocity ω 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 ω q = 0.42, resulting in r q ≈ 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 ω (0) 2 to the reference vorticity has a positive contribution to the total vorticity slope at r q .In parallel, we recall that μ (0) r can be directly interpreted as the sensibility at a time t of the transient gain to the addition of u (0) 2,θ to the reference flow.For this, μ (0) r computes the integrated effect of this addition between 0 and t (see formula (C10)).Therefore, the decreasing tendency of μ (0) r until t = 60 in figure 8, results from the negativity of the slope of ω (0) 2 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. (2000), Turner & Gilbert (2007) and Turner et al. (2008).On the contrary, from t = 60 onward the coefficient μ (0) r increases, for the presence of ω (0) 2 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 976 A10-30 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. (2001).In addition, the present amplitude equation leads to the conclusion that, when μ 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 × 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,θ + a 2 u (0) 2,θ , 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 û(r, t) .= ȗ(r) e σ t , and solving the subsequent eigenvalue problem.The azimuthal mean is linearly unstable if and only if there exists an eigenvalue σ such that Re(σ ) = σ 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.
For 10 ≤ t ≤ 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.phase speed, thus the reinforcement of the one by the other is preserved in time, eventually leading to instability.
From 20 ≤ 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 × 10 −2 in figure 9.This goes in the sense of Carton & Legras (1994) and Kossin et al. (2000), 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.

Summary and perspectives
In conclusion, we believe our work to have brought a twofold generalization to existing literature.The first lies on a purely methodological level.We have derived an amplitude equation for non-normal systems, describing the transient response to an initial condition, in a weakly nonlinear regime.Unlike in Ducimetière et al. (2022), the reference state of these systems can now depend arbitrarily on times, owing to the propagator formalism, without the need for this latter to take its particular operator exponential shape.This offers numerous possibilities of applications, and weak nonlinearities could be modelled, for instance, in pulsating pipe flows, which play a key role in the hemodynamic system of many species (Pier & Schmid 2017, 2021); it could also be applied to time-dependent stratified shear flows, which have revealed to support strong transient growth, for instance in Arratia, Caulfield & Chomaz (2013) and Parker et al. (2021).Not only could the weakly nonlinear evolution of the gain associated with the linear optimal perturbation be captured, but the amplitude equation could also be included in a Lagrangian optimization problem whose stationary conditions would constitute a weakly nonlinear optimal, parameterized by the amplitude of the initial condition.
The method does not rely on any modal ('eigenmodal') quantities, therefore the existence of a continuous spectrum and/or the absence of discrete eigenmode is not problematic.Corollary, the shape of the reference flow is not constrained, apart from the fact that it should lead to strong energy growth at some finite time.The equation is derived for the amplitude of the time-dependent (linear) optimal response, whose computation already encompassed the entirety of the spectrum, regardless of its precise nature.This was particularly convenient when applied to the two-dimensional Lamb-Oseen (Gaussian) vortex flow, the linear optimal response of which is characterized by vorticity filaments under constant shear by the reference flow.The work of Schecter et al. (2000) has shown this response to project well onto a very large number of eigenmodes constituting the continuum, all with different shapes and frequencies.Therefore it was extremely poorly described as a single eigenmode, which invalidates the use of classical weakly nonlinear methods.
Instead, by using the amplitude equation (3.33), we could correctly predict for different Re values the nonlinear evolution of the response for 0 ≤ t ≤ 130, and for small amplitudes of the initial condition.In this specific regime, at Re = 5000 as an example, nonlinearities have been found to reduce the transient gain for 0 ≤ t ≤ 87, but to enhance it for 87 ≤ t ≤ 130, as confirmed by DNS.Owing to the simplicity of the amplitude equation and its link to the sensitivity formula, this could be further related to the creation of an azimuthal mean flow distortion from the Reynolds stress of the linear response, which affects the Landau damping controlling the non-viscous dissipation of the response.The second harmonic effect has also been found to be important over the time interval where nonlinearities reinforce the gain.
However, for relatively large amplitudes of the initial condition, the predictions of the amplitude equation are found to remain accurate at small times only.After a well-captured episode of diminution of the gain, the DNS simulations with the largest considered initial condition amplitudes depart from the weakly nonlinear prediction and evidence a bifurcation towards a tripolar state.By performing a mean flow stability analysis, this departure can be related to the emergence of a shear instability.Specifically, the mean flow distortion by the Reynolds stress of the linear response, at the same time that it enhances the Landau damping by its effect at the radius r q , generates a vorticity ring closely around the central part of the vortex.A shear instability results from an interaction between the inner 'edge' (sharp slope) of this annular ring and the outer edge of the central region.The weakly nonlinear approach does not predict an equation for the unstable mode that emerges on the mean flow.Therefore, as soon as the unstable mode dominates the response, around t = 40 for the largest considered U 0 = 2.82 × 10 −2 and Re = 5000, the weakly nonlinear description rapidly degrades in terms of both energy and structure.
Nevertheless, at larger times, the amplitude equation can still give indirect information about the fully nonlinear response.Specifically, for Re = 5000, for a given and above a certain value for the amplitude of the initial condition, the amplitude equation has no solution over a finite time interval contained in 87 ≤ t ≤ 130, where the coefficient μ r is positive.The non-existence of any solution over a time interval may imply that, at least over the same time interval, the fully nonlinear solution must have reached another nonlinear state.This seems confirmed by the DNS.In this sense, the threshold initial condition amplitude predicted by the amplitude equation could be a reasonable approximation of the fully nonlinear one, even if a direct comparison between the former and the latter is found to be delicate.
For future research, the nonlinear self-sustaining mechanism(s) of the tripolar state remain to be clarified.In this perspective, a semi-linear approach such as the one deployed in Yim, Billant & Gallaire (2020) could be appropriate.This approach relies on the assumption that the dominant nonlinear mechanism is the Reynolds stress feedback onto the mean flow, thus neglecting the nonlinearity arising from the cross-coupling between different frequencies.In this sense, it is less rigorous than the weakly nonlinear expansion proposed here.Nevertheless, the semi-linear model does not assume the fluctuation over the mean flow to be small and typically retains its spatial degrees of freedom.Therefore the nonlinear structure it predicts is possibly considerably different from that of the linear regime and could evolve towards a tripole.
Eventually, we believe that the fact that the proposed method does not make assumptions on the shape of the base profile, nor on the values of external parameters, only that the reference flow should lead to some energy growth, could be exploited further.For instance, the proposed amplitude equation could be employed to assess and interpret weakly nonlinear effects on the optimal response on vorticity profiles from actual field measurements such as those reported in Kuo et al. (1999), Reasor et al. (2000), Kossin et al. (2000) and Kossin & Schubert (2001).Indeed, the mean vorticity profile predicted in figure 17 of this work is qualitatively very similar to the one of Hurricane Gilbert shown in figure 1 of Kossin et al. (2000).This raises the question of the relevance of non-normal mechanisms combined with nonlinear effects in tropical cyclones.δΨ (t, 0) with δL m (t) by taking the variation of (2.8), which leads to ∂ t (δΨ (t, 0)) = (δL m (t))Ψ (t, 0) + L m (t)δΨ (t, 0)

Figure 2 .Figure 3 .
Figure 2. The continuous line is the reproduction of the envelope shown in figure 1(a) for 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, . . ., 105, for which the corresponding vorticity structures are shown in figure 3.

Figure 4 .
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 φ(t) divided by 2π.(b) Amplitude |Q (2) (t)|.The black dashed lines correspond to a pure Landau damping Q(2) (t) = exp(−iω q t − γ t), where we use fitted values for ω q and γ .

Figure 6 .
Figure 6.Temporal evolution of the vorticity structure of the fully nonlinear perturbation for U 0 = 2.82 × 10 −2 , shown at the specific times highlighted by the black dots in figure 5.Each panel shows only [x, y] ∈ [−4, 4] × [−4, 4], the plus sign denotes the origin, and the dashed line is the unit circle.
) for Re = 5000 and t o = t o,max = 35.It is further split into the sum of a mean flow distortion contribution Figure8.(a) Real part of the weakly nonlinear coefficient μ 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 μ r > 0. In the white zone, it decreases monotonically.(b) The coefficient μ r is compared with its reconstruction from DNS data (magenta dashed line) according to (5.2).

Figure 10 .
Figure 10.Weakly and fully nonlinear gains for Re = 5000 and increasing amplitude of the initial condition U 0 ∈ [4.5, 7.1, 11.2, 17.8] × 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.

Figure 11
Figure11.(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 Re = [1.25,2.5, 5, 10] × 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 Re (also with a star symbol).The prediction from the weakly nonlinear model, U s 0 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 Re, with ∝Re −0.88 , whereas the thin dashed line is fitted on the weakly nonlinear data with ∝Re −0.66 .The inset shows the same in log-log scale.

Figure 12 .
Figure12.Same as in figure11, 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 ∝Re −0.69 , whereas the thin dashed line ∝Re −0.66 is similar to figure 11.

Figure 13
Figure 13.(a) For Re = 5000, the evolution of the energy * 2 of the linear response l (black dashed line), of the second harmonic û(2) 2 (blue dashed-dotted line) and, mainly, of the mean flow distortion u (0) 2 (red continuous line).The red dots correspond to the specific times t = 30, 40, . . ., 100, for which the vorticity structure of u (0) 2 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. ∂ r ω (0) 2 | r=rq , as a function of time.
), the mean flow distortion is forced by f Reynolds stress of the linear response.Under its action, the energy of u (0) 2 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 ≥ 85.It is striking to notice in figure13(a) that u (0) with the wavenumber m = 0, for which the equation for the radial perturbation is ∂ r (ru (0) 2,r ) = 0 from the continuity.This leads to u (0) 2,r = 0 and a forced diffusion equation for u (0) 2,θ , reading

Figure 14 .
Figure14.Temporal evolution of ω(0)  2 , the vorticity structure of u (0) 2 , the mean flow distortion induced by the Reynolds stress of the linear response shown at the specific times corresponding to the red dots in figure13(a).Each panel shows only[x, y] ∈ [−4, 4] × [−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 ω q = 0.42.
t ≤ 35.The Reynolds stress forcing f (0) , let us study now how it feeds back on the response.We show in figure13(b) the slope of ω

Figure 15 .
Figure 15.Temporal evolution of the forcing f (0) 2,θ (panels on the left) and of the velocity response u (0)

Figure 16
Figure 16.(a) Maximum growth rate of the eigenvalues of the Navier-Stokes operator at Re = 5000 and linearized around the (azimuthal) mean flow, in order to describe m = 2 perturbations.Mean flows corresponding to U 0 = 2.82 × 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.
The full linear response seeded by o ûo and such that l(t o ) = lo reads l(t) = Solving (3.24) is equivalent to solving its partial derivative with respect to the fast time scale t, reading subject to the initial condition (3.25).The system is written solely in terms of t by evaluating (3.25) and (3.26) at T = o t, and remembering that dt A = o d T A| T= o t , b|Ψ (0, t) ũ b| ûo , with A(0) = α ûo | ûh .(3.27)The amplitude A, previously undefined at t = 0, was prolonged by continuity there by stating A(0) = lim t→0 A(t).Note that such rewriting of the amplitude equation in terms of t was not done directly for (3.24), since it would have given an ordinary differential equation without an initial condition, (3.25) being intrinsically satisfied.
max Re 1/3 Figure 1.Linear transient growth in the two-dimensional, time-dependent Lamb-Oseen vortex flow for varying Re ∈ [1.25, 2.5, 5, 10] × 10 3 , larger 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 Re, that is, to G o (t o,max ) = 1/ o,max ; (b) o,max multiplied by Re 1/3 and plotted as a function of Re; (c) t o,max multiplied by Re −1/3 and shown as a function of Re.
(a) and for different Re.For all considered 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 Re.