# Weak nonlinearity for strong nonnormality

Yves-Marie Ducimetière<EMAIL_ADDRESS>Edouard Boujo  François
Gallaire  Laboratory of Fluid Mechanics and Instabilities, EPFL, CH1015
Lausanne, Switzerland

###### Abstract

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 two
application cases highlight two common nonnormal mechanisms in hydrodynamics:
the flow past a backward-facing step, subject to streamwise convective
nonnormal amplification, and the plane Poiseuille flow, subject to lift-up
nonnormality.

††preprint: APS/123-QED

##  I Introduction

Nonlinear dynamical systems can have one or several equilibrium solutions,
which form one of the building blocks of the phase space [1]. The linear
stability of an equilibrium can be deduced from the eigenvalues of the
linearized 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 generalized over the last decades to account for nonlinear
and nonmodal effects, although these two types of correction have generally
been opposed [2], culminating into a paper entitled “Nonlinear normality
versus non-normal linearity” [3]. 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.

###  I.1 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 ([4], §26), weakly nonlinear
analyses using a multiple-scale approach leading to an equation for the
bifurcated mode amplitude $A$ were performed in some pioneering works in the
context of thermal convection [5, 6], parallel shear flows [7, 8, 9] and non-
parallel [10] shear flows. In theses studies, a so-called Stuart-Landau
equation of the form $\mathrm{d}_{T}A=\lambda A-\kappa A\left|A\right|^{2}$ is
obtained as a condition for non-resonance. When the real part of the nonlinear
coefficient is strictly positive, $\Re(\kappa)>0$, the cubic term
$A\left|A\right|^{2}$ is sufficient to capture the saturation amplitude, and
the Stuart-Landau equation is an accurate model for supercritical
bifurcations; otherwise it can be extended to describe subcritical
bifurcations. Amplitude equations, which can also depend on space, are widely
used to describe the spatiotemporal pattern formation in physical systems near
threshold. Beyond hydrodynamics, this occurs in plasma physics, solidification
fronts, nonlinear optics, laser physics, oscillatory chemical reactions,
buckling of elastic rods, and many others fields of study (see [11] for a
comprehensive review). More generally, while the form of the amplitude
equation can often be deduced from symmetry considerations [12, 13], its
coefficients ($\lambda$ and $\kappa$ in the case of the Stuart-Landau
equation) are evaluated with scalar products of fields computed at the
bifurcation point.

Other approaches exist to deduce the normal form, i.e. the amplitude equation
which distillates the quintessence of the nonlinear behavior in the vicinity
of a bifurcation point [14, 15, 16]. Common to all these approaches is the
concept of center manifold, along which the dynamics are slow, while, under a
spectral gap assumption, an adiabatic elimination ensures the slaving of
damped modes.

Regardless of whether the system is forced or freely evolving, an amplitude
equation can only be constructed close to a bifurcation point. Indeed, only
linearized systems with a neutral or weakly damped eigenmode may experience
resonance, whose avoidance condition results in the amplitude equation. Not
all systems possess such eigenmodes, and systems with a significantly damped
spectrum are often encountered. Nevertheless, such systems are still of great
importance in practice, owing to so-called nonnormal amplifying mechanisms.

###  I.2 Strong nonnormality

Upon the choice of a scalar product, a linear operator is nonnormal if it does
not commute with its adjoint. Consequently, its eigenmodes do not form an
orthogonal set, and the response to an initial condition and to time-harmonic
forcing may be highly non-trivial (see [17] for an exhaustive presentation).
This response generally results from an intricate cooperation between a large
amount of eigenmodes. The leading (least stable or most unstable) eigenvalue
solely provides the asymptotic (long-time) linear behavior of the energy of
the unforced system. At finite time, restriction to the leading eigenmode is
generally irrelevant. In particular, a negative growth rate for all
eigenvalues is not a guarantee for the energy to decay monotonously for all
initial conditions: some small-amplitude perturbations may experience a large
transient amplification (Fig. 1a). The same is true for systems subject to
harmonic forcing: they may exhibit strong amplification, much larger than the
inverse of the smallest damping rate, and at forcing frequencies unpredictable
at the sight of the spectrum (Fig. 1b).

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 $\sigma_{1}$ of the eigenspectrum in (c) has
indeed a negative growth rate. (a) In the linear regime, the amplitude of the
solution eventually decays like $\exp(\sigma_{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 $\sigma_{1,i}$. Nonnormal systems can respond at
different frequencies, with an amplification much larger than predicted by
$\sigma_{1,r}$. Nonlinearity may be stabilizing or destabilizing.

Nonnormal operators are traditionally encountered in laser physics (see [17]
§60), which H. J. Landau described by developing the concept of
pseudospectrum, as a pertinent alternative to modal analysis [18, 19]. The
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 [20]. In astrophysics, pseudospectrum
analysis was used very recently to study the stability of black holes [21]. In
network science, a recent study [22] 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 (Fig. 1), and radically change
the behavior of dynamical systems. This is illustrated in [23], 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 linearization of the advective term
$(\bm{U}\cdot\nabla)\bm{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 [24, 25, 26, 2, 27, 28], pipe flow [29] and
parallel boundary layers [25, 30]. Nonnormality is also found in nonparallel
flows [31], for instance spatially developing boundary layers [32, 33, 34, 35,
36], jets [37, 38] and the flow past a backward-facing step [39, 40].
Exhaustive reviews of nonnormality in hydrodynamics can be found in [41, 42].
The crucial role played by nonnormality in the transition to turbulence has
become clear over the years [2, 43, 42]. As mentioned in the context of
nonnormal networks, if the flow is nonnormal, low-energy 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 [2]
is an excellent illustration of this so-called “bypass” scenario.

###  I.3 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
jeopardize a classical, single-mode amplitude equation.

Notwithstanding the relevance and usefulness of fully nonlinear solutions [44,
45], as well as the existence of a fully nonlinear nonnormal stability theory
able to compute nonlinear optimal initial conditions via Lagrangian
optimization [46, 47, 48, 49], 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 Fig. 1 illustrates
possible scenarios for the effects of nonlinearities subsequent to strong
nonnormal amplification. Finally, amplitude equations are useful for flow
control and optimization, as shown for instance in [50] in the more classical
context of a marginally stable flow, displaying little nonnormality, a well
isolated eigenvalue and a sufficiently large spectral gap.

For theses reasons, 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. 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}$. In both cases, we perform
an a priori weakly nonlinear prolongation of the gain, at very low numerical
cost. The applied 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.

A generic nonlinear dynamical system is considered,

$\partial_{t}\bm{U}=N(\bm{U})+\bm{F},\quad\quad\bm{U}(0)=\bm{U}_{0},$ (1)

where $N(*)$ is a nonlinear operator and $\bm{F}$ is a forcing term. An
appropriate and common place to begin the analysis of (1) is to linearize it
around an unforced equilibrium. The latter is denoted $\bm{U}_{e}$ and
satisfies $N(\bm{U}_{e})=\bm{0}$. Around this equilibrium are considered
small-amplitude perturbations in velocity $\epsilon\bm{u}$, forcing
$\epsilon\bm{f}$, and initial condition $\epsilon\bm{u}_{0}$, where
$\epsilon\ll 1$. An asymptotic expansion of (1) in terms of $\epsilon$ can
thus be performed, transforming the nonlinear equation into a succession of
linear ones. The fields $\bm{u}$, $\bm{f}$ and $\bm{u}_{0}$ are recovered at
order $\epsilon$ and linked trough the linear relation

$\partial_{t}\bm{u}=L\bm{u}+\bm{f},\quad\bm{u}(0)=\bm{u}_{0},$ (2)

where $L$ results from the linearization of $N$ around $\bm{U}_{e}$. For fluid
flows governed by the incompressible Navier-Stokes equations,
$L\bm{u}=-(\bm{U}_{e}\cdot\nabla)\bm{u}-(\bm{u}\cdot\nabla)\bm{U}_{e}+Re^{-1}\Delta\bm{u}-\nabla
p(\bm{u})$, where the pressure field $p$ is such that the velocity field
$\bm{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.

Figure 2: Sketch of the flow configurations. (a) Two-dimensional flow over a
backward-facing 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=\pm 1$, and invariant in the $x$
(streamwise) and $z$ (spanwise) directions

##  II Response to Harmonic Forcing

We first derive an amplitude equation for the weakly nonlinear amplification
of time-harmonic forcing
$\bm{f}(\bm{x},t)=\bm{\hat{f}}(\bm{x})e^{i\omega_{o}t}+c.c$ in a linearly
strictly stable system. In the long-time regime, only the same-frequency
harmonic response $\bm{u}(\bm{x},t)=\bm{\hat{u}}(\bm{x})e^{i\omega_{o}t}+c.c$
persists. Injecting the expressions of $\bm{f}$ and $\bm{u}$ in (2) leads to
$\bm{\hat{u}}=(i\omega_{o}I-L)^{-1}\bm{\hat{f}}\doteq
R(i\omega_{o})\bm{\hat{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

$\centering
G(i\omega_{o})=\max_{\bm{\hat{f}}}\frac{\left\|\bm{\hat{u}}\right\|}{\left\|\bm{\hat{f}}\right\|}=\left\|R(i\omega_{o})\right\|\doteq\frac{1}{\epsilon_{o}}.\@add@centering$
(3)

In the following, we choose the $L^{2}$ norm (or “energy” norm) induced by the
Hermitian inner product
$\left\langle\bm{\hat{u}}_{a},\bm{\hat{u}}_{b}\right\rangle=\int_{\Omega}\bm{\hat{u}}_{a}^{H}\bm{\hat{u}}_{b}\mathrm{d}\Omega$
(the superscript $H$ denotes the Hermitian transpose). The operator
$R(i\omega_{o})^{\dagger}$ denotes the adjoint of $R(i\omega_{o})$ under this
scalar product, such that $\left\langle
R(i\omega_{o})\bm{\hat{u}}_{a},\bm{\hat{u}}_{b}\right\rangle=\left\langle\bm{\hat{u}}_{a},R(i\omega_{o})^{\dagger}\bm{\hat{u}}_{b}\right\rangle$,
for any $\bm{\hat{u}}_{a},\bm{\hat{u}}_{b}$. Among all frequencies
$\omega_{o}$, the one leading to the maximum amplification is noted
$\omega_{o,m}$ and associated with an optimal gain
$G(i\omega_{o,m})=1/\epsilon_{o,m}$. The singular value decomposition of
$R(i\omega_{o})$ provides $G(i\omega_{o})=\epsilon_{o}^{-1}$ as the largest
singular value, and the associated pair of right singular vector
$\bm{\hat{f}}_{o}$ and left singular vector $\bm{\hat{u}}_{o}$. The former
represents the optimal forcing, whereas the latter characterizes the long-time
harmonic response reached, after the transients fade away:

$R(i\omega_{o})^{-1}\bm{\hat{u}}_{o}=\epsilon_{o}\bm{\hat{f}}_{o},\quad\left[R(i\omega_{o})^{\dagger}\right]^{-1}\bm{\hat{f}}_{o}=\epsilon_{o}\bm{\hat{u}}_{o},$
(4)

where $||\bm{\hat{f}}_{o}||=||\bm{\hat{u}}_{o}||=1$. Smaller singular values
of $R(i\omega_{o})$ constitute sub-optimal gains, and the associated right
singular vectors are sub-optimal forcing structures. Note that one can express
$\langle\bm{\hat{u}},\bm{\hat{u}}\rangle=\langle
R\bm{\hat{f}},R\bm{\hat{f}}\rangle$ as $\langle
R^{\dagger}R\bm{\hat{f}},\bm{\hat{f}}\rangle$, such that the singular values
of $R(i\omega_{o})$ are also the square root of the eigenvalues of the
symmetric operator $R(i\omega_{o})^{\dagger}R(i\omega_{o})$. An important
implication is that the singular vectors form an orthogonal set for the scalar
product $\langle*,*\rangle$. The practical computation of $\epsilon_{o}$,
$\bm{\hat{f}}_{o}$ and $\bm{\hat{u}}_{o}$ is detailed for the Navier-Stokes
equations in [38], for instance.

Since $L$ is strongly nonnormal, as assumed in the rest of the present study,
none of $\epsilon_{o}$, $\bm{\hat{u}}_{o}$ and $\bm{\hat{f}}_{o}$ are
immediately determined from its spectral (modal) properties. Strong
nonnormality implies $\epsilon_{o}\ll 1$, such that the inverse resolvent
$R(i\omega_{o})^{-1}$ appearing in (4) is almost singular. Perturbing it as

$\begin{split}\Phi\doteq
R(i\omega_{o})^{-1}-\epsilon_{o}P,\quad\text{where}\quad
P=\bm{\hat{f}}_{o}\left\langle\bm{\hat{u}}_{o},\ast\right\rangle,\end{split}$
(5)

leads to $\Phi\bm{\hat{u}}_{o}=\bm{0}$, such that $\Phi$ is exactly singular.
The norm of the perturbation operator is small since $||P||=1$. The field
$\bm{\hat{u}}_{o}$ constitutes the only non-trivial part of the kernel of
$\Phi$, and its associated adjoint mode is $\bm{\hat{f}}_{o}$ (see Appendix A)

We note that $\Phi$ can be rewritten as $\Phi=(i\omega_{o}I-L_{n})$ where
$L_{n}\doteq L+\epsilon_{o}P$, such that (5) 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\omega_{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\omega_{o}$. By construction,
$\epsilon_{o}$ is the smallest possible amplitude of the right-hand side of
(4) for a given $i\omega_{o}$, such that $\epsilon_{o}P$ is the smallest
perturbation of $L$ necessary to relocate an eigenvalue of $L$ on
$i\omega_{o}$. This fact can be formalized with the pseudospectrum theory
outlined in [17]. In the complex plane, $z\in\mathbb{C}$ belongs to the
$\epsilon$-pseudospectrum $\Lambda_{\epsilon}(L)$ if and only if
$\left\|R(z)\right\|\geq 1/\epsilon$. If $E$ is an operator with $||E||=1$,
eigenvalues of $L-\epsilon E$ can lie anywhere inside $\Lambda_{\epsilon}(L)$.
Eigenvalues of $L$ and singularities of $\left\|R(z)\right\|$ thus collide
with the $\epsilon$-pseudospectrum in the limit $\epsilon\rightarrow 0$. As
$\epsilon$ increases, the $\epsilon$-pseudospectrum may touch the imaginary
axis, such that any $z=i\omega_{o}$ can be an eigenvalue of $L-\epsilon E$ if
the amplitude of the perturbation is greater than or equal to
$\epsilon=\left\|R(i\omega_{o})\right\|^{-1}$. We recognize $\epsilon$ as the
inverse gain $\epsilon_{o}$ defined in (3), and thus $E$ as $P$. In
particular, if $\omega_{o}=\omega_{o,m}$, the associated $\epsilon_{o,m}$ is
referred to as the stability radius of $L$ since the
$\epsilon_{o,m}$-pseudospectrum is the first to touch the imaginary axis.

Figure 3:  Natural and perturbed spectra of the flow past a backward-facing
step (sketched in Fig. 2a) at $Re=500$. Blue circles: eigenvalues of the
linearized Navier-Stokes operator $L$. Red dots: eigenvalues of the linear
operator perturbed with
$\epsilon_{o}P=\epsilon_{o}\bm{\hat{f}}_{o}\left\langle\bm{\hat{u}}_{o},*\right\rangle$.
By construction, one eigenvalue of $L_{n}=L+\epsilon_{o}P$ lies on the
imaginary axis. Green isocontour: part of the $\epsilon_{o}$-pseudospectrum of
$L$, where $||R(z)||=1/\epsilon_{o}$. By construction, the
$\epsilon_{o}$-pseudospectrum is contained in the stable half-plane, except at
$i\omega_{o}$ where it touches the neutral axis.

As an illustration of the fact that a small-amplitude perturbation can easily
“neutralize” a nonnormal operator, we consider the Navier-Stokes operator
linearized around the steady flow past a backward-facing step (BFS), sketched
in Fig. 2, at $Re=500$. The most amplified frequency
$\omega_{o}=\omega_{o,m}\approx 0.47$ is associated with $\epsilon_{o}\approx
1.3\cdot 10^{-4}\ll 1$. The spectra of $L$ and $L_{n}$ are shown in Fig. 3,
together with part of the $\epsilon_{o}$-pseudospectrum of $L$. Clearly, the
very small perturbation $\epsilon_{o}P$ locates an eigenvalue exactly onto
$i\omega_{o}$, despite the strong stability of $L$. We stress that neither
$\omega_{o}$ nor $\epsilon_{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 matrix $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 $\epsilon_{o}\ll 1$ constitutes a natural choice of small
parameter. We choose the Navier-Stokes equations for their nonlinear term
$(\bm{U}\cdot\nabla)\bm{U}$, which yields both a nonnormal linearized operator
and a rich diversity of behaviors. The flow is weakly forced by
$\bm{F}=\phi\sqrt{\epsilon_{o}}^{3}\bm{\hat{f}}_{h}e^{i\omega_{o}t}+c.c$,
where $\bm{\hat{f}}_{h}$ is an arbitrary (not necessarily optimal) forcing
structure, and $\phi=O(1)$ is a real prefactor. Imposing
$||\bm{\hat{f}}_{h}||=1$, the forcing amplitude is
$F\doteq\phi\sqrt{\epsilon_{o}}^{3}$. A separation of time scales is invoked
for the flow response: its envelope is assumed to vary on a slow time scale
$T=\epsilon_{o}t$ (such that
$\mathrm{d}_{t}=\partial_{t}+\epsilon_{o}\partial_{T}$). This ensures a
comprehensive distinguished scaling and suggests the following multiple-scale
expansion:

$\displaystyle\bm{U}(t,T)=\bm{U}_{e}+\sqrt{\epsilon_{o}}\bm{u}_{1}(t,T)+\epsilon_{o}\bm{u}_{2}(t,T)+$
$\displaystyle\sqrt{\epsilon_{o}}^{3}\bm{u}_{3}(t,T)+O(\epsilon_{o}^{2}).$ (6)

The velocity field at each order $j$ is then Fourier-expanded as

$\bm{u}_{j}(t,T)=\bm{\overline{u}}_{j,0}(T)+\sum_{m}(\bm{\overline{u}}_{j,m}(T)e^{im\omega_{o}t}+c.c),$
(7)

$m=1,2,3\ldots$. This decomposition is certainly justified in the permanent
regime, of interest in this analysis. The proposed slow dynamics does not aim
to capture the transient regime but flow variations around the permanent
regime. Introducing (6)-(7) into the Navier-Stokes equations and using (5) to
perturb the operator $R(i\omega_{o})^{-1}$ appearing from time derivation
yields

$\displaystyle\begin{split}&\sqrt{\epsilon_{o}}\Big{[}\left(\Phi\bm{\overline{u}}_{1,1}e^{i\omega_{o}t}+c.c\right)+\bm{s}_{1}\Big{]}+\epsilon_{o}\Big{[}\left(\Phi\bm{\overline{u}}_{2,1}e^{i\omega_{o}t}+c.c\right)+\bm{s}_{2}+C(\bm{u}_{1},\bm{u}_{1})\Big{]}+\\\
&\sqrt{\epsilon_{o}}^{3}\Big{[}\left(\Phi\bm{\overline{u}}_{3,1}e^{i\omega_{o}t}+c.c\right)+\bm{s}_{3}+2C(\bm{u}_{1},\bm{u}_{2})+\partial_{T}\bm{u}_{1}+\left(P\bm{\overline{u}}_{1,1}e^{i\omega_{o}t}+c.c\right)\Big{]}+O(\epsilon_{o}^{2})\\\
&=\phi\sqrt{\epsilon_{o}}^{3}\bm{\hat{f}}_{h}e^{i\omega_{o}t}+c.c,\end{split}$
(8)

where

$\bm{s}_{j}\doteq-L\bm{\overline{u}}_{j,0}(T)+\left[\sum_{m}(im\omega_{o}-L)\bm{\overline{u}}_{j,m}(T)e^{im\omega_{o}t}+c.c\right]$

For $m=2,3,...,$ and
$C(\bm{a},\bm{b})\doteq\frac{1}{2}((\bm{a}\cdot\nabla)\bm{b}+(\bm{b}\cdot\nabla)\bm{a})$.
Note that the perturbation $\epsilon_{o}P$ modifying $R(i\omega_{0})^{-1}$
into $\Phi$ at leading order is compensated for at third order. Terms are then
collected at each order in $\sqrt{\epsilon_{o}}$, leading to a cascade of
linear problems, detailed hereafter.

At order $\sqrt{\epsilon_{o}}$, we collect
$(im\omega_{o}I-L)\bm{\overline{u}}_{1,m}=\bm{0}$ for $m=0,2,3\ldots$, and
$\Phi\bm{\overline{u}}_{1,1}=\bm{0}$. Since $L$ is strictly stable, the
unforced equation for $m\neq 1$ can only lead to
$\bm{\overline{u}}_{1,m}=\bm{0}$. Conversely, the kernel of $\Phi$ contains
the optimal response $\bm{\hat{u}}_{o}$, therefore
$\bm{\overline{u}}_{1,1}(T)=A(T)\bm{\hat{u}}_{o}$, where $A(T)\in\mathbb{C}$
is a slowly-varying scalar amplitude verifying $\partial_{t}A=0$. Finally, the
general solution at order $\sqrt{\epsilon_{o}}$ writes

$\bm{u}_{1}(t,T)=A(T)\bm{\hat{u}}_{o}e^{i\omega_{o}t}+c.c.$ (9)

At order $\epsilon_{o}$, we obtain the solution
$\bm{u}_{2}=\left|A\right|^{2}\bm{u}_{2,0}+\left(A^{2}e^{2i\omega_{o}t}\bm{\hat{u}}_{2,2}+c.c\right)$,
where

$\begin{split}-L\bm{u}_{2,0}&=-2C(\bm{\hat{u}}_{o},\bm{\hat{u}}^{*}_{o}),\\\
(2i\omega_{o}I-L)\bm{\hat{u}}_{2,2}&=-C(\bm{\hat{u}}_{o},\bm{\hat{u}}_{o}).\end{split}$
(10)

The homogeneous solution of the system $\Phi\bm{\overline{u}}_{2,1}=\bm{0}$ is
arbitrarily proportional to $\bm{\hat{u}}_{o}$, and written
$A_{2}(T)\bm{\hat{u}}_{o}$. It can be ignored
($\bm{\overline{u}}_{2,1}=\bm{0}$) without loss of generality. As shown in
[51], it could also be kept, provided it is included in the definition of the
amplitude, which would then become $A+\epsilon_{o}A_{2}$.

At order $\sqrt{\epsilon_{o}}^{3}$, we assemble two equations yielding the
Fourier components of the solution oscillating at $\omega_{o}$,

$\begin{split}\Phi\bm{\overline{u}}_{3,1}=-A|A|^{2}\left[2C(\bm{\hat{u}}_{o},\bm{u}_{2,0})+2C(\bm{\hat{u}}_{o}^{*},\bm{\hat{u}}_{2,2})\right]\\\
-\bm{\hat{u}}_{o}\frac{\mathrm{d}A}{\mathrm{d}T}-A\bm{\hat{f}}_{o}+\phi\bm{\hat{f}}_{h}\end{split}$
(11)

(recalling $P\bm{\hat{u}}_{o}=\bm{\hat{f}}_{o}$), and at $3\omega_{o}$,
$(3i\omega_{o}I-L)\bm{\overline{u}}_{3,3}=2A^{3}C(\bm{\hat{u}}_{o},\bm{\hat{u}}_{2,2})$.
The operator $\Phi$ being singular, the only way for $\bm{\overline{u}}_{3,1}$
to be non-diverging, and thus for the asymptotic expansion to make sense, is
that the right-hand side of (11) has a null scalar product with the kernel of
$\Phi^{\dagger}$, i.e. is orthogonal to the adjoint mode $\bm{\hat{f}}_{o}$
associated with $\bm{\hat{u}}_{o}$. This is known as the “Fredholm
alternative”. As a result, the amplitude $A(T)$ satisfies

$\frac{1}{\eta}\frac{\mathrm{d}A}{\mathrm{d}T}=\phi\gamma-A-\frac{\mu+\nu}{\eta}A\left|A\right|^{2},$
(12)

with the coefficients

$\displaystyle\eta=\frac{1}{\left\langle\bm{\hat{f}}_{o},\bm{\hat{u}}_{o}\right\rangle},\quad\gamma=\left\langle\bm{\hat{f}}_{o},\bm{\hat{f}}_{h}\right\rangle,$
$\displaystyle\frac{\mu}{\eta}=\left\langle\bm{\hat{f}}_{o},2C(\bm{\hat{u}}_{o},\bm{u}_{2,0})\right\rangle,\quad\frac{\nu}{\eta}=\left\langle\bm{\hat{f}}_{o},2C(\bm{\hat{u}}^{*}_{o},\bm{\hat{u}}_{2,2})\right\rangle.$
(13)

The coefficient $\gamma$ is the projection of the applied forcing on the
optimal forcing. The coefficient $\mu$ embeds the interaction between
$\bm{\hat{u}}_{o}$ and the static perturbation $\bm{u}_{2,0}$, i.e. it
corrects the gain according to the fact that $\bm{\hat{u}}_{o}$ extracts
energy from the time-averaged mean flow rather than from the original base
flow. We show in Appendix B that, in the regime of small variations around the
linear gain, the amplitude equation reduces to the standard sensitivity of the
gain [52] to a base flow modification induced by $\bm{u}_{2,0}$. In contrast,
the coefficient $\nu$ embeds the interaction between $\bm{\hat{u}}_{o}^{*}$
and the second harmonic $\bm{\hat{u}}_{2,2}$. Introducing the rescaled
quantities $a\doteq\sqrt{\epsilon_{o}}A$ and $F=\phi\sqrt{\epsilon_{o}}^{3}$,
such that the weakly nonlinear harmonic gain
$G=\left\|\sqrt{\epsilon_{o}}\bm{\overline{u}}_{1,1}\right\|/\left\|\phi\sqrt{\epsilon_{o}}^{3}\bm{\hat{f}}_{h}\right\|$
is simply $=\left|a\right|/F$, (12) becomes

$\frac{1}{\eta\epsilon_{o}}\frac{\mathrm{d}a}{\mathrm{d}t}=\frac{\gamma
F}{\epsilon_{o}}-a-\frac{\mu+\nu}{\eta\epsilon_{o}}a\left|a\right|^{2}.$ (14)

The gain associated to the linearized version of (14) is
$G=\left|\gamma\right|/\epsilon_{o}$, as expected for the linear prediction.
We recover $G=1/\epsilon_{o}$ when the optimal forcing is applied
($\gamma=1$). We also note that this expression predicts $G=0$ when
$\gamma=0$, which merely indicates that the linear response is orthogonal to
$\bm{\hat{u}}_{o}$, without stating anything on the gains associated with sub-
optimal forcings except that they should be at most $O(\epsilon_{o}^{-1/2})$,
assuming a sufficiently large “spectral” gap in the singular-value
decomposition of the resolvent operator.

In the nonlinear regime, equilibrium amplitudes $\left|\bar{a}\right|$ are
solution of a third-order polynomial in $Y=\left|\bar{a}\right|^{2}$,

$\displaystyle
Y^{3}\frac{\left|\mu+\nu\right|^{2}}{\epsilon_{o}^{2}\left|\eta\right|^{2}}+2Y^{2}\Re\left[\frac{\mu+\nu}{\epsilon_{o}\eta}\right]+Y=\left(\frac{F}{\epsilon_{o}}\right)^{2}.$
(15)

For a given forcing amplitude, equation (15) predicts either one or three
equilibrium amplitudes. When the coefficient of $Y^{2}$ has a positive real
part, $\Re\left[(\mu+\nu)/(\epsilon_{o}\eta)\right]>0$, (15) predicts a unique
solution and a monotonous decay of the harmonic gain with $F$; conversely,
when $\Re\left[(\mu+\nu)/(\epsilon_{o}\eta)\right]<0$, the nonlinear gain is
larger than the linear one at small $F$, and there might be three solutions.

Equation (14) 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 Fig. 2
is forced harmonically with the optimal structure $\bm{\hat{f}}_{o}$. At
$Re=500$ and the optimal forcing frequency, $\bm{\hat{f}}_{o}$ is shown in
Fig. 4(a) together with is associated response $\bm{\hat{u}}_{o}$ in Fig. 4(b)
(see Appendix C for details about the geometry and the numerical method).

(a)

(b)

Figure 4: (a) Streamwise ($x$) component of the optimal harmonic forcing
structure $\bm{\hat{f}}_{o}$ for the BFS (sketched in Fig. 2a) at $Re=500$ and
at the optimal forcing frequency
$\omega_{o}/(2\pi)=\omega_{o,m}/(2\pi)=0.075$. (b) Streamwise component of the
associated response $\bm{\hat{u}}_{o}$. Both structures are normalized as
$||\bm{\hat{f}}_{o}||=||\bm{\hat{u}}_{o}||=1$.

As shown in [39, 40], the BFS flow constitutes a striking illustration of
streamwise nonnormality. As seen in Fig. 4(a), 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 $\omega_{o}=2\pi\times 0.075$ close
to the most linearly amplified frequency $\omega_{o,m}$, which varies only
slightly with $Re$. The linear gain grows exponentially with $Re$ [40], as
seen in table 1. Since $\eta$ scales like $O(\epsilon_{o}^{-1/2})$, the term
in $\mathrm{d}A/\mathrm{d}T$ in (11) is asymptotically consistent only close
to equilibrium points where $\mathrm{d}A/\mathrm{d}T=0$, which is the regime
of primary interest in the context of harmonic forcing. In accordance, the
temporal derivative $\mathrm{d}A/\mathrm{d}T$ is kept in (12) to assess the
stability of such equilibria, determined by the analysis of the Jacobian
matrix obtained after separating modulus and phase.

Table 1:  WNNh coefficients for the backward-facing step flow, when the optimal forcing structure ($\gamma=1$) is applied at the optimal frequency $\omega_{o}/(2\pi)=\omega_{o,m}/(2\pi)=0.075$. $Re$ | $\epsilon_{0}$ | $\eta$ | $\mu/(\epsilon_{o}\eta)$ | $\nu/(\epsilon_{o}\eta)$
---|---|---|---|---
$200$ | $73.9^{-1}$ | $3.66+i\cdot 0.0163$ | $5.13+i\cdot 1.32$ | $0.137-i\cdot 1.13$
$500$ | $7456.6^{-1}$ | $117.1+i\cdot 0.653$ | $8.23+i\cdot 2.60$ | $0.364+i\cdot 0.396$
$700$ | $148080^{-1}$ | $1626.7+i\cdot 8.65$ | $9.06+i\cdot 4.38$ | $-0.729+i\cdot 1.39$

(a)

(b)

Figure 5:  Weakly and fully nonlinear harmonic gain in the backward-facing
step flow (sketched in Fig. 2a). At each frequency and each Reynolds number,
the optimal linear forcing structure $\bm{\hat{f}}_{o}$ is applied. (a) Fixed
frequency $\omega_{o}/(2\pi)=0.075$, varying Reynolds number $Re=200$ and
$700$ (larger $Re$ darker). Inset: log-log scale, $Re=200,300\ldots 700$. (b)
Fixed Reynolds number $Re=500$, varying forcing $rms$ amplitude
$F=\sqrt{2}^{-1}[1,2,4,10]\cdot 10^{-4}$ (larger amplitudes darker).

Predictions from the WNNh model are compared to fully nonlinear gains
extracted from direct numerical simulations (DNS) in Fig. 5(a). The DNS gains
are the ratio between the temporal rms of the kinetic energy of the
fluctuations at $\omega_{o}$ (extracted through a Fourier transform) and the
rms of the kinetic energy of the forcing (for instance, the forcing
$F\bm{\hat{f}}_{o}e^{i\omega_{o}t}+c.c$ with $||\bm{\hat{f}}_{o}||=1$
corresponds to an effective forcing rms amplitude of $\sqrt{2}F$). Since
$\Re\left[(\mu+\nu)/(\epsilon_{o}\eta)\right]>0$ 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 (14). 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 $\epsilon_{0}$.
Indeed, $F\sim\epsilon_{0}$ implies $\phi\sim 1/\sqrt{\epsilon_{o}}$, which
jeopardizes 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
stability of the finite amplitude limit cycle is for instance ensured by their
combination. The nonlinear coefficients remain of order one, which confirms
the validity of the chosen scalings. The real part of $\mu$ being larger than
that of $\nu$, the present analysis rationalizes a priori the predominance of
the mean flow distortion over the second harmonic in the saturation mechanism
reported a posteriori in [53].

Next, we select $Re=500$ and report in Fig. 5(b) harmonic gains as a function
of the frequency, for increasing forcing amplitudes. At each frequency, the
corresponding optimal forcing structure $\bm{\hat{f}}_{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 $\omega_{o}/(2\pi)=0.04$ where the second harmonics at frequency
$2\omega_{0}$ could in principle be amplified approximately four times more
than the fundamental. It happens, however, that the associated forcing
structure $-C(\bm{\hat{u}}_{o},\bm{\hat{u}}_{o})$ is located much farther
downstream than the optimal forcing at $2\omega_{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.

In should be mentioned that, in some situations, the amplitude equation (14)
may be in default. First, as just mentioned, when the optimal linear harmonic
gain at frequency $2\omega_{o}$ is $\sim 1/\sqrt{\epsilon_{o}}$ or larger and
projects well onto the optimal forcing, the asymptotic hierarchy is threatened
as $\bm{\hat{u}}_{2,2}$ may be substantial enough to reach order
$\sqrt{\epsilon_{o}}$ or above. It is thus important to assess that the norm
of $\bm{\hat{u}}_{2,2}$ remains of order one. A second delicate situation
arises, for the same reason, when a sub-optimal gain at the frequency
$\omega_{o}$ is $\sim 1/\epsilon_{o}$. In both cases, the model could be
extended by including in the kernel of $\Phi$ the optimal response at
frequency $2\omega_{o}$, or the sub-optimal response at frequency
$\omega_{o}$, respectively. For the sake of completeness, Appendix D shows how
to compute higher-order corrections of (14). It is worth mentioning, in
particular, that the action of $\Phi$ need not be computed explicitly and can
be replaced by the action of $(i\omega_{o}I-L)$ for all practical purposes.

##  III Response to Initial Condition

Next, we derive an amplitude equation for the weakly nonlinear transient
growth in an unforced ($\bm{f}=\bm{0}$) system, without restriction on its
linear stability. The solution to the linearized equation (2) is
$\bm{u}(t)=e^{Lt}\bm{u}(0)$, where $e^{Lt}$ is the matrix exponential of $Lt$.
In an unforced context, the propagator $e^{Lt}$ maps an initial structure at
time $t=0$ onto its evolution at $t\geq 0$. The largest linear amplification
at $t_{o}>0$ (subscript $o$ for “optimal”) is

$\centering
G(t_{o})=\max_{\bm{u}(0)}\frac{\left\|\bm{u}(t_{o})\right\|}{\left\|\bm{u}(0)\right\|}=\left\|e^{Lt_{o}}\right\|\doteq\frac{1}{\epsilon_{o}}.\@add@centering$
(16)

The singular value decomposition of the propagator $e^{Lt_{o}}$ provides the
transient gain $G(t_{o})$ as the largest singular value of $e^{Lt_{o}}$, as
well as the left and right singular pair $\bm{v}_{o}$ and $\bm{u}_{o}$,
respectively,

$e^{-Lt_{o}}\bm{v}_{o}=\epsilon_{o}\bm{u}_{o},\quad\left[\left(e^{Lt_{o}}\right)^{\dagger}\right]^{-1}\bm{u}_{o}=\epsilon_{o}\bm{v}_{o},$
(17)

where $||\bm{v}_{o}||=||\bm{u}_{o}||=1$. The field $\bm{u}_{o}$ is the optimal
initial structure for the propagation time $t=t_{o}$, and $\bm{v}_{o}$ is its
normalized evolution at $t_{o}$. The corresponding amplification is
$1/\epsilon_{o}$, as defined in (16). 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^{Lt_{o}}$ also are the eigenvectors of the symmetric operator
$(e^{Lt_{o}})^{\dagger}e^{Lt_{o}}$, 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_{t_{o}>0}G(t_{o})=G(t_{o,m})$.

By construction, the linear gain is independent of the amplitude of the
initial condition $\bm{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,
$\epsilon_{o}\ll 1$. While the previous section focused on the inverse
resolvent, it is now the inverse propagator $e^{-Lt_{o}}$ that appears close
to singular. The first equality of (17) can be rewritten as
$(e^{-Lt_{o}}-\epsilon_{o}\bm{u}_{o}\langle\bm{v}_{o},*\rangle)\bm{v}_{o}=\bm{0}$,
which shows that the operator
$(e^{-Lt_{o}}-\epsilon_{o}\bm{u}_{o}\langle\bm{v}_{o},*\rangle)$ is singular
since $\bm{v}_{o}\neq\bm{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

$\bm{l}(t)\doteq\epsilon_{o}e^{Lt}\bm{u}_{o}$ (18)

seeded by the optimal initial condition $\bm{u}_{o}$ and normalized in
$t=t_{o}$ such that $\bm{l}(t_{o})=\bm{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\geq 0$ as

$\Phi(t)=e^{-Lt}-\epsilon_{o}P(t),\quad\text{where}\quad P(t)\doteq
H(t)\frac{\bm{u}_{o}\langle\bm{l}(t),*\rangle}{\left\|\bm{l}(t)\right\|^{2}},$
(19)

and where the Heaviside distribution $H(t)$ satisfies $H(0)=0$ and $H(t>0)=1$.
The non-trivial kernel of $\Phi(t)$ is $\bm{l}(t)$ for all $t>0$; the kernel
reduces to $\bm{0}$ at $t=0$ since $\Phi(0)=I$. We show in Appendix E that,
for $t>0$, the non-trivial kernel of the adjoint operator $\Phi(t)^{\dagger}$
is

$\begin{split}\bm{b}(t)&\doteq\left(e^{Lt}\right)^{\dagger}\bm{l}(t).\end{split}$
(20)

Figure 6:  Restricted spectra (fifteen least stable eigenvalues) of the
natural and perturbed inverse propagators of the plane Poiseuille flow
(sketched in Fig. 2b) for $t=t_{o}=10$ and $(Re,k_{x},k_{z})=(3000,0.5,2)$.
Blue circles: eigenvalues of $e^{-Lt_{o}}$. Red dots: eigenvalues of
$\Phi(t_{o})$. By construction, one eigenvalue of $\Phi(t)$ lies at the
origin. Thin red lines: full locus of the eigenvalues of $\Phi(t)$ for $t\leq
t_{o}$. Green line: $\epsilon_{o}$-pseudospectrum of $e^{-Lt_{o}}$, such that
$\left\|(e^{-Lt_{o}}-zI)^{-1}\right\|=1/\epsilon_{o}$.

As an illustration of the singularization of $e^{-Lt_{o}}$, parts of the
spectra of $e^{-Lt_{o}}$ and $\Phi(t_{o})$ are shown in Fig. 6 for the plane
Poiseuille flow sketched in Fig. 2, a typical flow with component-wise
nonnormality [2, 42], associated more specifically for the chosen wavenumbers
to the so-called lift-up mechanism. The red dot at the origin is the null
singular eigenvalue of $\Phi(t_{o})$ associated with $\bm{l}(t_{o})$. Since
$||P(t_{o})||=1$, this singular eigenvalue lies on the
$\epsilon_{o}$-pseudospectrum of $e^{-Lt_{o}}$, meaning that a perturbation of
amplitude $\epsilon_{o}$ is sufficient to make the inverse propagator
singular.

Recalling that $L$ is assumed strongly nonnormal, we choose $\epsilon_{o}\ll
1$ as expansion parameter, introduce the slow time scale $T=\epsilon_{o}t$,
and propose the multiple-scale expansion

$\bm{U}(t,T)=\bm{U}_{e}+\epsilon_{o}\bm{u}_{1}(t,T)+\epsilon_{o}^{2}\bm{u}_{2}(t,T)+O(\epsilon_{o}^{3}).$
(21)

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 initialized with
$\bm{U}(0)=\alpha\epsilon_{o}^{2}\bm{u}_{o}$, where $\alpha=O(1)$ is a
prefactor. After injecting this expansion in the unforced Navier-Stokes
equations and perturbing this operator according to (19), the operator
$\Phi(t)$ appears at each order (see Appendix F for details),

$\displaystyle\epsilon_{o}e^{Lt}\partial_{t}(\Phi\bm{u}_{1})+$
$\displaystyle\epsilon_{o}^{2}\Big{[}e^{Lt}\partial_{t}(\Phi\bm{u}_{2})+C(\bm{u}_{1},\bm{u}_{1})+\partial_{T}\bm{u}_{1}+e^{Lt}\partial_{t}(P(t)\bm{u}_{1})\Big{]}+$
$\displaystyle+O(\epsilon_{o}^{3})=\bm{0},$ (22)

subject to $\bm{u}_{2}(0)=\alpha\bm{u}_{o}$, and $\bm{u}_{i}(0)=\bm{0}$ for
$i\neq 2$. Terms are then collected at each order in $\epsilon_{o}$, leading
to a cascade of linear problems, detailed hereafter.

At order $\epsilon_{o}$, we collect $\partial_{t}(\Phi\bm{u}_{1})=\bm{0}$,
subject to $\bm{u}_{1}(0)=\bm{0}$. We obtain
$\Phi\bm{u}_{1}=\Phi(0)\bm{u}_{1}(0)=\bm{0}$, therefore $\bm{u}_{1}(t,T)$ is
proportional to the kernel of $\Phi(t)$ for all $t\geq 0$. We choose the non-
trivial solution

$\bm{u}_{1}(t,T)=A(T)H(t)\bm{l}(t),$ (23)

where the initial condition $\bm{u}_{1}(0)=\bm{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 $\partial_{t}A=0$.

At order $\epsilon_{o}^{2}$, we collect

$\displaystyle\partial_{t}(\Phi\bm{u}_{2})+A^{2}He^{-Lt}C(\bm{l},\bm{l})+H\frac{\mathrm{d}A}{\mathrm{d}T}e^{-Lt}\bm{l}$
$\displaystyle+A\mathrm{d}_{t}(HP\bm{l})=\bm{0},$ (24)

where we used $H^{2}=H$. After a series of manipulations of (24) and time
integration detailed in Appendix G, we obtain for $t>0$

$\displaystyle\Phi\bm{u}_{2}$
$\displaystyle=\alpha\bm{u}_{o}+A^{2}e^{-Lt}\tilde{\bm{u}}_{2}-\frac{\mathrm{d}A}{\mathrm{d}T}\epsilon_{o}t\bm{u}_{o}-A\bm{u}_{o},$
(25)

where

$\frac{\mathrm{d}\tilde{\bm{u}}_{2}}{\mathrm{d}t}=L\tilde{\bm{u}}_{2}-C(\bm{l},\bm{l}),\quad\tilde{\bm{u}}_{2}(0)=\bm{0}.$
(26)

Invoking again the Fredholm alternative, (25) admits a solution if and only if
its right-hand side is orthogonal to $\bm{b}(t)$ for all $t>0$. A series of
manipulations (detailed in Appendix H) results in the amplitude equation

$\frac{\mathrm{d}A}{\mathrm{d}t}=\frac{\mathrm{d}\mu_{2}}{\mathrm{d}t}A^{2},\quad
A(0)=\alpha,$ (27)

where

$\mu_{2}(t)\doteq\frac{\left\langle
e^{-Lt}\tilde{\bm{u}}_{2}(t),\bm{b}(t)\right\rangle}{\left\langle\bm{u}_{o},\bm{b}(t)\right\rangle}=\epsilon_{o}\frac{\left\langle\tilde{\bm{u}}_{2}(t),\bm{l}(t)\right\rangle}{\left\langle\bm{l}(t),\bm{l}(t)\right\rangle}.$
(28)

We stress that the inverse propagator is not needed to solve the amplitude
equation (27). Just like the original problem considered in this section, (27)
is unforced and has a non-zero initial condition. In the linear regime,
$A=\alpha$ for all times, and the linear gain is
$\left\|\epsilon_{o}\alpha\bm{l}(t)\right\|/\left\|\alpha\epsilon_{o}^{2}\right\|=\left\|\bm{l}(t)\right\|/\epsilon_{o}$.
At $t=t_{o}$, in particular, we recover that it is equal to $1/\epsilon_{o}$
since $\left\|\bm{l}(t_{o})\right\|=\left\|\bm{v}_{o}\right\|=1$.

In the following, we call Equation (27) the Weakly Nonlinear Nonnormal
transient (WNNt) model. It can be corrected with higher-order terms, which
requires solving the linear singular system (25), as detailed in Appendix I.
We show in particular that singular higher-order solutions are orthogonal to
the first-order order solution $\bm{l}(t)$.

The WNNt is applied to the plane Poiseuille flow. While the base flow $U(y)$
has only one velocity component and depends only on one coordinate,
perturbations can be fully three dimensional. 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},t_{o})=(3000,0,2,t_{o,m}=230)$ is selected. See Appendix C
for details about the numerical method. The linear optimal gain is validated
with the result of [28]; the associated optimal initial condition and its
evolution at $t=t_{o}$ are shown in Fig. 7(a) and Fig. 7(b), 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 center
and high-velocity fluid from the center of the channel towards the walls, thus
generating alternated streamwise streaks.

(a)

(b)

Figure 7: (a) Optimal initial condition $\bm{u}_{o}$ for the plane Poiseuille
flow (sketched in Fig. 2b) for $(Re,k_{x},k_{z})=(3000,0,2)$ and
$t_{o}=t_{o,m}=230$. Arrows: cross-sectional velocity field
($u_{o,z}$,$u_{o,y}$). Contours: streamwise component $u_{o,x}$. (b) Evolution
$\bm{v}_{o}$ at $t=t_{o}$. Both fields are normalized as
$||\bm{u}_{o}||=||\bm{v}_{o}||=1$. Initial vortices have a null streamwise
component $u_{o,x}$, and streaks at $t=t_{o}$ have negligible cross-sectional
components($v_{o,z}$,$v_{o,y}$). Only one wavelength $-\pi\leq\beta z\leq\pi$
is shown.

Due to the spanwise periodicity of the optimal initial condition $\bm{u}_{o}$,
all the solutions at even orders $\epsilon_{o}^{2n}$ $(n=1,2,3,...)$ only
yield even spatial harmonics in $k_{z}$ and are orthogonal to $\bm{l}(t)$,
such that the coefficient $\mu_{2}(t)$ defined in (28) is null at all times.
In this flow, the WNNt model therefore needs to be extended to order
$\epsilon_{o}^{3}$ (see Appendix J), yielding

$\frac{\mathrm{d}A}{\mathrm{d}t}=\frac{\mathrm{d}\mu_{3}}{\mathrm{d}t}\epsilon_{o}A^{3},\quad
A(0)=\alpha,$ (29)

where

$\mu_{3}(t)\doteq\epsilon_{o}\frac{\left\langle\tilde{\bm{u}}_{3}(t),\bm{l}(t)\right\rangle}{\left\langle\bm{l}(t),\bm{l}(t)\right\rangle},$
(30)

and

$\frac{\mathrm{d}\tilde{\bm{u}}_{3}}{\mathrm{d}t}=L\tilde{\bm{u}}_{3}-2C(\tilde{\bm{u}}_{2},\bm{l}),\quad\tilde{\bm{u}}_{3}(0)=\bm{0},$
(31)

with $\tilde{\bm{u}}_{2}$ defined in (26). Equation (29) is similar to (27),
except that the quadratic term $A^{2}$ is replaced by a cubic term $A^{3}$. We
formulate the analytical solution in terms of the rescaled quantities
$a=\epsilon_{o}A$ and the amplitude of the initial condition
$U_{0}=||\bm{U}(0)||=\alpha\epsilon_{o}^{2}$:

$a(t)=\frac{U_{0}}{\epsilon_{o}}\left[1-\left(\frac{U_{0}}{\epsilon_{o}}\right)^{2}\frac{2\mu_{3}(t)}{\epsilon_{o}}\right]^{-1/2}.$
(32)

The weakly nonlinear transient gain becomes $G(t_{o})=a(t_{o})/U_{0}$. We show
in Appendix K that, at first order in the gain variation, (32) reduces to the
sensitivity of the transient gain to the base flow modification
$(U_{0}/\epsilon_{o})^{2}\tilde{\bm{u}}_{2}(t)$.

Predictions from equation (32) are shown in Fig. 8 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}$, or that perturbation
projected on $\bm{l}(t)$, for $G_{\bm{l}}$ (see Appendix L for explicit
expressions).

(a)

(b)

Figure 8:  Transient gain in the plane Poiseuille flow (sketched in Fig. 2b)
for $(Re,k_{x},k_{z})=(3000,0,2)$. (a) Gain squared $G(t_{o})^{2}$ for
$t_{o}=t_{o,m}=230$ as a function of the amplitude of the initial condition.
Streamwise invariance $k_{x}=0$ is enforced in the DNS as well. (b) History of
the gain squared for $0\leq t\leq t_{o,m}$ and for three amplitudes of initial
condition, $U_{0}/\epsilon_{o}=[0.088,0.18,0.37]$ (vertical dashed lines in
(a)); larger amplitudes darker. Inset: weakly nonlinear coefficient
$\mu_{3}(t)$ as a function of time.

The WNNt model predicts $G_{\bm{l}}$ accurately in the weakly nonlinear regime
for $t=t_{o,m}$, which supports our approach (Fig. 8(a)). In the strongly
nonlinear regime, beyond $U_{0}/\epsilon_{o}\approx 0.4$, the model
overestimates $G_{\bm{l}}$. This can be interpreted by noting that $G_{tot}$
is twice as large as $G_{\bm{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 (32) 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.

Figure 8(b) compares for $t\leq t_{o}$ the history of the approximated gain
$a(t)||\bm{l}(t)||/U_{0}$ with that of the DNS gain $G_{\bm{l}}$, and shows a
convincing overall agreement. The coefficient $\mu_{3}(t)$ is negative and
decays monotonously with time until $\mu_{3}(t_{o})=-0.0786$ (inset),
enhancing the saturation. This results in a reduction of the approximated
optimal time with forcing amplitude, as also observed in the DNS.

##  IV 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 in parallel and non-parallel two-dimensional flows. In both cases,
they predict accurately the nonlinear saturation 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.

It should be noted that the proposed method is not restricted to the Navier-
Stokes equations, but apply to all nonlinear systems whose linearized operator
exhibits a strong nonnormality (see [17] §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 [54]. 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 to other rich nonlinear phenomena such as sub-critical behavior,
currently under investigation, and three-dimensional flows. For the transient
growth aspect, the assumption of a time-independent state matrix and the
associated matrix 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. Their nonlinear interactions are particularly relevant when
distinct structures lead to comparable gains, for instance optimal and sub-
optimal initial conditions [25, 39], or perturbations of different spatial
wavenumbers like in jet flows forced with different azimuthal wavenumbers
[38]. 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 [27] and [55] with linear and self-consistent models,
respectively.

##  Appendix A Harmonic forcing: kernel of the perturbed inverse resolvent.

We show that
$\bm{\hat{f}}_{o}=\epsilon_{o}R(i\omega_{o})^{\dagger}\bm{\hat{u}}_{o}$, with
$\left\langle\bm{\hat{f}}_{o},\bm{\hat{f}}_{o}\right\rangle=1$, constitutes
the non-trivial kernel of $\Phi^{\dagger}$. Using that
$P^{\dagger}=\bm{\hat{u}}_{o}\langle\bm{\hat{f}}_{o},*\rangle$, we have

$\displaystyle\begin{split}\Phi^{\dagger}\bm{\hat{f}}_{o}&=\left[R(i\omega_{o})^{-1}\right]^{\dagger}\bm{\hat{f}}_{o}-\epsilon_{o}\bm{\hat{u}}_{o}\left\langle\bm{\hat{f}}_{o},\bm{\hat{f}}_{o}\right\rangle\\\
&=\left[R(i\omega_{o})^{\dagger}\right]^{-1}\bm{\hat{f}}_{o}-\epsilon_{o}\bm{\hat{u}}_{o}=\bm{0},\end{split}$

where we used the fact that the inverse of the adjoint is the adjoint of the
inverse.

##  Appendix B Harmonic forcing: sensitivity.

Let $G_{o}=1/\epsilon_{o}$ designate the linear harmonic gain. Then
$R^{\dagger}R\bm{\hat{f}}_{o}=G_{o}^{2}\bm{\hat{f}}_{o}$ holds by definition,
and implies $G_{o}^{2}=\left\langle
R^{\dagger}R\bm{\hat{f}}_{o},\bm{\hat{f}}_{o}\right\rangle$ thanks to the
chosen normalization $||\bm{\hat{f}}_{o}||=1$. We are interested in the
squared gain variation $\delta G_{o}^{2}$ (where this notation does not
designate the square of the gain variation) induced by a small perturbation
$\delta L$ of the state operator. The latter results in the following
perturbation $\delta R$ of the resolvent:

$\displaystyle\delta R$ $\displaystyle=(i\omega I-L-(\delta L))^{-1}-(i\omega
I-L)^{-1}$ $\displaystyle=\left[(i\omega I-L)(I-R(\delta L))\right]^{-1}-R$
$\displaystyle\approx(I+R\delta L)R-R$ $\displaystyle=R(\delta L)R.$

The gain variation is therefore

$\displaystyle\delta G_{o}^{2}$
$\displaystyle=\langle\delta(R^{\dagger}R)\bm{\hat{f}}_{o},\bm{\hat{f}}_{o}\rangle=\left\langle(\delta
R)^{\dagger}R\bm{\hat{f}}_{o}+R^{\dagger}(\delta
R)\bm{\hat{f}}_{o},\bm{\hat{f}}_{o}\right\rangle$ $\displaystyle=\left\langle
R^{\dagger}(\delta R)\bm{\hat{f}}_{o},\bm{\hat{f}}_{o}\right\rangle+c.c$
$\displaystyle=\left\langle R(\delta
L)R\bm{\hat{f}}_{o},R\bm{\hat{f}}_{o}\right\rangle+c.c$
$\displaystyle=\left\langle(\delta
L)G_{o}\bm{\hat{u}}_{o},G_{o}^{2}\bm{\hat{f}}_{o}\right\rangle+c.c.,$

so finally

$\displaystyle\delta G_{o}^{2}=2G_{o}^{3}\Re\left[\left\langle(\delta
L)\bm{\hat{u}}_{o},\bm{\hat{f}}_{o}\right\rangle\right].$ (33)

For instance, a base flow modification $\delta\bm{U}_{e}$ results in $(\delta
L)\bm{\hat{u}}_{o}=-(\bm{\hat{u}}_{o}\cdot\nabla)\delta\bm{U}_{e}-(\delta\bm{U}_{e}\cdot\nabla)\bm{\hat{u}}_{o}=-2C(\bm{\hat{u}}_{o},\delta\bm{U}_{e})$
yielding the same formula as in [52] with a different normalization.

On the other hand, the WNNh model predicts:

$\displaystyle
Y^{3}\frac{\left|\mu+\nu\right|^{2}}{\epsilon_{o}^{2}\left|\eta\right|^{2}}+2Y^{2}\Re\left[\frac{\mu+\nu}{\epsilon_{o}\eta}\right]+Y=\left(\frac{F}{\epsilon_{o}}\right)^{2}.$
(34)

where $Y=|\bar{a}|^{2}$. We identify the weakly nonlinear harmonic gain as
$G^{2}=Y/F^{2}$, and multiply (34) by $\epsilon_{o}^{2}/Y$:

$\displaystyle
Y^{2}\frac{\left|\mu+\nu\right|^{2}}{\left|\eta\right|^{2}}+\frac{2Y}{G_{o}}\Re\left[\frac{\mu+\nu}{\eta}\right]+\frac{1}{G_{o}^{2}}-\frac{1}{G^{2}}=0.$
(35)

Being interested in small variations around $G_{o}^{2}$ (that correspond to
the linear limit $Y=|\bar{a}|^{2}\rightarrow 0$), we write
$G^{2}=G_{o}^{2}+\delta G_{o}^{2}$ with $|\delta G_{o}^{2}/G_{o}^{2}|\ll 1$.
In this manner, $1/G_{o}^{2}-1/G^{2}=\delta G_{o}^{2}/G_{o}^{4}+h.o.t$,
eventually leading to

$\displaystyle\delta
G_{o}^{2}=-2G_{o}^{3}\Re\left[\frac{|\bar{a}|^{2}(\mu+\nu)}{\eta}\right]+h.o.t.$
(36)

We recognize at leading order equation (33) where $(\delta
L)\bm{\hat{u}}_{o}=-|\bar{a}|^{2}[2C\left(\bm{\hat{u}}_{o},\bm{u}_{2,0}\right)+2C\left(\bm{\hat{u}}_{o}^{*},\bm{\hat{u}}_{2,2}\right)]$.
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}\bm{u}_{2,0}$, and embeds the effect of the second harmonic
$\bm{\hat{u}}_{2,2}$ as well.

##  Appendix C Harmonic forcing/transient growth: the peculiar case of the
Navier–Stokes equations.

The incompressible Navier–Stokes equations write after linearizing around the
equilibrium velocity field $\bm{U}_{e}$

$\begin{split}B\frac{\mathrm{d}\bm{q}}{\mathrm{d}t}=L\bm{q}+\bm{d},\end{split}$
(37)

with the state vector $\bm{q}=[\bm{u},p]^{T}$, the forcing
$\bm{d}=[\bm{f},0]^{T}$, the singular mass matrix

$B=\begin{bmatrix}I&0\\\ 0&0\end{bmatrix},\quad$

and the linearized Navier–Stokes operator

$L=\begin{bmatrix}-(\bm{U}_{e}\cdot\nabla)*-(*\cdot\nabla)\bm{U}_{e}+Re^{-1}\Delta(*)\
&\nabla(*)\\\ \nabla\cdot(*)&0\end{bmatrix}.$

Several subtleties arise from the peculiarity of the pressure variable, that
ensures the instantaneous satisfaction of the incompressibility condition: (i)
the absence of time-derivative 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 generalized as $R(i\omega_{o})=\left(i\omega_{o}B-L\right)^{-1}$,
and the gain is measured according to

$G^{2}(i\omega_{o})=\frac{\langle\hat{\bm{q}},\hat{\bm{q}}\rangle_{B}}{\langle\hat{\bm{d}},\hat{\bm{d}}\rangle},$

where we used the following scalar products

$\begin{split}\langle\bm{\hat{q}}_{a},\bm{\hat{q}}_{b}\rangle_{B}&=\int_{\Omega}\left(\hat{u}_{a,x}^{*}\hat{u}_{b,x}+\hat{u}_{a,y}^{*}\hat{u}_{b,y}+\hat{u}_{a,z}^{*}\hat{u}_{b,z}\right)\mathrm{d}\Omega,\quad\text{and}\\\
\langle\bm{\hat{q}}_{a},\bm{\hat{q}}_{b}\rangle&=\int_{\Omega}\left(\hat{u}_{a,x}^{*}\hat{u}_{b,x}+\hat{u}_{a,y}^{*}\hat{u}_{b,y}+\hat{u}_{a,z}^{*}\hat{u}_{b,z}+\hat{p}_{a}^{*}\hat{p}_{b}\right)\mathrm{d}\Omega.\\\
\end{split}$

The $B$-scalar product excludes pressure, such that the pseudonorm
$\langle\bm{\hat{q}},\bm{\hat{q}}\rangle_{B}=||\bm{\hat{q}}||_{B}^{2}$ is the
total kinetic energy of the response. The scalar product at the denominator
includes pressure, although this will not change the norm of $\bm{\hat{d}}$,
namely $\langle\bm{\hat{d}},\bm{\hat{d}}\rangle=||\bm{\hat{d}}||^{2}$, since
we have no source/sink of mass. The weakly nonlinear coefficients must be
considered under these scalar products. Let
$\bm{\hat{q}}_{o}=[\bm{\hat{u}}_{o},p_{o}]^{T}$ with
$||\bm{\hat{q}}_{o}||_{B}=1$, $\bm{\hat{d}}_{o}=[\bm{\hat{f}}_{o},0]^{T}$ with
$||\bm{\hat{d}}_{o}||=1$, and $\bm{\hat{d}}_{h}=[\bm{\hat{f}}_{h},0]^{T}$ with
$||\bm{\hat{d}}_{h}||=1$, then:

$\begin{split}&\gamma=\langle\bm{\hat{d}}_{o},\bm{\hat{d}}_{h}\rangle,\quad
1/\eta=\langle\bm{\hat{d}}_{o},B\bm{\hat{q}}_{o}\rangle=\langle\bm{\hat{d}}_{o},\bm{\hat{q}}_{o}\rangle_{B},\\\
&\mu/\eta=\langle\bm{\hat{d}}_{o},\bm{\hat{d}}_{3,1}^{(3)}\rangle,\quad\nu/\eta=\langle\bm{\hat{d}}_{o},\bm{\hat{d}}_{3,1}^{(4)}\rangle,\end{split}$

where $\bm{\hat{d}}_{3,1}^{(3)}=[\bm{\hat{f}}_{3,1}^{(3)},0]^{T}$,
$\bm{\hat{d}}_{3,1}^{(4)}=[\bm{\hat{f}}_{3,1}^{(4)},0]^{T}$,
$\bm{\hat{f}}_{3,1}^{(3)}=2C(\bm{\hat{u}}_{o},\bm{u}_{2,0})$ and
$\bm{\hat{f}}_{3,1}^{(4)}=2C(\bm{\hat{u}}^{*}_{o},\bm{\hat{u}}_{2,2})$. The
pressure field has no influence on the weakly nonlinear coefficients. For
instance

$\begin{split}1/\eta&=\int_{\Omega}\left(\hat{f}_{o,x}^{*}\hat{u}_{o,x}+\hat{f}_{o,y}^{*}\hat{u}_{o,y}+\hat{f}_{o,z}^{*}\hat{u}_{o,z}\right)\mathrm{d}\Omega,\quad\text{and}\\\
\mu/\eta&=\int_{\Omega}\left(\hat{f}_{o,x}^{*}\hat{f}_{(3,1),x}^{(3)}+\hat{f}_{o,y}^{*}\hat{f}_{(3,1),y}^{(3)}+\hat{f}_{o,z}^{*}\hat{f}_{(3,1),z}^{(3)}\right)\mathrm{d}\Omega.\end{split}$

For the two-dimensional flow past a BFS, 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 [38]. Finally, DNS are performed by applying a
time scheme based on the characteristic–Galerkin method as described in [56].
For the validation of the codes with existing literature and the mesh
convergence, we refer to [53] where the same codes have been used. The length
of the outlet channel is chosen as $L_{out}=50$ for $Re\leq 500$ [53],
$L_{out}=65$ for $Re=600$, and $L_{out}=80$ for $Re=700$. This ensures the
convergence of the linear gain and the weakly nonlinear coefficients.

For the transient growth model, the linearized problem writes

$B\frac{\mathrm{d}\bm{q}}{\mathrm{d}t}=L\bm{q}\quad\text{subject
to}\quad\bm{q}(0)=\begin{bmatrix}\bm{u}(0)\\\ 0\end{bmatrix},$ (38)

and the gain is measured as

$G(t_{o})^{2}=\frac{\langle\bm{q}(t_{o}),\bm{q}(t_{o})\rangle_{B}}{\langle\bm{q}(0),\bm{q}(0)\rangle_{B}},$
(39)

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 $\mu_{2}(t)$ writes

$\mu_{2}(t)=\epsilon_{o}\frac{\left\langle\tilde{\bm{q}}_{2}(t),\bm{q}_{\bm{l}}(t)\right\rangle_{B}}{\left\langle\bm{q}_{\bm{l}}(t),\bm{q}_{\bm{l}}(t)\right\rangle_{B}},$
(40)

where $\bm{q}_{\bm{l}}(t)=[\bm{l}(t),p_{\bm{l}}(t)]^{T}$, and where
$\tilde{\bm{q}}_{2}=[\tilde{\bm{u}}_{2}(t),\tilde{p}_{2}(t)]^{T}$ is solution
of

$B\frac{\mathrm{d}\tilde{\bm{q}}_{2}}{\mathrm{d}t}=L\tilde{\bm{q}}_{2}-\begin{bmatrix}C(\bm{l},\bm{l})\\\
0\end{bmatrix},\quad\tilde{\bm{q}}_{2}(0)=\bm{0}.$ (41)

Again, pressure does not influence the weakly nonlinear coefficient since only
velocity fields are involved in the scalar product. In particular at
$t=t_{o}$, $\bm{q}_{\bm{l}}(t_{o})=[\bm{v}_{o},p_{\bm{l}}(t_{o})]^{T}$ thus
$\left\langle\bm{q}_{\bm{l}}(t_{o}),\bm{q}_{\bm{l}}(t_{o})\right\rangle_{B}=1$
by construction, and

$\mu_{2}(t_{o})=\epsilon_{o}\int_{\Omega}\tilde{u}_{2,x}v_{o,x}+\tilde{u}_{2,y}v_{o,y}+\tilde{u}_{2,z}v_{o,z}\mathrm{d}\Omega.$
(42)

For the plane Poiseuille flow, the software FreeFem++ is used to solve for
($u_{x}$,$u_{y}$,$u_{z}$,$p$) by means of the Finite Element Method with
Taylor-Hood (P2, P2, P2, P1) elements, respectively. In all 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. The computations are
systematically performed in the spanwise-periodic box
$(y,z)\in[-1,1]\times[-\pi/k_{z},\pi/k_{z}]\equiv\Omega$. All the scalar
products are taken upon integration inside this periodic box, in particular
for the normalization
$\langle\bm{u}_{o},\bm{u}_{o}\rangle=\langle\bm{v}_{o},\bm{v}_{o}\rangle=1$
and the evaluation of the weakly nonlinear coefficient. The practical
computation of the gain (39) proposed in [37] is followed. The application of
the propagator $e^{Lt}$ (resp. its adjoint $(e^{Lt})^{\dagger}$) are performed
by integrating in time the linearized problem (38) (resp. the adjoint problem)
with the Crank-Nicolson method. The application of the inverse propagator
$e^{-Lt}$ is never needed. Our linear optimization codes are validated upon
comparison with the open-source results of [28], obtained with a Chebyshev
polynomials discretization, and where the singular value decomposition of the
matrix exponential $e^{Lt}$ is performed directly. For the chosen set of
parameter $(Re,k_{x},k_{z},t_{o})=(3000,0,2,230)$, convergence was achieved
for a squared linear gain of $G(t_{o})^{2}=1761.8$, against
$G(t_{o})^{2}=1761.9$ in [28].

##  Appendix D Harmonic forcing: higher-order corrections of the amplitude
equation.

$\displaystyle\Phi\bm{\overline{u}}_{3,1}=-A|A|^{2}\left[2C(\bm{\hat{u}}_{o},\bm{u}_{2,0})+2C(\bm{\hat{u}}_{o}^{*},\bm{\hat{u}}_{2,2})\right]$
$\displaystyle-\bm{\hat{u}}_{o}\frac{\mathrm{d}A}{\mathrm{d}T}-A\bm{\hat{f}}_{o}+\phi\bm{\hat{f}}_{h},$
(43)

For higher-order corrections of the WNNh model, the field
$\bm{\overline{u}}_{3,1}$ is needed and is solution of (43) where the operator
$\Phi$ is singular since $\bm{\hat{u}}_{o}\neq\bm{0}$ belongs to its kernel.
Since by construction
$\bm{\hat{f}}_{o}=\epsilon_{o}R(i\omega_{o})^{\dagger}\bm{\hat{u}}_{o}$, it
follows immediately that
$\langle\mbox{RHS},\bm{\hat{f}}_{o}\rangle=\epsilon_{o}\langle
R(i\omega_{o})\mbox{RHS},\bm{\hat{u}}_{o}\rangle$. Thus, thanks to the
(imposed) orthogonality of the RHS with $\bm{\hat{f}}_{o}$, solving the
equation replacing $\Phi$ by $(i\omega_{o}I-L)$ leads directly to
$\bm{\overline{u}}_{3,1}$ being orthogonal to $\bm{\hat{u}}_{o}$. Note that
$\bm{\overline{u}}_{3,1}$ appears as a true correction to $\bm{\hat{u}}_{o}$
in the sense of the scalar product. In addition, $P\bm{\overline{u}}_{3,1}=0$
and $(i\omega_{o}I-L)\bm{\overline{u}}_{3,1}=\Phi\bm{\overline{u}}_{3,1}$,
which implies that the field $\bm{\overline{u}}_{3,1}$ computed with
$(i\omega_{o}I-L)$ instead of $\Phi$ is directly the particular solution of
(43). This property has the striking and important consequence that the
operator $\Phi$ never needs to be constructed explicitly, whatever the order
of the amplitude equation. The homogeneous part of the solution of (43) is
arbitrarily proportional to $\bm{\hat{u}}_{o}$. They can be ignored without
loss of generality [51]. Eventually, the term
$P\bm{\overline{u}}_{j,1}e^{i\omega_{o}t}+c.c$ collected at order
$O(\sqrt{\epsilon_{o}}^{j+2})$ disappears if $j\geq 2$. This is due to the
nullity of $\bm{\overline{u}}_{j,1}$ for even $j$, and to the nullity of
$P\bm{\overline{u}}_{j,1}$ for odd $j$.

##  Appendix E Transient growth: kernel of the perturbed inverse propagator.

We show that $\bm{b}(t)=\left(e^{Lt}\right)^{\dagger}\bm{l}(t)$ constitutes
the non-trivial kernel of $\Phi(t)^{\dagger}$ for $t>0$. Using that
$P^{\dagger}=\bm{l}(t)\langle\bm{u}_{o},*\rangle/\left\langle\bm{l}(t),\bm{l}(t)\right\rangle$
for $t>0$, we have

$\begin{split}\Phi(t)^{\dagger}\bm{b}(t)&=\left(e^{-Lt}\right)^{\dagger}\bm{b}(t)-\epsilon_{o}\bm{l}(t)\frac{\left\langle\bm{u}_{o},\bm{b}(t)\right\rangle}{\left\langle\bm{l}(t),\bm{l}(t)\right\rangle}\\\
&=\left(e^{-Lt}\right)^{\dagger}\bm{b}(t)-\epsilon_{o}\bm{l}(t)\frac{\left\langle
e^{Lt}\bm{u}_{o},\bm{l}(t)\right\rangle}{\left\langle\bm{l}(t),\bm{l}(t)\right\rangle}\\\
&=\left(e^{-Lt}\right)^{\dagger}\bm{b}(t)-\bm{l}(t)\\\
&=\left[\left(e^{Lt}\right)^{\dagger}\right]^{-1}\bm{b}(t)-\bm{l}(t)\\\
&=\bm{0},\end{split}$

where we used the fact that the inverse of the adjoint is the adjoint of the
inverse.

##  Appendix F Transient growth: re-writing the multiple scale expansion.

After injecting the flow field expansion in the unforced Navier-Stokes
equations, we obtain

$\epsilon_{o}(\partial_{t}-L)\bm{u}_{1}+\epsilon_{o}^{2}\left[(\partial_{t}-L)\bm{u}_{2}+C(\bm{u}_{1},\bm{u}_{1})+\partial_{T}\bm{u}_{1}\right]+O(\epsilon_{o}^{3})=\bm{0},$
(44)

In its primary quality of inverse propagator, the following property holds for
$e^{-Lt}$: $\partial_{t}(e^{-Lt})=-e^{-Lt}L$, where the commutation of
$e^{-Lt}$ and $L$ has not been used. Thanks to this relation, we write
$(\partial_{t}-L)\bm{u}_{i}=e^{Lt}\partial_{t}(e^{-Lt}\bm{u}_{i})$. As a
result, $L$ disappears from the asymptotic expansion but $e^{-Lt}$ appears.
The latter is perturbed according to $\Phi(t)=e^{-Lt}-\epsilon_{o}P(t)$,
leading to
$e^{Lt}\partial_{t}(e^{-Lt}\bm{u}_{i})=e^{Lt}\partial_{t}(\Phi(t)\bm{u}_{i})+\epsilon_{o}e^{Lt}\partial_{t}(P(t)\bm{u}_{i})$
for $i=1,2,...$. The asymptotic expansion (44) becomes

$\epsilon_{o}e^{Lt}\partial_{t}(\Phi\bm{u}_{1})+\epsilon_{o}^{2}\Big{[}e^{Lt}\partial_{t}(\Phi\bm{u}_{2})+C(\bm{u}_{1},\bm{u}_{1})+\partial_{T}\bm{u}_{1}+e^{Lt}\partial_{t}(P(t)\bm{u}_{1})\Big{]}+O(\epsilon_{o}^{3})=\bm{0}.$
(45)

Note that the transformation performed from (44) to (45) is not restricted to
time-independent base flows, as the property
$\partial_{t}(\Psi(t)^{-1})=-\Psi(t)^{-1}L(t)$ holds for a time-varying
operator $L(t)$ and the associated propagator $\Psi(t)$. This can be shown
easily by taking the time derivative of $\Psi(t)^{-1}\bm{u}(t)=\bm{u}(0)$.

##  Appendix G Transient growth: algebraic manipulations.

At order $\epsilon_{o}^{2}$, we collect

$\partial_{t}(\Phi\bm{u}_{2})+A^{2}He^{-Lt}C(\bm{l},\bm{l})+H\frac{\mathrm{d}A}{\mathrm{d}T}e^{-Lt}\bm{l}+A\mathrm{d}_{t}(HP\bm{l})=\bm{0},$
(46)

subject to $\bm{u}_{2}(0)=\alpha\bm{u}_{o}$. We used the property
$H(t)^{2}=H(t)$, which will henceforth be understood. The particular solution
of (46) yields

$\displaystyle\bm{u}_{2}(t,T)=$
$\displaystyle\bm{u}_{2}^{(a)}(t)+A(T)^{2}\bm{u}_{2}^{(b)}(t)+\frac{\mathrm{d}A(T)}{\mathrm{d}T}\bm{u}_{2}^{(c)}(t)+A(T)\bm{u}_{2}^{(d)}(t),$
(47)

where

$\displaystyle\mathrm{d}_{t}\left[\Phi\bm{u}_{2}^{(a)}\right]=\bm{0},\quad\mathrm{d}_{t}\left[\Phi\bm{u}_{2}^{(b)}\right]=-He^{-Lt}C(\bm{l},\bm{l}),$
$\displaystyle\mathrm{d}_{t}\left[\Phi\bm{u}_{2}^{(c)}\right]=-He^{-Lt}\bm{l},\quad\text{and}\quad\mathrm{d}_{t}\left[\Phi\bm{u}_{2}^{(d)}\right]=-\mathrm{d}_{t}\left[HP\bm{l}\right],$

subject to the initial conditions $\bm{u}_{2}^{(a)}(0)=\alpha\bm{u}_{o}$ and
$\bm{u}_{2}^{(b)}(0)=\bm{u}_{2}^{(c)}(0)=\bm{u}_{2}^{(d)}(0)=\bm{0}$. Time
integration can now be performed without ambiguity as all the partial
derivatives ($\partial_{t}...$) have been replaced by total derivatives
($\mathrm{d}_{t}...$). After time integration between $t=0$ and $t>0$, we
obtain a series of problems for $\bm{u}_{2}^{(a)}$, $\bm{u}_{2}^{(b)}$,
$\bm{u}_{2}^{(c)}$ and $\bm{u}_{2}^{(d)}$:

$\displaystyle\Phi(t)\bm{u}_{2}^{(a)}(t)=\Phi(0)\bm{u}_{2}^{(a)}(0)=\alpha\bm{u}_{o},$

since $\Phi(0)=I$;

$\displaystyle\Phi(t)\bm{u}_{2}^{(b)}(t)=-\int_{0}^{t}H(s)e^{-Ls}C[\bm{l}(s),\bm{l}(s)]\mathrm{d}s=e^{-Lt}\tilde{\bm{u}}_{2}(t),$

where

$\frac{\mathrm{d}\tilde{\bm{u}}_{2}}{\mathrm{d}t}=L\tilde{\bm{u}}_{2}-C(\bm{l},\bm{l}),\quad\tilde{\bm{u}}_{2}(0)=\bm{0},$
(48)

and where we used that the general solution of
$\mathrm{d}_{t}\bm{x}=L\bm{x}+\bm{F}$ is
$\bm{x}(t)=e^{Lt}[\bm{x}(0)+\int_{0}^{t}e^{-Ls}\bm{F}(s)\mathrm{d}s]$;

$\displaystyle\Phi(t)\bm{u}_{2}^{(c)}(t)$
$\displaystyle=-\int_{0}^{t}H(s)e^{-Ls}\bm{l}(s)\mathrm{d}s$
$\displaystyle=-\int_{0}^{t}H(s)\epsilon_{o}\bm{u}_{o}\mathrm{d}s=-\epsilon_{o}\bm{u}_{o}t,$

since $e^{-Lt}\bm{l}(t)=\epsilon_{o}\bm{u}_{o}$ holds by construction; and

$\displaystyle\Phi(t)\bm{u}_{2}^{(d)}(t)=-\left[H(t)P(t)\bm{l}(t)-H(0)P(0)\bm{l}(0)\right]=-\bm{u}_{o},$

since, by construction, $H(t)P(t)\bm{l}(t)=H(t)\bm{u}_{o}$. Note that the
presence of the Heaviside distribution inside the integral is unimportant.
Eventually,

$\displaystyle\Phi\bm{u}_{2}$
$\displaystyle=\alpha\bm{u}_{o}+A^{2}e^{-Lt}\tilde{\bm{u}}_{2}-\frac{\mathrm{d}A}{\mathrm{d}T}\epsilon_{o}\bm{u}_{o}t-A\bm{u}_{o},$
(49)

for $t>0$.

##  Appendix H Transient growth: derivation of the amplitude equation.

The following relation has been derived:

$\Phi\bm{u}_{2}=\alpha\bm{u}_{o}+A^{2}e^{-Lt}\tilde{\bm{u}}_{2}-\frac{\mathrm{d}A}{\mathrm{d}T}\epsilon_{o}\bm{u}_{o}t-A\bm{u}_{o},\quad
0<t.$ (50)

Invoking the Fredholm alternative, (50) admits a solution if and only if the
right-hand side is orthogonal to $\bm{b}(t)$ for all $t>0$. This leads to:

$\displaystyle\left\langle\bm{u}_{o},\bm{b}(t)\right\rangle\left(\alpha-A\right)+A^{2}\left\langle
e^{-Lt}\tilde{\bm{u}}_{2}(t),\bm{b}(t)\right\rangle$
$\displaystyle-\epsilon_{o}t\frac{\mathrm{d}A}{\mathrm{d}T}\left\langle\bm{u}_{o},\bm{b}(t)\right\rangle=0,\quad
0<t.$ (51)

Dividing (51) by $\left\langle\bm{u}_{o},\bm{b}(t)\right\rangle$ leads to

$\begin{split}\left(\alpha-A\right)+A^{2}\mu_{2}(t)-\epsilon_{o}t\frac{\mathrm{d}A}{\mathrm{d}T}=0,\quad
0<t,\end{split}$ (52)

where

$\mu_{2}(t)=\frac{\left\langle
e^{-Lt}\tilde{\bm{u}}_{2}(t),\bm{b}(t)\right\rangle}{\left\langle\bm{u}_{o},\bm{b}(t)\right\rangle}.$
(53)

Equation (52) is re-expressed as $E(t,T)=0$ for $0<t$, where
$E(t,T)=\left(\alpha-A\right)+A^{2}\mu_{2}(t)-\epsilon_{o}t\mathrm{d}_{T}A$.
Since $\int_{t\rightarrow
0}^{t}\partial_{s}E(s,T)\mathrm{d}s=E(t,T)-E(t\rightarrow 0,T)$, solving
$E(t,T)=0$ is equivalent to solving $\partial_{t}E(t,T)=0$ for $0<t$ subject
to $E(t\rightarrow 0,T)=0$. Thereby, the partial derivative of (52) with
respect to the short time scale $t$ is taken, leading to

$\begin{split}A^{2}\frac{\mathrm{d}\mu_{2}(t)}{\mathrm{d}t}-\epsilon_{o}\frac{\mathrm{d}A}{\mathrm{d}T}=0,\quad
0<t,\end{split}$ (54)

where we have used that $\partial_{t}A=0$ by construction since $A=A(T)$ does
not explicitly depend on $t$. Furthermore, the relation (54) is subject to
$E(t\rightarrow 0,T)=\lim_{t\rightarrow 0}(\alpha-A)=0$ where we have used
that $\tilde{\bm{u}}_{2}(t\rightarrow 0)=\tilde{\bm{u}}_{2}(0)=0$. To be
meaningful, equation (54) and its initial condition must be re-written solely
in terms of $t$, which is done by evaluating $T$ along $T=\epsilon_{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,
$d_{t}A=\partial_{t}A+\epsilon_{o}\partial_{T}A=\epsilon_{o}d_{T}A$, such that
the final amplitude equation reads

$\frac{\mathrm{d}A}{\mathrm{d}t}=A^{2}\frac{\mathrm{d}\mu_{2}(t)}{\mathrm{d}t},\quad\text{with}\quad
A(0)=\alpha,$ (55)

as $\lim_{t\rightarrow 0}\left(\alpha-A(\epsilon_{o}t)\right)=0$ implies
$A(t\rightarrow 0)=\alpha$, and the amplitude $A$ is extended by continuity in
$t=0$ so as to eventually impose $A(0)=\alpha$. Note that the evaluation in
$T=\epsilon_{o}t$ and the passage to the total derivative would lead to
indeterminacy in its solution if performed directly in (52), since that
equation is not subject to any initial condition. Indeed, at linear level for
instance, it would yield $\mathrm{d}_{t}A=(\alpha-A)/t$, which admits the
family of solutions $A(t)=\alpha+Ct$, with $C$ an undetermined constant.

##  Appendix I Transient growth: orthogonality properties and higher-order
corrections of the amplitude equation.

After satisfaction of the Fredholm alternative, the equation at order
$O(\epsilon_{o}^{2})$ is

$\begin{split}\Phi\bm{u}_{2}=A^{2}\left(e^{-Lt}\tilde{\bm{u}}_{2}(t)-\mu_{2}(t)\bm{u}_{o}\right),\quad
0<t.\end{split}$ (56)

The orthogonality of the right-hand side (RHS) with $\bm{b}(t)$ is ensured by
construction of $\mu_{2}(t)$ in (53). The general solution reads
$\bm{u}_{2}=\bm{u}_{2}^{(\perp\bm{l})}+A_{2}(T)\bm{l}(t)$. The particular
solution of the system, $\bm{u}_{2}^{(\perp\bm{l})}$, is obtained by solving
(56) after replacing $\Phi(t)$ by $e^{-Lt}$. Indeed, such
$\bm{u}_{2}^{(\perp\bm{l})}$ must be orthogonal to $\bm{l}(t)$, since
$0=\left\langle\mbox{RHS},\bm{b}(t)\right\rangle$ = $\left\langle
e^{Lt}(\mbox{RHS}),\bm{l}(t)\right\rangle$ =
$\left\langle\bm{u}_{2}^{(\perp\bm{l})},\bm{l}(t)\right\rangle$. On the other
hand, the term $A_{2}(T)\bm{l}(t)$ constitutes the homogeneous part of the
solution, where $A_{2}(T)$ is a scalar amplitude. It can be kept in further
calculations, provided it is included in the final amplitude for $\bm{l}(t)$,
which would then become
$\epsilon_{o}A+\epsilon_{o}^{2}A_{2}+O(\epsilon_{o}^{3})$. Instead, and
without loss of generality [51], we propose to set $A_{2}=0$. In particular,
this implies that the term $\partial_{t}(P\bm{u}_{2})$ that appears at order
$O(\epsilon_{o}^{3})$ actually vanishes since
$P\bm{u}_{2}=P\bm{u}_{2}^{(\perp\bm{l})}=\bm{0}$. If this is performed at each
order $j\geq 3$, all the terms $\partial_{t}(P\bm{u}_{j})$ vanish. In this
way, the “retroaction” forcing due to the operator perturbation only appears
at order $O(\epsilon_{o}^{2})$.

##  Appendix J Transient growth: higher-order amplitude equation for spatially
periodic flows.

Deriving a higher-order amplitude equation for transient growth first requires
introducing a very long time scale $\tau=\epsilon_{o}^{2}t$, such that
$A=A(T,\tau)$. As mentioned in the main text, the $z$-periodicity of the
initial condition implies that the coefficient $\mu_{2}(t)$ defined in (53) is
null. Thus, at order $O(\epsilon_{o}^{2})$, the amplitude equation is simply
$\partial_{T}A=0$, subject to $\lim_{t\rightarrow 0}(\alpha-A)=0$, and
$\Phi\bm{u}_{2}=A^{2}e^{-Lt}\tilde{\bm{u}}_{2}$ for $t>0$. The particular
solution is $\bm{u}_{2}=A^{2}\tilde{\bm{u}}_{2}$ (see S9), and the homogeneous
solution is not needed. At order $O(\epsilon_{o}^{3})$, we gather

$\displaystyle\partial_{t}(\Phi\bm{u}_{3})=$
$\displaystyle-2A^{3}He^{-Lt}C(\bm{l},\tilde{\bm{u}}_{2})-e^{-Lt}\partial_{T}\bm{u}_{2}-e^{-Lt}\partial_{\tau}\bm{u}_{1}-\partial_{t}(P\bm{u}_{2})$
(57) $\displaystyle=$
$\displaystyle-2A^{3}He^{-Lt}C(\bm{l},\tilde{\bm{u}}_{2})-H\epsilon_{o}\bm{u}_{o}\partial_{\tau}A,$
(58)

since $\partial_{T}\bm{u}_{2}=2A(\partial_{T}A)\tilde{\bm{u}}_{2}=\bm{0}$,
$P\bm{u}_{2}=\bm{0}$, and
$e^{-Lt}\partial_{\tau}\bm{u}_{1}=He^{-Lt}\bm{l}\partial_{\tau}A=H\epsilon_{o}\bm{u}_{o}\partial_{\tau}A$.
This equation is subject to $\bm{u}_{3}(0)=\bm{0}$. The particular solution of
(58) yields

$\bm{u}_{3}(t,T,\tau)=A(T,\tau)^{3}\bm{u}_{3}^{(a)}(t)+\frac{\partial
A(T,\tau)}{\partial\tau}\bm{u}_{3}^{(b)}(t),$ (59)

where

$\displaystyle\mathrm{d}_{t}\left[\Phi\bm{u}_{3}^{(a)}\right]=-2He^{-Lt}C(\bm{l},\tilde{\bm{u}}_{2}),\quad\text{and}\quad$
$\displaystyle\mathrm{d}_{t}\left[\Phi\bm{u}_{3}^{(b)}\right]=-H\epsilon_{o}\bm{u}_{o}$

subject to the initial conditions
$\bm{u}_{3}^{(a)}(0)=\bm{u}_{3}^{(b)}(0)=\bm{0}$. After time integration
between $0$ and $t>0$, we obtain

$\displaystyle\Phi\bm{u}_{3}$
$\displaystyle=A^{3}e^{-Lt}\tilde{\bm{u}}_{3}-\epsilon_{o}\bm{u}_{o}t\partial_{\tau}A,\quad
0<t,$ (60)

where $\tilde{\bm{u}}_{3}$ is solution of

$\frac{\mathrm{d}\tilde{\bm{u}}_{3}}{\mathrm{d}t}=L\tilde{\bm{u}}_{3}-2C(\bm{l},\tilde{\bm{u}}_{2}),\quad\tilde{\bm{u}}_{3}(0)=\bm{0}.$

Canceling the projection of the RHS of (60) on $\bm{b}(t)$, dividing the
ensuing relation by $\left\langle\bm{b}(t),\bm{u}_{o}\right\rangle$, and
taking the partial derivative with respect to $t$ leads to

$A^{3}\frac{\mathrm{d}\mu_{3}(t)}{\mathrm{d}t}=\epsilon_{o}\frac{\partial
A}{\partial\tau},\quad\quad 0<t,$ (61)

where

$\mu_{3}(t)=\frac{\left\langle
e^{-Lt}\tilde{\bm{u}}_{3}(t),\bm{b}(t)\right\rangle}{\left\langle\bm{u}_{o},\bm{b}(t)\right\rangle}.$

To be meaningful, equation (61) and its initial condition must be re-written
solely in terms of $t$, which is done by evaluating $T=\epsilon_{o}t$ and
$\tau=\epsilon_{o}^{2}t$. The total derivative of $A$, denoted $d_{t}A$, is
now needed, as it takes into account the implicit dependence of $A$; it reads
$\mathrm{d}_{t}A=\partial_{t}A+\epsilon_{o}\partial_{T}A+\epsilon_{o}^{2}\partial_{\tau}A=\epsilon_{o}^{2}\partial_{\tau}A$,
such that

$\epsilon_{o}A^{3}\frac{\mathrm{d}\mu_{3}(t)}{\mathrm{d}t}=\frac{\mathrm{d}A}{\mathrm{d}t},\quad
0<t,$

subject to $\lim_{t\rightarrow
0}(\alpha-A(\epsilon_{o}t,\epsilon_{o}^{2}t))=0$ so $A(t\rightarrow 0)=\alpha$
and the amplitude $A$ is extended by continuity in $t=0$ so as to impose
$A(0)=\alpha$.

##  Appendix K Transient growth: sensitivity.

We consider a linear system $\partial_{t}\bm{u}=L\bm{u}$ subject to the
initial condition $\bm{u}(0)$ with $||\bm{u}(0)||=1$. The linear transient
gain at $t=t_{o}$ writes $G_{o}=||\bm{u}(t_{o})||$. A variational method is
used to derive the variation of the optimal transient gain induced by a small
perturbation $\delta L$ of operator. Let us introduce the Lagrangian

$\mathcal{L}=G_{o}^{2}-\int_{0}^{t_{o}}\left\langle\partial_{t}\bm{u}-L\bm{u},\bm{u}^{\dagger}\right\rangle\mathrm{d}t-\beta\left(1-||\bm{u}(0)||^{2}\right),$

where the Lagrange multipliers $\bm{u}^{\dagger}$ and $\beta$ enforce the
constraints on the state equation and on the norm of the initial condition,
respectively. Imposing
$\left\langle\partial_{\bm{u}}\mathcal{L},\delta\bm{u}\right\rangle=0$ for all
$\delta\bm{u}$ leads to the adjoint equation
$\partial_{t}\bm{u}^{\dagger}=-L^{\dagger}\bm{u}^{\dagger}$, to be integrated
backward in time from the terminal condition
$\bm{u}^{\dagger}(t_{o})=2\bm{u}(t_{o})$. Eventually, the gain variation
induced by $\delta L$ is

$\delta(G_{o}^{2})=\left\langle\frac{\partial\mathcal{L}}{\partial L},\delta
L\right\rangle=\int_{0}^{t_{o}}\left\langle(\delta
L)\bm{u},\bm{u}^{\dagger}\right\rangle\mathrm{d}t.$ (62)

On the other hand, we derived in the main text for the WNNt model:

$\displaystyle a(t)$
$\displaystyle=\frac{U_{0}}{\epsilon_{o}}\left[1-\left(\frac{U_{0}}{\epsilon_{o}}\right)^{2}\frac{2\mu_{3}(t)}{\epsilon_{o}}\right]^{-1/2}.$
(63)

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_{o}^{2}=(1/\epsilon_{o})^{2}$, such that

$\displaystyle\frac{1}{G(t_{o})^{2}}-\frac{1}{G_{o}^{2}}$
$\displaystyle=-U_{0}^{2}\frac{2\mu_{3}(t_{o})}{\epsilon_{o}}.$

We are interested in small variations around $G_{o}^{2}$, thus we write
$G(t_{o})^{2}=G_{o}^{2}+\delta G_{o}^{2}$ with $|\delta
G_{o}^{2}/G_{o}^{2}|\ll 1$. In this manner,
$1/G_{o}^{2}-1/G(t_{o})^{2}=\delta(G_{o}^{2})/G_{o}^{4}+h.o.t$, eventually
leading to

$\displaystyle\delta(G_{o}^{2})$
$\displaystyle=2\mu_{3}(t_{o})\frac{U_{0}^{2}}{\epsilon_{o}^{5}}.$

In addition,

$\displaystyle\mu_{3}(t_{o})=\frac{\left\langle
e^{-Lt_{o}}\tilde{\bm{u}}_{3}(t_{o}),\bm{b}(t_{o})\right\rangle}{\left\langle\bm{u}_{o},\bm{b}(t_{o})\right\rangle},$

with
$\tilde{\bm{u}}_{3}(t)=-e^{Lt}\int_{0}^{t}2e^{-Ls}C[\tilde{\bm{u}}_{2}(s),\bm{l}(s)]\mathrm{d}s$.
Therefore

$\displaystyle\mu_{3}(t_{o})$
$\displaystyle=-\frac{\left\langle\int_{0}^{t_{o}}2e^{-Ls}C[\tilde{\bm{u}}_{2}(s),\bm{l}(s)]\mathrm{d}s,\bm{b}(t_{o})\right\rangle}{\left\langle\bm{u}_{o},\bm{b}(t_{o})\right\rangle}$
$\displaystyle=-\frac{\int_{0}^{t_{o}}\left\langle
2e^{-Ls}C[\tilde{\bm{u}}_{2}(s),\bm{l}(s)],\bm{b}(t_{o})\right\rangle\mathrm{d}s}{\left\langle\bm{u}_{o},\bm{b}(t_{o})\right\rangle}.$

By definition, $\bm{b}(t_{o})=(e^{Lt_{o}})^{\dagger}\bm{l}(t_{o})$ and
$\left\langle\bm{u}_{o},\bm{b}(t_{o})\right\rangle=\left\langle
e^{Lt_{o}}\bm{u}_{o},\bm{l}(t_{o})\right\rangle=1/\epsilon_{o}$. In addition,
$\left(e^{-Ls}\right)^{\dagger}(e^{Lt_{o}})^{\dagger}=\left(e^{Lt_{o}}e^{-Ls}\right)^{\dagger}=\left(e^{-L(s-t_{o})}\right)^{\dagger}$,
such that

$\displaystyle\mu_{3}(t_{o})$
$\displaystyle=-\epsilon_{o}\int_{0}^{t_{o}}\left\langle
2C[\tilde{\bm{u}}_{2}(s),\bm{l}(s)],\left(e^{-L(s-t_{o})}\right)^{\dagger}\bm{l}(t_{o})\right\rangle\mathrm{d}s.$

In terms of our previous notations, we have the direct correspondence
$\bm{u}(t)=\bm{l}(t)/\epsilon_{o}$ and
$\bm{u}^{\dagger}(s)=\left(e^{-L(s-t_{o})}\right)^{\dagger}2\bm{u}(t_{o})$, so
we can write

$\displaystyle\delta(G_{o}^{2})$
$\displaystyle=\left(\frac{2U_{0}^{2}}{\epsilon_{o}^{5}}\right)\left(-\frac{\epsilon_{o}^{3}}{2}\int_{0}^{t_{o}}\left\langle
2C[\tilde{\bm{u}}_{2}(s),\bm{u}(s)],\left(e^{-L(s-t_{o})}\right)^{\dagger}2\bm{u}(t_{o})\right\rangle\mathrm{d}s\right)$
$\displaystyle=-\int_{0}^{t_{o}}\left\langle
2C\left[(U_{0}/\epsilon_{o})^{2}\tilde{\bm{u}}_{2}(s),\bm{u}(s)\right],\bm{u}^{\dagger}(s)\right\rangle\mathrm{d}s.$

The sensitivity relation (62) is immediately recognized, where $\delta L$ is
here induced by the addition of $(U_{0}/\epsilon_{o})^{2}\tilde{\bm{u}}_{2}$
to the base flow. Indeed, $U_{0}/\epsilon_{o}=a_{lin}$ is the linear solution
of (63) corresponding the limit $U_{0}\rightarrow 0$, such that the flow field
is described in this limit by
$\bm{U}(t)=\bm{U}_{e}+a_{lin}\bm{l}(t)+a_{lin}^{2}\tilde{\bm{u}}_{2}(t)+O(\epsilon_{o}^{3})$.

##  Appendix L Transient growth: formal definition of the DNS gains $G_{tot}$
and $G_{l}$.

We recall the asymptotic expansion of the solution
$\bm{U}=\bm{U}_{e}+a\bm{l}+\epsilon_{o}^{2}\bm{u}_{2}+O(\epsilon_{o}^{3})$ for
$0<t$, with $\left\langle\bm{l},\bm{u}_{i}\right\rangle=0$ for all $i\geq 2$
as a consequence of the Fredholm alternative (see Appendix I). This gives:

$a(t)=\frac{\left\langle\bm{U}-\bm{U}_{e},\bm{l}(t)\right\rangle}{||\bm{l}(t)||^{2}}.$

Thus, from the knowledge of $\bm{U}-\bm{U}_{e}$ computed by DNS, the gain that
should be compared with the weakly nonlinear prediction $a||\bm{l}(t)||/U_{0}$
is

$G_{\bm{l}}(t)=\frac{\left\langle\bm{U}(t)-\bm{U}_{e},\bm{l}(t)\right\rangle}{U_{0}||\bm{l}(t)||},$

which, evaluated at $t=t_{o}$ (with $\bm{l}(t_{o})=\bm{v}_{o}$ is of unit
norm), gives
$G_{\bm{l}}(t_{o})=\left\langle\bm{U}(t_{o})-\bm{U}_{e},\bm{v}_{o}\right\rangle/U_{0}$.
By contrast, one can also consider a gain that includes all higher-order
fields, namely:

$G_{tot}(t)=\frac{||\bm{U}(t)-\bm{U}_{e}||}{U_{0}}.$

## References

  * Strogatz [2015] S. H. Strogatz,  _Nonlinear Dynamics and Chaos_ , 2nd ed. (CRC Press, 2015).
  * Trefethen  _et al._ [1993] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll, Hydrodynamic stability without eigenvalues, Science 261, 578 (1993).
  * Waleffe [1995] F. Waleffe, Transition in shear flows. nonlinear normality versus non‐normal linearity, Phys. Fluids 7, 3060 (1995).
  * Landau and Lifshitz [1987] L. Landau and E. Lifshitz,  _Fluid Mechanics, 2nd edn_ (Oxford: Butterworth-Heinemann, 1987).
  * Gor’kov [1957] L. P. Gor’kov, Stationary convection in a plane liquid layer near the critical heat transfer point, Zh. Eksp. Teor. Fiz. 6, 311–15 (1957).
  * Malkus and Veronis [1958] W. V. R. Malkus and G. Veronis, Finite amplitude cellular convection, J. Fluid Mech. 4, 225–60 (1958).
  * Stuart [1958] J. T. Stuart, On the nonlinear mechanics of hydrodynamic stability, J. Fluid Mech. 4, 1 (1958).
  * Stuart [1960] J. T. Stuart, On the nonlinear mechanics of wave disturbances in stable and unstable parallel flows, part 1. the basic behavior in plane Poiseuille flow, J. Fluid Mech. 9, 353 (1960).
  * Watson [1960] J. Watson, On the nonlinear mechanics of wave disturbances in stable and unstable parallel flows, part 2. the development of a solution for plane Poiseuille and for plane Couette flow, J. Fluid Mech. 9, 371–89 (1960).
  * Sipp and Lebedev [2007] D. Sipp and A. Lebedev, Global stability of base and mean flows : a general approach and its applications to cylinder and open cavity flows, J. Fluid Mech. 593, 333 (2007).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).
  * Fauve [1998] S. Fauve, Pattern forming instabilities, in  _Hydrodynamics and Nonlinear Instabilities_ , edited by C. Godréche and P. Manneville (Cambridge University Press, 1998) pp. 387–492.
  * Crawford and Knobloch [1991] J. D. Crawford and E. Knobloch, Symmetry and symmetry-breaking bifurcations in fluid dynamics, Annu. Rev. Fluid Mech. 23, 341 (1991).
  * Manneville [2004] P. Manneville,  _Instabilities, Chaos and Turbulence_ (Imperial College Press, 2004).
  * Guckenheimer and Holmes [1983] J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of vector fields_ (Springer, New York, NY, 1983).
  * Haragus and Iooss [2011] M. Haragus and G. Iooss,  _Local bifurcations, center manifolds, and normal forms in infinite-dimensional dynamical systems_ (Springer-Verlag London, 2011).
  * Trefethen and Embree [2005] L. N. Trefethen and M. Embree,  _Spectra and pseudo-spectra_ (Princeton University Press, 2005).
  * Landau [1976] H. J. Landau, Loss in unstable resonators, J. Opt. Soc. Am. 66, 525 (1976).
  * Landau [1977] H. J. Landau, The notion of approximate eigenvalues applied to an integral equation of laser theory, Quart. Appl. Math. 35, 165 (1977).
  * Petermann [1979] K. Petermann, Calculated spontaneous emission factor for double-heterostructure injection lasers with gain-induced waveguiding, IEEE J. Quant. Elect. 15, 566 (1979).
  * Jaramillo  _et al._ [2021] J. L. Jaramillo, R. P. Macedo, and L. A. Sheikh, Pseudospectrum and black hole quasinormal mode instability, Phys. Rev. X 11, 031003 (2021).
  * Asllani  _et al._ [2018] M. Asllani, R. Lambiotte, and T. Carletti, Structure and dynamical behavior of non-normal networks, Science Advances 4, eaau9403 (2018).
  * Asllani and Carletti [2018] M. Asllani and T. Carletti, Topological resilience in non-normal networked systems, Phys. Rev. E 97, 042302 (2018).
  * Gustavsson [1991] L. H. Gustavsson, Energy growth of three-dimensional disturbances in plane Poiseuille flow, J. Fluid Mech. 224, 241–260 (1991).
  * Butler and Farrell [1992] K. M. Butler and B. F. Farrell, Three‐dimensional optimal perturbations in viscous shear flow, Phys. Fluids A 4, 1637 (1992).
  * Reddy and Henningson [1993] S. C. Reddy and D. S. Henningson, Energy growth in viscous channel flows, J. Fluid Mech. 252, 209–238 (1993).
  * Farrell and Ioannou [1993] B. F. Farrell and P. J. Ioannou, Stochastic forcing of the linearized Navier–Stokes equations, Phys. Fluids A 5, 2600 (1993).
  * Schmid and Henningson [2001] P. Schmid and D. Henningson,  _Stability and transition in shear flows_ (Springer, 2001).
  * Schmid and Henningson [1994] P. J. Schmid and D. S. Henningson, Optimal energy density growth in Hagen–Poiseuille flow, J. Fluid Mech. 277, 197–225 (1994).
  * Corbett and Bottaro [2000] P. Corbett and A. Bottaro, Optimal perturbations for boundary layers subject to stream-wise pressure gradient, Phys. Fluid 12, 120 (2000).
  * Cossu and Chomaz [1997] C. Cossu and J.-M. Chomaz, Global measures of local convective instabilities, Phys. Rev. Lett. 78, 4387 (1997).
  * Ehrenstein and Gallaire [2005] U. Ehrenstein and F. Gallaire, On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer, J. Fluid Mech. 536, 209 (2005).
  * Åkervik  _et al._ [2008] E. Åkervik, U. Ehrenstein, F. Gallaire, and D. S. Henningson, Global two-dimensional stability measures of the flat plate boundary-layer flow, Eur. J. Mech. B/Fluids 27, 501 (2008).
  * Ehrenstein and Gallaire [2008] U. Ehrenstein and F. Gallaire, Two-dimensional global low-frequency oscillations in a separating boundary-layer flow, J. Fluid Mech. 614, 315 (2008).
  * Alizard  _et al._ [2009] F. Alizard, S. Cherubini, and J. Robinet, Sensitivity and optimal forcing response in separated boundary layer flows, Phys. Fluids 21, 064108 (2009).
  * Monokrousos  _et al._ [2010] A. Monokrousos, E. Åkervik, L. Brandt, and D. S. Henningson, Global three-dimensional optimal disturbances in the Blasius boundary-layer flow using time-steppers, J. Fluid Mech. 650, 181 (2010).
  * Garnaud  _et al._ [2013a] X. Garnaud, L. Lesshafft, P. J. Schmid, and P. Huerre, Modal and transient dynamics of jet flows, Phys. Fluids 25 (2013a).
  * Garnaud  _et al._ [2013b] X. Garnaud, L. Lesshafft, P. J. Schmid, and P. Huerre, The preferred mode of incompressible jets: linear frequency response analysis, J. Fluid Mech. 716, 189–202 (2013b).
  * Blackburn  _et al._ [2008] H. M. Blackburn, D. Barkley, and S. J. Sherwin, Convective instability and transient growth in flow over a backward-facing step, J. Fluid Mech. 603, 271–304 (2008).
  * Boujo and Gallaire [2015] E. Boujo and F. Gallaire, Sensitivity and open-loop control of stochastic response in a noise amplifier flow: the backward-facing step, Journal of Fluid Mechanics 762, 361–392 (2015).
  * Chomaz [2005] J.-M. Chomaz, Global instabilities in spatially developing flows: non-normality and nonlinearity, Annu. Rev. Fluid Mech. 37, 357 (2005).
  * Schmid [2007] P. J. Schmid, Nonmodal stability theory, Annu. Rev. Fluid Mech. 39, 129 (2007).
  * Baggett and Trefethen [1997] J. S. Baggett and L. N. Trefethen, Low-dimensional models of subcritical transition to turbulence, Phys. Fluids 9, 1043 (1997).
  * Hof  _et al._ [2004] B. Hof, C. W. H. van Doorne, J. Westerweel, F. T. M. Nieuwstadt, H. Faisst, B. Eckhardt, H. Wedin, R. R. Kerswell, and F. Waleffe, Experimental observation of nonlinear traveling waves in turbulent pipe flow, Science 305, 1594 (2004).
  * Schneider  _et al._ [2010] T. M. Schneider, J. F. Gibson, and J. Burke, Snakes and ladders: localized solutions of plane Couette flow, Phys. Rev. Lett. 104, 104501 (2010).
  * Cherubini  _et al._ [2010] S. Cherubini, P. De Palma, J.-C. Robinet, and A. Bottaro, Rapid path to transition via nonlinear localized optimal perturbations in a boundary-layer flow, Phys. Rev. E 82, 066302 (2010).
  * Cherubini  _et al._ [2011] S. Cherubini, P. De Palma, J.-C. Robinet, and A. Bottaro, The minimal seed of turbulent transition in the boundary layer, J. Fluid Mech. 689, 221–253 (2011).
  * Pringle and Kerswell [2010] C. C. T. Pringle and R. R. Kerswell, Using nonlinear transient growth to construct the minimal seed for shear flow turbulence, Phys. Rev. Lett. 105, 154502 (2010).
  * Kerswell [2018] R. R. Kerswell, Nonlinear nonmodal stability theory, Annu. Rev. Fluid Mech. 50, 319 (2018).
  * Sipp [2012] D. Sipp, Open-loop control of cavity oscillations with harmonic forcings, J. Fluid Mech. 708, 439 (2012).
  * Fujimura [1991] K. Fujimura, Methods of centre manifold and multiple scales in the theory of weakly nonlinear stability for fluid motions, Proc. R. Soc. Lond. A 434, 719 (1991).
  * Brandt  _et al._ [2011] L. Brandt, D. Sipp, J. O. Pralits, and O. Marquet, Effect of base-flow variation in noise amplifiers: the flat-plate boundary layer, J. Fluid Mech. 687, 503 (2011).
  * Mantic-Lugo and Gallaire [2016a] V. Mantic-Lugo and F. Gallaire, Self-consistent model for the saturation mechanism of the response to harmonic forcing in the backward-facing step flow, J. Fluid Mech. 793, 777–97 (2016a).
  * Neubert and Caswell [1997] M. G. Neubert and H. Caswell, Alternatives to resilience for measuring the responses of ecological systems to perturbations, Ecology 78, 653 (1997).
  * Mantic-Lugo and Gallaire [2016b] V. Mantic-Lugo and F. Gallaire, Saturation of the response to stochastic forcing in two-dimensional backward-facing step flow: a self-consistent approximation, Phys. Fluids 1 (2016b).
  * Benitez and Bermudez [2011] M. Benitez and A. Bermudez, A second order characteristics finite element scheme for natural convection problems, J. Comput. Appl. Maths 235, 3270 (2011).