1. Introduction
Pressure-driven two-dimensional (2-D) plane Poiseuille flow (or channel flow) is a widely studied canonical flow in fluid mechanics (Orszag & Kells Reference Orszag and Kells1980). It consists of an incompressible fluid with kinematic viscosity
$\nu$
governed by the Navier–Stokes equations flowing between two parallel walls, with half-channel height
$h_c$
, and driven by a constant streamwise pressure gradient producing a parabolic laminar velocity profile with centreline velocity
$U_c$
, so that the Reynolds number is defined by
$\textit{Re}=U_c h_c/\nu$
. Despite its fundamental status, there is a lot that remains unknown about its dynamics.
To study how this flow transitions, we analyse the stability of its steady laminar state, i.e. whether or not initial velocity perturbations about this state vanish in time. Its linear stability limit has been computed to be
$\textit{Re}_L\approx 5772$
(Orszag Reference Orszag1971), meaning that just above
$\textit{Re}_L$
, there exists an initial (small) unstable perturbation that will not decay in time. Having said that, this number is not as informative as it seems, because transition to turbulence in this flow is subcritical, being observed well below
$\textit{Re}_L$
. Indeed, travelling waves (i.e. self-sustaining non-laminar periodic exact solutions of the flow) have been accurately computed as early as
$\textit{Re}\approx 2939$
(Casas & Jorba Reference Casas and Jorba2012). Thus
$\textit{Re}\approx 2939$
constitutes a clear upper bound to the global stability limit of the flow,
$\textit{Re}_G$
, defined as the largest Reynolds number such that every initial perturbation (no matter the magnitude) can be proved to vanish in time, resulting in the global asymptotic (nonlinear) stability of its laminar flow (Schmid & Henningson Reference Schmid and Henningson2001). Meanwhile, the only known lower bound to
$\textit{Re}_G$
is given by its energy stability limit, i.e. the largest Reynolds number such that the kinetic energy of any perturbation monotonically decays in time. It was originally computed by Orr (Reference Orr1907) to be
$\textit{Re}_E\approx 88$
, with a more accurate update being
$\textit{Re}_E\approx 87.59$
, attained by a critical periodic perturbation of non-dimensional length
$L_E= {2\pi }/{\alpha _E}\approx 2.99$
. (The three-dimensional doubly periodic pressure-driven version of plane Poiseuille flow shares the same linear stability limit, but has travelling waves detected at
$\textit{Re}\approx 331$
(Zammert & Eckhardt Reference Zammert and Eckhardt2016), and its energy stability limit is
$\textit{Re}_E\approx 49.6$
(Joseph & Carmi Reference Joseph and Carmi1969).)
Thus much remains unknown about the flow’s global stability for
$87.59\leq \textit{Re} \leq 2939$
. The purpose of this paper is to establish the global (nonlinear) stability of 2-D plane Poiseuille flow beyond the energy stability limit, resulting in improved lower bounds on the global stability of the flow. These results, settled for channels up to a certain streamwise length, represent the first global stability certificates in over a century since Orr (Reference Orr1907) first computed
$\textit{Re}_E$
. Proving global stability of a fluid flow beyond
$\textit{Re}_E$
is usually a difficult endeavour that has historically been confined to seeking ‘quadratic’ Lyapunov functionals. Instead, we obtained these new bounds by constructing ‘quartic’ Lyapunov functionals using a computer, following the original ideas of Goulart & Chernyshenko (Reference Goulart and Chernyshenko2012), which were later implemented and refined by Huang et al. (Reference Huang, Chernyshenko, Goulart, Lasagna, Tutty and Fuentes2015) and Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022). The computational ‘SOS–Lyapunov framework’ searches for functionals that satisfy certain polynomial sum-of-squares (SOS) inequalities. These are recast as a semidefinite program (SDP), whose feasibility certifies that the flow is globally stable. (We remark that these Lyapunov functionals apply to the genuine nonlinear dynamics of the Navier–Stokes equations, whereas the results of Fraternale et al. (Reference Fraternale, Domenicale, Staffilani and Tordella2018), upon inspection, apply only to the linearised dynamics, because they linearised the so-called
$C^1$
nonlinear compatibility conditions described by Temam (Reference Temam2006) in the admissible set of velocities.)
The same general approach was used by Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022) to establish the global stability of 2-D plane Couette flow, but here, due to changing the flow, we had to use different numerical methods, based on finite elements, to solve the auxiliary partial differential equations necessary to accurately set up and apply the methodology. Moreover, a new expression was derived to more quickly and accurately compute some of the bounds involved in the method (see Appendix C). We note that compared to 2-D plane Couette flow, where travelling waves have not been found (Ehrenstein et al. Reference Ehrenstein, Nagata and Rincon2008), 2-D plane Poiseuille flow does have known upper bounds on
$\textit{Re}_G$
, which makes it a richer 2-D dynamical system to study. The Lyapunov functionals have as building blocks a small set of velocity ‘modes’, which we first need to identify. We make relevant observations for this specific flow, complemented with experimentation, to successfully construct several of these ‘mode sets’ (see § 3). Importantly, even though the computed expressions for Lyapunov functionals are often not particularly informative of the underlying physics, the selection of mode sets is insightful and is likely to be a key ingredient for more analytical approaches to proving global stability. Indeed, the smallest mode set that we found in this work was used in the recent work of Darrow et al. (Reference Darrow, Carlson and Goluskin2026) to analytically construct structured quartic Lyapunov functionals for this flow, though these are only valid up to certain values of
$\textit{Re}$
, which are below the best bounds showcased in this work for the larger mode sets.
2. Formulation and numerical implementation
For a fixed non-dimensional streamwise period
$L$
and Reynolds number
$\textit{Re}$
, our goal is to certify global asymptotic stability of the laminar plane Poiseuille flow with respect to 2-D periodic perturbations. We consider velocity and pressure perturbations,
$\boldsymbol{u}=(u,v)$
and
$p$
, about the non-dimensional laminar plane Poiseuille flow
$\boldsymbol{U}(y)=(1-y^2,0)$
and constant pressure gradient
$-\partial _x P= {2}/{\textit{Re}}$
in the domain
$\varOmega =(0,L)\times (-1,1)$
, with periodic boundary conditions in
$x$
, and no-slip conditions at the walls,
$y=\pm 1$
. According to the incompressible Navier–Stokes equations, the perturbations are governed by
\begin{align} \partial _t \boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} +(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} +(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{U} =-\boldsymbol{\nabla }\!p+\frac {1}{\textit{Re}}{\nabla} ^2\boldsymbol{u}, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0. \end{align}
To certify global stability, we use the SOS–Lyapunov framework developed by Goulart & Chernyshenko (Reference Goulart and Chernyshenko2012) and refined by Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022). Specifically, we seek a Lyapunov functional
$V(\boldsymbol{u})$
for the perturbation dynamics: if one can construct a functional that is positive away from the laminar state (i.e.
$\boldsymbol{u}=\boldsymbol{0}$
) and whose value monotonically decreases in time along every non-zero solution of (2.1), then the perturbation must decay to zero as
$t\to \infty$
, implying that the laminar state is globally asymptotically stable. Therefore, the existence of such a functional for a given
$\textit{Re}$
provides a certified lower bound on
$\textit{Re}_G$
at
$L$
.
2.1. The SOS–Lyapunov framework
We begin by decomposing the velocity field as
\begin{align} \boldsymbol{u}(x,y,t)= {{\sum}_{i=1}^m}a_i(t)\,\boldsymbol{u}_i(x,y)+\boldsymbol{v}(x,y,t), \end{align}
where
$\mathscr{U}=\{\boldsymbol{u}_1,\ldots ,\boldsymbol{u}_m\}$
is a finite orthonormal set of solenoidal modes, and
$\boldsymbol{v}$
is an infinite-dimensional ‘tail’ orthogonal to each
$\boldsymbol{u}_i$
. That is,
$\langle \boldsymbol{u}_i,\boldsymbol{u}_{\!j}\rangle = \delta _{ij}$
,
$\langle \boldsymbol v,\boldsymbol u_i\rangle =0$
and
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_i=0$
for all
$i,j=1,\ldots ,m$
, where
$\delta _{ij}$
is the Kronecker delta,
\begin{align} \langle \boldsymbol{u}, \boldsymbol{v}\rangle = \textstyle {\int _\varOmega }\, \boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{v}\, \mathrm d\varOmega \quad \text{and}\quad \|\boldsymbol u\|^2=\langle \boldsymbol u,\boldsymbol u \rangle . \end{align}
Thus the time dependence of
$\boldsymbol{u}$
is determined by the functions
$a_i$
and
$\boldsymbol{v}$
. Indeed, note that the time-dependent perturbation energy
$E(t)=\tfrac{1}{2}\|\boldsymbol{u}\|^2$
can be written as
\begin{align} E=\frac {1}{2}\big (\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{a}+q^2\big ),\quad \text{with}\quad q^2=\|\boldsymbol{v}\|^2\ \text{and}\ \boldsymbol{a}=(a_1,\ldots ,a_m)^{\rm T}. \end{align}
We consider quartic functionals of the form
\begin{align} V(\boldsymbol{a},q)=E(\boldsymbol{a},q)^2+P(\boldsymbol{a},q), \end{align}
where
$P$
is a cubic polynomial, with no constant or linear terms, and only even powers of
$q$
. (One could consider more general ansätze of the form
$V=E^d+P$
for
$d\gt 2$
and
$P$
of at most degree
$2d-1$
, but the resulting computational cost is significantly higher, so we focus on the
$d=2$
case.) We call
$V$
a Lyapunov functional if
\begin{align} V(\boldsymbol{a},q)\ge 0\ \text{and}\ \mathcal{L}V(\boldsymbol{a},q)\le 0, \quad \text{with equality if and only if }(\boldsymbol{a},q)=(\boldsymbol{0}, 0). \end{align}
Here,
$\mathcal LV(\boldsymbol a(t),q(t))=({\mathrm dV}/{\mathrm dt})(\boldsymbol a(t),q(t))$
denotes the Lie derivative of
$V$
along solutions of (2.1) written in the form (2.2). The existence of such a (radially unbounded) functional implies global asymptotic stability of the laminar state (Mironchenko & Wirth Reference Mironchenko and Wirth2019, Prop. 3.2). Importantly, these conditions are sufficient to establish the global stability of the infinite-dimensional system (2.1), and not merely that of a truncated finite-dimensional approximation. This is because the full perturbation velocity field
$\boldsymbol{u}$
is properly encoded in
$(\boldsymbol a,q)$
, and, when taking into account the dynamics,
$\mathcal LV$
includes the time-evolution of both
$\boldsymbol a(t)$
and
$q(t)$
.
For example, consider the energy
$E$
, which is positive definite. If additionally, for any possible perturbation governed by the Navier–Stokes equations (2.1), its energy satisfies
$ {\mathrm dE}/{\mathrm dt}\leq -cE$
for some
$c\gt 0$
(which occurs precisely when
$\textit{Re}\leq \textit{Re}_E$
), Grönwall’s inequality yields
$E(t)\leq E(0)\,{\rm e}^{-ct}$
, so any perturbation, regardless of its initial magnitude or form, eventually vanishes in time, resulting in a globally stable flow. In these cases, the energy is a Lyapunov functional. This can be reformulated as polynomials, since
$E(\boldsymbol{a},q)=\tfrac{1}{2}(\boldsymbol{a}\,\!\boldsymbol{\cdot }\,\!\boldsymbol{a}+q^2)$
is a quadratic polynomial in
$(\boldsymbol{a},q)$
, which will be a Lyapunov functional at the Reynolds numbers for which
$\mathcal{L}E(\boldsymbol{a},q)\lt 0$
at every non-zero
$(\boldsymbol{a},q)$
.
The idea is to find more general quartic Lyapunov functionals for
$\textit{Re}\gt \textit{Re}_E$
. In the ansatz (2.5), the leading term
$E^2$
guarantees radial unboundedness (meaning that
$V\to \infty$
when
$\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{a}+q^2\to \infty$
) and that
$\mathcal{L}V$
is of even degree, while the lower-degree polynomial
$P$
provides the flexibility to enforce the inequalities in (2.6) with the aid of a computer (see Goulart & Chernyshenko (Reference Goulart and Chernyshenko2012) and Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022) for more details).
Projecting (2.1) onto
$\operatorname {span}(\mathscr U)$
and its orthogonal complement yields equations for
${\mathrm d \boldsymbol{a}}/{\mathrm d t}$
and
${\mathrm d q^2}/{\mathrm d t}$
. These equations are not closed in
$(\boldsymbol a,q)$
, since some terms depend on the tail
$\boldsymbol v$
beyond its energy
$\tfrac{1}{2}q^2$
. Following Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022), such terms can nevertheless be conveniently bounded when
$\mathscr U$
is chosen as a finite set of ‘energy eigenmodes’ (described in § 2.2). This leads to explicit sufficient conditions for (2.5) to define a Lyapunov functional satisfying (2.6). Specifically, global stability of the flow is verified by finding a cubic polynomial
$P$
with no constant or linear terms, quartic polynomials
$r_i$
with no constant or linear terms, and quadratic polynomials
$s_i$
, all involving only even powers of
$q$
, such that
\begin{gather} E(\boldsymbol{a},q)^2 + P(\boldsymbol{a},q) - \varepsilon\, E(\boldsymbol{a},q) \ge 0, \nonumber \\-\left(\tilde{G}(\boldsymbol{a},q)+\sum\nolimits _{i=1}^{m} (r_i(\boldsymbol{a},q)+C_i q^2\, s_i(\boldsymbol{a},q))+\varepsilon E(\boldsymbol{a},q)\right) \ge 0, \nonumber \\\frac {\partial V}{\partial q^2} \ge 0,\quad s_i(\boldsymbol{a},q)-M_i(\boldsymbol{a},q) \ge 0, \quad s_i(\boldsymbol{a},q)+M_i(\boldsymbol{a},q) \ge 0, \quad i=1,\ldots ,m, \\ \begin{bmatrix} w_1 & w_2 \end{bmatrix} \begin{bmatrix} q^2\tilde {\boldsymbol{a}}^{\rm T}{{\unicode{x1D642}}}_i\tilde {\boldsymbol{a}}\,r_i(\boldsymbol{a},q) & q^2\tilde {\boldsymbol{a}}^{\rm T}{{\unicode{x1D642}}}_i\tilde {\boldsymbol{a}}\,M_i(\boldsymbol{a},q) \\[5pt] q^2\tilde {\boldsymbol{a}}^{\rm T}{{\unicode{x1D642}}}_i\tilde {\boldsymbol{a}}\,M_i(\boldsymbol{a},q) & r_i(\boldsymbol{a},q) \end{bmatrix} \begin{bmatrix} w_1 \\[5pt] w_2 \end{bmatrix} \ge 0,\quad i=1,\ldots ,m, \nonumber \end{gather}
for all
$(\boldsymbol{a},q,w_1,w_2)\in \mathbb{R}^{m+3}$
, where
$\tilde {\boldsymbol{a}} = (1, a_1,\ldots ,a_m)^{{\rm T}}$
and
$\varepsilon \gt 0$
. Here,
$\tilde{G}$
and
$M_i$
are polynomials,
${{\unicode{x1D642}}}_i$
are constant tensors, and
$C_i$
are constants defined in Appendix A.
Note that all the terms in (2.7) are polynomials. For a fixed choice of
$\textit{Re}$
,
$L$
and
$\mathscr U$
, to reliably solve (2.7) with a computer, the inequalities are strengthened to SOS polynomial constraints and recast as an SDP that can be solved numerically. If the SDP is found to be feasible, then its numerical solution includes the explicit coefficients of monomials in
$P(\boldsymbol{a},q)$
, with the resulting
$V(\boldsymbol{a},q)$
in (2.5) being a Lyapunov functional satisfying (2.6) (see (D1) for an explicit example).
2.2. Energy-rate eigenvalue problem
The modes in
$\mathscr U$
are chosen from the solutions to the 2-D second-order energy-rate eigenvalue problem in
$\varOmega$
(Doering & Gibbon Reference Doering and Gibbon1995, § 2.2),
\begin{align} \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u} -\frac 12\big (\boldsymbol{\nabla }\boldsymbol{U} + (\boldsymbol{\nabla }\boldsymbol U)^{\rm T}\big )\boldsymbol{u}-\boldsymbol{\nabla } \zeta = \lambda \boldsymbol u, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \quad \|\boldsymbol{u}\|^2 = 1\,. \end{align}
Here,
$\zeta$
is an
$x$
-periodic Lagrange multiplier enforcing incompressibility,
$\boldsymbol{u}$
satisfies periodic boundary conditions in
$x$
and no-slip at
$y=\pm 1$
, and
$\lambda =({\mathrm dE}/{\mathrm dt})\big |_{t=0}$
is the instantaneous energy growth rate of (2.1) when the ‘energy eigenmode’
$\boldsymbol{u}$
is chosen as the initial condition. Since the problem is 2-D, we represent the velocity in terms of a streamfunction
$\psi$
, writing
$\boldsymbol u=(\partial _y \psi ,-\partial _x\psi )$
, which automatically satisfies the incompressibility constraint in (2.8), while the third equation is enforced by a posteriori rescaling of the eigenmodes.
Letting
$\psi (x,y) = \sum _{n\in \mathbb Z} \varphi _n(y)\,{\rm e}^{{i}\alpha _n x}$
with
$\alpha _n ={2\pi n}/{L}$
, and taking the curl of the first equation in (2.8), yields a one-dimensional
$\mathbb{C}$
-valued fourth-order generalised Hermitian eigenvalue problem for each
$n$
. For
$\alpha \neq 0$
(dropping the
$n$
for convenience),
\begin{align} {\frac {1}{\textit{Re}}}\bigl (\alpha ^4\varphi -2\alpha ^2\varphi ''+\varphi ''''\bigr )-\text{i}\alpha \bigl (2y\varphi '+\varphi \bigr ) \!=\! \lambda \bigl (\varphi ''-\alpha ^2\varphi \bigr ),\quad\text{with} \quad \varphi (\pm 1)=\varphi '(\pm 1)=0. \end{align}
For each eigenvalue
$\lambda$
, there are two linearly independent real-valued velocity eigenmodes, chosen to be orthogonal via
$\psi _{A}(x,y)=\textrm{Re} ({\rm e}^{{\rm i}\theta }\,\varphi (y)\,{\rm e}^{{\rm i}\alpha x} )$
and
$\psi _{B}(x,y)=\textrm{Re} (\text{i}\,{\rm e}^{{\rm i}\theta }\,\varphi (y)\,{\rm e}^{{\rm i}\alpha x} )$
(thus mutually streamwise shifted by
$ {L}/{4n}$
). Here,
$|z|$
,
$\textrm{Re}(z)$
and
$\bar z$
denote the modulus, real part and complex conjugate of
$z\in \mathbb{C}$
, respectively. Any fixed
$\theta$
, representing a streamwise shift, will work, but we select it such that
${\rm e}^{{\rm i}\theta }\, (\varphi (0) + \text{i}\,\varphi '(0) )=|\varphi (0)| + |\varphi '(0)|$
, which ensures that
$\psi _{\!A}$
has a stagnation point at the origin.
For
$\alpha = 0$
, the analytical solutions are the eigenvalues
$\lambda =-({1}/({4\,{Re}}))(k+1)^2\pi ^2$
, each having a single velocity eigenmode
$\psi (y)=\textrm{Re} ({\rm e}^{{\rm i}(\pi /2)(k+1)}\,{\rm e}^{{\rm i}(\pi /2)(k+1)y} )$
, where
$k$
is a non-negative integer. Indeed, this solution comes from a simpler second-order eigenvalue problem for the eigenmodes
$\boldsymbol{u}=(u_0,0)=(\partial _y\psi ,0)$
, directly derived from (2.8) when
$\alpha =0$
, namely
\begin{align} \frac {1}{\textit{Re}}u_0'' = \lambda u_0,\quad \text{with} \quad u_0(\pm 1)=0. \end{align}
2.3. Numerical methods and computational implementation
To implement the inequalities in (2.7), we must first numerically solve (2.8), select a finite subset of energy eigenmodes comprising
$\mathscr{U} = \{\boldsymbol{u}_1,\ldots , \boldsymbol{u}_m\}$
(selection criteria are discussed in § 3), and then compute the relevant tensors in the polynomials
$\tilde{G}$
, and the terms
${{\unicode{x1D642}}}_i$
and
$C_i$
for
$i=1,\ldots ,m$
, which are specified in Appendix A.
To numerically solve (2.9), we use finite element methods (FEMs). With this in mind, we multiply by a sufficiently regular test function and integrate by parts accordingly, yielding a variational formulation of the form
\begin{align} \text{find}\ \varphi \in \mathscr{V}\ \text{such that}\, \mathscr{a}(\varphi , \eta ) = \lambda\, \mathscr{b}(\varphi ,\eta )\quad \forall \eta \in \mathscr{V}. \end{align}
Here,
$\mathscr{V}$
is chosen as the Sobolev space
$H_0^2(-1,1)$
, which is well approximated by piecewise polynomial functions
$\varphi$
with continuous derivative satisfying
$\varphi (\pm 1)=\varphi '(\pm 1)=0$
, and
\begin{align} \begin{aligned} \mathscr{a}(\varphi , \eta ) &= {\frac {1}{\textit{Re}}\int _{-1}^1} \big ( \varphi ''\overline {\eta ''} +2\alpha ^2\varphi '\overline {\eta '} +\alpha ^4\varphi \overline {\eta } \big )\,\mathrm{d}y -\text{i}\alpha \int _{-1}^1 (2y\varphi '+\varphi )\overline {\eta }\,\mathrm{d}y,\\ \mathscr{b}(\varphi ,\eta ) &= -{\int _{-1}^1} \big ( \varphi '\overline {\eta '} +\alpha ^2\varphi \overline {\eta } \big )\,\mathrm{d}y. \end{aligned} \end{align}
To discretise
$\mathscr{V}$
, we choose an
$H^2$
-conforming discretisation subordinate to a uniform mesh of
$\varOmega _y = (-1,1)$
with element size
$h$
, comprised of continuously differentiable locally supported cubic Hermite finite elements. This leads to a generalised eigenvalue problem of the form
${\unicode{x1D63C}}\mathsf{z} = \lambda {\unicode{x1D63D}}\mathsf{z}$
, where
${\unicode{x1D63C}}$
and
${\unicode{x1D63D}}$
are Hermitian matrices. For each
$\alpha _n$
, we assemble these matrices and solve the problem using an in-house MATLAB (The MathWorks, Inc. 2021) implementation with mesh size
$h=0.001$
, and then construct the velocity eigenmodes as described previously in § 2.2. The largest energy eigenvalue not associated with the eigenmodes in
$\mathscr{U}$
is labelled
$\kappa$
and recorded for future use (see (A1)).
The relevant tensors in
$\tilde{G}$
(namely,
${\unicode{x1D647}}$
and
${\unicode{x1D649}}$
in (A2)) involve integrals whose values are computed numerically using quadrature as they only entail piecewise polynomial finite element integrands. The matrices
${{\unicode{x1D642}}}_i$
are computed as described in Appendix B, while the bounds
$C_i$
are computed using a new accurate procedure detailed in Appendix C. Equipped with these quantities, and setting
$\varepsilon =10^{-5}$
, we are able to set up the SDP associated with (2.7) in MATLAB using YALMIP (Löfberg Reference Löfberg2004). We then determine its feasibility by solving with MOSEK v8.0.0.81 (MOSEK ApS 2017). Given
$\mathscr{U}$
, for each
$L$
we recorded the largest
$\textit{Re}$
, to within
$0.025$
, with a valid quartic Lyapunov functional. For other details of the numerical implementation and its verification, see Appendix D.
3. Results and discussion
3.1. Selection of mode sets
The physical reasoning behind this whole framework relies on selecting a relatively small explicit ‘mode set’
$\mathscr U$
, which should at least capture the nonlinear dynamics of any existing energy growth in the flow and, left to its own means (i.e. ignoring the tail
$\boldsymbol{v}$
in (2.2)), be able to globally stabilise such energy growth, eventually dissipating it entirely. In other words, the dynamics of (2.1) truncated to
$\mathscr U$
(with reference to (2.2), given by
${\mathrm d a_i}/{\mathrm d t}=\sum _{j=1}^m{\unicode{x1D647}}_{ij}a_{\!j}+\sum _{j,k=1}^m{\unicode{x1D649}}_{ijk}a_{\!j} a_k$
for
$\unicode{x1D647}$
and
$\unicode{x1D649}$
in (A2)) ought to be globally stable. That said, the residual, represented by the tail
$\boldsymbol{v}$
in (2.2) and whose isolated dynamics should have monotonically decreasing energy (i.e.
$({\mathrm d}/{\mathrm d t})(\tfrac{1}{2}q^2)=({1}/{\textit{Re}})\langle \boldsymbol{v},{\nabla} ^2 \boldsymbol{v}\rangle - \langle \boldsymbol{v},(\boldsymbol{\nabla }\boldsymbol{U}+(\boldsymbol{\nabla }\boldsymbol U)^{\rm T})\boldsymbol{v}\rangle \lt 0$
), is essential to yield rigorous global stability results for the complete dynamics described by (2.1). Thus a careful selection of
$\mathscr U$
is fundamental for the methodology to work.
In line with the reasoning above, two quick preliminary checks are essential when selecting the energy eigenmodes comprising
$\mathscr U$
. First, the largest energy eigenvalue among the omitted modes, which we refer to as
$\kappa$
, must be negative. That is, all the energy eigenmodes with instantaneous energy growth (i.e.
$\lambda \gt 0$
in (2.8)) at given
$L$
and
$\textit{Re}$
should be present in
$\mathscr U$
. Second, the truncated dynamics of (2.1) to
$\operatorname {span}(\mathscr U)$
must be linearly stable (i.e. eigenvalues of
$\unicode{x1D647}$
in (A2) should all have negative real part). Having said that, these conditions are not sufficient for (2.7) to hold, so careful experimentation is required to settle on mode sets
$\mathscr U$
with non-trivial interactions resulting in global stability.
(a) The energy eigenvalue branches of 2-D Poiseuille flow at
$\textit{Re}=100$
as a function of the streamwise wavenumber
$\alpha$
, labelled by the parity of the corresponding streamfunction. Eigenvalues consistent with
$L=3$
are labelled as
$(n,k)$
, which when
$n\neq 0$
correspond to the
$k$
th largest eigenvalues at the fixed wavenumber
$\alpha _n=( {2\pi }/{L})n$
, and when
$n=0$
correspond to
$\lambda =-( {1}/({4\,\textit{Re}}))(k+1)^2\pi ^2$
for
$k\geq 0$
. (b) The corresponding eigenmodes, including their multiplicity, along with a list of the mode sets
$\mathscr{U}_m$
to which they belong.
Figure 1. Long description
Panel A: A line graph shows the energy eigenvalue branches of 2-D Poiseuille flow at a Reynolds number of 100 as a function of the streamwise wavenumber. The x-axis is labeled with the streamwise wavenumber (alpha) ranging from 0 to 3*(2*pi/L), and the y-axis is labeled with the eigenvalue (lambda) ranging from -0.8 to 0. The graph includes solid lines for even eigenvalues and dashed lines for odd eigenvalues. Specific points on the graph are labeled with coordinates (m, n), indicating the mode sets. Panel B: Multiple contour plots display the corresponding eigenmodes, including their multiplicity. Each subplot is labeled with coordinates (m, n) and shows the spatial distribution of the eigenmodes with contour lines representing different values.
In this work, we selected several mode sets, labelled for simplicity as
$\mathscr U_m$
, with
$m$
coinciding with cardinality. In what follows, we provide some insight into the choice of the
$\mathscr U_m$
. To facilitate the discussion, it is useful to look at the energy spectrum in detail, and try to understand the eigenmode interactions. As a guide, see figure 1, which illustrates the energy eigenvalues and eigenmodes when
$L=3$
and
$\textit{Re}=100$
. Eigenvalue branches are clearly observed, so at the wavenumber
$\alpha _n=2\pi n/L$
for
$n\neq 0$
, we label the energy eigenvalue as
$(n,k)$
if it sits in the
$k$
th branch, for
$k=1,2,\ldots$
, and recall that there are two linearly independent eigenmodes associated with it (see § 2.2), clearly depicted in figure 1. Meanwhile, the
$(0,k)$
modes are a special case, labelled from
$k=0$
, with the eigenvalues being
$\lambda =-({1}/({4\,\textit{Re}}))(k+1)^2\pi ^2$
and having a single streamwise-independent energy eigenmode. Moreover, the eigenvalue branches have been classified by parity, since along them, the streamfunctions associated with the eigenmodes are either odd or even about
$y=0$
(e.g. the first branch has even streamfunctions satisfying
$\psi (x,-y)=\psi (x,y)$
). This parity is caused by the symmetries of Poiseuille flow, and was first pointed out by Orr (Reference Orr1907, p. 76).
To hone in on a minimal mode set leading to an improvement on energy stability, let
$L=3$
and consider
$\textit{Re}\gt \textit{Re}_E\approx 87.6$
. As seen in figure 1, the
$(1,1)$
mode is the only one leading to energy growth (this is true up to
$\textit{Re}\approx 124$
), so it must be in
$\mathscr U_m$
(and will lead to
$\kappa \lt 0$
). Now, to have linear stability within
$\mathscr U_m$
, there needs to be another mode that linearly interacts and stabilises the
$(1,1)$
modes, so that at the very least, the corresponding non-diagonal entry
${\unicode{x1D647}}_{ij}$
is non-zero (see (A2)). All modes of the form
$(n,k)$
for
$n\neq 1$
are naturally orthogonal with the
$(1,1)$
modes, so they do not linearly couple, and because of different parity, the same is true for the
$(1,k)$
modes for even
$k$
. Thus the first modes that linearly couple with the
$(1,1)$
modes are the
$(1,3)$
modes, which are fortunately enough to linearly stabilise them. Lastly, we add the
$(0,0)$
mode as well, yielding the minimal
$5$
-mode set
$\mathscr U_5$
. The
$(0,k)$
modes, which always have decreasing energy, do not interact linearly with the other modes, but they do, however, couple nonlinearly (through the quadratic tensor
${\unicode{x1D649}}$
in (A2)), representing a genuinely nonlinear effect, typically stabilising, which appears to be necessary to obtain non-trivial quartic Lyapunov global stability certificates. This claim is due to methodical experimentation with both the full system and the truncated system obtained when the tail is ignored. The discussion of the larger mode sets is continued below after looking at the initial results.
3.2. New global stability bounds
The energy stability limit as a function of
$L$
, denoted by
$\textit{Re}_E(L)$
, which attains its minimum
$\textit{Re}_E\approx 87.59$
at
$L_E\approx 2.99$
, is plotted in figure 2 along with the global stability results using
$\mathscr U=\mathscr U_m$
for the different mode sets. Focusing on
$\mathscr U_5$
first, we see that our construction yields a moderate yet visible improvement on
$\textit{Re}_E(L)$
for
$1.5\lt L\lt 4.3$
, where only the
$(1,1)$
modes produce energy growth. Indeed, the regions where the
$(n,k)$
modes lead to energy growth are very informative, so are shaded and labelled as
$A_{(n,k)}$
in figure 2. For
$L\gt 4.3$
, the results of
$\mathscr U_5$
coincide with those of
$\textit{Re}_E(L)$
, which is to be expected: energy growth at these
$L$
is due to the
$(2,1)$
modes, so to get a non-trivial improvement, the
$(2,1)$
modes (which are absent from
$\mathscr U_5$
) should be included in
$\mathscr U$
.
Global stability curves resulting from using the SOS–Lyapunov framework with the mode sets
$\mathscr{U}_m$
for
$m=5,6,7,9,11,13$
(see figure 1 for details), together with the energy stability limit
$\textit{Re}_E(L)$
. For each
$L$
, plotted are the largest Reynolds numbers for which a quartic Lyapunov functional was found certifying global stability of 2-D plane Poiseuille flow. Regions where the
$(n,k)$
energy eigenvalues are positive (see figure 1) are shaded and labelled as
$A_{(n,k)}$
accordingly.
Adding the
$(0,1)$
and
$(0,2)$
modes successively to
$\mathscr{U}_5$
leads to
$\mathscr{U}_6$
and
$\mathscr{U}_7$
, respectively. As can be seen from figure 2, these
$n=0$
modes provide further small improvements in stability, but to see more notable improvements, other
$n\neq 0$
modes have to be added to the mix. The addition of the
$(2,1)$
modes in
$\mathscr{U}_9$
, for example, produces much better results, which, as expected, stop abruptly at the boundary of
$A_{(2,1)}$
, when those modes yield energy growth and require other modes to linearly stabilise their effects. Inclusion of the
$(2,3)$
modes in
$\mathscr{U}_{11}$
does precisely that, and results in global stability even in
$A_{(1,1)}\cap A_{(2,1)}$
, where four modes grow in energy. In truth, reasonable and patient manual testing of mode sets by trial and error is important: just the right combination of modes might result in a significant improvement. This is precisely what happens with
$\mathscr{U}_{13}$
, where the
$(1,2)$
and
$(2,2)$
modes are added to
$\mathscr{U}_9$
, showing that the branches with odd parity may also play an important stabilising role.
The methodology is limited in that it will not yield
$L$
-independent improvements in global stability, but the results in figure 2 are stronger than they appear: for
$L\gt 4.3$
, a different
$\tilde {\mathscr{U}}_5$
comprised of the modes from
$(0,0)$
,
$(2,1)$
and
$(2,3)$
, is found to be globally stable in
$A_{(2,1)}$
, and the same with the analogous
$\tilde {\mathscr{U}}_6$
. Despite not sharing some of the shift-invariant properties of
$\textit{Re}_E(L)$
(such as
$\textit{Re}_E(nL_E)=\textit{Re}_E(L_E)$
for every positive integer
$n$
), this still hints at a parametrised family of modes, which could be studied with more analytical approaches, with the ultimate aim of achieving a result holding at a much broader range of
$L$
.
Beyond 13 modes, we run into computational limitations produced by the large memory footprint and long optimisation solve times involved in the construction of the Lyapunov functionals. These computational requirements grow very quickly as a function of
$m$
, taking only a couple of seconds and less than 1 GB when
$m=5$
, and scaling tenfold for every two modes added, reaching approximately 15 hours and 40 GB, respectively, when
$m=13$
. Note that for each
$L$
, several solves are required before finding the largest globally stable
$\textit{Re}$
at that
$L$
. A single solve for a 15-mode system was tested and took 6 days and over 512 GB of memory, which is why constructing a full curve is not currently viable. Other than the development of specialised and efficient distributed-memory parallelisable optimisation algorithms (note that MOSEK is limited to shared-memory systems), alternative computational reformulations of the methodology are desirable and likely needed for the computation of Lyapunov functionals involving a larger number of modes.
4. Concluding remarks
Using the SOS–Lyapunov framework originally proposed by Goulart & Chernyshenko (Reference Goulart and Chernyshenko2012) and refined by Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022), we find 2-D plane Poiseuille flow to be globally stable beyond the energy stability limit. For example, at the critical energy-stable streamwise length
$L_E\approx 2.99$
, where
$\textit{Re}_E\approx 87.59$
, the flow is certified to be globally stable up to
$\textit{Re} \approx 106.8$
, representing a
$22\,\%$
improvement. This is the first enhancement of a lower bound for the global stability limit of this flow since Orr (Reference Orr1907) originally computed
$\textit{Re}_E$
more than a hundred years ago, and shows, for the first time, that this laminar flow is certified to be globally stable even at
$\textit{Re}$
where transient energy growth is observed.
To achieve this, finitely many eigenmodes are carefully selected from the energy-rate eigenvalue problem (which is solved numerically using a finite element method). These small mode sets are capable of capturing sufficient features of the nonlinear dynamics of energy growth and subsequent decay. Several mode sets resulting in successful probes for quartic Lyapunov functionals at
$\textit{Re}\gt \textit{Re}_E$
are proposed, with the smallest having five modes. These mode sets could be an important ingredient in being able to prove the existence of Lyapunov functionals using more analytical approaches. In fact, in the recent work of Darrow et al. (Reference Darrow, Carlson and Goluskin2026), the same five-mode architecture, which they identified independently, undergirds their analytical construction of quartic Lyapunov functionals for this flow. It remains to study whether the larger mode sets found in this work may be used for similar analytical constructions, or whether the analytical ansätze shed light into further reducing computational costs within the SOS–Lyapunov framework.
The computational methodology is limited to certifying global stability in the range of streamwise periods,
$L$
, that are tested. Stronger global stability results in
$\textit{Re}$
or for a larger range of
$L$
values would require larger mode sets, but computational bottlenecks currently prevent us from achieving that goal. Similar statements are expected for analogous three-dimensional flows, which are of obvious interest for future research. Thus as it relates to their computational cost and memory footprint, reformulations of this methodology or development of specialised optimisation solvers resulting in better scaling properties are extremely desirable and are being actively explored.
Acknowledgements
The authors thank D. Goluskin for helpful discussions and suggestions on the original manuscript.
Funding
F.F. acknowledges the partial support of the National Center for Artificial Intelligence CENIA FB210017, Basal ANID based in Chile, and the Fondecyt Grant N. 11261732 from ANID in Chile. All the authors also gratefully acknowledge the partial support of the Office of Naval Research (ONR) award N629092312098.
Declaration of interests
The authors report no conflict of interest.
Appendix A. The SOS formulation
For completeness, we collect the definitions of the operators and constants entering the SOS formulation (2.7):
\begin{align} M_i(\boldsymbol{a},q) &= \frac {\partial V}{\partial a_i}-2\frac {\partial V}{\partial q^2}a_i, \quad \tilde{G}(\boldsymbol{a},q) = {\sum}_{i,j,k=1}^m \frac {\partial V}{\partial a_i}\big ({\unicode{x1D647}}_{ij}a_{\!j}+{\unicode{x1D649}}_{ijk}a_{\!j} a_k\big ) +2\kappa q^2\frac {\partial V}{\partial q^2}, \\[-28pt]\nonumber\end{align}
\begin{align} {\unicode{x1D647}}_{ij} &= \left \langle \boldsymbol{u}_i,\frac {1}{\textit{Re}}{\nabla} ^2\boldsymbol{u}_{\!j} -(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_{\!j} -(\boldsymbol{u}_{\!j}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{U}\right \rangle , \quad {\unicode{x1D649}}_{ijk} = -\big \langle \boldsymbol{u}_i,(\boldsymbol{u}_{\!j}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_k\big \rangle , \\[-28pt]\nonumber\end{align}
\begin{align} [{\unicode{x1D642}}_i]_{jk} &= \big \langle \tilde {\boldsymbol{h}}_{ij},\tilde {\boldsymbol{h}}_{ik}\big \rangle , \quad C_i = \max _{(x,y)\in \varOmega } \rho \left (\frac 12\big (\boldsymbol{\nabla }\boldsymbol{u}_i+(\boldsymbol{\nabla }\boldsymbol{u}_i)^{\rm T}\big )\right )\!, \\[-28pt]\nonumber \end{align}
\begin{align} \boldsymbol{h}_{i0} &= \frac {1}{\textit{Re}}\nabla^2 \boldsymbol{u}_i +(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_i -(\boldsymbol{\nabla }\boldsymbol{U})^{\rm T}\boldsymbol{u}_i, \quad \boldsymbol{h}_{ij} = (\boldsymbol{u}_{\!j}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_i -(\boldsymbol{\nabla }\boldsymbol{u}_{\!j})^{\rm T}\boldsymbol{u}_i, \\[0pt] \nonumber \end{align}
where
$\kappa$
is the largest energy eigenvalue not included in
$\mathscr U$
,
$\rho ({\cdot })$
denotes the spectral radius, and
$\tilde {\boldsymbol{h}}_{ij}$
are the solenoidal projections of
$\boldsymbol{h}_{ij}$
onto the orthogonal complement of
$\mathscr U$
. See Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022) for derivation.
Appendix B. Gram matrix calculation
To construct the Gram matrices
${{\unicode{x1D642}}}_i\in \mathbb{R}^{(m+1)\times (m+1)}$
, first note that each
$\boldsymbol h_{ij}$
in (A4) has a Helmholtz decomposition of the form
${\boldsymbol h}_{ij}={\boldsymbol h}_{ij}^{{div}}+\boldsymbol{\nabla }\phi$
, where
${\boldsymbol h}_{ij}^{{div}}$
is an
$x$
-periodic divergence-free solenoidal projection (i.e.
$\boldsymbol{\nabla }\boldsymbol{\cdot } {\boldsymbol h}_{ij}^{{div}} = 0$
) of
$\boldsymbol{h}_{ij}$
, in this case uniquely determined by the boundary conditions
${\boldsymbol h}_{ij}^{{div}}\boldsymbol{\cdot } \boldsymbol n=0$
at
$y=\pm 1$
. To compute it, note that the
$x$
-periodic scalar potential
$\phi$
satisfies the Poisson equation
\begin{align} {\nabla} ^2 \phi = \boldsymbol{\nabla }\boldsymbol{\cdot } {\boldsymbol h}_{ij}, \quad \text{with} \quad \boldsymbol{\nabla }\phi \boldsymbol{\cdot } \boldsymbol n = {\boldsymbol h}_{ij}\boldsymbol{\cdot } \boldsymbol n \quad \text{on } \quad y=\pm 1. \end{align}
Letting
$\phi (x,y) = \sum _{n\in \mathbb Z} \varphi _n(y)\,{\rm e}^{{\rm i}\alpha _n x}$
with
$\alpha _n = {2\pi n}/{L}$
allows (B1) to be projected for each
$n$
, and recast as
$\varphi _n''-\alpha _n^2\varphi _n=f_{ij,n}$
, with Neumann boundary conditions
$\varphi _n'(\pm 1)=g_{ij,n}^{\pm }$
. Here,
\begin{align} f_{ij,n}(y) = \frac {1}{L}\int _0^L \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol h}_{ij}(x,y)\,{\rm e}^{-{\rm i}\alpha _n x}\,\mathrm dx \quad \text{and}\quad g_{ij,n}^{\pm } = \frac {1}{L}\int _0^L ({\boldsymbol h}_{ij})_2(x,\pm 1)\,{\rm e}^{-{\rm i}\alpha _n x}\,\mathrm dx. \end{align}
When
$n\neq 0$
, this equation has a unique solution, but when
$n=0$
, the Neumann problem determines
$\varphi _0$
only up to an additive constant, so we fix this by imposing
$\varphi _0(-1) = 0$
(this choice is immaterial, since only
$\boldsymbol{\nabla }\phi$
enters the construction). In either case, we choose to multiply by a test function and integrate by parts to produce variational formulations of the form
\begin{equation} \begin{gathered} \text{find}\ \ \varphi \in \mathscr{W}\ \ \text{such that}\ \ \mathscr a_n^{\boldsymbol h}(\varphi , \eta ) = \ell _n^{\boldsymbol h}(\eta )\ \forall \eta \in \mathscr{ W}, \quad \text{where} \nonumber \\ \mathscr a_n^{\boldsymbol h}(\varphi , \eta ) = \int _{-1}^1 ( \varphi '\overline {\eta }' +\alpha ^2_n\varphi \overline {\eta } )\,\mathrm dy \quad \text{and}\\ \ell _n^{\boldsymbol h}(\eta ) = -\int _{-1}^1 f_{ij,n}\overline {\eta }\,\mathrm dy +g_{ij,n}^{+}\overline {\eta }(1) -(1-\delta _{0n})g_{ij,n}^{-}\overline {\eta }(-1). \nonumber \end{gathered} \end{equation}
When
$n\neq 0$
,
$\mathscr{W}$
is selected as the Sobolev space
$H^1(-1, 1)$
, whereas when
$n=0$
, we select
$\mathscr{W}$
as the subset of
$H^1(-1,1)$
whose functions vanish at
$y=-1$
. Note that (B3) must only be solved for the values of
$n$
present in
$\boldsymbol{h}_{ij}$
(which typically involves only very few modes), after which we reconstruct
$\phi$
and
${\boldsymbol h}_{ij}^{{div}} = \boldsymbol{h}_{ij} - \boldsymbol{\nabla }\phi$
.
Finally, we project
${\boldsymbol h}_{ij}^{{div}}$
onto the space orthogonal to
$\operatorname {span}(\mathscr U)$
, i.e.
\begin{align} \tilde {\boldsymbol h}_{ij} = {\boldsymbol h}_{ij}^{{div}} - \sum\nolimits _{k=1}^m \langle {\boldsymbol h}_{ij}^{{div}},\boldsymbol u_k\rangle \,\boldsymbol u_k, \end{align}
since
$\mathscr U = \{\boldsymbol{u}_1,\ldots ,\boldsymbol{u}_m \}$
is a solenoidal orthonormal set, and then compute
$[{{\unicode{x1D642}}}_i]_{jk}=\langle \tilde {\boldsymbol h}_{ij},\tilde {\boldsymbol h}_{ik}\rangle$
.
To numerically solve (B3), we use an FEM by discretising
$\mathscr W$
using Hermite finite elements as described in § 2.3, while the remaining integrals in
$\tilde {\boldsymbol{h}}_{ij}$
and
${\unicode{x1D642}}_i$
are exactly computed numerically.
Appendix C. Spectral radii calculation
To compute the
$C_i$
associated with
$\boldsymbol{u}_i=(\partial _y \psi _i,-\partial _x\psi _i)$
and specified in (A3), we use that
$\boldsymbol{\nabla }\boldsymbol{u}_i$
has zero trace (due to its incompressibility) to derive
\begin{align} C_i= \max\nolimits _{(x,y)\in \varOmega } \rho \left (\tfrac {1}{2}\big (\boldsymbol{\nabla }\boldsymbol{u}_i+ \left (\boldsymbol{\nabla }\boldsymbol{u}_i\right )^{\rm T}\big )\right) =\max\nolimits _{(x,y)\in \varOmega } \sqrt {(\partial _{xy}^2\psi _i)^2 +\tfrac {1}{4}(\partial _y^2\psi _i-\partial _x^2\psi _i)^2}. \end{align}
(Note that there is a typo in (S16) of Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022), but the formula above, which underlies the actual computations in Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022) and in this work, is the correct one.) To avoid a 2-D optimisation over
$\varOmega$
, the idea is to fix
$y$
and analytically optimise over
$x$
first, exploiting the periodicity in this direction. Assuming that the
$\boldsymbol{u}_i$
are energy eigenmodes,
$\psi _i(x,y) = \textrm{Re}(\varphi (y)\,{\rm e}^{{\rm i}\alpha x})$
for some
$\varphi$
and
$\alpha$
. We leave as an exercise to the reader to show that
\begin{gather} (\partial _{xy}^2\psi _i)^2 +\frac {1}{4}(\partial _y^2\psi _i - \partial _x^2\psi _i)^2=A(y) + B(y)\cos (2\alpha x - \delta (y)), \end{gather}
with
\begin{gather}\begin{gathered} A(y)=\tfrac{1}{2}\big(\beta(y)^2+\gamma(y)^2\big),\\ B(y)=\tfrac{1}{2}\sqrt{\beta(y)^4+\gamma(y)^4 +2\beta(y)^2\gamma(y)^2\cos\big(2(\theta_1(y)-\theta_2(y))\big)}, \end{gathered}\end{gather}
\begin{gather} \beta (y)=\alpha\, |\varphi '(y)|,\quad \gamma (y)=\frac {1}{2}|\varphi ''(y)+\alpha ^2\,\varphi (y)|, \end{gather}
\begin{gather} \theta _1(y) = \arg \big (\text{i}\alpha\, \varphi '(y)\big ), \quad \theta _2(y) = \arg \big (\varphi ''(y)+\alpha ^2\,\varphi (y)\big ), \end{gather}
where
$\delta (y)$
is some angle. The maximum of (C2) for
$x\in (0,L)$
is attained when
$\cos (2\alpha x - \delta (y))=1$
, so that
\begin{align} C_i^2=\max\nolimits_{-1\leq y\leq 1}A(y)+B(y). \end{align}
Appendix D. Verification and other implementation details
To reduce the computational cost, we restrict the ansätze of
$P$
,
$r_i$
and
$s_i$
in (2.7) by forcing them to respect the symmetries of the truncated dynamics governed by
$ {\mathrm d a_i}/{\mathrm d t}=\sum _{j=1}^m{\unicode{x1D647}}_{ij}a_{\!j}+\sum _{j,k=1}^m{\unicode{x1D649}}_{ijk}a_{\!j} a_k=\boldsymbol{f}_i(\boldsymbol{a})$
(for
$\unicode{x1D647}$
and
$\unicode{x1D649}$
in (A2)). This imposes linear relations between the coefficients of the monomials in the ansätze (such as sign symmetries), which are precisely the optimisation parameters, often resulting in some coefficients outright vanishing, and in a significant reduction in the number of free parameters. The symmetries themselves are dictated by the flow and choice of mode set
$\mathscr U$
, but ultimately, the idea is that for each equivariant symmetry
$\varPsi$
of the dynamical system (meaning a linear isometry
$\varPsi :\mathbb{R}^m\to \mathbb{R}^m$
with
$\varPsi ^K={\unicode{x1D644}}$
for some
$K\in \mathbb{N}$
satisfying
$\varPsi (\boldsymbol{f}(\boldsymbol{a}))=\boldsymbol{f}(\varPsi (\boldsymbol{a}))$
), we enforce that
$P(\varPsi (\boldsymbol{a}),q)=P(\boldsymbol{a},q)$
, and the same with the
$r_i$
and
$s_i$
. For this flow, if
$(a_{\!j},a_{j+1})$
is the pair of coefficients of the eigenmodes associated with
$(n,k)$
for
$n\neq 0$
(see figure 1), then the dynamics are invariant under
$\mathcal{R}^n$
, where
$\mathcal{R}(a_{\!j})=a_{j+1}$
,
$\mathcal{R}(a_{j+1})=-a_{\!j}$
and
$\mathcal{R}(a_k)=a_k$
for
$k\notin \{j,j+1\}$
(note that
$\mathcal{R}^4={\unicode{x1D644}}$
).
For example, for the 5-mode set
$\mathscr{U}_5$
described in § 3, one such equivariant symmetry (of ‘order’
$K=4$
) is described by the isometry
$(a_1, a_2, a_3, a_4, a_5) \mapsto (a_1, -a_3, a_2, -a_5, a_4)$
. Notice that the monomial
$a_2a_3$
is mapped to
$-a_2a_3$
under this transformation, so if
$P(\varPsi (\boldsymbol{a}),q)=P(\boldsymbol{a},q)$
, then clearly the monomial
$a_2a_3$
cannot be present in
$P$
(i.e. its monomial coefficient is zero). Similarly,
$a_3a_5$
is mapped to
$a_2a_4$
, and
$a_2a_4$
is mapped back to
$a_3a_5$
, so the invariance implies that both monomials must have the same coefficient, resulting in an invariant term of the form
$a_3a_5+a_2a_4$
in
$P$
accompanied by the same coefficient. Likewise, the combination
$a_3a_4-a_2a_5$
is invariant, as is
$\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{a}$
. As an illustration, at
$L=2.99$
and
$\textit{Re}=92.3$
, a feasible SDP solution results in the following Lyapunov functional
\begin{equation} \begin{aligned} V(\boldsymbol a, q)= E^2 &+ 0.215 a_1 E+ 0.0283 E \\ &\quad-0.00131 (a_3a_5+a_2a_4) - 0.000271 (a_3a_4-a_2a_5), \end{aligned}\end{equation}
where
$E=\tfrac{1}{2}(\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{a}+q^2)$
, and only three significant digits of each coefficient are being displayed. In this case, the Lyapunov functional shares the same structure as the one proposed by Darrow et al. (Reference Darrow, Carlson and Goluskin2026) for
$\mathscr{U}_5$
, but in principle, our ansatz can vary more freely and is applicable to any
$\mathscr{U}_m$
(which, in practice, may come at the cost of obfuscating a possible physical interpretation).
The energy eigenvalue problem (2.8) (in fact, both (2.9) and (2.10)), along with the auxiliary Poisson equations in Appendix B, are solved using FEMs with element size
$h=0.001$
, as noted in the main text. To verify that this mesh size was sufficiently small, we performed a mesh convergence analysis on specific extreme test cases, complemented by visual inspection to ensure that all solution details were resolved. Moreover, we solved the entire global stability problem with
$\mathscr{U}_9$
(see figure 1 for its definition) for
$h=0.01$
,
$h=0.001$
and
$h=0.0005$
, achieving essentially the same results in all three cases, with a maximum discrepancy in
$\textit{Re}$
of 0.025. In truth,
$h=0.01$
is probably sufficiently accurate, and an adaptive mesh with fewer elements could also attain the same level of accuracy, but for simplicity we use uniform
$h=0.001$
throughout. Finally, we used the results of 2-D plane Couette flow from Fuentes et al. (Reference Fuentes, Goluskin and Chernyshenko2022) with a 6-mode
$\mathscr{U}$
as a benchmark to compare against the independently coded global stability computations via the FEM described in this work, with the results, once again, being nearly indistinguishable.
As a completely separate endeavour, we coded and solved (2.8) and the bounds in Appendix B using ultra-spherical pseudo-spectral methods. The results were then verified against those of the FEM in both 2-D plane Couette and plane Poiseuille flows, also getting nearly indiscernible results in the aforementioned test cases. That said, sometimes, for a fixed
$L$
, we occasionally observed small isolated intervals of apparent infeasibility bracketed by feasible Reynolds numbers, which we interpreted as false negatives. These events were nearly non-existent when using the FEM implementation compared to the pseudo-spectral one. Thus all results reported in the main text involve the FEM, since it yielded much more robust and consistent results.
Near the stability boundary, in rare cases, repeated runs of the SDP solver occasionally returned different feasibility outcomes for the exact same inputs. We mitigated this sensitivity by checking neighbouring Reynolds numbers. This behaviour may be related to near-degeneracies in the SDP, and the possible use of random initial guesses in the optimisation algorithms. The use of a fixed streamwise shift convention for the
$n\neq 0$
energy eigenmodes (see § 2.2, where we ensure a stagnation point at the origin for one of the eigenmodes) may provide some numerical stability, because when this was left unspecified, we observed minor discrepancies, typically below
$\pm 0.1\,\textit{Re}$
. Finally, the use of an ‘old’ version of MOSEK, namely v8.0.0.81, is due to a tuning change in v8.1 and above, which results in a memory footprint that is simply too large for modern machines.
Feasibility (or infeasibility) of the SDP was determined using the stopping criteria of the conic interior-point solver in MOSEK, with a solution being classified as feasible only when MOSEK returned the status PRIMAL_AND_DUAL_FEASIBLE. We refer to § 9 of the manual for the MOSEK v8.0 toolbox for MATLAB (2017) for more details (or to § 13.3 for the version v11.2.2). The default termination (relative) tolerances would likely work well, but we used slightly stricter tolerances given by
$\varepsilon _p = 10^{-11}$
for primal feasibility,
$\varepsilon _d = 10^{-11}$
for dual feasibility,
$\varepsilon _g = 10^{-11}$
for the relative gap,
$\varepsilon _i = 10^{-15}$
for infeasibility,
$10^{-11}$
for the complementarity gap, and
$10^5$
for the parameter MSK_DPAR_INTPNT_CO_TOL_NEAR_REL.
Re=100
α
L=3
(n,k)
n≠0
k
αn=(2π/L)n
n=0
λ=−(1/(4Re))(k+1)2π2
k≥0
Um
Um
m=5,6,7,9,11,13
ReE(L)
L
(n,k)
A(n,k)
Open access
