Weak nonlinearity for strong nonnormality

We propose a theoretical approach to derive amplitude equations governing the weakly nonlinear evolution of nonnormal dynamical systems when they experience transient growth or respond to harmonic forcing. This approach reconciles the nonmodal nature of these growth mechanisms and the need for a center manifold to project the leading-order dynamics. Under the hypothesis of strong nonnormality, we take advantage of the fact that small operator perturbations suffice to make the inverse resolvent and the inverse propagator singular, which we encompass in a multiple-scale asymptotic expansion. The methodology is outlined for a generic nonlinear dynamical system, and four application cases highlight common nonnormal mechanisms in hydrodynamics: the streamwise convective nonnormal amplification in the flow past a backward-facing step, and the Orr and lift-up mechanisms in the plane Poiseuille flow.


Introduction
Nonlinear dynamical systems can have one or several equilibrium solutions, which form one of the building blocks of the phase space Strogatz (2015). The linear stability of an equilibrium can be deduced from the eigenvalues of the linearised operator: linear modal analysis thus helps to distinguish between linearly unstable, neutral (marginally stable) and strictly stable equilibria (when the largest growth rate is positive, null and negative, respectively), and to detect bifurcations. It sometimes remains too simplistic, however, and has therefore been generalised over the last decades to account for nonlinear and nonmodal effects, although these two types of correction have generally been opposed Trefethen et al. (1993), culminating into a paper entitled "Nonlinear normality versus non-normal linearity" Waleffe (1995). The objective of the present study is precisely to contribute to reconcile nonlinearity and nonnormality, and to rigorously derive weakly nonlinear amplitude equations ruling nonnormal systems.

Weak nonlinearity
While the most unstable eigenmode eventually dominates the linear, small-amplitude dynamics, its shape and frequency may differ significantly from those of the nonlinear state when moving away from the bifurcation point. Fundamentally, the saturation amplitude can only be determined trough nonlinear considerations, and one must resort to a weakly or fully nonlinear analysis. Following the insight of Lev Landau, who introduced amplitude equations in analogy to phase transitions (Landau & Lifshitz (1987), §26), weakly nonlinear analyses using a multiple-scale approach leading to an linearity and nonnormality, illustrated in the time domain (a) and ystem (in the eigenspectrum (c), the least stable eigenvalue σ 1 has regime, the energy eventually decays like exp(2σ 1,r t). Nonnormal nsient growth. Nonlinearity may be stabilizing or destabilizing. , for linearly stable system (in the eigenspectrum (c), the least stable eigenvalue σ 1 has a negative growth rate). (a) In the linear regime, the energy eventually decays like exp(2σ 1,r t). Nonnormal systems can experience a very large transient growth. Nonlinearity may be stabilizing or destabilizing. (b) Normal systems subject to external forcing respond preferentially at frequency σ 1,i . Nonnormal systems can respond at different frequencies, with an amplification much larger than predicted by 1/σ r . Nonlinearity may be stabilizing or destabilizing. Figure 1: Cartoon representation of nonlinearity and nonnormality, illustrated in the time domain (a) and frequency domain (b), for a linearly stable system; the least stable eigenvalue σ 1 of the eigenspectrum in (c) has indeed a negative growth rate. (a) In the linear regime, the amplitude of the perturbations eventually decays like exp(σ 1,r t). Nonnormal systems can experience a very large transient growth. Nonlinearity may be stabilising or destabilising. (b) Normal systems subject to external forcing respond preferentially at frequency σ 1,i . Nonnormal systems can respond at different frequencies, with an amplification much larger than predicted by σ 1,r . Nonlinearity may be stabilising or destabilising.
physical implication of nonnormality in the unstable laser cavity is profound, as it results in a substantial increase in the linewidth of the laser beam signal compared to a perfect resonator Petermann (1979). In astrophysics, pseudospectrum analysis was used very recently to study the stability of black holes Jaramillo et al. (2021). In network science, a recent study  reports a systematic analysis of a large set of (directed) empirical networks from a variety of disciplines including "biology, sociology, communication, transport, and many more", and gives evidence that all of them present a strong nonnormality. Nonnormality may shrink the basin of attraction of a linearly strictly stable equilibrium, as strong amplification may trigger nonlinearities (figure 1), and radically change the behaviour of dynamical systems. This is illustrated in , where the nonnormality of the London Tube network results in an hypothetical outbreak of measles epidemic, although the linear stability theory predicts an asymptotic decay of the number of contagions. In hydrodynamics, nonnormality is frequent and inherited from the linearisation of the advective term (U · ∇)U . This term gives a preferential direction to the fluid flow, which breaks the normality of the linear operator. In the context of parallel flows, nonnormality is found for instance in the canonical plane Couette and Poiseuille flows Gustavsson (1991); Butler & Farrell (1992); Reddy & Henningson (1993); Trefethen et al. (1993); Farrell & Ioannou (1993); Schmid & Henningson (2001), pipe flow Schmid & Henningson (1994) and parallel boundary layers Butler & Farrell (1992); Corbett & Bottaro (2000). Nonnormality is also found in nonparallel flows Cossu & Chomaz (1997), for instance spatially developing boundary layers Ehrenstein & Gallaire (2005);Åkervik et al. (2008); ; Alizard et al. (2009); Monokrousos et al. (2010), jets Garnaud et al. (2013a,b) and the flow past a backward-facing step Blackburn et al. (2008); Boujo & Gallaire (2015). Exhaustive reviews of nonnormality in hydrodynamics can be found in Chomaz (2005); Schmid (2007). The crucial role played by nonnormality in the transition to turbulence has become clear over the years Trefethen et al. (1993); Baggett & Trefethen (1997); Schmid (2007). As mentioned in the context of nonnormal networks, if the flow is nonnormal, lowenergy perturbations such as free-stream turbulence or wall roughness can be amplified strongly enough to lead to a regime where nonlinearities come into play, which may lead to turbulence through a sub-critical bifurcation. The toy system presented in Trefethen et al. (1993) is an excellent illustration of this so-called "bypass" scenario.

Amplitude equations without eigenvalues
Through nonnormality, systems with strongly damped spectra can bear strong amplification of specific structures at relatively selective frequencies and/or temporal horizons. To the best of the authors' knowledge, it is currently impossible to construct an amplitude equation for such systems, again because no neutral bifurcation point exists. Furthermore, a systems may well have a weakly damped mode and still exhibit nonnormality, which would jeopardise a classical, single-mode amplitude equation.
Notwithstanding the relevance and usefulness of fully nonlinear solutions Hof et al. (2004);Schneider et al. (2010), as well as the existence of a fully nonlinear nonnormal stability theory able to compute nonlinear optimal initial conditions via Lagrangian optimisation Cherubini et al. (2010Cherubini et al. ( , 2011Pringle & Kerswell (2010); Kerswell (2018), we believe that establishing a rigorous reduced-order model for weak nonlinearities is relevant. In specific regimes, such a model would quantify the respective contribution of each dominant nonlinear interaction, thus bringing insight on the saturation mechanisms of harmonic and transient amplification. It would also predict efficiently if such a saturation actually exists or if, on the contrary, nonlinearities tend to yield even stronger amplification than in the linear regime, thus leading to sub-critical or non-monotonous behaviors; the sketch in figure 1 illustrates possible scenarios for the effects of nonlinearities subsequent to strong nonnormal amplification. Finally, amplitude equations are useful for flow control and optimisation, as shown for instance in Sipp (2012) in the more classical context of a marginally stable flow, displaying little nonnormality, a well isolated eigenvalue and a sufficiently large spectral gap.
The present work proposes to reconcile amplitude equations and nonnormality. Specifically, a method is advanced to derive amplitude equations in the context of (i) harmonic forcing and (ii) transient growth. In case (i), we vary the amplitude of a given harmonic forcing at a prescribed frequency and predict the gain (energy growth) of the asymptotic response ( §2). In case (ii), we vary the amplitude of a given initial condition and predict the gain of the response at a selected time t = t o ( §3). In both cases, we perform an a priori weakly nonlinear prolongation of the gain, at very low numerical cost. The applied Figure 2: Sketch of the flow configurations. (a) Two-dimensional flow over a backwardfacing step, with fully developed parabolic profile of unit maximum centerline velocity at the inlet. (b) Three-dimensional plane Poiseuille flow, confined between two solid walls at y = ±1, and invariant in the x (streamwise) and z (spanwise) directions harmonic forcing and initial condition are allowed to be arbitrarily different from any eigenmode. The method does not rely on the presence of an eigenvalue close to the neutral axis; instead, it applies to any sufficiently nonnormal operator. If such an eigenvalue was nevertheless present on the neutral axis, we recover a classical, modal amplitude equation. The method is illustrated with two flows, the nonparallel flow past a backward-facing step (sketched in figure 2a) and the parallel plane Poiseuille flow (figure 2b). These two nonnormal flows exhibit large gains, both in the context of harmonic forcing ( §2.1-2.2) and transient growth ( §3.1-3.2).
In both contexts, a generic nonlinear dynamical system is considered, where N ( * ) is a nonlinear operator and F is a forcing term. An appropriate and common place to begin the analysis of (1.1) is to linearise it around an unforced equilibrium. The latter is denoted U e and satisfies N (U e ) = 0. Around this equilibrium are considered small-amplitude perturbations in velocity u, forcing f , and initial condition u 0 , where 1. An asymptotic expansion of (1.1) in terms of can thus be performed, transforming the nonlinear equation into a succession of linear ones. The fields u, f and u 0 are recovered at order and linked trough the linear relation where L results from the linearisation of N around U e . For fluid flows governed by the incompressible Navier-Stokes equations, where the pressure field p is such that the velocity field u is divergence-free. Both fields are linked trough a linear Poisson equation. In practice, pressure is included in the state variable, resulting in a singular mass matrix; it is omitted here, for the sake of clarity.

Response to Harmonic Forcing
We first derive an amplitude equation for the weakly nonlinear amplification of timeharmonic forcing f (x, t) =f (x)e iωot + c.c in a linearly strictly stable system. In the long-time regime, only the same-frequency harmonic response u(x, t) =û(x)e iωot + c.c persists. Injecting the expressions of f and u in (1.2) leads toû = (iω o I − L) −1f . = R(iω o )f , where R(z) = (zI − L) −1 is the resolvent operator. In the current context, it maps a harmonic forcing structure onto its asymptotic linear response at the same frequency. A measure of the maximum gain is (2.1) In the following, we choose the L 2 norm (or "energy" norm) induced by the Hermitian inner product û a ,û b = Ωû H aûb dΩ (the superscript H denotes the Hermitian transpose). The operator R(iω o ) † denotes the adjoint of R(iω o ) under this scalar product, such that R(iω o )û a ,û b = û a , R(iω o ) †û b , for anyû a ,û b . Among all frequencies ω o , the one leading to the maximum amplification is noted ω o,m and associated with an optimal gain G(iω o,m ) = 1/ o,m . The singular value decomposition of R(iω o ) provides G(iω o ) = −1 o as the largest singular value, and the associated pair of right singular vector f o and left singular vectorû o . The former represents the optimal forcing, whereas the latter characterises the long-time harmonic response reached, after the transients fade away: leads to Φû o = 0, such that Φ is exactly singular. The norm of the perturbation operator is small since ||P || = 1. The fieldû o constitutes the only non-trivial part of the kernel of Φ, and its associated adjoint mode isf o . Indeed, using that where we used the fact that the inverse of the adjoint is the adjoint of the inverse. We note that Φ can be rewritten as Φ = (iω o I − L n ) where L n . = L + o P , such that (2.3) seems to imply that the state operator L has been perturbed. In this process, the operator L n has acquired an eigenvalue equal to iω o , and therefore has become neutral. However, it has also lost its reality and therefore does not, in general, possess an eigenvalue equal to −iω o . By construction, o is the smallest possible amplitude of the right-hand side of (2.2) for a given iω o , such that o P is the smallest perturbation of L necessary to relocate an eigenvalue of L on iω o . This fact can be formalised with the pseudospectrum theory outlined in Trefethen & Embree (2005). In the complex plane, z ∈ C belongs to the -pseudospectrum Λ (L) if and only if R(z) 1/ . If E is an operator with ||E|| = 1, eigenvalues of L − E can lie anywhere inside Λ (L). Eigenvalues of L and Weak nonlinearity for strong nonnormality singularities of R(z) thus collide with the -pseudospectrum in the limit → 0. As increases, the -pseudospectrum may touch the imaginary axis, such that any z = iω o can be an eigenvalue of L − E if the amplitude of the perturbation is greater than or equal to = R(iω o ) −1 . We recognise as the inverse gain o defined in (2.1), and thus E as P . In particular, if ω o = ω o,m , the associated o,m is referred to as the stability radius of L since the o,m -pseudospectrum is the first to touch the imaginary axis.
As an illustration of the fact that a small-amplitude perturbation can easily "neutralise" a nonnormal operator, we consider the Navier-Stokes operator linearised around the steady flow past a backward-facing step (BFS), sketched in figure 2, at Re = 500. The most amplified frequency ω o = ω o,m ≈ 0.47 is associated with o ≈ 1.3 · 10 −4 1. The spectra of L and L n are shown in figure 3, together with part of the o -pseudospectrum of L. Clearly, the very small perturbation o P locates an eigenvalue exactly onto iω o , despite the strong stability of L. We stress that neither ω o nor o can be deduced only by inspecting the spectrum of L.
Nevertheless, in what follows, it is really the inverse resolvent and not the state operator L that we propose to perturb. Indeed, L is generally a real operator whereas L n is necessarily a complex one, and only one side of the spectrum of L n can generally be made neutral at a time, depending on whether L is perturbed with P or its complex conjugate P * .
The inverse gain o 1 constitutes a natural choice of small parameter. We choose the Navier-Stokes equations for their nonlinear term (U · ∇)U , which yields both a nonnormal linearised operator and a rich diversity of behaviours. The flow is weakly h e iωot + c.c, wheref h is an arbitrary (not necessarily optimal) forcing structure, and φ = O(1) is a real prefactor. Imposing ||f h || = 1, the forcing amplitude is F . = φ √ o 3 . A separation of time scales is invoked for the flow response: its envelope is assumed to vary on a slow time scale The velocity field at each order j is then Fourier-expanded as For m = 2, 3, ..., and C(a, b) .
. Note that the perturbation o P modifying R(iω 0 ) −1 into Φ at leading order is compensated for at third order. Terms are then collected at each order in √ o , leading to a cascade of linear problems, detailed hereafter.
At order √ o , we collect (imω o I−L)u 1,m = 0 for m = 0, 2, 3 . . ., and Φu 1,1 = 0. Since L is strictly stable, the unforced equation for m = 1 can only lead to u 1,m = 0. Conversely, the kernel of Φ contains the optimal responseû o , therefore u 1,1 (T ) = A(T )û o , where A(T ) ∈ C is a slowly-varying scalar amplitude verifying ∂ t A = 0. Finally, the general solution at order At order o , we obtain the solution u 2 = |A| 2 u 2,0 + A 2 e 2iωotû 2,2 + c.c , where (2.8) The homogeneous solution of the system Φu 2,1 = 0 is arbitrarily proportional toû o , and written A 2 (T )û o . It can be ignored (u 2,1 = 0) without loss of generality. As shown in Fujimura (1991), it could also be kept, provided it is included in the definition of the amplitude, which would then become A + o A 2 . At order √ o 3 , we assemble two equations yielding the Fourier components of the solution oscillating at ω o , (recalling Pû o =f o ), and at 3ω o , (3iω o I − L)u 3,3 = 2A 3 C(û o ,û 2,2 ). The operator Φ being singular, the only way for u 3,1 to be non-diverging, and thus for the asymptotic Weak nonlinearity for strong nonnormality 9 expansion to make sense, is that the right-hand side of (2.9) has a null scalar product with the kernel of Φ † , i.e. is orthogonal to the adjoint modef o associated withû o . This is known as the "Fredholm alternative". As a result, the amplitude A(T ) satisfies The coefficient γ is the projection of the applied forcing on the optimal forcing. The coefficient µ embeds the interaction betweenû o and the static perturbation u 2,0 , i.e. it corrects the gain according to the fact thatû o extracts energy from the time-averaged mean flow rather than from the original base flow. We show in Appendix A that, in the regime of small variations around the linear gain, the amplitude equation reduces to the standard sensitivity of the gain Brandt et al. (2011) to a base flow modification induced by u 2,0 . In contrast, the coefficient ν embeds the interaction betweenû * o and the second harmonicû 2,2 . Introducing the rescaled quantities a . (2.12) The gain associated with the linearised version of (2.12) is G = |γ| / o , as expected for the linear prediction. We recover G = 1/ o when the optimal forcing is applied (γ = 1). We also note that this expression predicts G = 0 when γ = 0, which merely indicates that the linear response is orthogonal toû o , without stating anything on the gains associated with sub-optimal forcings except that they should be at most O( −1/2 o ), assuming a sufficiently large "spectral" gap in the singular-value decomposition of the resolvent operator. For the rest of the paper, we set γ = 1. Expressing a in terms of an amplitude |a| ∈ R + and a phase ρ ∈ R such that a(t) = |a(t)| e iρ(t) , the time-independent equilibrium solutions, or fixed points, of equation (2.12), named (|a e | , ρ e ), solve: Squaring and adding the real and imaginary parts of (2.13) leads to a third-order polynomial for the equilibrium amplitude of (2.12): and where Y = |a e | 2 > 0. Let p(Y ) = DY 3 + 2BY 2 + Y be the left-hand side of (2.14). We further distinguish two cases: (i) if B 0 , p(Y ) is increasing monotonously with Y and can only cross the constant line (F/ o ) 2 once. We have in addition p(Y ) > Y , thus the gain smaller than the linear prediction and monotonously decaying while increasing F . Conversely, if (ii) B < 0, we have p(Y ) < Y in the interval 0 < Y < −2B/D, and the gain should then be greater than the linear one in the corresponding range of forcing 0 < (F/ o ) 2 < −2B/D. Furthermore, p(Y ) may vary non-monotonously over this interval and cross the constant line (F/ o ) 2 three times (leading to three solutions for Y ); namely, p(Y ) may be decreasing on a certain interval of Y while dominated by the negative term ∝ Y 2 , bridging two other intervals where p(Y ) is increasing due to the respective positive terms ∝ Y and ∝ Y 3 . A necessary and sufficient condition for such a case to occur is that the equation dP/dY = 3DY 2 + 4BY + 1 = 0 possesses two real and and positive solutions. This is guaranteed if and only if the determinant ∆ . = 16B 2 − 12D is strictly positive. Finally, for −2B/D Y , p(Y ) must be monotonously increasing again with p(Y ) Y , resulting in a gain smaller than the linear one and monotonously decreasing while increasing F .
The stability of the limit cycle associated with the equilibrium solution(s) (|a e | , ρ e ) can be established from the amplitude equation (2.12). Although in a different context, this was demonstrated for instance in Tuckerman & Barkley (1990), where the bifurcation diagram of the Eckhaus instability is determined directly from the Ginzburg-Landau equation for the envelope of the critical eigenfunction. Equation (2.12) can be expressed as a two-by-two amplitude/phase nonlinear dynamical system: (2.16) Perturbing this system around the equilibrium solution (|a e | , ρ e ) + (|a| (t), ρ (t)) and neglecting nonlinear terms leads to the following equation for the perturbation If at least one of the two eigenvalues of J has a positive real part, the associated equilibrium is linearly unstable. Note that equations (2.15) and (2.16) for the amplitude and the phase of the oscillating linear response, are similar to those that would be obtained for a classical Duffing-Van der Pol oscillator with appropriate parameters and harmonically forced around its natural frequency. If the latter is set to one, η r o and η i o are respectively proportional to the damping ratio and the detuning parameter. The coefficient (µ i +ν i ) is proportional to the the cubic stiffness parameter (Duffing nonlinearity ∝ x 3 ), and (µ r + ν r ) to the nonlinear damping parameter (Van der Pol nonlinearity ∝ẋx 2 ).
For the sake of completeness, Appendix C shows how to compute higher-order corrections of (2.12). It is worth mentioning, in particular, that the action of Φ need not be computed explicitly and can be replaced by the action of (iω o I − L) for all practical purposes.
2.1. Application case: the flow past a backward-facing step Equation (2.12) is the first main result of this study and will be further referred to as the Weakly Nonlinear Nonnormal harmonic (WNNh) model. We discuss its performance when the stationary flow past a BFS sketched in figure 2 is forced harmonically with the optimal structuref o . At Re = 500 and the optimal forcing frequency,f o is shown in figure 11a together with is associated responseû o in figure 11b (see Appendix B for details about the geometry and the numerical method). As shown in Blackburn et al. (2008); Boujo & Gallaire (2015), the BFS flow constitutes a striking illustration of streamwise 200 73.9 −1 3.66 + i · 0.0163 5.13 + i · 1.32 0.137 − i · 1.13 5.27 500 7456.6 −1 117.1 + i · 0.653 8.23 + i · 2.60 0.364 + i · 0.396 8.59 700 148080 −1 1626.7 + i · 8.65 9.06 + i · 4.38 −0.729 + i · 1.39 8.33 Table 1: WNNh coefficients for the backward-facing step flow, when the optimal forcing structure (γ = 1) is applied at the optimal frequency ω o /(2π) = ω o,m /(2π) = 0.075.
nonnormality. As seen in figure 11a, the optimal forcing structure is located upstream and triggers a spatially growing response along the shear layer adjoining the recirculation region, as the result of the convectively unstable nature of the shear layer. We first set the Reynolds number Re between 200 and 700, and the frequency ω o = 2π × 0.075 close to the most linearly amplified frequency ω o,m , which varies only slightly with Re. The linear gain grows exponentially with Re Boujo & Gallaire (2015), as seen in table 1. Since η scales like O( −1/2 o ), the term in dA/dT in (2.9) is asymptotically consistent only close to equilibrium points where dA/dT = 0, which is the regime of primary interest in the context of harmonic forcing. In accordance, the temporal derivative dA/dT is kept in (2.10) to assess the stability of such equilibria, determined by the analysis of the Jacobian matrix (2.17).
Predictions from the WNNh model are compared to fully nonlinear gains extracted from direct numerical simulations (DNS) in figure 5a. The DNS gains are the ratio between the temporal rms of the kinetic energy of the fluctuations at ω o (extracted through a Fourier transform) and the rms of the kinetic energy of the forcing (for instance, the forcing Ff o e iωot +c.c with ||f o || = 1 corresponds to an effective forcing rms amplitude of √ 2F ). Since the coefficient B defined in (2.14) is strictly positive for all Re, the WNNh model predicts nonlinearities to saturate the energy of the response, and thus the gain to decrease monotonously with the forcing amplitude. This is confirmed by the comparison with DNS, displaying an excellent overall agreement. As shown in the inset (in logarithmic scale), the nonlinear gain transitions from a constant value in the linear regime to a −2/3 power-law decay when nonlinearities prevail, as predicted from (2.12). This transition is delayed when the Reynolds number (and therefore the linear gain)  decreases, and compares well with DNS data. The main plot (in linear scale) confirms the agreement with the DNS, and the improvement over the linear model. Re-scaled WNNh curves appear similar for Re = 200 and Re = 700, and a slight overestimate is observed as the forcing amplitude approaches 0 . Indeed, F ∼ 0 implies φ ∼ 1/ √ o , which jeopardises the asymptotic hierarchy. Nonetheless, the error remains small for this flow in the considered range of forcing amplitudes. Further physical insight is gained from the WNNh coefficients gathered in table 1. The nonlinear coefficients remain of order one, which confirms the validity of the chosen scalings. The real part of µ being larger than that of ν, the present analysis rationalises a priori the predominance of the mean flow distortion over the second harmonic in the saturation mechanism reported a posteriori in Mantic-Lugo & Gallaire (2016b).
Next, we select Re = 500 and report in figure 5b harmonic gains as a function of the frequency, for increasing forcing amplitudes. At each frequency, the corresponding optimal forcing structuref o is applied. The comparison between DNS and WNNh is conclusive over the whole range of frequencies. The saturating character of nonlinearities is well captured. Such a good agreement may appear surprising in the low-frequency regime, for instance at ω o /(2π) = 0.04 where the second harmonics at frequency 2ω 0 could in principle be amplified approximately four times more than the fundamental. It happens, however, that the associated forcing structure −C(û o ,û o ) is located much farther downstream than the optimal forcing at 2ω o , with a weak overlap region which results in a poor projection. Therefore, the second-order contribution does not reach amplitudes of concern in this flow, as a consequence of its streamwise nonnormality.

Application case: Orr mechanism in the plane Poiseuille flow
The weakly nonlinear evolution of the harmonic gain is now sought for the plane Poiseuille flow sketched in figure 2, a typical flow with component-wise nonnormality Trefethen et al. (1993); Schmid (2007). Periodicity is imposed in the streamwise and spanwise directions with wavenumbers k x and k z , respectively. The set of parameters (Re, k x , k z ) = (3000, 1.2, 0) is selected. According to the classical work of Orszag (1971), the base flow at this Re number is linearly stable since instability first occurs at Re cr ≈ 5772 and k x,cr ≈ 1.02. In both the linear and nonlinear computations, the spanwise invariance k z = 0 is systematically maintained. While the base flow U (y) has only one velocity component and depends only on one coordinate, the perturbations are here two-dimensional (i.e, u = (u x (x, y), u y (x, y))). The computations are performed in the streamwise-periodic box (x, y) ∈ [0, 2π/k x ]×[−1, 1] ≡ Ω. All the scalar products are taken upon integration inside this periodic box, in particular for the normalisation û o ,û o = f o ,f o = 1, and latter for the evaluation of the weakly nonlinear coefficients. The linear optimal gain (2.1) is computed in the frequency interval 0 ω o 0.8 (figure 6a), together with the associated optimal forcing and responses structures. Results are validated with the 1D results of Schmid & Henningson (2001) based on a Fourier expansion of wavenumbers k x = 1.2 and k x = 0 in the streamwise direction. Eigenspectra are also reported in figure 6(b). The SVD algorithm applied to the periodic box automatically selects the most amplified wavenumber among all spatial harmonics nk x with n = 0, 1, 2, ... Below ω o ≈ 0.12, the harmonic 0 · k x = 0 is dominant due to the concentration of weakly damped eigenvalues along the imaginary axis. The gain G(ω o = 0) = 1216 is equal to the inverse of the smallest damping rate among all these spatially invariant modes. The large value of the gain associated with those modes is understood considering that the small pressure gradient (2/Re, 0) T = (2/3000, 0) T is sufficient to induce the Poiseuille base flow (equal to unity in the centerline). Above ω o ≈ 0.12, the fundamental wavenumber 1·k x = 1.2 prevails. The corresponding harmonic gain presents a local and selective maximum for ω o = 0.38, certainly linked to the presence of the weakly damped eigenvalue σ 1 = −0.0103 + 0.380i. Nevertheless, G(ω o = 0.38) = 416 is significantly bigger than 1/0.0103 ≈ 97. This is a direct consequence of the nonnormality of the plane Poiseuille flow. Unlike the backward-facing step flow, the nonnormality at play here is not due to the presence of a convectively unstable region but to the Orr mechanism suggested for the first time in Orr (1907). Namely, an initial condition or forcing field constituted of spanwise vortices tilted towards the upstream direction ( fig. 7a), tilts downstream under the action of the mean shear ( fig. 7b), which leads to a significant gain in the kinetic energy of the perturbation. The coefficient B is shown in figure 8a, and the associated WNNh prolongation of the harmonic gain in figure 8b. The coefficient B is negative in the interval 0.378 ω o 0.486, and B and A are such that three equilibrium amplitudes |a e | exist for some values of F in the sub-interval 0.389 ω o 0.428. Among them, none or only one is found to be stable. Consequently, as the forcing amplitude is increased, the harmonic gain curve leans toward the higher frequencies in figure 8b; in the meantime, a frequency interval where no stable solution is predicted appears and grows larger. Note that in the absence of a stable equilibrium, it is natural to consider completing (2.12) up to O( √ o 5 ). It is shown in Appendix C, however, that such an approach is problematic in the present case, because the non-oscillating forcing terms appearing   1.2, 0), and when the optimal forcing structure (γ = 1) is applied. The corresponding WNNh prolongation of the linear gain as a function of the forcing amplitude is shown in figure 9, together with DNS results. For comparison, the prediction of a "classical" (modal) amplitude equation constructed around the weakly damped eigenvalue σ 1 and its associated direct and adjoint modes is also added. Its derivation is briefly recalled in Appendix D.
For ω o = 0.3810 (figure 9a), the WNNh gain initially increases with F due to the negativity of B. As visible in table 2, this is mostly due to the contribution of [µ/( o η)] which is ten times larger than that of [ν/( o η)]. Thus, at this frequency, the principal factor for the initial increase of the WNNh gain is the Reynolds stress of the response aû o . The latter creates a mean flow that amplifies the linear forcingf o more than the base flow does. This may be interpreted considering the displacement of the eigenvalue σ 1 . Letq 1 (resp.â 1 ) denote the eigenmode (resp. adjoint mode) associated with the eigenvalue σ 1 . The sensibility of the latter to the base flow deformation δU b due to the Reynolds stress of aû o writes: where δU b = |a| 2 u 2,0 . For ω o = 0.3810, we obtain δσ 1 = |a| 2 (1.2+i·3.9). Since [δσ 1 ] > 0, the eigenvalue is moving towards the unstable part of the complex plane under the action of the Reynolds stress. This is in accordance with the fact that the plane Poiseuille flow is subcritical, and may explain the initial increase in the gain with F . Meanwhile, [δσ 1 ] > 0 and σ 1 is shifting toward higher frequencies.
Thus ω o ceases to be the least damped frequency, which could shed light on the fact that increasing F further leads to a monotonous decay in the WNNh gain at ω o . Because of the flow nonnormality, however, this explanation based solely on the location of σ 1 remains qualitative.
The overall agreement with the DNS results is excellent. Nevertheless, the WNNh model slightly underestimates the threshold in F above which a stable equilibrium does not exist any more. It stands at F/ o = 0.087 against F/ o = 0.11 for the DNS. This loss of a proper harmonic response may be symptomatic of the fact that σ 1 eventually crosses the neutral line and becomes unstable. Indeed for F/ o = 0.11 (blue circle in the grey zone in figure 9a), the FFT of the flow in its stationary regime presents two dominant neighbouring frequencies: the forcing one at ω = ω o and a second "natural" one at ω ≈ 0.404. As these two frequencies are very close, a beating behaviour is visible in the inset of figure 9a at a frequency consistent with ∆ω = 0.023.
The classical modal amplitude equation leads to a prediction that is only qualitative. Even for F = 0 the linear harmonic gain | f o ,â 1 /( q 1 ,â 1 σ 1,r )| (see Appendix D for its derivation) is overestimated, as it is deduced from the modal quantities linked to σ 1 only. As mentioned earlier, in nonnormal flows a high number of eigenmode is generally necessary to describe its harmonic response, even in the presence of a weakly damped eigenvalue. Thus, relying on a single mode constitutes a poor description of the response to forcing.
We now consider ω o = 0.4025, and the associated results in figure 9b. The WNNh model yields multiple equilibrium solutions in the range 0 < F/ o < 0.0264. Only the one represented by a thick continuous line is stable, and corresponds to a monotonous growth of the gain with F . The DNS results validate the existence of this solution. The two other solutions, depicted by the dash-dotted and dashed lines, are unstable in one eigendirection and two eigendirections, respectively. Above F/ o = 0.0264 the WNNh models predicts the loss of the stable equilibrium solution, which is accurately confirmed by the DNS whose threshold is located around F/ o = 0.0286. Slightly above, the signal of u y (0, 0) in the inset suggests again the presence of a "natural" frequency due to the subcritical destabilisation of σ 1 . Indeed, u y (0, 0) alternates between an algebraic growth typical of a true resonance (both natural and forcing frequencies collapse), and a beatinglike behaviour whose period is very long (the natural frequency drifts slightly from the forcing one).
Across this threshold, the evolution of the average kinetic energy of the response appears discontinuous. This loss of a stable equilibrium is to be distinguished with its destabilisation encountered for ω = 0.3810. Overall, the difference of behaviours between figures 9a and 9b may be explained by the difference of proximity between ω o and [σ 1 ] of the mean flow. As the forcing is progressively increased above F/ o = 0.0286, the flow response quickly becomes chaotic, and then turbulent.
It should be mentioned that, in some situations, the amplitude equation (2.12) may be in default. First, as just mentioned, when the optimal linear harmonic gain at frequency 2ω o is ∼ 1/ √ o or larger and projects well onto the optimal forcing, the asymptotic hierarchy is threatened asû 2,2 may be substantial enough to reach order √ o or above. It is thus important to assess that the norm ofû 2,2 remains of order one. A second delicate situation arises, for the same reason, when a sub-optimal gain at the frequency ω o is ∼ 1/ o . In both cases, the model could be extended by including in the kernel of Φ the optimal response at frequency 2ω o , or the sub-optimal response at frequency ω o , respectively.

Transient Growth
Next, we derive an amplitude equation for the weakly nonlinear transient growth in an unforced (f = 0) system, without restriction on its linear stability. The solution to the linearised equation (1.2) is u(t) = e Lt u(0), where e Lt is the operator exponential of Lt. In an unforced context, the propagator e Lt maps an initial structure at time t = 0 onto its evolution at t 0. The largest linear amplification at t o > 0 (subscript o for "optimal") is (3.1) The singular value decomposition of the propagator e Lto provides the transient gain G(t o ) as the largest singular value of e Lto , as well as the left and right singular pair v o and u o , respectively, where ||v o || = ||u o || = 1. The field u o is the optimal initial structure for the propagation time t = t o , and v o is its normalised evolution at t o . The corresponding amplification is 1/ o , as defined in (3.1). Smaller singular values are sub-optimal gains, associated with orthogonal sub-optimal initial conditions. Their orthogonality is ensured by the fact that singular vectors of the operator e Lto also are the eigenvectors of the symmetric operator (e Lto ) † e Lto , the singular values of the former being the square root of the eigenvalues of the latter. Of all the t o , the time leading to the largest optimal gain will be highlighted with the subscript m (for "maximum") such that max to>0 G(t o ) = G(t o,m ). By construction, the linear gain is independent of the amplitude of the initial condition u(0). As this amplitude increases, however, nonlinearities may come into play and the nonlinear gain may depart from the linear gain G. Similar to the previous section on harmonic gain, we propose a method for capturing the effect of weak nonlinearities on the transient gain.
Due to the assumed nonnormality of L, the inverse gain is small, o 1. While the previous section focused on the inverse resolvent, it is now the inverse propagator e −Lto that appears close to singular. The first equality of (3.2) can be rewritten as is singular since v o = 0 belongs to its kernel. Mirroring our previous reasoning for the WNNh model, we now wish to construct a perturbed inverse propagator whose kernel is the linear trajectory seeded by the optimal initial condition u o and normalised in t = t o such that l(t o ) = v o . One conceptual difficulty lies in that the linear response is not a fixed vector field, but a time-dependent trajectory; therefore, the perturbed inverse propagator too should depend on time. We propose to perturb the inverse propagator for all t 0 as and where the Heaviside distribution H(t) satisfies H(0) = 0 and H(t > 0) = 1. As the time t → t o , the perturbation operator P → u o v o , * such that ||P || → 1 and the expansion (3.4) is certainly justified. The non-trivial kernel of Φ(t) is l(t) for all t > 0; the kernel reduces to 0 at t = 0 since Φ(0) = I. We show in addition that, for t > 0, the non-trivial kernel of the adjoint operator Φ(t) † is b(t) . = e Lt † l(t). (3.5) Indeed, using that P † = l(t) u o , * / l(t), l(t) for t > 0, we have As an illustration of the singularisation of e −Lto , parts of the spectra of e −Lto and Φ(t o ) are shown in figure 10 for the plane Poiseuille flow sketched in figure 2. The red dot at the origin is the null singular eigenvalue of Φ(t o ) associated with l(t o ). Since ||P (t o )|| = 1, this singular eigenvalue lies on the o -pseudospectrum of e −Lto , meaning that a perturbation of amplitude o is sufficient to make the inverse propagator singular.
Recalling that L is assumed strongly nonnormal, we choose o 1 as expansion parameter, introduce the slow time scale T = o t, and propose the multiple-scale expansion The square root scaling of the previous section is not made here, as resonance at second order cannot be excluded a priori. The flow is initialised with U (0) = α 2 o u o , where α = O(1) is a prefactor. After injecting this expansion in the unforced Navier-Stokes equations, we obtain subject to u 2 (0) = αu o , and u i (0) = 0 for i = 2. In its primary quality of inverse propagator, the following property holds for e −Lt : ∂ t (e −Lt ) = −e −Lt L, where the commutation of e −Lt and L has not been used. Thanks to this relation, we write (∂ t − L)u i = e Lt ∂ t (e −Lt u i ). As a result, L disappears from the asymptotic expansion but e −Lt appears. The latter is perturbed according to (  varying operator L(t) and the associated propagator Ψ (t). This can be shown easily by taking the time derivative of Ψ (t) −1 u(t) = u(0). Terms of (3.8) are then collected at each order in o , leading to a succession of linear problems, detailed hereafter. At order o , we collect ∂ t (Φu 1 ) = 0, subject to u 1 (0) = 0. We obtain Φu 1 = Φ(0)u 1 (0) = 0, therefore u 1 (t, T ) is proportional to the kernel of Φ(t) for all t 0. We choose the non-trivial solution where the initial condition u 1 (0) = 0 is enforced by H(t), while the slowly-varying scalar amplitude A(T ) is continuous in T and modulates the linear trajectory. This choice is motivated by the observation that, since A must be constant in time in the linear regime, we expect it to be weakly time-dependent in the weakly nonlinear regime. We stress that A(T ) does not depend explicitly on t, such that ∂ t A = 0. Note that the choice u 1 (t) = A(t)H(t)l(t) would also have been possible, and the assumption of the amplitude depending on a slow time scale is made solely to simplify the ensuing calculations. At order 2 o , we collect ∂ t (Φu 2 ) + A 2 He −Lt C(l, l) + H dA dT e −Lt l + Ad t (HP l) = 0, (3.10) subject to u 2 (0) = αu o . We used the property H(t) 2 = H(t), which will henceforth be understood. The particular solution of (3.10) yields subject to the initial conditions u and where we used that the general solution of d t Invoking again the Fredholm alternative, (3.13) admits a non-diverging particular solution if and only if its right-hand side is orthogonal to b(t) for all t > 0. This leads to: (3.14) Dividing (3.14) by u o , b(t) leads to (3.16) is equivalent to solving ∂ t E(t, T ) = 0 for t > 0 subject to E(t → 0, T ) = 0. Thereby, the partial derivative of (3.15) with respect to the short time scale t is taken, leading to where we have used that ∂ t A = 0 by construction since A = A(T ) does not explicitly depend on t. Furthermore, the relation (3.17) is subject to E(t → 0, T ) = lim t→0 (α−A) = 0 where we have used thatũ 2 (t → 0) =ũ 2 (0) = 0. To be meaningful, equation (3.17) and its initial condition must be re-written solely in terms of t, which is done by evaluating T along T = o t. The total derivative of A, denoted d t A, is now needed, as it takes into account the implicit dependence of A on t. By definition, such that the final amplitude equation reads as lim t→0 (α − A( o t)) = 0 implies A(t → 0) = α, and the amplitude A is extended by continuity in t = 0 so as to eventually impose A(0) = α. Note that the evaluation in T = o t and the passage to the total derivative would lead to indeterminacy in its solution if performed directly in (3.15), since that equation is not subject to any initial condition. Indeed, at linear level for instance, it would yield d t A = (α − A)/t, which admits the family of solutions A(t) = α + Ct, with C an undetermined constant. We stress that the inverse propagator is not needed to solve the amplitude equation (3.18). Just like the original problem considered in this section, (3.18) is unforced and has a non-zero initial condition. In the linear regime, A = α for all times, and the linear gain In the following, we call Equation (3.18) the Weakly Nonlinear Nonnormal transient (WNNt) model. It can be corrected with higher-order terms, which requires solving the linear singular system (3.13), as detailed in Appendix E. We show in particular that singular higher-order solutions are orthogonal to the first-order order solution l(t).  happens to be insufficient to capture the nonlinear saturation of the transient gain for this particular flow, in particular because of the weak value of the coefficient µ 2 (t). Indeed, l(t o ) = v o appears to be dominated by a specific spatial wavenumber (see figure 11b), thus the fieldũ 2 , being generated by the nonlinear interaction of l(t) with itself, it is dominated by spatial harmonics and its projection on l(t) is close to zero. For this flow the WNNt model therefore needs to be extended to order 3 o (see Appendix E), yielding dA dt (3.20) and Equation ( (3.21) In this manner, the weakly nonlinear transient gain becomes G(t o ) = a(t o )/U 0 . Note that in (3.21) the amplitude a(t) does not depend on U 0 nor on o independently, but on their ratio U 0 / o . Thus, as expected, increasingly nonlinear regimes are found when the amplitude of the initial condition increases with respect to the linear gain. Predictions from equation (3.21) are shown in figure 12 together with the linear and fully nonlinear DNS gains evaluated in two ways: using either the total perturbation around the base flow, for G tot with or using that perturbation projected on l(t), for G l . Namely, we recall the asymptotic expansion of the solution U = U e + al + 2 o u 2 + O( 3 o ) for t > 0, with l, u i = 0 for all i 2 as a consequence of the Fredholm alternative (see Appendix E). This gives: Thus, from the knowledge of U −U e computed by DNS, the gain that should be compared with the weakly nonlinear prediction a||l(t)||/U 0 is figure 12a, the WNNt model extended to O( 3 o ) appears to capture the weakly evolution of the transient gain with precision, in particular G l . When G l and G tot depart from each other, the higher-order fieldsũ 2 ,ũ 3 , . . . are expected to have a significant amplitude, and thus the WNNt prediction deteriorates since it is based on an asymptotic hierarchy. For this specific flow, however, the error remains small and the prediction satisfactory even in the fully nonlinear regime. In figure 12b, the gain history of G l for all times 0 t t o is successfully compared to a(t)||l(t)||/U 0 . The coefficient µ 3 (t) is much larger than µ 2 (t) (inset), and is largely dominated by the part ofũ 3 generated by the forcing term C(l,ũ 2 ). Since µ 3 (t) is monotonously decreasing toward µ 3 (t o ) = −7.77, larger times are subject to a stronger saturation. This leads to a decrease of the time for which the specific initial condition u o leads to a maximum transient gain, consistently with the DNS results.

Application case: Lift-up in the plane Poiseuille flow
The WNNt is now applied to the plane Poiseuille flow. The set of parameters (Re, k x , k z , t o ) = (3000, 0, 2, t o,m = 230) is selected. In both the linear and nonlinear computations, the wavenumber k x = 0 is maintained such that the fields are constant in x, and only the dependence in y and z is computed. Contrarily to the application case §2.2, perturbations can now be fully three dimensional (i.e. u = (u x (y, z), u y (y, z), u z (y, z)). The computations are performed in the spanwiseperiodic box (y, z) ∈ [−1, 1] × [−π/k z , π/k z ] ≡ Ω. All the scalar products are taken upon integration inside this periodic box, in particular for the normalisation u o , u o = v o , v o = 1, and for the evaluation of the weakly nonlinear coefficients. The linear optimal gain is validated with the result of Schmid & Henningson (2001); the associated optimal initial condition and its evolution at t = t o are shown in figures 13a and 13b, respectively. The optimal initial condition consists of vortices aligned in the streamwise direction; as these streamwise vortices are superimposed on the parabolic base flow, they bring low-velocity fluid from the wall towards the channel centre and high-velocity fluid from the centre of the channel towards the walls, thus generating alternated streamwise streaks. Due to the spanwise periodicity of the optimal initial condition u o , all the solutions at even orders 2n o (n = 1, 2, 3, ...) only yield even spatial harmonics in k z and are orthogonal to l(t), such that the coefficient µ 2 (t) defined in  (3.16) is null at all times. Therefore, (3.21) reduces to andũ 3 solves the simplified equation The analytical solution of (3.22) writes We show in Appendix F that, at first order in the gain variation, (3.24) reduces to the sensitivity of the transient gain to the base flow modification (U 0 / o ) 2ũ 2 (t). Predictions from equation (3.24) are shown in figure 14 together with the linear and fully nonlinear DNS gains. The WNNt model predicts G l accurately in the weakly nonlinear regime for t = t o,m , which supports our approach (figure 14a). In the strongly nonlinear regime, beyond U 0 / o ≈ 0.4, the model overestimates G l . Again, this can be interpreted by noting that G tot is twice as large as G l , i.e. more energy is contained in the higher-order terms generated by the linear response than in the linear response itself. Therefore, this is not the amplitude equation (3.24) that breaks down, but the very idea of an asymptotic expansion. Whether higher-order terms remain smaller than the fundamental is certainly flow-dependent, and the WNNt model is expected to be even more accurate when this is the case, as shown in §3.1 for the flow past a backward-facing step, which generated rather weak higher-order fields. Figure 14b compares for t t o the history of the approximated gain a(t)||l(t)||/U 0 with that of the DNS gain G l , and shows a convincing overall agreement. The coefficient µ 3 (t) is negative and decays monotonously with time until µ 3 (t o ) = −3.30 (inset), enhancing the saturation. This results in a reduction of the approximated optimal time with forcing amplitude and associated with the initial condition u o , as also observed in the DNS.

Conclusions
In summary, we have derived two weakly nonlinear amplitude equations for nonnormal systems, describing the asymptotic response to harmonic forcing and the transient response to initial condition. In both cases, the presence of a neutral or weakly damped mode is unnecessary. Both approaches are based on the same observation: in nonnormal systems, the resolvent and propagator operators may lead to a notable amplification, so their inverses may in contrast lead to a notable mitigation of the response and are close to singular. A small perturbation is then sufficient to fill their kernel with the response and render them singular. This can be encompassed in a multiple-scale expansion, closed by means of classical compatibility conditions.
The resulting amplitude equations have been validated with fully nonlinear simulations, both in parallel and non-parallel two-dimensional flows. In all cases, they predict accurately the supercritical or subcritical nonlinear evolution of the response, and bring insight on the weakly nonlinear mechanisms that modify the gains as the amplitude of the harmonic forcing or the initial condition varies. In particular, the efficiency of the WNNh model to capture a subcritical behaviour may prove useful in the search of optimal paths to chaos or turbulence. Indeed, equation (2.14) could be included as an additional constraint in a Lagrangian optimisation problem, whose stationary point would constitute a weakly-nonlinear optimal. Such an approach could complement fully nonlinear optimisations (Pringle & Kerswell (2010)), by providing physical understanding at a numerical cost close to the linear one.
It should be noted that the proposed method is not restricted to the Navier-Stokes equations, but apply to all nonlinear systems whose linearised operator exhibits a strong nonnormality (see Trefethen & Embree (2005) §55-60 for a comprehensive discussion, as well as the situations discussed in the introduction). For instance in ecological models describing the temporal evolution of a population, such as the canonical Lotka-Volterra predator-prey equations, the so-called resilience of a community (spectral abscissa of the Jacobian of the system) is known to be sometimes a misleading or incomplete measure Neubert & Caswell (1997). The conjunction of nonnormality and nonlinearity is then key to predict the extinction of a population.
The amplitude equations proposed in the present study are expected to be relevant for three-dimensional flows. For the transient growth aspect, the assumption of a timeindependent state operator and the associated operator exponential formalism are unnecessary, and we believe that the model can be extended to time-varying base flows. Another interesting extension of the models, apart from higher-order corrections, is the inclusion of multiple forcing structures or trajectories, originating for instance from additional singular vectors in the asymptotic expansions. The nonlinear interaction of multiple harmonic forcings or initial conditions is particularly relevant when distinct structures lead to comparable gains, for instance optimal and sub-optimal initial conditions Butler & Farrell (1992); Blackburn et al. (2008), or perturbations of different spatial wavenumbers like in jet flows forced with different azimuthal wavenumbers Garnaud et al. (2013b). The ensuing system of coupled amplitude equations may bear rich dynamics, such as hysteresis and chaos. Our current efforts also involve deriving an amplitude equation for the response to stochastic forcing, as investigated in Farrell & Ioannou (1993) and Mantic-Lugo & Gallaire (2016a) with linear and self-consistent models, respectively. We are interested in the squared gain variation δG 2 o (where this notation does not designate the square of the gain variation) induced by a small perturbation δL of the state operator. The latter results in the following perturbation δR of the resolvent: The gain variation is therefore On the other hand, the WNNh model predicts: where Y = |ā| 2 . We identify the weakly nonlinear harmonic gain as G 2 = Y /F 2 , and multiply (A 2) by 2 o /Y : Being interested in small variations around G 2 o (that correspond to the linear limit Y = |ā| 2 → 0), we write We recognise at leading order equation Thus, in the small gain variation limit, the WNNh model both contains the sensitivity formula of the harmonic gain to the base flow static perturbation |a| 2 u 2,0 , and embeds the effect of the second harmonicû 2,2 as well.
Appendix B. Applying the WNN models to the Navier-Stokes equations.
The incompressible Navier-Stokes equations write after linearising around the equilibrium velocity field U e Several subtleties arise from the peculiarity of the pressure variable, that ensures the instantaneous satisfaction of the incompressibility condition: (i) the absence of timederivative of the pressure results in a singular mass matrix, (ii) forcing terms remain restricted to the momentum equations as we choose to have no source/sink of mass and (iii) the pressure is not included in the energy norm of the response. This complicates slightly the practical computation of the gain. For the harmonic response model, the resolvent operator is generalised as R(iω o ) = (iω o B − L) −1 , and the gain is measured according to where we used the following scalar products q a ,q b B = Ω û * a,xûb,x +û * a,yûb,y +û * a,zûb,z dΩ, and q a ,q b = Ω û * a,xûb,x +û * a,yûb,y +û * a,zûb,z +p * apb dΩ.

28
Y.-M. Ducimetière, E. Boujo and F. Gallaire The B-scalar product excludes pressure, such that the pseudonorm q,q B = ||q|| 2 B is the total kinetic energy of the response. The scalar product at the denominator includes pressure, although this will not change the norm ofd, namely d ,d = ||d|| 2 , since we have no source/sink of mass. The weakly nonlinear coefficients must be considered under these scalar products.
The linear and nonlinear Navier-Stokes equations are solved for (u x ,u y ,p) by means of the Finite Element Method with Taylor-Hood (P2, P2, P1) elements, respectively, after implementation of the weak form in the software FreeFem++. The steady solutions of the Navier-Stokes equations are solved using the iterative Newton-Raphson method, and the linear operators are built thanks to a sparse solver implemented in FreeFem++. The singular value decomposition is performed in Matlab following directly Garnaud et al. (2013b). Finally, DNS are performed by applying a time scheme based on the characteristic-Galerkin method as described in Benitez & Bermudez (2011). For the two-dimensional flow past a BFS presented in §2.1, we refer to Mantic-Lugo & Gallaire (2016b) for the validation of the codes with existing literature and the mesh convergence, since the same codes have been used. The length of the outlet channel is chosen as L out = 50 for Re 500 (Mantic-Lugo & Gallaire 2016b), L out = 65 for Re = 600, and L out = 80 for Re = 700. This ensures the convergence of the linear gain and weakly nonlinear coefficients. For the plane Poiseuille studied in §2.2, the validation is proposed in the main text.
For the transient growth model, the linearised problem writes and the gain is measured as where the pressure component of the initial condition can be chosen as p(0) = 0. The orthogonality properties holds under the B-scalar product, and the weakly nonlinear coefficient µ 2 (t) writes where q l (t) = [l(t), p l (t)] T , and whereq 2 = [ũ 2 (t),p 2 (t)] T is solution of Again, pressure does not influence the weakly nonlinear coefficient since only velocity fields are involved in the scalar product. In particular at , q l (t o ) B = 1 by construction, and The software FreeFem++ is again used to solve for the velocity and pressure by means of the Finite Element Method with Taylor-Hood elements, (P2 for velocity and P1 for pressure). The practical computation of the gain (B 3) proposed in Garnaud et al. (2013a) is followed. The application of the propagator e Lt (resp. its adjoint (e Lt ) † ) are performed by integrating in time the linearised problem (B 2) (resp. the adjoint problem) with the Crank-Nicolson method. The application of the inverse propagator e −Lt is never needed. Appendix C. Higher-order corrections of the WNNh equation.
Recall the equation (2.9) obtained at order √ o 3 : After imposition of the Fredholm alternative, leading to the equation (2.10) for dA/dT , the relation (C 1) becomes: where ζ = (µ + ν). For higher-order corrections of the WNNh model, the field u 3,1 is needed and is solution of (C 2) where the operator Φ is singular sinceû o = 0 belongs to its kernel. Since by Thus, thanks to the (imposed) orthogonality of the RHS withf o ( RHS,f o = 0), solving the equation replacing Φ by (iω o I − L) leads directly to u 3,1 being orthogonal toû o . Therefore, P u 3,1 = 0 and (iω o I − L)u 3,1 = Φu 3,1 , which implies that the field u 3,1 computed with (iω o I −L) instead of Φ is directly the particular solution of (C 2). Note that u 3,1 appears as a true correction toû o in the sense of the scalar product. This property has the striking and important consequence that the operator Φ never needs to be constructed explicitly, whatever the order of the amplitude equation. The homogeneous part of the solution of (C 2) is arbitrarily proportional toû o . It can be ignored without loss of generality Fujimura (1991). Eventually, the term P u j,1 e iωot + c.c collected at order O( √ o j+2 ) disappears if j 2. This is due to the nullity of u j,1 for even j, and to the nullity of P u j,1 for odd j. Overall, the particular solution at order √ o 3 writes: The equation at order 2 o is assembled as: As mentioned, P u 2,1 = 0 since u 2,1 = 0, and the forcing terms are −2C(u 1 , u 3 ), −C(u 2 , u 2 ) and −∂ T u 2 . We first develop C(u 1 , u 3 ) as: then C(u 2 , u 2 ) as: C(u 2 , u 2 ) = |A| 4 [C(û 2,2 ,û * 2,2 ) + c.c] + |A| 4 C(u 2,0 , u 2,0 ) + 2A 2 |A| 2 e 2iωot C(u 2,0 ,û 2,2 ) + c.c + A 4 e 4iωot C(û 2,2 ,û 2,2 ) + c.c .
Eventually, collecting all terms leads to the following particular solution for u 4 : 10.5 · 10 4 4.2 · 10 4 5.3 · 10 4 1.6 · 10 2 4.0 · 10 2 1.6 · 10 2   table 3 for the plane Poiseuille flow at (Re, k x , k z ) = (3000, 1.2, 0) considered in §2.2 and forced at ω o = 0.3810. Despite a large harmonic gain for ω = 0 as visible in figure 6, the stationary field u 2,0 remains of reasonable amplitude as the associated Reynolds stress forcing C(û o ,û * o ) projects poorly on the most amplified singular mode for ω = 0. However, the same does not hold for the stationary fieldsû 3,1 ||, ||û 3,3 ||) such that until order √ o 3 each order appears as a true correction of the previous one, this does not maintain for order 2 o . As 2 o ||û (a) 4,0 || is of order unity, it cannot be considered as a correction of the order √ o 3 but appears directly at the base flow level, which is asymptotically ill-posed.
Note that for ω o = ω 1 (i.e, the detuning parameter β = 0) the amplitude in the linear regime, A l , reads: which corresponds to the following linear harmonic gain: Canceling the projection of the RHS of (E 5) on b(t), dividing the ensuing relation by b(t), u o , and taking the partial derivative with respect to t leads to = ũ 3 (t), l(t) l(t), l(t) .
To be meaningful, equation (E 6) and its initial condition must be re-written solely in terms of t, which is done by evaluating T = o t and τ = 2 o t. The total derivative of A, denoted d t A, is now needed, as it takes into account the implicit dependence of A; it subject to lim t→0 (α − A( o t, 2 o t)) = 0 so A(t → 0) = α and the amplitude A is extended by continuity in t = 0 so as to impose A(0) = α.
Appendix F. Transient gain sensitivity and comparison with the WNNt model.
We consider a linear system ∂ t u = Lu subject to the initial condition u(0) with ||u(0)|| = 1. The linear transient gain at t = t o writes G o = ||u(t o )||. A variational method is used to derive the variation of the optimal transient gain induced by a small perturbation δL of operator. Let us introduce the Lagrangian where the Lagrange multipliers u † and β enforce the constraints on the state equation and on the norm of the initial condition, respectively. Imposing ∂ u L, δu = 0 for all δu leads to the adjoint equation ∂ t u † = −L † u † , to be integrated backward in time from the terminal condition u † (t o ) = 2u(t o ). Eventually, the gain variation induced by δL is On the other hand, we derived in the main text for the WNNt model: The weakly-nonlinear transient gain squared can be expressed as G(t o ) 2 = (a(t o )/U 0 ) 2 , while the linear gain squared is G 2 o = (1/ o ) 2 , such that 1 We  The sensitivity relation (F 1) is immediately recognised, where δL is here induced by the addition of (U 0 / o ) 2ũ 2 to the base flow. Indeed, U 0 / o = a lin is the linear solution of (F 2) corresponding the limit U 0 → 0, such that the flow field is described in this limit by U (t) = U e + a lin l(t) + a 2 linũ 2 (t) + O( 3 o ).