# The background method: Theory and computations

Giovanni Fantuzzi Corresponding author<EMAIL_ADDRESS>Department of Aeronautics, Imperial College London, London, SW7 2AZ, UK Ali
Arslan  Department of Aeronautics, Imperial College London, London, SW7 2AZ,
UK Andrew Wynn  Department of Aeronautics, Imperial College London, London,
SW7 2AZ, UK

###### Abstract

The background method is a widely used technique to bound mean properties of
turbulent flows rigorously. This work reviews recent advances in the
theoretical formulation and numerical implementation of the method. First, we
describe how the background method can be formulated systematically within a
broader “auxiliary function” framework for bounding mean quantities, and
explain how symmetries of the flow and constraints such as maximum principles
can be exploited. All ideas are presented in a general setting and are
illustrated on Rayleigh–Bénard convection between stress-free isothermal
plates. Second, we review a semidefinite programming approach and a
timestepping approach to optimizing bounds computationally, revealing that
they are related to each other through convex duality and low-rank matrix
factorization. Open questions and promising directions for further numerical
analysis of the background method are also outlined.

MSC 2020 subject classification: 76M30, 76F25, 90C22
Keywords: Background method, Bounds, Auxiliary functions, Variational methods,
Semidefinite programming

## 1 Introduction

Making quantitative predictions for key properties of incompressible turbulent
flows, such as the mean energy dissipation in a pipe or the average amount of
heat transported by natural convection, is a fundamental problem in fluid
mechanics. A particular challenge is to determine the functional dependence of
these mean quantities on relevant nondimensional parameters, such as the
Reynolds number, directly from the equations governing the flow and without
introducing physically reasonable, but unproven, assumptions.

One way to study mean quantities rigorously, pioneered by Malkus, Howard and
Busse [1, 2, 3, 4, 5, 6], is to bound their values using variational
techniques. The idea is simple: rather than optimizing the mean quantity of
interest over solutions to the flow’s governing equations, one optimizes it
over a larger set of incompressible flow fields that satisfy only integral
constraints in the form of energy and flux balances, which are weaker but more
tractable. The latter optimum is often still hard to evaluate, but
conservative one-sided estimates (lower bounds for minima and upper bounds for
maxima) can be obtained with the background method of Doering&Constantin [7,
8, 9, 10].

To illustrate the method, suppose we seek bounds on the mean energy
dissipation for the flow of an incompressible fluid in a domain $\Omega$ with
no-slip boundaries driven by an externally applied force. The velocity
$\mathbf{u}$ of the fluid is governed by the nondimensional Navier–Stokes
equations

$\displaystyle\partial_{t}\mathbf{u}=\operatorname{Re}^{-1}\Delta\mathbf{u}-\mathbf{u}\cdot\nabla\mathbf{u}-\nabla
p+\mathbf{f},$ (1.1a) $\displaystyle\nabla\cdot\mathbf{u}=0,$ (1.1b)

where $\operatorname{Re}$ is the Reynolds number, $p$ is the pressure, and
$\mathbf{f}$ is the nondimensional external forcing. The mean energy
dissipation is

$\overline{\Phi}=\overline{\int_{\Omega}\left|\nabla\mathbf{u}\right|^{2}{\rm
d}\mathbf{x}},$ (1.2)

where we use overbars to indicate averaging over infinite time. The background
method begins by considering the energy
$\frac{1}{2}\int_{\Omega}\left|\mathbf{u}-\mathbf{U}\right|^{2}{\rm
d}\mathbf{x}$ of perturbations from a background field $\mathbf{U}$, which is
incompressible and satisfies no-slip boundary conditions, but can otherwise be
chosen arbitrarily. Assuming that this perturbation energy remains uniformly
bounded in time, as is reasonable to expect for a broad class of force vectors
$\mathbf{f}$, its time derivative has zero infinite-time average:

$\displaystyle 0$ $\displaystyle=\overline{\frac{{\rm d}}{{\rm
d}t}\frac{1}{2}\int_{\Omega}\left|\mathbf{u}-\mathbf{U}\right|^{2}{\rm
d}\mathbf{x}}$
$\displaystyle=\overline{\int_{\Omega}\left(\mathbf{u}-\mathbf{U}\right)\cdot\left(\operatorname{Re}^{-1}\Delta\mathbf{u}-\mathbf{u}\cdot\nabla\mathbf{u}-\nabla
p+\mathbf{f}\right){\rm d}\mathbf{x}}$
$\displaystyle=\overline{\int_{\Omega}\left(-\operatorname{Re}^{-1}\left|\nabla\mathbf{u}\right|^{2}+\operatorname{Re}^{-1}\nabla\mathbf{U}\cdot\nabla\mathbf{u}-\mathbf{u}\cdot\nabla\mathbf{U}\cdot\mathbf{u}+\mathbf{f}\cdot\mathbf{u}-\mathbf{f}\cdot\mathbf{U}\right){\rm
d}\mathbf{x}}.$ (1.3)

Here, the second line follows from differentiation in time using the momentum
equation 1.1a, while the last line is obtained after integration by parts
using the no-slip boundary conditions and the incompressibility of
$\mathbf{u}$ and $\mathbf{U}$. This identity enables us to express the mean
energy dissipation as

$\overline{\Phi}=\overline{\int_{\Omega}\left((1-\tfrac{\lambda}{\operatorname{Re}})\left|\nabla\mathbf{u}\right|^{2}+\tfrac{\lambda}{\operatorname{Re}}\nabla\mathbf{U}\cdot\nabla\mathbf{u}-\lambda\mathbf{u}\cdot\nabla\mathbf{U}\cdot\mathbf{u}+\lambda\mathbf{f}\cdot\mathbf{u}-\lambda\mathbf{f}\cdot\mathbf{U}\right){\rm
d}\mathbf{x}},$ (1.4)

where $\lambda$ is an arbitrary scalar often called the balance parameter.
Suppose now that, for a given value of $\operatorname{Re}$, the background
field $\mathbf{U}$ and the balance parameter $\lambda$ can be chosen such that

$\int_{\Omega}\left((1-\tfrac{\lambda}{\operatorname{Re}})\left|\nabla\mathbf{u}\right|^{2}+\tfrac{\lambda}{\operatorname{Re}}\nabla\mathbf{U}\cdot\nabla\mathbf{u}-\lambda\mathbf{u}\cdot\nabla\mathbf{U}\cdot\mathbf{u}+\lambda\mathbf{f}\cdot\mathbf{u}-\lambda\mathbf{f}\cdot\mathbf{U}\right)\leq\gamma$
(1.5)

for some constant $\gamma$ and all incompressible vector fields $\mathbf{u}$
satisfying the no-slip boundary conditions. Then, the inequality holds also
for all possible solutions of the Navier–Stokes equations (1.1a,b) and
identity 1.4 implies the upper bound $\overline{\Phi}\leq\gamma$ at the given
Reynolds number. To obtain a lower bound $\overline{\Phi}\geq\gamma$, it
suffices to reverse inequality 1.5. One is of course especially interested in
obtaining bounds $\gamma=\gamma(\operatorname{Re})$ for all or a range of
Reynolds numbers, so as to understand how the mean energy dissipation (the
mean quantity of interest) scales with $\operatorname{Re}$ (the control
parameter). As demonstrated in section 3, these ideas can be extended to a
very broad class of flows by considering generalized perturbation energies and
generalized background fields.

In the above example and in general, optimizing the background field and the
balance parameter to produce the best possible bound is a dual problem (in the
sense of convex duality) to the Malkus–Howard–Busse approach described above
[11, 12, 13] and is generally difficult. Suboptimal background fields,
however, can usually be constructed using only elementary calculus and
functional inequalities, and often yield useful bounds. For these reasons, the
background method has enjoyed tremendous success since its introduction in the
1990s (see Table 1 for a nonexhaustive list of flows to which it has been
applied) and, to this date, it remains one of the key tools for rigorous flow
analysis.

Table 1: A nonexhaustive summary of flows to which the background method has
been applied.

| Flow | Analysis | Computations |
---|---|---|---|---
| Rayleigh–Bénard convection | [10, 14, 15, 16, 17, 13, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] |  [33, 20, 21, 34, 35, 36, 30, 37, 31, 38], [39]† |
| Bénard–Marangoni convection | [40, 41] | [42] |
| Porous-media convection | [43] | [44, 45, 46, 47] |
| Internally heated convection | [48, 49, 26, 50, 51, 52, 53] | [51, 52] |
| Double-diffusive convection | [54] | none |
| Horizontal convection | [55] | none |
| Parallel shear flows | [7, 8, 56, 15, 57, 58, 11, 59, 60, 61] | [62, 63, 39, 64, 65, 66] |
| Taylor–Couette flow | [67, 68, 69] | [69] |
| Pressure-driven channel flow | [9, 70] | none |
| Precessing flow | [71] | none |
| Flows in unbounded domains | [72, 73] | none |

†Computations for Couette flow in [39] imply optimal bounds for
Rayleigh–Bénard convection [13, §4].

This work reviews recent advances in the theoretical formulation and
computational implementation of the background method. The latter are of
particular interest because, even though no algorithm can produce rigorous
bounds for all possible values of a flow’s governing parameters, numerical
approximations of the best bound available to the background method can reveal
whether suboptimal bounds proven analytically exhibit optimal parameter
dependence. If they do not, moreover, computations can guide improved
analysis. This was recently demonstrated for Bénard–Marangoni convection at
infinite Prandtl number [42, 41] and for internally heated convection [51,
53].

On the theoretical side, we show how to formulate the background method
systematically within a more general framework for bounding infinite-time
averages, which is based on auxiliary functions [74, 75, 76] and is known to
be sharp for well-posed ordinary and partial differential equations [77, 78].
This interpretation of the background method, discussed in [74, 79, 30, 76,
51, 52] and already used in the example given above, will be described in a
general setting in section 3. Its main advantages compared to the original
formulation of the method by Doering & Constantin [7, 8, 9, 10] are that (i)
it allows for the consideration of generalized perturbation energies and
generalized background fields, (ii) it can be used easily to bound mean
quantities not equivalent to the dissipation rate, and (iii) it always reduces
the search for a bound to a convex variational principle. Indeed, observe that
inequality 1.5 is jointly convex in the balance parameter $\lambda$, the
scaled background field $\lambda\mathbf{U}$, and the bound $\gamma$. We also
explain how this convex variational principle can be simplified using
symmetries and, sometimes, improved by incorporating constraints such as
maximum principles.

On the computational side (section 4), we first review two recent approaches
to optimizing background fields—one based on semidefinite programming [80, 65,
42, 36] and one based on timestepping [46, 35]—that stand out from other
strategies employed in the literature (cf. Table 1) for their generality,
robustness, and simplicity. We then reveal a previously unrecognised
connection between these two approaches, which opens new avenues for the
numerical analysis of the background method and may enable the efficient
optimization of bounds for complex flows.

To highlight the broad applicability of the background method (including its
numerical implementation strategies) and separate its key ingredients from
specific aspects pertaining to particular examples, we work in an abstract
setting similar to that in [81], which generalizes the basic Navier–Stokes
equations in 1.1 but preserves their key properties (see section 2). To
minimize the technicalities, however, we omit a rigorous functional-analytic
setup. The abstract ideas are illustrated in the context of Rayleigh–Bénard
convection between stress-free isothermal plates.

## 2 Governing equations in abstract form

We consider an incompressible fluid in a domain $\Omega\subset\mathbb{R}^{d}$
(usually, $d=2$ or $3$), whose state at time $t$ is described by a pressure
field $p(t):\Omega\to\mathbb{R}$ and a vector field
$\mathbf{u}(t):\Omega\to\mathbb{R}^{n}$. The components of $\mathbf{u}$
represent $n$ physical variables of interest, such as velocity, temperature,
and vorticity. At each time $t$, the state $\mathbf{u}(t)$ and pressure $p(t)$
are functions of the spatial coordinate $\mathbf{x}\in\Omega$, and we assume
that:

  1. (A1)

$\mathbf{u}(t)$ belongs to a linear subspace $\mathbb{U}$ of a Hilbert space
$\mathbb{V}$ with norm $\left\|\cdot\right\|_{\mathbb{V}}$, dual space
$\mathbb{V}^{\prime}$, and which is continuously embedded into the Lebesgue
space $L^{2}(\Omega;\mathbb{R}^{n})$;

  2. (A2)

$p(t)$ belongs to a Hilbert space $\mathbb{P}$ continuously embedded into
$L^{2}(\Omega)$.

For typical flows, $\mathbb{V}$ encodes regularity and boundary conditions,
while $\mathbb{U}$ includes further constraints such as incompressibility. In
the body-forced flow considered in the introduction, for instance,
$\mathbf{u}$ is simply the velocity field,
$\mathbb{V}=H^{1}_{0}(\Omega;\mathbb{R}^{n})$ is the Sobolev space of (weakly)
differentiable square-integrable vector fields that vanish on the domain
boundary and have square-integrable derivatives, $\mathbb{U}$ is its
divergence-free subspace, and $\mathbb{P}=L^{2}(\Omega)$. In the
Rayleigh–Bénard problem discussed in Example 1 below, instead, $\mathbf{u}$
represents both the fluid’s velocity and its temperature. Note that assumption
1 requires all boundary conditions defining $\mathbb{V}$ and all other
constraints defining $\mathbb{U}$ to be homogeneous; this is usually achieved
by letting $\mathbf{u}(t)$ and $p(t)$ be perturbations from a reference flow
state (say, a laminar flow).

The state of the fluid evolves according to a differential equation in the
form

$\frac{{\rm d}\mathbf{u}}{{\rm
d}t}=\mathscr{A}\mathbf{u}+\mathscr{B}(\mathbf{u},\mathbf{u})+\mathscr{C}p+\mathbf{f},\qquad\mathbf{u}(0)=\mathbf{u}_{0},$
(2.1)

where $\mathscr{A}:\mathbb{V}\to\mathbb{V}^{\prime}$ is a linear differential
operator that typically represents the effects of diffusion, advection by a
mean flow, and buoyancy forces proportional to the temperature of the fluid;
$\mathscr{B}:\mathbb{V}\times\mathbb{V}\to\mathbb{V}^{\prime}$ is a bilinear
differential operator representing nonlinear advection;
$\mathscr{C}:\mathbb{P}\to\mathbb{V}^{\prime}$ is linear differential operator
representing pressure forces; $\mathbf{f}\in L^{2}(\Omega;\mathbb{R}^{n})$ is
an externally applied steady force111The methods discussed in this work can be
extended to the case of time-dependent external forces if an explicit bound on
the spatial $L^{2}$ norm $\left\|\mathbf{f}(t)\right\|_{2}$ is available
pointwise in time, but we do not consider this for brevity.; and
$\mathbf{u}_{0}\in\mathbb{U}$ is the initial state of the fluid. Equation 2.1
is essentially an abstract version of the usual momentum equation 1.1a, where
$\mathscr{A}\mathbf{u}=\operatorname{Re}^{-1}\Delta\mathbf{u}$,
$\mathscr{B}(\mathbf{u},\mathbf{u})=-\mathbf{u}\cdot\nabla\mathbf{u}$ and
$\mathscr{C}p=-\nabla p$. Such an abstraction enables us to consider flows
that require evolution equations for quantities beyond the fluid’s velocity,
such as its temperature or vorticity, whilst using a compact notation.

The operators $\mathscr{A}$, $\mathscr{B}$ and $\mathscr{C}$ are associated to
a bounded222A $k$-linear form
$f:\mathbb{V}\times\cdots\times\mathbb{V}\to\mathbb{R}$ is bounded if
$\left|f(\mathbf{u}_{1},,\ldots,\mathbf{u}_{k})\right|\leq
C\left\|\mathbf{u}_{1}\right\|_{\mathbb{V}}\cdots\left\|\mathbf{u}_{k}\right\|_{\mathbb{V}}$
for some constant $C>0$ and all
$\mathbf{u}_{1},\ldots,\mathbf{u}_{k}\in\mathbb{V}$. bilinear form
$a:\mathbb{V}\times\mathbb{V}\to\mathbb{R}$ (not necessarily symmetric), a
bounded trilinear form
$b:\mathbb{V}\times\mathbb{V}\times\mathbb{V}\to\mathbb{R}$, and a bounded
bilinear form $c:\mathbb{V}\times\mathbb{P}\to\mathbb{R}$ such that

$\displaystyle\langle\mathbf{u}_{2},\mathscr{A}\mathbf{u}_{1}\rangle$
$\displaystyle=a(\mathbf{u}_{2},\mathbf{u}_{1})$
$\displaystyle\forall\mathbf{u}_{1},\mathbf{u}_{2}\in\mathbb{V},$ (2.2a)
$\displaystyle\langle\mathbf{u}_{2},\mathscr{B}(\mathbf{u}_{1},\mathbf{u}_{1})\rangle$
$\displaystyle=b(\mathbf{u}_{2},\mathbf{u}_{1},\mathbf{u}_{1})$
$\displaystyle\forall\mathbf{u}_{1}\in\mathbb{U},\mathbf{u}_{2}\in\mathbb{V},$
(2.2b) $\displaystyle\langle\mathbf{u},\mathscr{C}p\rangle$
$\displaystyle=c(\mathbf{u},p)$
$\displaystyle\forall\mathbf{u}\in\mathbb{V},\,p\in\mathbb{P},$ (2.2c)

where angled brackets denote the inner product in
$L^{2}(\Omega;\mathbb{R}^{n})$:

$\langle\mathbf{u}_{1},\mathbf{u}_{2}\rangle:=\int_{\Omega}\mathbf{u}_{1}\cdot\mathbf{u}_{2}\,{\rm
d}\mathbf{x}\qquad\forall\mathbf{u}_{1},\mathbf{u}_{2}\in
L^{2}(\Omega;\mathbb{R}^{n}).$ (2.3)

For the body-forced flow considered in the introduction, for example,
$a(\mathbf{u}_{2},\mathbf{u}_{1})=\operatorname{Re}^{-1}\int_{\Omega}\nabla\mathbf{u}_{1}\cdot\nabla\mathbf{u}_{2}{\rm
d}\mathbf{x}$,
$b(\mathbf{u}_{2},\mathbf{u}_{1},\mathbf{u}_{1})=\int_{\Omega}\mathbf{u}_{1}\cdot\nabla\mathbf{u}_{2}\cdot\mathbf{u}_{1}{\rm
d}\mathbf{x}$ and $c(\mathbf{u},p)=\int_{\Omega}(\nabla\cdot\mathbf{u})p\,{\rm
d}\mathbf{x}$. Recalling that for this example
$\mathbb{V}=H^{1}_{0}(\Omega;\mathbb{R}^{n})$, $\mathbb{U}$ is the divergence-
free subspace of $\mathbb{V}$, and $\mathbb{P}=L^{2}(\Omega)$, the identities
(2.2a–c) can be verified using integration by parts, the no-slip boundary
conditions of functions in $\mathbb{V}$ and, for 2.2b, the incompressibility
of $\mathbf{u}_{1}\in\mathbb{U}$.

It is also useful to consider an alternative representation for the trilinear
form $b$. For this, we define a linear operator $\mathscr{H}$ mapping a vector
field $\mathbf{u}_{1}\in\mathbb{V}$ into a linear operator
$\mathscr{H}(\mathbf{u}_{1}):\mathbb{V}\to\mathbb{V}^{\prime}$ via the
identity

$\langle\mathbf{u}_{2},\mathscr{H}(\mathbf{u}_{1})\mathbf{u}_{3}\rangle=b(\mathbf{u}_{1},\mathbf{u}_{2},\mathbf{u}_{3})\qquad\forall\mathbf{u}_{2},\mathbf{u}_{3}\in\mathbb{V}.$
(2.4)

For the body-forced flow we may take $\mathscr{H}$ to satisfy
$\mathscr{H}(\mathbf{u}_{1})\mathbf{u}_{3}=\nabla\mathbf{u}_{1}\cdot\mathbf{u}_{3}$,
and Example 1 provides the relevant definition for Rayleigh–Bénard convection.

Finally, to model the energy conservation properties of the nonlinear
advection and pressure terms in Navier–Stokes equations, we make two technical
assumptions that are usually satisfied by incompressible flows subject to no-
slip, periodic, or stress-free boundary conditions:

  1. (A3)

There exists a linear space $\mathbb{M}$ of symmetric $n\times n$ matrices
such that $b(\Lambda\mathbf{u},\mathbf{u},\mathbf{u})=0$ for all
$\Lambda\in\mathbb{M}$ and all $\mathbf{u}\in\mathbb{U}$.

  2. (A4)

There exists a linear subspace $\mathbb{M}_{0}\subseteq\mathbb{M}$ such that
$c(\Lambda\mathbf{u},p)=0$ for all $\Lambda\in\mathbb{M}_{0}$,
$\mathbf{u}\in\mathbb{U}$ and $p\in\mathbb{P}$.

In the case of body-forced flows governed by 1.1, for example,
$\mathbb{M}\equiv\mathbb{M}_{0}$ contains all multiples of the identity
matrix, since
$b(\lambda\mathbf{u},\mathbf{u},\mathbf{u})=\lambda\int_{\Omega}\mathbf{u}\cdot\nabla\mathbf{u}\cdot\mathbf{u}{\rm
d}\mathbf{x}=0$ and
$c(\lambda\mathbf{u},p)=\lambda\int_{\Omega}(\nabla\cdot\mathbf{u})p{\rm
d}\mathbf{x}=0$ by virtue of incompressibility and the no-slip boundary
conditions. For the Rayleigh–Bénard problem discussed in Example 1 below,
instead, the spaces $\mathbb{M}$ and $\mathbb{M}_{0}$ are more general and do
not coincide.

Within this abstract framework, we seek a priori upper and lower bounds on the
infinite-time average of a quantity that can be represented by a bounded
bilinear form $\Phi:\mathbb{V}\times\mathbb{V}\to\mathbb{R}$ (not necessarily
symmetric) and is associated to a linear operator
$\mathscr{P}:\mathbb{V}\to\mathbb{V}^{\prime}$ via

$\Phi(\mathbf{u},\mathbf{u}):=\langle\mathbf{u},\mathscr{P}\mathbf{u}\rangle.$
(2.5)

One example is the energy dissipation for the body-forced flow considered in
the introduction, which is given by
$\Phi(\mathbf{u},\mathbf{u})=\int_{\Omega}\left|\nabla\mathbf{u}\right|^{2}{\rm
d}\mathbf{x}$ and can be expressed in the inner-product form LABEL:{e:general-
phi-def} with $\mathscr{P}\mathbf{u}=-\Delta\mathbf{u}$ upon integration by
parts. Since the infinite-time average of $\Phi(\mathbf{u}(t),\mathbf{u}(t))$
may depend on the initial state $\mathbf{u}(0)$, we focus on bounding from
above the maximal average

$\overline{\Phi}^{*}:=\sup_{\mathbf{u}(0)\in\mathbb{U}}\overline{\Phi(\mathbf{u}(t),\mathbf{u}(t))},$
(2.6)

where overbars denote averaging over infinite time:

$\overline{\Phi(\mathbf{u}(t),\mathbf{u}(t))}:=\limsup_{T\to\infty}\frac{1}{T}\int_{0}^{T}\Phi(\mathbf{u}(t),\mathbf{u}(t))\,{\rm
d}t.$ (2.7)

Lower bounds on the minimal average are of equal interest, but can be deduced
from upper bounds on the maximal average of $-\Phi$ and will therefore not be
discussed.

######  Example 1.

Consider the problem of bounding the mean vertical heat flux in
Rayleigh–Bénard (RB) convection between stress-free plates. Using standard
nondimensional variables (see, e.g., [82]) and writing $\mathbf{x}=(x,y,z)$
for the position vector, the fluid occupies the domain
$\Omega=[-L_{x},L_{x}]\times[-L_{y},L_{y}]\times[0,1]$, which we assume to be
periodic in the horizontal directions ($x$ and $y$). The bottom plate ($z=0$)
is held at a constant temperature $T=1$, while the top one ($z=1$) is held at
$T=0$.

To lift the inhomogeneous boundary conditions on the temperature, we consider
the state vector $\mathbf{u}(t)=(\mathbf{w}(t),\vartheta(t))$ that describes
incompressible velocity perturbations $\mathbf{w}=(w_{1},w_{2},w_{3})$ and
temperature perturbations $\vartheta$ from a purely conductive state, with no
fluid motion and temperature distribution $1-z$. The flow is then described by
the nondimensional Boussinesq equations in perturbation form,

$\frac{{\rm d}}{{\rm d}t}\begin{pmatrix}\mathbf{w}\\\
\vartheta\end{pmatrix}=\underbrace{\begin{pmatrix}\operatorname{Pr}\,\Delta\mathbf{w}+\operatorname{Pr}\operatorname{Ra}\,\vartheta\,\mathbf{e}_{3}\\\
\Delta\vartheta+\mathbf{w}\cdot\mathbf{e}_{3}\end{pmatrix}}_{\mathscr{A}\mathbf{u}}+\underbrace{\begin{pmatrix}-\mathbf{w}\cdot\nabla\mathbf{w}\\\
-\mathbf{w}\cdot\nabla\vartheta\end{pmatrix}}_{\mathscr{B}(\mathbf{u},\mathbf{u})}+\underbrace{\begin{pmatrix}-\nabla
p\\\ 0\end{pmatrix}}_{\mathscr{C}p},$ (2.8)

where the Prandtl number $\Pr$ is the ratio of kinematic viscosity and thermal
diffusivity, the Rayleigh number $\operatorname{Ra}$ describes the strength of
the thermal forcing, $\mathbf{e}_{3}$ is the unit vector in the vertical
direction, and there is no external forcing ($\mathbf{f}=0$). The
instantaneous vertical heat flux is

$\Phi(\mathbf{u}(t),\mathbf{u}(t))=\int_{\Omega}\vartheta(t,\mathbf{x})\,w_{3}(t,\mathbf{x})\,{\rm
d}\mathbf{x}=\int_{\Omega}\mathbf{u}(t,\mathbf{x})\cdot\underbrace{\begin{pmatrix}0&0\\\
\mathbf{e}_{3}^{\scriptscriptstyle\top}&0\end{pmatrix}\mathbf{u}(t,\mathbf{x})}_{\mathscr{P}\mathbf{u}}\,{\rm
d}\mathbf{x}.$ (2.9)

Suitable pressure and velocity spaces are $\mathbb{P}=L^{2}(\Omega)$,
$\mathbb{V}=\mathbb{W}\times H^{1}_{p,0}(\Omega)$ with
$\mathbb{W}=H^{1}_{p}(\Omega)\times H^{1}_{p}(\Omega)\times
H^{1}_{p,0}(\Omega)$, and
$\mathbb{U}=\\{(\mathbf{w},\vartheta)\in\mathbb{V}:\nabla\cdot\mathbf{w}=0\\}$,
where $H^{1}_{p}(\Omega)$ is the space of horizontally-periodic square-
integrable functions with square-integrable weak derivatives and $H^{1}_{p,0}$
is its subspace of functions that vanish at $z=0$ and $z=1$. Writing
$\mathbf{v}=(\boldsymbol{\varphi},\tau)$ with
$\boldsymbol{\varphi}=(\varphi_{1},\varphi_{2},\varphi_{3})$, the forms $a$,
$b$ and $c$ in 2.2 are

$\displaystyle
a(\mathbf{v},\mathbf{u})=\int_{\Omega}\left(-\Pr\,\nabla\boldsymbol{\varphi}\cdot\nabla\mathbf{w}+\Pr\operatorname{Ra}\,\vartheta\varphi_{3}+\tau
w_{3}-\nabla\vartheta\cdot\nabla\tau\right){\rm d}\mathbf{x},$ (2.10a)
$\displaystyle
b(\mathbf{v},\mathbf{u},\mathbf{u})=\int_{\Omega}\left(\mathbf{w}\cdot\nabla\boldsymbol{\varphi}\cdot\mathbf{w}+\nabla\tau\cdot\mathbf{w}\vartheta\right){\rm
d}\mathbf{x},$ (2.10b) $\displaystyle
c(\mathbf{v},p)=\int_{\Omega}(\nabla\cdot\boldsymbol{\varphi})\,p\,{\rm
d}\mathbf{x}.$ (2.10c)

These expression are obtained by taking the inner product of
$\mathbf{v}\in\mathbb{V}$ with the right-hand side of 2.8, integrating some
terms by parts using the divergence-free and boundary conditions under the
assumption that $(\mathbf{w},\vartheta)\in\mathbb{U}$, and recognizing that
the result is well defined for all $(\mathbf{w},\vartheta)\in\mathbb{V}$. From
2.10b, one sees that the action of the operator $\mathscr{H}$ defined via 2.4
is
$\mathscr{H}(\mathbf{v})\mathbf{u}=\left(\nabla\boldsymbol{\varphi}\cdot\mathbf{w},\,\nabla\tau\cdot\mathbf{w}\right)$.
Finally, assumptions 3 and 4 hold with

$\mathbb{M}=\left\\{\Lambda=\begin{pmatrix}\frac{\lambda_{1}}{\operatorname{Pr}\operatorname{Ra}}I&\lambda_{3}\mathbf{e}_{3}\\\
\lambda_{3}\mathbf{e}_{3}^{\scriptscriptstyle\top}&\lambda_{2}\end{pmatrix}:\lambda_{i}\in\mathbb{R}\right\\}\quad\text{and}\quad\mathbb{M}_{0}=\\{\Lambda\in\mathbb{M}:\,\lambda_{3}=0\\},$
(2.11)

where we divide $\lambda_{1}$ by $\operatorname{Pr}\operatorname{Ra}$ to ease
the notation in subsequent examples. Indeed, integration by parts using the
divergence-free and the boundary conditions imply that
$\int_{\Omega}\mathbf{w}\cdot\nabla\mathbf{w}\cdot\mathbf{w}\,{\rm
d}\mathbf{x}=0$, $\int\nabla\vartheta\cdot\mathbf{w}\vartheta\,{\rm
d}\mathbf{x}=0$, and $\int_{\Omega}\nabla w_{3}\cdot\mathbf{w}\vartheta\,{\rm
d}\mathbf{x}=-\int_{\Omega}\mathbf{w}\cdot\nabla(\vartheta\mathbf{e}_{3})\cdot\mathbf{w}\,{\rm
d}\mathbf{x}$. Thus, for any $\Lambda$ in the space $\mathbb{M}$ specified
above and any $(\mathbf{w},\vartheta)\in\mathbb{U}$, the quantity

$b(\Lambda\mathbf{u},\mathbf{u},\mathbf{u})=\int_{\Omega}\left(\mathbf{w}\cdot\nabla\left(\tfrac{\lambda_{1}\mathbf{w}}{\operatorname{Pr}\operatorname{Ra}}+\lambda_{3}\vartheta\mathbf{e}_{3}\right)\cdot\mathbf{w}+\nabla(\lambda_{3}w_{3}+\lambda_{2}\vartheta)\cdot\mathbf{w}\vartheta\right){\rm
d}\mathbf{x}$ (2.12)

vanishes identically as required. Similarly, since $\mathbf{w}$ is divergence-
free, the quantity

$c(\Lambda\mathbf{u},p)=\int_{\Omega}\nabla\cdot\left(\tfrac{\lambda_{1}\mathbf{w}}{\operatorname{Pr}\operatorname{Ra}}+\lambda_{3}\vartheta\mathbf{e}_{3}\right)\,p\,{\rm
d}\mathbf{x}=\lambda_{3}\int_{\Omega}p\,\partial_{z}\vartheta\,{\rm
d}\mathbf{x}$ (2.13)

vanishes identically for all $(\mathbf{w},\vartheta)\in\mathbb{U}$ and all
$p\in\mathbb{P}$ only if $\lambda_{3}=0$, so $\Lambda\in\mathbb{M}_{0}$.

## 3 The background method via auxiliary functions

Upper bounds on the maximal infinite-time average $\overline{\Phi}^{*}$ can be
derived using a general strategy based on the following simple observation:
the time derivative of any bounded and differentiable function
$V:\mathbb{U}\to\mathbb{R}$, which we refer to as an auxiliary function, has
zero infinite-time average. Specifically, let
$\mathcal{L}V:\mathbb{U}\times\mathbb{P}\to\mathbb{R}$ be the Lie derivative
of $V$ along solutions of 2.1, meaning that
$\mathcal{L}V(\mathbf{u}(t),p(t))=\frac{{\rm d}}{{\rm d}t}V(\mathbf{u}(t))$.
By the fundamental theorem of calculus, the infinite-time average of
$\mathcal{L}V(\mathbf{u}(t),p(t))$ vanishes if $V(\mathbf{u}(t))$ is bounded
uniformly in time, so

$\overline{\Phi(\mathbf{u}(t),\mathbf{u}(t))}=\overline{\Phi(\mathbf{u}(t),\mathbf{u}(t))+\mathcal{L}V(\mathbf{u}(t),p(t))}$
(3.1)

irrespective of the initial state $\mathbf{u}(0)$. If $V$ can be chosen such
that

$\Phi(\mathbf{u},\mathbf{u})+\mathcal{L}V(\mathbf{u},p)\leq\gamma\qquad\forall(\mathbf{u},p)\in\mathbb{U}\times\mathbb{P}$
(3.2)

for some constant $\gamma$, then the right-hand side of 3.1 is also bounded
above by $\gamma$ irrespective of the initial state, and we conclude that
$\overline{\Phi}^{*}\leq\gamma$.

The background method amounts to searching for the smallest possible bound
$\gamma$ over quadratic auxiliary functions $V$ constructed to ensure that
inequality 3.2 is quadratic in $\mathbf{u}$ and $p$ [74, 79]. The body-forced
flow example given in the Introduction, for instance, uses the auxiliary
function
$V(\mathbf{u})=\frac{\lambda}{2}\int_{\Omega}\left|\mathbf{u}-\mathbf{U}\right|{\rm
d}\mathbf{x}$, which represents the energy of perturbations from an
incompressible background field $\mathbf{U}$, scaled by the balance parameter
$\lambda$. In this case, the general inequality 3.2 reduces exactly to 1.5 and
the pressure drops out by virtue of the incompressibility of $\mathbf{U}$.
Section 3.1 extends this example, showing how generalized background fields
and generalized perturbation energies can be used to derive bounds for any
flow whose governing equations fit within the abstract setting introduced in
section 2. Section 3.2, instead, discusses how symmetries in the governing
equations and in $\Phi$ can be exploited when optimizing $\gamma$. Finally,
section 3.3 outlines how additional constraints on the asymptotic behaviour of
the flow can be taken into account to try and improve the optimal bound.

### 3.1 General formulation

The background method considers quadratic auxiliary functions of the form

$V(\mathbf{u})=\tfrac{1}{2}\langle\Lambda\mathbf{u},\mathbf{u}\rangle-\langle\mathbf{v},\mathbf{u}\rangle,$
(3.3)

where $\mathbf{v}\in\mathbb{V}$ is a time-independent vector field and
$\Lambda\in\mathbb{M}$ is a symmetric matrix satisfying assumption 3.
Traditional formulations of the method [8, 81, 10], including that given in
the Introduction for body-forced flows, additionally take
$\mathbf{v}=\Lambda\mathbf{U}$ for some $\mathbf{U}\in\mathbb{U}$. In this
case, one may rewrite the $V$ in 3.3 as

$V(\mathbf{u})=\tfrac{1}{2}\langle\Lambda(\mathbf{u}-\mathbf{U}),\mathbf{u}-\mathbf{U}\rangle+\tfrac{1}{2}\langle\Lambda\mathbf{U},\mathbf{U}\rangle$
(3.4)

and interpret this as the “energy” of perturbations from the background field
$\mathbf{U}$, weighted by the balance parameter $\Lambda$, plus a constant
term that could be dropped if desired because it gives no contribution to the
Lie derivative $\mathcal{L}V$ appearing in 3.2. Here, however, we remove the
restriction that $\mathbf{v}=\Lambda\mathbf{U}$ and consider arbitrary
$\mathbf{v}\in\mathbb{V}$, which will be called the generalized background
field. This generalization is convenient for two reasons. First, it may
improve the best upper bound on $\overline{\Phi}^{*}$ that can be proven,
because the vectors $\mathbf{v}=\Lambda\mathbf{U}$ generated by classical
background fields usually span only a strict subspace of $\mathbb{V}$. Second,
as shown below, minimizing the upper bound $\gamma$ on $\overline{\Phi}^{*}$
over $\Lambda\in\mathbb{M}$ and $\mathbf{v}\in\mathbb{V}$ is a convex problem,
while doing so over $\Lambda\in\mathbb{M}$ and $\mathbf{U}\in\mathbb{U}$ is
not. This convexity allows for the efficient computational optimization of
$\gamma$ using the approaches described in section 4.

Differentiating 3.3 along solutions of 2.1, using the identities (2.2a–c), and
invoking assumption 3 yields

$\displaystyle\Phi(\mathbf{u},\mathbf{u})+\mathcal{L}V(\mathbf{u},p)$
$\displaystyle=\Phi(\mathbf{u},\mathbf{u})+\langle\Lambda\mathbf{u}-\mathbf{v},\,\mathscr{A}\mathbf{u}+\mathscr{B}(\mathbf{u},\mathbf{u})+\mathscr{C}p+\mathbf{f}\rangle$
$\displaystyle=c(\Lambda\mathbf{u}-\mathbf{v},p)-\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u})-q(\mathbf{u},\mathbf{u}),$
(3.5)

where $q$ and $\ell$ are quadratic and linear forms, respectively, defined as

$\displaystyle q(\mathbf{u},\mathbf{u})$
$\displaystyle:=b(\mathbf{v},\mathbf{u},\mathbf{u})-a(\Lambda\mathbf{u},\mathbf{u})-\Phi(\mathbf{u},\mathbf{u}),$
(3.6a) $\displaystyle\ell(\mathbf{u})$
$\displaystyle:=\tfrac{1}{2}a(\mathbf{v},\mathbf{u})-\tfrac{1}{2}\langle\mathbf{f},\Lambda\mathbf{u}\rangle.$
(3.6b)

Inequality 3.2 then reduces to

$c(\Lambda\mathbf{u}-\mathbf{v},p)-\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u})-q(\mathbf{u},\mathbf{u})\leq\gamma\quad\forall(\mathbf{u},p)\in\mathbb{U}\times\mathbb{P},$
(3.7)

and the smallest value of $\gamma$ for which this condition holds is the best
upper bound on $\overline{\Phi}^{*}$ that can be proven with the background
method as formulated in this section. As anticipated above, minimizing
$\gamma$ over the generalized background field $\mathbf{v}$ and balance
parameter $\Lambda$ is a convex problem because 3.7 is linear in the
optimization variables.

For fixed $\Lambda$ and $\mathbf{v}$, the best choice of $\gamma$ clearly
coincides with the supremum of the left-hand side of 3.7 over
$\mathbb{U}\times\mathbb{P}$. Since the dependence on $p$ is linear, this
supremum is finite only if $c(\Lambda\mathbf{u}-\mathbf{v},p)$ vanishes for
all $\mathbf{u}\in\mathbb{U}$ and $p\in\mathbb{P}$. In particular, we must
have
$c(\mathbf{v},p)=\langle\mathbf{v},\mathscr{C}p\rangle=\langle\mathscr{C}^{\star}\mathbf{v},p\rangle=0$
for all $p\in\mathbb{P}$, where $\mathscr{C}^{\star}$ is the adjoint of
$\mathscr{C}$, and $c(\Lambda\mathbf{u},p)=0$ for all $p\in\mathbb{P}$ and all
$\mathbf{u}\in\mathbb{U}$. The former condition requires
$\mathscr{C}^{\star}\mathbf{v}=0$, which in typical applications constrains
the velocity component of $\mathbf{v}$ to be incompressible (see Example 2
below for details in the context of RB convection). The latter, instead, can
be ensured by taking $\Lambda\in\mathbb{M}_{0}$ according to assumption 4.
These observations can be summarized as follows.

######  Theorem 1 (General background method).

The best upper bound $\gamma^{*}$ on $\overline{\Phi}^{*}$ provable with the
background method is

$\displaystyle\gamma^{*}$
$\displaystyle=\mathop{\vphantom{p}\inf}\limits_{\begin{subarray}{c}(\Lambda,\mathbf{v})\in\mathbb{M}_{0}\times\mathbb{V}\\\
\mathscr{C}^{\star}\mathbf{v}=0\end{subarray}}\;\left\\{\gamma:\;-\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u})-q(\mathbf{u},\mathbf{u})\leq\gamma\;\;\forall\mathbf{u}\in\mathbb{U}\right\\}$
(3.8a)
$\displaystyle=\mathop{\vphantom{p}\inf}\limits_{\begin{subarray}{c}(\Lambda,\mathbf{v})\in\mathbb{M}_{0}\times\mathbb{V}\\\
\mathscr{C}^{\star}\mathbf{v}=0\end{subarray}}\;\sup_{\begin{subarray}{c}\mathbf{u}\in\mathbb{U}\end{subarray}}\;\left\\{-\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u})-q(\mathbf{u},\mathbf{u})\right\\}.$
(3.8b)

Next, we state necessary and sufficient condition for this upper bound to be
finite.

######  Proposition 1 (Spectral constraint).

The upper bound on $\overline{\Phi}^{*}$ in Theorem 1 is finite if and only if
there exist $\Lambda\in\mathbb{M}_{0}$ and $\mathbf{v}\in\mathbb{V}$
satisfying $\mathscr{C}^{\star}\mathbf{v}=0$ and such that:

  1. (i)

$\ell(\mathbf{u})=0$ whenever $q(\mathbf{u},\mathbf{u})=0$;

  2. (ii)

$q(\mathbf{u},\mathbf{u})$ satisfies

$\inf_{\begin{subarray}{c}\mathbf{u}\in\mathbb{U}:\;\left\|\mathbf{u}\right\|_{2}=1\end{subarray}}q(\mathbf{u},\mathbf{u})\geq
0.$ (3.9)

######  Remark 1.

Condition ii in Proposition 1 is a spectral constraint, so called because it
requires the nonnegativity of all (real) eigenvalues $\mu$ of the self-adjoint
linear eigenvalue problem

$\left(\mathscr{H}(\mathbf{v})+\mathscr{H}(\mathbf{v})^{\star}-\Lambda\mathscr{A}-\mathscr{A}^{\star}\Lambda-\mathscr{P}-\mathscr{P}^{\star}\right)\mathbf{u}=\mu\mathbf{u},$
(3.10)

where the linear operators appearing on the left-hand side are defined in
section 2 and stars indicate adjoints (see, e.g., [8, 81, 10]). The
corresponding eigenfunctions are subject to all constraints embedded in the
definition of the space $\mathbb{U}$, plus any additional “natural” boundary
conditions that arise when deriving 3.10 as a necessary optimality condition
for $\mathbf{u}$ in 3.9. Condition i in Proposition 1, instead, is a
solvability condition for the Euler–Lagrange equation characterizing the
optimal $\mathbf{u}$ for problem 3.8b,

$\left(\mathscr{H}(\mathbf{v})+\mathscr{H}(\mathbf{v})^{\star}-\Lambda\mathscr{A}-\mathscr{A}^{\star}\Lambda-\mathscr{P}-\mathscr{P}^{\star}\right)\mathbf{u}=\Lambda\mathbf{f}-\mathscr{A}^{\star}\mathbf{v}.$
(3.11)

###### Proof of Proposition 1.

To see that conditions i, LABEL: and ii suffice to obtain a finite bound on
$\overline{\Phi}^{*}$, observe that either all $\mathbf{u}$-dependent terms in
3.8b vanish, or the nonnegative quadratic term dominates as
$\|\mathbf{u}\|_{2}$ increases. For the necessity of ii, note that the optimal
bound on $\overline{\Phi}^{*}$ stated in Theorem 1 cannot be finite unless
$q(\mathbf{u},\mathbf{u})$ is bounded below. Recalling from 3.6a, 2.5, 2.2a
and 2.2b that $q$ is a homogeneous bilinear form, the only possible finite
lower bound is zero, and it suffices to establish it for $\mathbf{u}$ with
unit $L^{2}$ norm. Condition i, instead, is necessary because if
$q(\mathbf{u}_{0},\mathbf{u}_{0})=0$ but $\ell(\mathbf{u}_{0})\neq 0$, then
one may set $\mathbf{u}=k\mathbf{u}_{0}$ and select the constant $k$ to make
$\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u})-q(\mathbf{u},\mathbf{u})=\langle\mathbf{v},\mathbf{f}\rangle-2k\ell(\mathbf{u}_{0})$
arbitrarily large in 3.8b. ∎

######  Example 2.

Let us apply the approach described above to bound the mean vertical heat
transport in the RB convection problem introduced in Example 1. With a balance
parameter $\Lambda$ in the admissible set $\mathbb{M}$ in 2.11 and a
generalized background field
$\mathbf{v}=(\boldsymbol{\varphi},\tau)\in\mathbb{V}$, the quadratic auxiliary
function in 3.3 becomes

$V(\mathbf{u})=\frac{1}{2}\underbrace{\int_{\Omega}\left(\frac{\lambda_{1}}{\operatorname{Pr}\operatorname{Ra}}\left|\mathbf{w}\right|^{2}+\lambda_{2}\left|\vartheta\right|^{2}+2\lambda_{3}\vartheta
w_{3}\right){\rm
d}\mathbf{x}}_{\langle\Lambda\mathbf{u},\mathbf{u}\rangle}-\underbrace{\int_{\Omega}\left(\boldsymbol{\varphi}\cdot\mathbf{w}+\tau\vartheta\right){\rm
d}\mathbf{x}}_{\langle\mathbf{v},\mathbf{u}\rangle}.$ (3.12)

This can be rewritten as

$\displaystyle V(\mathbf{u})=$
$\displaystyle\;\frac{\lambda_{1}}{2\operatorname{Pr}\operatorname{Ra}}\int_{\Omega}\Big{|}\mathbf{w}-\frac{\operatorname{Pr}\operatorname{Ra}}{\lambda_{1}}\boldsymbol{\varphi}\Big{|}^{2}{\rm
d}\mathbf{x}+\frac{\lambda_{2}}{2}\int_{\Omega}\Big{|}\vartheta-\frac{1}{\lambda_{2}}\tau\Big{|}^{2}{\rm
d}\mathbf{x}$ $\displaystyle+\lambda_{3}\int_{\Omega}w_{3}\vartheta\,{\rm
d}\mathbf{x}-\frac{(\operatorname{Pr}\operatorname{Ra})^{2}}{2\lambda_{1}}\int_{\Omega}\left|\boldsymbol{\varphi}\right|^{2}{\rm
d}\mathbf{x}-\frac{1}{2\lambda_{2}}\int_{\Omega}\tau^{2}\,{\rm d}\mathbf{x},$
(3.13)

where the first three terms represent the energy of deviations from the
background velocity
$(\operatorname{Pr}\operatorname{Ra}/\lambda_{1})\boldsymbol{\varphi}$, the
energy of perturbations from the background temperature $\tau/\lambda_{2}$,
and the vertical convective heat flux. Using the definition of the
instantaneous vertical heat transport $\Phi$ in 2.9 and the forms $a$, $b$ and
$c$ in 2.10a–2.10c, inequality 3.7 becomes

$\displaystyle\int_{\Omega}\left(\lambda_{3}\left|w_{3}\right|^{2}+\lambda_{3}\operatorname{Pr}\operatorname{Ra}\left|\vartheta\right|^{2}-\tfrac{\lambda_{1}}{\operatorname{Ra}}\left|\nabla\mathbf{w}\right|^{2}-\lambda_{2}\left|\nabla\vartheta\right|^{2}-\mathbf{w}\cdot\nabla\boldsymbol{\varphi}\cdot\mathbf{w}\right)\,$
$\displaystyle{\rm d}\mathbf{x}$
$\displaystyle+\int_{\Omega}\left([(1+\lambda_{1}+\lambda_{2})\mathbf{e}_{3}-\nabla\tau]\cdot\mathbf{w}\vartheta-\lambda_{3}(\operatorname{Pr}+1)\nabla\vartheta\cdot\nabla
w_{3}\right)\,$ $\displaystyle{\rm d}\mathbf{x}$
$\displaystyle+\int_{\Omega}\left(\Pr\nabla\boldsymbol{\varphi}\cdot\nabla\mathbf{w}-\Pr\operatorname{Ra}\vartheta\varphi_{3}-\tau
w_{3}+\nabla\tau\cdot\nabla\vartheta\right)\,$ $\displaystyle{\rm
d}\mathbf{x}$
$\displaystyle+\int_{\Omega}\left(\lambda_{3}\partial_{z}\vartheta-\nabla\cdot\boldsymbol{\varphi}\right)p\;$
$\displaystyle{\rm d}\mathbf{x}\;\leq\;\gamma.$ (3.14)

This must hold, by a suitable choice of $\gamma$, $\lambda_{1}$,
$\lambda_{2}$, $\lambda_{3}$, $\tau$ and $\boldsymbol{\varphi}$, for all
velocity and temperature perturbations $(\mathbf{w},\vartheta)\in\mathbb{U}$
and all pressure fields $p\in\mathbb{P}$.

Forcing the pressure terms to vanish yields
$\mathscr{C}^{\star}\mathbf{v}:=\nabla\cdot\boldsymbol{\varphi}=0$, meaning
that the velocity component of the generalized background field must be
incompressible, and $\lambda_{3}=0$, so the balance parameter $\Lambda$
belongs to the space $\mathbb{M}_{0}$ in 2.11. With these simplifications, the
spectral constraint 3.9 requires

$\int_{\Omega}\\!\left(\\!\tfrac{\lambda_{1}}{\operatorname{Ra}}\left|\nabla\mathbf{w}\right|^{2}+\lambda_{2}\left|\nabla\vartheta\right|^{2}+[\nabla\tau-(1+\lambda_{1}+\lambda_{2})\mathbf{e}_{3}]\cdot\mathbf{w}\vartheta+\mathbf{w}\cdot\nabla\boldsymbol{\varphi}\cdot\mathbf{w}\\!\right)\\!{\rm
d}\mathbf{x}\geq 0$ (3.15)

for all $(\mathbf{w},\vartheta)\in\mathbb{U}$ with unit $L^{2}$ norm, i.e.,
$\int_{\Omega}(\left|\mathbf{w}\right|^{2}+\left|\vartheta\right|^{2})\,{\rm
d}\mathbf{x}=1$. This is true if and only if the eigenvalue problem

$\displaystyle-2\lambda_{2}\,\Delta\vartheta+[\nabla\tau-(1+\lambda_{1}+\lambda_{2})\mathbf{e}_{3}]\cdot\mathbf{w}=\mu\vartheta$
(3.16a)
$\displaystyle-\frac{2\lambda_{1}}{\operatorname{Ra}}\Delta\mathbf{w}+(\nabla\boldsymbol{\varphi}+\nabla\boldsymbol{\varphi}^{\scriptscriptstyle\top})\mathbf{w}+[\nabla\tau-(1+\lambda_{1}+\lambda_{2})\mathbf{e}_{3}]\vartheta+\nabla
p=\mu\mathbf{w}$ (3.16b) $\displaystyle\nabla\cdot\mathbf{w}=0$ (3.16c)

has only nonnegative eigenvalues $\mu$, where the “pressure” $p$ is a Lagrange
multiplier that enforces the incompressibility of $\mathbf{w}$.

######  Remark 2.

Requiring $\lambda_{3}=0$ in Example 2 means that including the flux
$\int_{\Omega}w_{3}\vartheta\,{\rm d}\mathbf{x}$ in 3.13 does not help when
trying to bound the mean vertical convective heat flux in RB convection with
the approach described in this section. Similar flux terms, however, are
necessary in double-diffusive convection [54]. Moreover, they can be exploited
even in RB convection by enforcing additional constraints on the admissible
flow states [37]. This is discussed in section 3.3.

### 3.2 Symmetry reduction

Analysis and numerical treatment of the optimization problem 3.8a for the best
bound available to the background method can be considerably simplified in the
presence of symmetries. Such symmetries may include, but are not limited to,
reflection of the flow variables about a symmetry plane and translation
invariance for flows on periodic domains. Abstractly, a symmetry can be
described using a group $\mathbb{G}$ of orthogonal linear transformations
$(\mathbf{u},p)\mapsto\mathscr{G}(\mathbf{u},p)=(\mathscr{G}_{1}\mathbf{u},\mathscr{G}_{2}p)$
acting on $L^{2}(\Omega;\mathbb{R}^{n})\times\mathbb{P}$. We assume that the
subspaces $\mathbb{U}$ and $\mathbb{V}$ are closed under $\mathbb{G}$, meaning
that $\mathbf{u}\in\mathbb{U}$ (resp. $\mathbb{V}$) implies
$\mathscr{G}_{1}\mathbf{u}\in\mathbb{U}$ (resp. $\mathbb{V}$). If the
governing equation 2.1 and the quantity of interest
$\Phi(\mathbf{u},\mathbf{u})$ are invariant under $\mathbb{G}$ in a sense made
precise by Proposition 2 below, then a symmetrization argument similar to that
in [76, Appendix A] shows that auxiliary functions used to bound
$\overline{\Phi}^{*}$ may be taken to be invariant without loss of generality.

######  Proposition 2.

Let $\mathbb{G}$ be a group of orthogonal linear transformations
$(\mathbf{u},p)\mapsto\mathscr{G}(\mathbf{u},p)=(\mathscr{G}_{1}\mathbf{u},\mathscr{G}_{2}p)$.
Suppose that:

  1. 1)

$\Phi$ is invariant under $\mathbb{G}$, meaning that
$\Phi(\mathscr{G}_{1}\mathbf{u},\mathscr{G}_{1}\mathbf{u})=\Phi(\mathbf{u},\mathbf{u})$
for all $\mathscr{G}\in\mathbb{G}$.

  2. 2)

The operators $\mathscr{A}$, $\mathscr{B}$ and $\mathscr{C}$ in 2.1 are
equivariant under $\mathbb{G}$, meaning that
$\mathscr{A}\mathscr{G}_{1}\mathbf{u}=\mathscr{G}_{1}\mathscr{A}\mathbf{u}$,
$\mathscr{B}(\mathscr{G}_{1}\mathbf{u},\mathscr{G}_{1}\mathbf{u})=\mathscr{G}_{1}\mathscr{B}(\mathbf{u},\mathbf{u})$,
$\mathscr{C}\mathscr{G}_{2}p=\mathscr{G}_{1}\mathscr{C}p$ for all
$\mathscr{G}\in\mathbb{G}$.

  3. 3)

The forcing $\mathbf{f}$ is invariant under $\mathbb{G}$, meaning that
$\mathscr{G}_{1}\mathbf{f}=\mathbf{f}$ for all $\mathscr{G}\in\mathbb{G}$.

Then, any upper bound $\overline{\Phi}^{*}\leq\gamma$ provable using an
auxiliary function satisfying 3.2 can also be proven using an auxiliary
function that is invariant under $\mathbb{G}$.

This result is true for general auxiliary functions, not just the quadratic
one in 3.3 used in the background method. When applied to the latter, it
imposes structure to the generalized background field $\mathbf{v}$ that
reduces the number of optimization variables in 3.8a and, more importantly,
reveals which dynamical information the background method extracts from the
governing equation 2.1. The next example, for instance, proves that the
generalized background velocity field $\boldsymbol{\varphi}$ for the RB
convection problem in Example 2 may be set to zero without worsening the bound
on the heat transport.333This fact was known to Charlie Doering, who described
it to one of the authors as due to a “convexity argument”. Unfortunately, we
could not find a published reference to this argument. This confirms
rigorously the recent observation [38] that the background method as described
above cannot exploit the flow’s momentum equation beyond the “energy” balances
encoded by the quadratic term in the auxiliary function 3.3.

######  Example 3.

The Boussinesq equations 2.8 and the associated spaces $\mathbb{U}$,
$\mathbb{V}$ and $\mathbb{P}$ defined in Example 1 are invariant under
horizontal translations and under the “flow reversal” operation

$\begin{pmatrix}\mathbf{w}(\mathbf{x},t)\\\ \vartheta(\mathbf{x},t)\\\
p(\mathbf{x},t)\end{pmatrix}\mapsto\begin{pmatrix}G\mathbf{w}(G\mathbf{x},t)\\\
\vartheta(G\mathbf{x},t)\\\ p(G\mathbf{x},t)\end{pmatrix},\qquad
G=\begin{pmatrix}-1&0&0\\\ 0&-1&0\\\ 0&0&1\end{pmatrix}.$ (3.17)

Requiring the quadratic auxiliary function in 3.13 to be invariant under
horizontal translation forces the generalized background velocity and
temperature fields $\boldsymbol{\varphi}$ and $\tau$ to depend only on the
vertical coordinate $z$. Since $\boldsymbol{\varphi}$ must be incompressible
(cf. Example 2) and vanish at the top and bottom boundaries, we must therefore
have $\boldsymbol{\varphi}=(\varphi_{1}(z),\varphi_{2}(z),0)$. Invariance
under the transformation in (3.17) requires
$\boldsymbol{\varphi}(\mathbf{x})=G\boldsymbol{\varphi}(G\mathbf{x})$, i.e.,
$\varphi_{1}(z)=-\varphi_{1}(z)$ and $\varphi_{2}(z)=-\varphi_{2}(z)$. This
can be true only if $\boldsymbol{\varphi}=0$, proving that the generalized
background velocity field may be taken to vanish identically.

### 3.3 Imposing additional constraints

The bound on $\overline{\Phi}^{*}$ obtained in Theorem 1 can sometimes be
improved by considering further information about the asymptotic behaviour of
the flow state $\mathbf{u}$ and pressure $p$. Indeed, suppose there exists an
absorbing subset $\mathbb{A}\subsetneq\mathbb{U}\times\mathbb{P}$ in which all
solutions $\mathbf{u}(t),p(t)$ of 2.1 remain after an initial transient, or
which is at least approached exponentially quickly. Then,
$\overline{\Phi}^{*}$ is determined by trajectories that start and remain in
$\mathbb{A}$ for all times, and to derive an upper bound $\gamma$ it suffices
to replace constraint 3.7 with

$c(\Lambda\mathbf{u}-\mathbf{v},p)-\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u})-q(\mathbf{u},\mathbf{u})\leq\gamma\qquad\forall(\mathbf{u},p)\in\mathbb{A}.$
(3.18)

Clearly, imposing the inequality on $\mathbb{A}$ rather than on the full space
$\mathbb{U}\times\mathbb{P}$ cannot worsen the optimal bound
$\overline{\Phi}^{*}\leq\gamma$ available to the background method.

To ensure that the bound remains rigorous, the set $\mathbb{A}$ must be
defined using constraints that are derived from the governing equations 2.1.
When studying convective flows, for example, one can often invoke extremum
principles to enforce uniform pointwise bounds for the fluid’s temperature,
and one can sometimes consider integral estimates involving the pressure; see
[51, 37] and Example 4 below for more details. If no rigorous constraints are
available, one may also consider sets $\mathbb{A}$ defined by reasonable but
unproven conditions, such as bounds on the energy spectrum of the flow
variables [14]. In this case, one obtains conditional bounds that apply only
to solutions of 2.1 satisfying the imposed conditions, which may or may not
exist.

Irrespective of whether the absorbing set $\mathbb{A}$ is defined rigorously
or not, we now assume for concreteness that it can be described using one
pointwise linear inequality (modelling, say, an extremum principle) and one
inequality involving quadratic forms (modelling, say, an integral estimate).
Specifically, we let

$\displaystyle\mathbb{A}:=\big{\\{}(\mathbf{u},p)\in\mathbb{U}\times\mathbb{P}:\quad$
$\displaystyle h+\mathscr{K}\mathbf{u}+\mathscr{L}p\geq 0\text{ a.e. on
}\Omega,$ $\displaystyle
d_{0}+d_{1}(\mathbf{u},\mathbf{u})+d_{2}(\mathbf{u},p)+d_{3}(p,p)\geq
0\big{\\}},$ (3.19)

where the function $h\in L^{2}(\Omega)$, the linear operators
$\mathscr{K}:\mathbb{U}\to L^{2}(\Omega)$ and $\mathscr{L}:\mathbb{P}\to
L^{2}(\Omega)$, the constant $d_{0}\in\mathbb{R}$, and the bilinear forms
$d_{1}:\mathbb{U}\times\mathbb{U}\to\mathbb{R}$,
$d_{2}:\mathbb{U}\times\mathbb{P}\to\mathbb{R}$ and
$d_{3}:\mathbb{P}\times\mathbb{P}\to\mathbb{R}$ are given. Our discussion can
be easily generalized to include pointwise quadratic inequalities, equality
constraints, and multiple conditions of each type.

The constraints that define the absorbing set $\mathbb{A}$ can be enforced in
3.18 using Lagrange multipliers. Indeed, suppose that if $s\in L^{2}(\Omega)$
and $\alpha\in\mathbb{R}$ satisfy $\alpha\geq 0$ and $s(\mathbf{x})\geq 0$
almost everywhere on the fluid’s domain $\Omega$. Then, inequality 3.18 holds
if

$\displaystyle\gamma\geq\,c(\Lambda\mathbf{u}-\mathbf{v},p)-\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u})-q(\mathbf{u},\mathbf{u})\qquad$
$\displaystyle+\alpha[d_{0}+d_{1}(\mathbf{u},\mathbf{u})+d_{1}(\mathbf{u},p)+d_{3}(p,p)]\quad$
$\displaystyle+\langle s,h+\mathscr{K}\mathbf{u}+\mathscr{L}p\rangle$
$\displaystyle\quad\forall(\mathbf{u},p)\in\mathbb{U}\times\mathbb{P}.$ (3.20)

We therefore arrive at the following result.

######  Theorem 2 (Background method with extra constraints).

The best upper bound $\gamma^{*}$ on $\overline{\Phi}^{*}$ provable by
restricting the background method to the absorbing set $\mathbb{A}$ in 3.19 is

$\gamma^{*}=\inf_{\begin{subarray}{c}\alpha\in\mathbb{R},\;\Lambda\in\mathbb{M}\\\
\mathbf{v}\in\mathbb{V},\;s\in
L^{2}(\Omega)\end{subarray}}\left\\{\gamma:\;\text{\lx@cref{creftype~refnum}{e:s-procedure}},\;\alpha\geq
0,\;s(\mathbf{x})\geq 0\text{ a.e. on }\Omega\right\\}.$ (3.21)

######  Remark 3.

Since 3.21 reduces to (3.8a) when $s(\mathbf{x})=0$ and $\alpha=0$, the upper
bound on $\overline{\Phi}^{*}$ in Theorem 2 is always at least as good as that
in Theorem 1. However, it can be strictly better; this happens, for instance,
in convection driven by uniform internal heating [51] and in RB convection
between stress-free boundaries [37].

######  Remark 4.

Proposition 2 on symmetry reduction extends to the case in which the full
space $\mathbb{U}\times\mathbb{P}$ is replaced by the absorbing set
$\mathbb{A}$ if the latter is closed under the symmetry group $\mathbb{G}$,
i.e., $(\mathbf{u},p)\in\mathbb{A}$ implies that
$\mathscr{G}(\mathbf{u},p)\in\mathbb{A}$ for all $\mathscr{G}\in\mathbb{G}$. A
sufficient condition for this is that the operators $\mathscr{K},\mathscr{L}$
and the bilinear forms $d_{1},d_{2},d_{3}$ in 3.19 be invariant under
$\mathbb{G}$.

######  Example 4.

Let us reconsider the problem of RB convection between stress-free boundaries
discussed in Examples 2 and 3. By the maximum principle, the temperature
perturbation $\vartheta$ satisfies $z-\vartheta\geq 0$ and $1-z+\vartheta\geq
0$ at all times if the initial condition also does so, and it approaches such
a state exponentially quickly otherwise [83]. Moreover, analysis in [37] shows
that

$\underbrace{\int_{\Omega}\left(\tfrac{1}{4}\left|\nabla\mathbf{w}\right|^{2}\pm\left|w_{3}\right|^{2}\pm\Pr\operatorname{Ra}\left|\partial_{z}\nabla\Delta^{-1}\vartheta\right|^{2}\right){\rm
d}\mathbf{x}}_{=d_{1}(\mathbf{u},\mathbf{u})}\underbrace{\pm\int_{\Omega}p\,\partial_{z}\vartheta\,{\rm
d}\mathbf{x}}_{=d_{2}(\mathbf{u},p)}\geq 0$ (3.22)

pointwise in time. The absorbing set $\mathbb{A}$ defined by these
constraints, which are specific to stress-free boundaries and have the form
considered in 3.19, is invariant under horizontal translation and the “flow
reversal” operation 3.17. Thus, an extension of Proposition 2 guarantees that
there is no loss of generality in taking a zero generalized background
velocity field ($\boldsymbol{\varphi}=0$) and a generalized background
temperature field $\tau=\tau(z)$ depending only on the vertical coordinate.
With nonnegative Lagrange multipliers $s_{1},s_{2}\in L^{2}(\Omega)$ for the
two inequalities on $\vartheta$ and nonnegative scalar Lagrange multipliers
$\alpha_{1},\alpha_{2}$ for the two inequalities in 3.22, condition 3.20 for
the RB problem then reads

$\displaystyle\gamma\geq\int_{\Omega}\left(\left(\tfrac{\alpha_{1}}{4}+\tfrac{\alpha_{2}}{4}-\tfrac{\lambda_{1}}{\operatorname{Ra}}\right)\left|\nabla\mathbf{w}\right|^{2}+(\alpha_{1}-\alpha_{2}+\lambda_{3})\left|w_{3}\right|^{2}-\lambda_{2}\left|\nabla\vartheta\right|^{2}\right)\,$
$\displaystyle{\rm d}\mathbf{x}$
$\displaystyle+\int_{\Omega}\left((1+\lambda_{1}+\lambda_{2}-\tau^{\prime})w_{3}\vartheta-\lambda_{3}(\operatorname{Pr}+1)\nabla\vartheta\cdot\nabla
w_{3}\right)\,$ $\displaystyle{\rm d}\mathbf{x}$
$\displaystyle+\int_{\Omega}\left((\alpha_{1}-\alpha_{2})\operatorname{Pr}\operatorname{Ra}\left|\partial_{z}\nabla\Delta^{-1}\vartheta\right|^{2}+\lambda_{3}\operatorname{Pr}\operatorname{Ra}\left|\vartheta\right|^{2}\right)\,$
$\displaystyle{\rm d}\mathbf{x}$
$\displaystyle+\int_{\Omega}\left(\tau^{\prime}\partial_{z}\vartheta-\tau
w_{3}+(z-\vartheta)s_{1}+(1-z+\vartheta)s_{2}\right)\,$ $\displaystyle{\rm
d}\mathbf{x}$
$\displaystyle+\int_{\Omega}(\alpha_{1}-\alpha_{2}+\lambda_{3})\,p\,\partial_{z}\vartheta\;$
$\displaystyle{\rm d}\mathbf{x},$ (3.23)

where primes denote total derivatives in the $z$ coordinate. Since this
inequality depends linearly on the pressure $p$ and must hold for all $p$ in
the linear space $\mathbb{P}$, one is forced to choose
$\lambda_{3}=\alpha_{2}-\alpha_{1}$ so the pressure drops out of the problem.
In contrast to the basic formulation in Example 2, however, one is not forced
to set $\lambda_{3}=0$, which in principle could result in better
$\operatorname{Pr}$-dependent bounds on the mean vertical heat flux.
Unfortunately, this appears not to be the case in practice: computations using
the methods described in section 4.1 for a two-dimensional version of the
problem, obtained by dropping the $y$ direction, yield the bounds and
generalized background temperature fields shown in Figure 1 at all values of
$\operatorname{Pr}$ we tested. However, the same approach yields
$\operatorname{Pr}$-dependent bounds on the mean poloidal kinetic energy [37].

Figure 1:  Left: Optimal bounds $\gamma^{*}$ for a two-dimensional version of
the RB convection problem in Example 4 with horizontal half-period
$L_{x}=\sqrt{2}$. These bounds were obtained for fixed combinations of
Rayleigh numbers $\operatorname{Ra}$ and Prandtl numbers $\operatorname{Pr}$
by minimizing the value of $\gamma$ for which inequality Example 4 holds using
the methods described in section 4.1. The results are shifted by 1 to ease the
visualization and are independent of $\operatorname{Pr}$, as demonstrated in
the inset for $\operatorname{Ra}=10^{5}$. Right: Optimal generalized
background temperature fields at $\operatorname{Ra}=10^{3}$ (orange), $10^{4}$
(red), $10^{5}$ (brown) and $10^{6}$ (black). Only half of the domain is shown
because all profiles are antisymmetric about $z=\frac{1}{2}$.

## 4 Computational implementation

Calculating the best upper bound on $\overline{\Phi}^{*}$ that the background
method has to offer requires solving the minimization problem 3.8a or its
improved version 3.21. Usually, these problems are analytically intractable
and must be tackled numerically.

The computational implementation of the background method is complicated by
the fact that the Euler–Lagrange equations for the optimal $\Lambda$ and
$\mathbf{v}$ in 3.8a and for the optimal $\mathbf{u}$ in 3.8b are nonlinear
and, often, admit multiple solutions. Only solutions satisfying the spectral
constraint 3.9 correspond to the saddle point of problem 3.8b and give the
optimal bound on $\overline{\Phi}^{*}$; the others are “spurious” stationary
points for which the eigenvalue problem 3.10 has negative eigenvalues (see
[35, Figure 2] for a schematic illustration), and which therefore do not
produce valid bounds.

This issue has historically been addressed using delicate numerical schemes
based on continuation and bifurcation analysis [33, 62, 63, 39, 44, 64, 20,
34]. Recently, two simpler alternatives have emerged that do not require
continuation and have been applied successfully to a variety of flows. One
approach [80, 84, 65, 85, 36, 42, 37, 51] discretizes the minimization problem
3.8a into a semidefinite program (SDP)—a convex optimization problem where
matrices with affine dependence on the optimization variables are constrained
to be positive semidefinite. This SDP can then be solved using algorithms with
polynomial-time complexity [86, 87, 88, 89, 90]. The other approach [91, 46,
35, 47, 69, 66] applies the steepest gradient method to the saddle-point
formulation 3.8b to derive time-dependent versions of the Euler–Lagrange
equations, which can be timestepped until convergence to a stationary
solution. The next subsections review these two approaches and reveal a new
connection between them.

### 4.1 Optimal bounds with semidefinite programming

To discretize problem 3.8a into an SDP, suppose that all balance parameters in
the linear space $\mathbb{M}_{0}$ can be represented as

$\Lambda=\sum_{i=1}^{s}\lambda_{i}\Lambda_{i}$ (4.1)

for fixed symmetric matrices $\Lambda_{1},\ldots,\Lambda_{s}$ and some vector
of coefficients $\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{s})$.
Further, replace the spaces $\mathbb{V}$, $\mathbb{U}$ and $\mathbb{P}$ with
finite-dimensional subspaces

$\displaystyle\mathbb{V}_{h}$
$\displaystyle:=\operatorname{span}(\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{m_{v}}),$
(4.2a) $\displaystyle\mathbb{U}_{h}$
$\displaystyle:=\operatorname{span}(\boldsymbol{\eta}_{1},\ldots,\boldsymbol{\eta}_{m_{u}}),$
(4.2b) $\displaystyle\mathbb{P}_{h}$
$\displaystyle:=\operatorname{span}(\zeta_{1},\ldots,\zeta_{m_{p}}),$ (4.2c)

where $\boldsymbol{\xi}_{i}$, $\boldsymbol{\eta}_{i}$ and $\zeta_{i}$ are
prescribed basis functions and $h$ is a discretization parameter inversely
proportional to the dimension of the discrete spaces (e.g., $h$ represents the
size of a mesh for the fluid’s domain). We assume for simplicity that
$\boldsymbol{\xi}_{i}$ and $\boldsymbol{\eta}_{i}$ satisfy all constraints
used to define $\mathbb{V}$ and $\mathbb{U}$, such as boundary conditions and
incompressibility; other choices may be more convenient in practice, but
require enforcing these conditions using projection operators or Lagrange
multipliers that complicate the exposition.

Expanding

$\displaystyle\mathbf{v}(\mathbf{x})$
$\displaystyle=\sum_{i=1}^{m_{v}}\hat{v}_{i}\,\boldsymbol{\xi}_{i}(\mathbf{x}),$
$\displaystyle\mathbf{u}(\mathbf{x})$
$\displaystyle=\sum_{i=1}^{m_{u}}\hat{u}_{i}\boldsymbol{\eta}_{i}(\mathbf{x}),$
$\displaystyle p(\mathbf{x})$
$\displaystyle=\sum_{i=1}^{m_{p}}\hat{p}_{i}\,\zeta_{i}(\mathbf{x}),$ (4.3)

and writing $\hat{\mathbf{v}}=(\hat{v}_{1},\ldots,\hat{v}_{m_{v}})$ and
$\hat{\mathbf{u}}=(\hat{u}_{1},\ldots,\hat{u}_{m_{u}})$ for the vectors of
expansion coefficients, the inequality constraint in 3.8a discretizes into

$\begin{pmatrix}1\\\
\hat{\mathbf{u}}\end{pmatrix}^{\scriptscriptstyle\top}\underbrace{\begin{pmatrix}\gamma+\mathbf{f}^{h}\cdot\hat{\mathbf{v}}&\mathbf{l}(\boldsymbol{\lambda},\hat{\mathbf{v}})^{\scriptscriptstyle\top}\\\
\mathbf{l}(\boldsymbol{\lambda},\hat{\mathbf{v}})&\tfrac{1}{2}Q(\boldsymbol{\lambda},\hat{\mathbf{v}})+\tfrac{1}{2}Q(\boldsymbol{\lambda},\hat{\mathbf{v}})^{\scriptscriptstyle\top}\end{pmatrix}}_{=:S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})}\begin{pmatrix}1\\\
\hat{\mathbf{u}}\end{pmatrix}\geq
0\quad\forall\hat{\mathbf{u}}\in\mathbb{R}^{m_{u}},$ (4.4)

where the (column) vectors $\mathbf{f}^{h}$ and $\mathbf{l}$ and the matrix
$Q$ have entries

$\displaystyle f^{h}_{i}:=\langle\mathbf{f},\boldsymbol{\xi}_{i}\rangle,\quad
i=1,\ldots,m_{v},$ (4.5a)
$\displaystyle{l}_{i}:=\frac{1}{2}\sum_{j=1}^{m_{v}}a(\boldsymbol{\xi}_{j},\boldsymbol{\eta}_{i})\,\hat{v}_{j}-\frac{1}{2}\sum_{j=1}^{s}\langle\mathbf{f},\Lambda_{j}\boldsymbol{\eta}_{i}\rangle\lambda_{j},\quad
i=1,\ldots,m_{u},$ (4.5b) $\displaystyle
Q_{ij}:=\sum_{k=1}^{m_{v}}b(\boldsymbol{\xi}_{k},\boldsymbol{\eta}_{i},\boldsymbol{\eta}_{j})\hat{v}_{k}-\sum_{k=1}^{s}a(\Lambda_{k}\boldsymbol{\eta}_{i},\boldsymbol{\eta}_{j})\lambda_{k}-\Phi(\boldsymbol{\eta}_{i},\boldsymbol{\eta}_{j}),\quad
i,j=1,\ldots,m_{u}.$ (4.5c)

It is well known that 4.4 holds if and only if the matrix
$S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})$ is positive semidefinite, a
linear matrix inequality (LMI) constraint which we denoted by
$S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})\succeq 0$. The constraint
$\mathscr{C}^{\star}\mathbf{v}=0$, instead, can be imposed in a weak sense by
requiring $c(\mathbf{v},p)=0$ for all $p\in\mathbb{P}_{h}$, leading to the
linear constraints $C^{\scriptscriptstyle\top}\hat{\mathbf{v}}=0$ with matrix
coefficients $C_{ij}:=c(\boldsymbol{\xi}_{i},\zeta_{j})$ for
$i=1,\ldots,m_{v}$ and $j=1,\ldots,m_{p}$. Thus, one obtains the SDP

$\gamma_{h}^{*}=\inf_{\begin{subarray}{c}\lambda\in\mathbb{R}^{s}\\\
\hat{\mathbf{v}}\in\mathbb{R}^{m_{v}}\end{subarray}}\;\left\\{\gamma:\;S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})\succeq
0\text{ and }C^{\scriptscriptstyle\top}\hat{\mathbf{v}}=0\right\\}.$ (4.6)

If the discrete spaces $\mathbb{V}_{h}$, $\mathbb{U}_{h}$ and $\mathbb{P}_{h}$
enjoy suitable approximation properties in $\mathbb{V}$, $\mathbb{U}$ and
$\mathbb{P}$, it is reasonable to expect that the optimal solution
$\gamma_{h}^{*}$ of this SDP converges to the optimal bound $\gamma^{*}$ in
Theorem 1 as the discretization is refined ($h\to 0$). To the best of our
knowledge, however, rigorous convergence results are yet to be proven.

The SDP 4.6 can be solved using a variety of general-purpose algorithms [86,
87, 88, 89, 90], which can be divided into the two categories of interior-
point methods and first-order methods. The former are versions of Newton’s
method in which positive semidefinite constraints are imposed via so-called
barrier functions, and converge to $\varepsilon$-suboptimal solutions in
$O(\left|\log\varepsilon\right|)$ iterations [92, §4.3.3] (often, no more than
10–50). However, for an SDP with an $n\times n$ LMI and $m$ optimization
variables, each iteration of a standard interior-point method requires
$O(n^{2}+m^{2})$ memory and $O(m^{3}+n^{2}m^{2}+n^{3}m)$ floating-point
operations [92, §4.3.3], both of which increase quickly as $n$ or $m$ are
raised. First-order methods, instead, have a lower complexity because they use
only gradient information, but also converge more slowly. For example, the
algorithm in [90] finds $\varepsilon$-suboptimal solutions in
$O(\frac{1}{\varepsilon})$ iterations that require $O(m^{2}+n^{3})$ floating-
point operations each.

It must be observed that the solution $\gamma_{h}^{*}$ of 4.6 is not a
rigorous upper bound on $\overline{\Phi}^{*}$, because replacing the spaces
$\mathbb{U}$ and $\mathbb{P}$ with their subspaces $\mathbb{U}_{h}$ and
$\mathbb{P}_{h}$ relaxes the constraint in the original problem 3.8a. On the
other hand, one can always optimize $\mathbf{v}$ over the finite-dimensional
subspace $\mathbb{V}_{h}$ of $\mathbb{V}$ at the expense of worsening the
bound on $\overline{\Phi}^{*}$. Consequently, $\gamma_{h}^{*}$ estimates from
below the best (suboptimal) upper bound on $\overline{\Phi}^{*}$ that can be
proven with $\mathbf{v}\in\mathbb{V}_{h}$. To obtain an estimate from above,
one must modify the constraints of 4.6 to account for the difference between
$\mathbb{U}$ and $\mathbb{P}$ and their discrete counterparts. This is
achieved in [65, 84, 85] for discretizations based on Legendre expansions in
one coordinate and Fourier expansions in the remaining ones: using estimates
that rely on the special properties of the Fourier and Legendre bases, such as
orthogonality and differentiation rules, the infinite-dimensional constraint
in 3.8a can be strengthened into sufficient finite-dimensional LMIs, which
depend on a finite set of expansion coefficients and on the $L^{2}$ norm of
the expansion ‘tail’. Generalizing these results to other expansion bases
remains an open problem.

Another crucial observation is that the choice of basis functions for
$\mathbb{V}_{h}$, $\mathbb{U}_{h}$ and $\mathbb{P}_{h}$ strongly influences
the structure of the LMI constraint in 4.6 and, therefore, its computational
complexity. For instance, Figure 2 illustrates different LMI structures
obtained when a two-dimensional version of the RB problem considered in
Example 2, simplified as described in Example 3, is discretized using
sinusoidal functions with wavenumber ${k}$ in the horizontal direction and
either Legendre polynomials or compactly-supported piecewise-cubic functions
in the vertical direction. The block-diagonal structure corresponds to the
decoupling of different wavenumbers, and the positive semidefiniteness of each
block can be imposed separately to obtain an SDP with multiple smaller LMIs.
This is convenient because SDPs of this type can currently be solved more
efficiently than SDPs with a single large LMI; for instance, standard interior
point methods require $O(m^{3}+ks^{2}m^{2}+ks^{3}m)$ floating-point operations
per iteration to handle a $ks\times ks$ block-diagonal LMI with $m$
optimization variable and $k$ diagonal blocks of size $s\times s$, which is
usually much smaller than the $O(m^{3}+k^{2}s^{2}m^{2}+k^{3}s^{3}m)$
operations required when the block-diagonal structure is not exploited.

Figure 2:  Top: Structure of the LMI 4.4 for a two-dimensional version of the
RB convection problem in Example 2 with $\lambda_{3}=0$ and after symmetry
reduction. The basis functions are products of sinusoidal functions with
wavenumber $k$ in the horizontal direction, and either Legendre polynomials or
piecewise (PW) cubic basis functions with compact support in the vertical
direction. Bottom: Chordal decomposition of a block with nonzero wavenumber
$k$ for the PW cubic case.

The structure of each block can be exploited in a similar way using chordal
decomposition techniques for SDPs [93, 94, 95, 96, 97], which decompose sparse
LMIs into smaller ones by considering their dense principal submatrices (see
the bottom panel in Figure 2 for an illustration). This requires introducing
additional optimization variables to account for the overlap between dense
submatrices, but, if these are small and do not overlap significantly, then
the added cost is negligible compared to the savings associated by the
reduction in LMI dimension. As illustrated in Figure 3, these savings can be
significant for general LMIs with the structure shown in the bottom panel of
Figure 2 ($s$ is the number of partially overlapping $8\times 8$ blocks, equal
to $4$ in Figure 2, while $m$ is the number of optimization variables). One
can therefore optimize generalized background fields with very sharp and
possibly nonsmooth boundary layers [42, 51]. However, it must be kept in mind
that a decomposition into dense submatrices does not strengthen the LMI in 4.6
only when its structure satisfies the following technical condition, which
holds for the piecewise-cubic case in Figure 2 but not for the case of
Legendre polynomials: for every integer $\alpha\geq 3$, if the entries of
$S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})$ in position $(i_{1},i_{2})$,
$(i_{2},i_{3})$, $\ldots$, $(i_{\alpha-1},i_{\alpha})$ and
$(i_{\alpha},i_{1})$ are nonzero, then there exist
$\beta,\gamma\in\\{1,\ldots,\alpha\\}$ with $\gamma>\beta+1$ such that the
entry in position $(i_{\beta},i_{\gamma})$ is also nonzero. Employing matrix
decomposition when this condition does not hold is restrictive and, in
general, leads to conservative bounds on $\overline{\Phi}^{*}$. It remains to
be determined whether the restriction is or is not mild enough for these
conservative bounds to be useful.

Figure 3:  Ratio between the floating-point operations required by standard
interior point methods before and after the chordal decomposition of an LMI
with the structure illustrated in the bottom panel of Figure 2. Here, $m$ is
the number of optimization variable, and $s$ is the number of $8\times 8$
partially overlapping blocks ($s=4$ in Figure 2). The size of the LMI is
$(4s+8)\times(4s\times 8)$. Large contour values indicate large computational
savings due to chordal decomposition. To the left of the solid line (small $s$
or $m$), chordal decomposition brings no savings. The dashed line marks values
of $m$ for which, given $s$, chordal decomposition brings maximum advantage.

Finally, the improved background method formulation outlined in section 3.3
can be implemented using semidefinite programming following a very similar
discretization strategy. However, care must be taken to enforce the
nonnegativity of any space-dependent Lagrange multipliers exactly. The task
can be simplified considerably by a careful choice of the basis functions used
to discretize these Lagrange multipliers; piecewise-linear basis functions on
a finite-element mesh, for instance, lead to simple linear inequalities on the
values of the Lagrange multipliers on the nodes of the mesh [42, 51].

### 4.2 Optimal bounds with timestepping

An alternative approach, used to optimize bounds for the Kuramoto–Sivashinsky
equation in [91] and extended to fluid flows in [46, 35], is to view the
optimization variables in the saddle-point problem 3.8b as time-dependent and
derive evolution equations ensuring the steepest ascent in $\mathbf{u}$ and
the steepest descent in $\mathbf{v}$ and $\Lambda$. Well-known and efficient
numerical schemes for differential equations can then be employed to solve
such equations until convergence to a steady solution, which is a critical
point for 3.8b and, very often, is the correct saddle point. This approach has
been shown to work robustly for Rayleigh–Bénard convection [35], porous media
convection [47], plane Couette flow with injection and suction [66], and
Taylor–Couette flow [69].

A subtle but essential observation (see, e.g., [98, Ch. 3] and [38]) is that
tracking the time evolution of a single field $\mathbf{u}\in\mathbb{U}$ may
not be sufficient to solve 3.8b correctly. This is because, when the balance
parameter $\Lambda$ and the generalized background field $\mathbf{v}$ are
optimal, the quadratic form $q$ may vanish for a number $\rho$ of fields
$\mathbf{u}_{1},\ldots,\mathbf{u}_{\rho}\in\mathbb{U}$, which are often called
critical modes [39, 46, 35] and are eigenfunctions with zero eigenvalue for
the eigenvalue problem 3.10 associated to the spectral constraint. For reasons
that will become clear in section 4.3 these critical fields must be tracked
explicitly to impose the spectral constraint, even though they give no
contribution to the optimal value of 3.8b by Proposition 1i. One must also
track a generic field, denoted by $\mathbf{u}_{0}$, for which $q$ need not
vanish and which plays the role of $\mathbf{u}$ in 3.8b. In other words,
rather than considering the objective functional in 3.8b, one must look for a
saddle point of the modified Lagrangian

$L(\mathbf{u}_{0},\ldots,\mathbf{u}_{\rho},p,\mathbf{v},\Lambda):=-c(\mathbf{v},p)-\langle\mathbf{v},\mathbf{f}\rangle-2\ell(\mathbf{u}_{0})-\sum_{r=0}^{\rho}q(\mathbf{u}_{r},\mathbf{u}_{r}),$
(4.7)

where $p$ is a Lagrange multiplier enforcing the constraint
$\mathscr{C}^{\star}\mathbf{v}=0$.

Writing $\Lambda=\sum_{j=1}^{s}\lambda_{j}\Lambda_{j}$ as in section 4.1, the
time derivative of this Lagrangian is

$\frac{{\rm d}L}{{\rm d}t}=\left\langle\frac{\delta
L}{\delta\mathbf{v}},\frac{{\rm d}\mathbf{v}}{{\rm
d}t}\right\rangle+\sum_{j=1}^{s}\left\langle\frac{\partial
L}{\partial\lambda_{j}},\frac{{\rm d}\lambda_{j}}{{\rm
d}t}\right\rangle+\sum_{r=0}^{\rho}\left\langle\frac{\delta
L}{\delta\mathbf{u}_{r}},\frac{{\rm d}\mathbf{u}_{r}}{{\rm
d}t}\right\rangle-\left\langle\mathscr{C}^{\star}\mathbf{v},\frac{{\rm
d}p}{{\rm d}t}\right\rangle,$ (4.8)

where $\frac{\delta L}{\delta\psi}$ denotes the variation (Fréchet derivative)
of $L$ with respect to a field $\psi$. Insisting that
$\mathscr{C}^{\star}\mathbf{v}=0$ pointwise in time, so the last term
vanishes, the evolution equations ensuring the fastest growth in $L$ with
respect to each $\mathbf{u}_{r}$ and the fastest reduction with respect to
$\mathbf{v}$ and each $\lambda_{j}$ are

$\displaystyle\frac{{\rm d}\mathbf{v}}{{\rm d}t}$ $\displaystyle=-\frac{\delta
L}{\delta\mathbf{v}}\equiv\mathscr{A}\mathbf{u}_{0}+\sum_{r=0}^{\rho}\mathscr{B}(\mathbf{u}_{r},\mathbf{u}_{r})+\mathscr{C}p+\mathbf{f},$
(4.9a) $\displaystyle\frac{{\rm d}\lambda_{j}}{{\rm d}t}$
$\displaystyle=-\frac{\partial
L}{\partial\lambda_{j}}\equiv-\langle\mathbf{f},\Lambda_{j}\mathbf{u}_{0}\rangle-\sum_{r=0}^{\rho}a(\Lambda_{j}\mathbf{u}_{r},\mathbf{u}_{r})$
$\displaystyle\text{for }j=1,\ldots,s,$ (4.9b) $\displaystyle\frac{{\rm
d}\mathbf{u}_{r}}{{\rm d}t}$ $\displaystyle=\;\frac{\delta
L}{\delta\mathbf{u}_{r}}\;\equiv\bigg{(}\sum_{j=1}^{s}\lambda_{j}\Lambda_{j}\mathbf{f}-\mathscr{A}^{\star}\mathbf{v}\bigg{)}\delta_{0r}-\mathscr{F}(\mathbf{v})\mathbf{u}_{r}$
$\displaystyle\text{for }r=0,\ldots,\rho,$ (4.9c)

where $\delta_{0r}$ is the usual Kronecker delta and we have introduced the
linear operator
$\mathscr{F}(\mathbf{v}):=\mathscr{H}(\mathbf{v})+\mathscr{H}(\mathbf{v})^{\star}-\Lambda\mathscr{A}-\mathscr{A}^{\star}\Lambda-\mathscr{P}-\mathscr{P}^{\star}$
for notational convenience. Stationary solutions to these equations are
stationary points of the Lagrangian $L$, and those satisfying the spectral
constraint 3.9 provide optimal solutions to 3.8b.

Projecting 4.9a, LABEL:, 4.9b, LABEL: and 4.9c onto the finite-dimensional
spaces $\mathbb{U}_{h}$, $\mathbb{V}_{h}$ and $\mathbb{P}_{h}$ defined in 4.2
yields the set of ODEs

$\displaystyle\frac{{\rm d}(V\hat{\mathbf{v}})}{{\rm d}t}$
$\displaystyle=A\hat{\mathbf{u}}_{0}+\sum_{r=0}^{\rho}\mathbf{b}(\hat{\mathbf{u}}_{r},\hat{\mathbf{u}}_{r})+C\hat{\mathbf{p}}+\mathbf{f}^{h},$
(4.10a) $\displaystyle\frac{{\rm d}\lambda_{j}}{{\rm d}t}$
$\displaystyle=-\sum_{i=1}^{m_{u}}\langle\mathbf{f},\Lambda_{j}\boldsymbol{\eta}_{i}\rangle\hat{u}_{0i}-\sum_{r=0}^{\rho}\,\sum_{i,k=1}^{m_{u}}\\!a(\Lambda_{j}\boldsymbol{\eta}_{i},\boldsymbol{\eta}_{k})\hat{u}_{ri}\hat{u}_{rk}$
$\displaystyle\text{for }j=1,\ldots,s,$ (4.10b) $\displaystyle\frac{{\rm
d}(U\hat{\mathbf{u}}_{r})}{{\rm d}t}$
$\displaystyle=-2\mathbf{l}(\boldsymbol{\lambda},\hat{\mathbf{v}})\delta_{0r}-\left[Q(\boldsymbol{\lambda},\hat{\mathbf{v}})+Q(\boldsymbol{\lambda},\hat{\mathbf{v}})^{\scriptscriptstyle\top}\right]\hat{\mathbf{u}}_{r},\phantom{\sum_{r=0}^{\rho}}$
$\displaystyle\text{for }r=0,\ldots,\rho.$ (4.10c)

Here, $\hat{\mathbf{u}}_{r}=(\hat{u}_{r1},\ldots,\hat{u}_{rm_{u}})$ is the
vector of expansion coefficients of $\mathbf{u}_{r}$ in the chosen basis for
$\mathbb{U}_{h}$, the vector
$\mathbf{b}(\hat{\mathbf{u}}_{r},\hat{\mathbf{u}}_{r})$ has entries

$b_{i}(\hat{\mathbf{u}}_{r},\hat{\mathbf{u}}_{r}):=\sum_{j,k=1}^{m_{u}}b(\boldsymbol{\xi}_{i},\boldsymbol{\eta}_{j},\boldsymbol{\eta}_{k})\hat{u}_{rj}\,\hat{u}_{rk},\qquad
i=1,\ldots,m_{v},$ (4.11)

the matrices $V$, $A$ and $U$ have entries

$\displaystyle V_{ij}$
$\displaystyle=\langle\boldsymbol{\xi}_{i},\boldsymbol{\xi}_{j}\rangle,$
$\displaystyle i,j$ $\displaystyle=1,\ldots,m_{v},$ (4.12a) $\displaystyle
A_{ij}$ $\displaystyle=a(\boldsymbol{\xi}_{i},\boldsymbol{\eta}_{j}),$
$\displaystyle i$ $\displaystyle=1,\ldots,m_{v},\;j=1,\ldots,m_{u},$ (4.12b)
$\displaystyle U_{ij}$
$\displaystyle=\langle\boldsymbol{\eta}_{i},\boldsymbol{\eta}_{j}\rangle,$
$\displaystyle i,j$ $\displaystyle=1,\ldots,m_{u},$ (4.12c)

and all other quantities are defined in section 4.1. The ODEs in 4.10,
complemented by the algebraic condition
$C^{\scriptscriptstyle\top}\hat{\mathbf{v}}=0$ obtained upon discretizing the
constraint $\mathscr{C}^{\star}\mathbf{v}=0$ as explained in section 4.1, can
be integrated in time until convergence to a stationary solution using
standard numerical schemes, starting from any randomly generated nonzero
initial condition.

######  Remark 5.

The value of $\rho$ in 4.7 must be chosen a priori and must be large enough to
ensure that the saddle point of $L$ is the only stationary attractor for 4.9,
so the timestepping process cannot converge to a “spurious” solution. The
discussion in section 4.3 below and the numerical results in [46, 35, 47, 69,
66] and [98, Ch. 3] lead us to conjecture that $\rho$ should be equal to the
largest possible multiplicity of the principal eigenvalue of the eigenvalue
problem 3.10 associated to the spectral constraint 3.9. Analysis in [35]
establishes this conjecture when the background method is applied with fixed
balance parameter $\Lambda$ (i.e., fixed $\lambda_{1},\ldots,\lambda_{s}$) to
RB convection with stress-free isothermal boundaries, plane Couette flow, and
porous media convection. It remains an open theoretical problem to extend
these specialized arguments to the general setting considered in the present
work.

######  Remark 6.

In highly turbulent regimes (e.g., convection at very high Rayleigh number),
the numerical solution of the ODE system 4.10 typically requires small
timesteps to ensure numerical stability. In such regimes, therefore,
convergence to a stationary solution can be slow. One way to mitigate this
issue is to employ a two-step approach [35] in which the timestepping is
stopped as soon as the magnitude of the time derivatives drops below a
moderate tolerance, and the resulting approximately stationary solution is
used to initialize Newton–Kantorovich iterations for the Euler–Lagrange
equations of the saddle-point problem 4.7. As discussed in [46, 35],
efficiency can be further improved by adaptively dropping from the computation
any vectors $\hat{\mathbf{u}}_{r}$ with $r\geq 1$ that appear to converge to
zero, because they have no effect on the long-time behaviour of solutions to
4.10.

### 4.3 Connecting the two approaches

The timestepping strategy of section 4.2 is closely related to the
semidefinite programming approach of section 4.1, in the sense that it
provides a particular numerical scheme for solving the SDP 4.6. Precisely, the
ODEs in 4.10 can be viewed as a steepest-gradient algorithm applied to the SDP
4.6 when the matrix Lagrange multiplier for its LMI constraint is represented
using a low-rank factorization.

To make this interpretation evident, let us assume that the SDP 4.6 is
strictly feasible, meaning that there exist $\boldsymbol{\lambda}$ and
$\hat{\mathbf{v}}$ such that $C^{\scriptscriptstyle\top}\hat{\mathbf{v}}=0$
and that the matrix $S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})$ is
strictly positive definite. With this additional (mild) assumption, standard
duality theory for SDPs [99, §5.9.1] guarantees that

$\gamma_{h}^{*}=\mathop{\vphantom{p}\inf}\limits_{\begin{subarray}{c}\gamma\in\mathbb{R},\;\hat{\mathbf{v}}\in\mathbb{R}^{m_{v}}\\\
\boldsymbol{\lambda}\in\mathbb{R}^{s}\end{subarray}}\;\sup_{\begin{subarray}{c}\hat{p}\in\mathbb{R}^{m_{p}}\\\
Z\succeq 0\end{subarray}}\;\left\\{\gamma-\hat{\mathbf{v}}\cdot
C\hat{\mathbf{p}}-Z\cdot
S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})\right\\},$ (4.13)

where $\hat{\mathbf{p}}$ and $Z$ are Lagrange multipliers for the constraints
of 4.6. Moreover, by a condition known as complementary slackness, the ranks
of the optimal $Z$ and of the optimal
$S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})$ must sum to the matrix size
$m_{u}+1$. Therefore, if one can anticipate that the optimal solution of the
SDP 4.6 satisfies
$\operatorname{rank}(S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}}))\geq
m_{u}-\rho$ for some integer $\rho\leq m_{u}$, then the corresponding optimal
multiplier must have $\operatorname{rank}(Z)\leq\rho+1$.

Now, any positive semidefinite matrix $Z$ of dimension $m_{u}+1$ and rank no
larger than $\rho+1$ can be expressed as the sum of at most $\rho+1$ rank-1
matrices. In particular, we can write

$Z=\begin{pmatrix}t^{2}&t\hat{\mathbf{u}}_{0}^{\scriptscriptstyle\top}\\\
t\hat{\mathbf{u}}_{0}&\hat{\mathbf{u}}_{0}\hat{\mathbf{u}}_{0}^{\scriptscriptstyle\top}\end{pmatrix}+\sum_{r=1}^{\rho}\begin{pmatrix}0&\mathbf{0}^{\scriptscriptstyle\top}\\\
\mathbf{0}&\hat{\mathbf{u}}_{r}\hat{\mathbf{u}}_{r}^{\scriptscriptstyle\top}\end{pmatrix}$
(4.14)

for some vectors
$\hat{\mathbf{u}}_{0},\ldots,\hat{\mathbf{u}}_{\rho}\in\mathbb{R}^{m_{u}}$ and
some nonnegative scalar $t$. Substituting this decomposition and the
definition of $S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})$ from 4.4 into
4.13, and unfolding matrix inner products involving rank-1 matrices into
standard matrix-vector products, leads to

$\gamma_{h}^{*}=\mathop{\vphantom{p}\inf}\limits_{\begin{subarray}{c}\gamma\in\mathbb{R}\\\
\hat{\mathbf{v}}\in\mathbb{R}^{m_{v}}\\\
\boldsymbol{\lambda}\in\mathbb{R}^{s}\end{subarray}}\;\sup_{\begin{subarray}{c}\hat{\mathbf{u}}_{0},\ldots,\hat{\mathbf{u}}_{\rho}\in\mathbb{R}^{m_{u}}\\\
\hat{p}\in\mathbb{R}^{m_{p}}\\\ t\geq
0\end{subarray}}\bigg{\\{}(1-t^{2})\gamma-\hat{\mathbf{v}}\cdot
C\hat{\mathbf{p}}-t^{2}\,\mathbf{f}^{h}\cdot\hat{\mathbf{v}}\\\\[-23.24991pt]
-2t\,\mathbf{l}(\boldsymbol{\lambda},\hat{\mathbf{v}})\cdot\hat{\mathbf{u}}_{0}-\sum_{r=0}^{\rho}\hat{\mathbf{u}}_{r}\cdot
Q(\boldsymbol{\lambda},\hat{\mathbf{v}})\hat{\mathbf{u}}_{r}\bigg{\\}}.$
(4.15)

Optimizing over $\gamma$ yields $t=1$, in which case the argument of this inf-
sup problem reduces exactly to the discrete version of the Lagrangian $L$ in
4.7. Further optimizing over $\hat{\mathbf{p}}$ returns the constraint
$C^{\scriptscriptstyle\top}\hat{\mathbf{v}}=0$, while the ODEs in 4.10 are
recovered upon fixing $t=1$ and applying steepest ascent in each
$\hat{\mathbf{u}}_{r}$ and steepest descent in $\hat{\mathbf{v}}$ and
$\boldsymbol{\lambda}$ with respect to the discrete norms
$\|\hat{\mathbf{u}}_{r}\|_{\mathbb{U}_{h}}=(\hat{\mathbf{u}}_{r}\cdot
U\hat{\mathbf{u}}_{r})^{1/2}$,
$\|\hat{\mathbf{v}}\|_{\mathbb{V}_{h}}=(\hat{\mathbf{v}}\cdot
V\hat{\mathbf{v}})^{1/2}$ and
$\|\boldsymbol{\lambda}\|_{2}=\sqrt{\boldsymbol{\lambda}\cdot\boldsymbol{\lambda}}$.
At the discrete level, therefore, the timestepping approach of section 4.2 is
a particular example of general low-rank factorization strategies for SDPs
that have been studied extensively [100, 101, 102, 103, 104, 105].

This realization enables us to make interesting observations regarding the
theoretical convergence properties of the timestepping approach—specifically,
about whether one can guarantee the avoidance of spurious solutions in
general, as suggested by the numerical evidence in [46, 35, 47, 69, 66]. On
the one hand, the timestepping approach should not be expected to work
correctly unless $\rho+1$ is at least as large as the rank of the optimal
matrix $Z$ in 4.13. This is why it is usually necessary to include the fields
$\mathbf{u}_{1},\ldots,\mathbf{u}_{\rho}$ in 4.7. By complementary slackness,
the rank of the optimal $Z$ coincides with the number of zero eigenvalues of
the optimal $S(\gamma,\boldsymbol{\lambda},\hat{\mathbf{v}})$. Since the
bottom-right block of this matrix is the discrete counterpart of the self-
adjoint operator on the left-hand side of the eigenvalue problem 3.10, we
conclude that $\rho$ should be no smaller than the maximum possible
multiplicity of the principal eigenvalue of that problem as $\mathbf{v}$ and
$\Lambda=\sum_{j=1}^{s}\lambda_{j}\Lambda_{j}$ are varied.

On the other hand, if we assume that the SDP 4.6 and its dual have smooth
constraint sets in the sense of [104], then general results for low-rank
factorization methods for SDPs [104] guarantee that the timestepping approach
avoids spurious solutions when $\rho+1$ is larger than an upper bound on the
rank of the optimal $Z$ derived in [106], which depends only on the size of
the optimization variables in 4.13. This bound, however, is usually much
larger than the maximum eigenvalue multiplicity mentioned above, so its use
results in significantly more expensive computations. Unfortunately, it is
also sharp [105]: there exist SDPs (not necessarily arising from applications
of the background method) for which setting $\rho+1$ to be larger than the
rank of the optimal $Z$ but smaller than the upper bound from [106] can result
in convergence to locally stable spurious solutions. This difficulty,
moreover, cannot be resolved by replacing gradient-based timestepping with
another algorithm for saddle-point problems because the results in [105] are
independent of the numerical scheme used to solve 4.15. It would be
interesting to investigate whether locally stable spurious solutions can arise
for SDPs coming from the background method, and one possible way forward is to
try and combine the general analysis in [104, 105] with the convergence proofs
for specific flows in [35]. Based on the available numerical evidence in [46,
35, 47, 69, 66] and in [98, Ch. 3], we make the following conjecture.

######  Conjecture.

Taking $\rho$ equal to the maximum possible multiplicity of the eigenvalue
problem 3.10 suffices to optimize background fields correctly.

While proving or disproving this statement remains an open theoretical
problem, we stress that the timestepping approach works robustly in practice
and we are not aware of any examples where locally stable spurious solutions
arise once $\rho$ is chosen according to our conjecture. (An example that
smaller values can result in convergence to spurious solutions, instead, is
given in [98, Ch. 3]). Furthermore, in particular implementations one can
simply choose $\rho$ as conjectured, check if timestepping returns a spurious
solution by solving the eigenvalue problem 3.10 and, in that case, repeat the
calculation with a larger $\rho$.

## 5 Conclusions

Despite being established over 30 years ago, the background method remains one
of the main tools to place rigorous bounds on time-averaged properties of
turbulent flows. Here, we have reviewed a recent interpretation of the method
as an application of the auxiliary function framework [74, 75, 79] with
quadratic auxiliary functions in the form 3.3, which are parametrized by
generalized background fields and balance parameters. This interpretation
allows one to systematically formulate convex variational principles to search
for bounds on mean quantities that, contrary to the traditional scope of the
background method [81], need not be equivalent to the mean dissipation in the
flow. We have also shown that symmetries can be used to simplify the choice of
(generalized) background fields—proving, in particular, that nonzero
background velocity fields cannot improve existing bounds on the heat
transport in Rayleigh–Bénard convection—and that additional constraints, such
as maximum principles, can be enforced using Lagrange multipliers. We hope
that our general presentation using the abstract flow description of section 2
can provide a “recipe” for future applications of the background method to a
wider variety of flows, and offer a basis for future refinements and
extensions.

We have also provided a general description of two approaches for the
computational implementation of the background method, which have been applied
successfully to a variety of shear and convective flows in horizontally
periodic domains (see Table 1 for references). The first approach, based on
semidefinite programming, is extremely flexible, has general convergence
guarantees, and can be implemented using general-purpose interior-point SDP
solvers that are available open-source. To use such solvers in turbulent
regimes that require accurate discretization of thin boundary layers, however,
one must employ advanced decomposition techniques that have only recently been
developed by the optimization community (see [97] for a review). The second
approach to optimizing background fields we have reviewed, instead, relies on
timestepping. This is generally more familiar to the fluid mechanics community
and can more easily reach highly turbulent regimes, but can converge slowly
and presently lacks general theoretical convergence guarantees beyond the
particular ones in [35]. Table 2 summarizes the main advantages and drawbacks
of the two methods, in the hope of helping readers choose which one to use in
particular applications.

Table 2: Comparison between the two approaches to numerically optimize background fields reviewed in this work. For each property, checkmarks indicate which approach should currently be preferred. |  | SDP (interior-point solvers) | Timestepping
---|---|---|---
| Core algorithms known in fluid mechanics |  | ✓
| Scalability to highly turbulent regimes |  | ✓
| Good practical convergence | ✓ | ✓
| Good theoretical convergence | ✓ |
| Easy inclusion of additional constraints | ✓ |

Given the poor scalability of current interior-point SDP solvers and the
possibly slow convergence of the timestepping approach in dynamically
complicated flow regimes (see remark 6), further progress is needed before
complex flows of industrial relevance can be tackled robustly and at a
reasonable computational cost. We expect that significant efficiency gains may
be achieved by exploiting sparsity in SDPs via matrix decomposition, as
outlined in section 4.1, as well as the newly established connection between
the timestepping approach of section 4.2 and low-rank factorization strategies
for SDPs. Particular open questions that warrant deeper investigation include
whether the (restrictive) technical assumptions required by chordal matrix
decomposition can be weakened, and whether sophisticated analysis techniques
and optimization algorithms for SDPs with low-rank solutions can guarantee
fast convergence to the optimal bounds whilst provably avoiding spurious
solutions.

Despite the many recent advances, therefore, the numerical analysis of the
background method remains in our opinion far from complete and offers numerous
opportunities for interesting research. Pursuing these opportunities promises
to extend the range of flows for which bounds on mean quantities can be
computed efficiently. Moreover, and perhaps more excitingly, any lessons
learnt in the process may reveal practical ways to optimize auxiliary
functions more general than the quadratic ones underpinning the background
method. These can in principle produce arbitrarily sharp bounds on the
properties of turbulent flows [78], but currently remain beyond reach.

Acknowledgements. We are indebted to Charlie Doering for encouraging our
numerical investigation of the background method with his characteristic
enthusiasm, generously hosting AW in 2013 at the University of Michigan and
through fruitful discussions with GF at the Woods Hole Oceanographic
Institution in 2015, 2017 and 2018. GF is also grateful to Antonis
Papachristodoulou and Paul Goulart for the introduction to the world of
structured SDPs, without which the connections made in section 4.3 would not
have been possible. Finally, we thank two anonymous reviewers and the
Associate Editor for providing comments that helped us improve the original
manuscript. GF was supported by an Imperial College Research Fellowship. AA is
funded by the EPSRC Centre for Doctoral Training in Fluid Dynamics across
Scales (award number EP/L016230/1).

## References

  * Malkus [1954] W. V. R. Malkus.  The heat transport and spectrum of thermal turbulence.  _Proc. Royal Society A_ , 225(1161):196–212, 1954\.  (doi:10.1098/rspa.1954.0197).
  * Howard [1963] L. N. Howard.  Heat transport by turbulent convection.  _J. Fluid Mech._ , 17(3):405–432, 1963.  (doi:10.1017/S0022112063001427).
  * Busse [1969] F. H. Busse.  On Howard’s upper bound for heat transport by turbulent convection.  _J. Fluid Mech._ , 37(3):457–477, 1969.  (doi:10.1017/S0022112069000668).
  * Busse [1970] F. H. Busse.  Bounds for turbulent shear flow.  _J. Fluid Mech._ , 41(1):219–240, 1970.  (doi:10.1017/S0022112070000599).
  * Howard [1972] L. N. Howard.  Bounds on flow quantities.  _Annu. Rev. Fluid Mech._ , 4(1):473–494, 1972\.  (doi:10.1146/annurev.fl.04.010172.002353).
  * Busse [1979] F. H. Busse.  The optimum theory of turbulence.  _Adv. Appl. Mech._ , 18:77–121, 1979.  (doi:10.1016/S0065-2156(08)70265-5).
  * Doering and Constantin [1992] C. R. Doering and P. Constantin.  Energy dissipation in shear driven turbulence.  _Phys. Rev. Lett._ , 69(11):1648–1651, 1992.  (doi:10.1103/PhysRevLett.69.1648).
  * Doering and Constantin [1994] C. R. Doering and P. Constantin.  Variational bounds on energy dissipation in incompressible flows: Shear flow.  _Phys. Rev. E_ , 49(5):4087–4099, 1994.  (doi:10.1103/PhysRevE.49.4087).
  * Constantin and Doering [1995a] P. Constantin and C. R. Doering.  Variational bounds on energy dissipation in incompressible flows. II. Channel flow.  _Phys. Rev. E_ , 51(4):3192–3198, 1995a.  (doi:10.1103/PhysRevE.51.3192).
  * Doering and Constantin [1996] C. R. Doering and P. Constantin.  Variational bounds on energy dissipation in incompressible flows. III. Convection.  _Phys. Rev. E_ , 53(6):5957–5981, 1996.  (doi:10.1103/PhysRevE.53.5957).
  * Kerswell [1998] R. R. Kerswell.  Unification of variational principles for turbulent shear flows: the background method of Doering–Constantin and the mean-fluctuation formulation of Howard–Busse.  _Phys. D_ , 121(1-2):175–192, 1998.  (doi:10.1016/S0167-2789(98)00104-3).
  * Kerswell [1999] R. R. Kerswell.  Variational principle for the Navier–Stokes equations.  _Phys. Rev. E_ , 59(5):5482–5494, 1999.  (doi:10.1103/PhysRevE.59.5482).
  * Kerswell [2001] R. R. Kerswell.  New results in the variational approach to turbulent Boussinesq convection.  _Phys. Fluids_ , 13(1):192–209, 2001.  (doi:10.1063/1.1327295).
  * Constantin and Doering [1996] P. Constantin and C. R. Doering.  Heat transfer in convective turbulence.  _Nonlinearity_ , 9(4):1049–1060, 1996.  (doi:10.1088/0951-7715/9/4/013).
  * Kerswell [1997] R. R. Kerswell.  Variational bounds on shear-driven turbulence and turbulent Boussinesq convection.  _Phys. D_ , 100(3–4):355–376, 1997.  (doi:10.1016/S0167-2789(96)00227-8).
  * Constantin and Doering [1999] P. Constantin and C. R. Doering.  Infinite Prandtl number convection.  _J. Stat. Phys._ , 94(1-2):159–172, 1999.  (doi:10.1023/A:1004511312885).
  * Doering and Constantin [2001] C. R. Doering and P. Constantin.  On upper bounds for infinite Prandtl number convection with or without rotation.  _J. Math. Phys._ , 42(2):784–795, 2001.  (doi:10.1063/1.1336157).
  * Otero et al. [2002] J. Otero, R. W. Wittenberg, R. A. Worthing, and C. R. Doering.  Bounds on Rayleigh–Bénard convection with an imposed heat flux.  _J. Fluid Mech._ , 473:191–199, 2002.  (doi:10.1017/S0022112002002410).
  * Yan [2004] X. Yan.  On limits to convective heat transport at infinite Prandtl number with or without rotation.  _J. Math. Phys._ , 45(7):2718–2743, 2004.  (doi:10.1063/1.1763246).
  * Plasting and Ierley [2005] S. C. Plasting and G. R. Ierley.  Infinite-Prandtl-number convection. Part 1. Conservative bounds.  _J. Fluid Mech._ , 542(2005):343–363, 2005.  (doi:10.1017/S0022112005006555).
  * Ierley et al. [2006] G. R. Ierley, R. R. Kerswell, and S. C. Plasting.  Infinite-Prandtl-number convection. Part 2. A singular limit of upper bound theory.  _J. Fluid Mech._ , 560:159–227, 2006.  (doi:10.1017/S0022112006000450).
  * Doering et al. [2006] C. R. Doering, F. Otto, and M. G. Reznikoff.  Bounds on vertical heat transport for infinite Prandtl number Rayleigh–Bénard convection.  _J. Fluid Mech._ , 560:229–241, 2006.  (doi:10.1017/S0022112006000097).
  * Wittenberg [2010] R. W. Wittenberg.  Bounds on Rayleigh–Bénard convection with imperfectly conducting plates.  _J. Fluid Mech._ , 665:158–198, 2010.  (doi:10.1017/S0022112010003897).
  * Whitehead and Doering [2011a] J. P. Whitehead and C. R. Doering.  Ultimate state of two-dimensional Rayleigh–Bénard convection between free-slip fixed-temperature boundaries.  _Phys. Rev. Lett._ , 106(24):244501, 2011a.  (doi:10.1103/PhysRevLett.106.244501).
  * Otto and Seis [2011] F. Otto and C. Seis.  Rayleigh–Bénard convection: Improved bounds on the Nusselt number.  _J. Math. Phys._ , 52(8):083702, 2011.  (doi:10.1063/1.3623417).
  * Whitehead and Doering [2012] J. P. Whitehead and C. R. Doering.  Rigid bounds on heat transport by a fluid between slippery boundaries.  _J. Fluid Mech._ , 707:241–259, 2012.  (doi:10.1017/jfm.2012.274).
  * Whitehead and Wittenberg [2014] J. P. Whitehead and R. W. Wittenberg.  A rigorous bound on the vertical transport of heat in Rayleigh–Bénard convection at infinite Prandtl number with mixed thermal boundary conditions.  _J. Math. Phys._ , 55(9):093104, 2014.  (doi:10.1063/1.4896223).
  * Goluskin and Doering [2016] D. Goluskin and C. R. Doering.  Bounds for convection between rough boundaries.  _J. Fluid Mech._ , 804:370–386, 2016.  (doi:10.1017/jfm.2016.528).
  * Nobili and Otto [2017] C. Nobili and F. Otto.  Limitations of the background field method applied to Rayleigh–Bénard convection.  _J. Math. Phys._ , 58(9):093102, 2017.  (doi:10.1063/1.5002559).
  * Fantuzzi [2018a] G. Fantuzzi.  Bounds for Rayleigh–Bénard convection between free-slip boundaries with an imposed heat flux.  _J. Fluid Mech._ , 837:R5, 2018a.  (doi:10.1017/jfm.2017.907).
  * Pachev et al. [2020] B. Pachev, J. P. Whitehead, G. Fantuzzi, and I. Grooms.  Rigorous bounds on the heat transport of rotating convection with Ekman pumping.  _J. Math. Phys._ , 61(2):023101, 2020.  (doi:10.1063/1.5134054).
  * Christopher and Llewellyn Smith [2021] T. W. Christopher and S. G. Llewellyn Smith.  Bounding temperature dissipation in time-modulated Rayleigh–Bénard convection.  _Phys. Rev. Fluids_ , 6(5):1–9, 2021.  (doi:10.1103/PhysRevFluids.6.L051501).
  * Doering and Hyman [1997] C. R. Doering and J. M. Hyman.  Energy stability bounds on convective heat transport: Numerical study.  _Phys. Rev. E_ , 55(6):7775–7778, 1997.  (doi:10.1103/PhysRevE.55.7775).
  * Wittenberg and Gao [2010] R. W. Wittenberg and J. Gao.  Conservative bounds on Rayleigh–Bénard convection with mixed thermal boundary conditions.  _Eur. Phys. J. B_ , 76(4):565–580, 2010.  (doi:10.1140/epjb/e2010-00227-x).
  * Wen et al. [2015] B. Wen, G. P. Chini, R. R. Kerswell, and C. R. Doering.  Time-stepping approach for solving upper-bound problems: Application to two-dimensional Rayleigh–Bénard convection.  _Phys. Rev. E_ , 92(4):043012, 2015.  (doi:10.1103/PhysRevE.92.043012).
  * Tilgner [2017] A. Tilgner.  Bounds on poloidal kinetic energy in plane layer convection.  _Phys. Rev. Fluids_ , 2(12):123502, 2017.  (doi:10.1103/PhysRevFluids.2.123502).
  * Tilgner [2019] A. Tilgner.  Time evolution equation for advective heat transport as a constraint for optimal bounds in Rayleigh–Bénard convection.  _Phys. Rev. Fluids_ , 4(1):1–11, 2019.  (doi:10.1103/PhysRevFluids.4.014601).
  * Ding and Kerswell [2020] Z. Ding and R. R. Kerswell.  Exhausting the background approach for bounding the heat transport in Rayleigh–Bénard convection.  _J. Fluid Mech._ , 889:A33(1–33), 2020.  (doi:10.1017/jfm.2020.41).
  * Plasting and Kerswell [2003] S. C. Plasting and R. R. Kerswell.  Improved upper bound on the energy dissipation rate in plane Couette flow: the full solution to Busse’s problem and the Constantin–Doering–Hopf problem with one-dimensional background field.  _J. Fluid Mech._ , 477:363–379, 2003.  (doi:10.1017/S0022112002003361).
  * Hagstrom and Doering [2010] G. I. Hagstrom and C. R. Doering.  Bounds on heat transport in Bénard–Marangoni convection.  _Phys. Rev. E_ , 81(4):047301, 2010.  (doi:10.1103/PhysRevE.81.047301).
  * Fantuzzi et al. [2020] G. Fantuzzi, C. Nobili, and A. Wynn.  New bounds on the vertical heat transport for Bénard–Marangoni convection at infinite Prandtl number.  _J. Fluid Mech._ , 885:R4(1—-12), 2020.  (doi:10.1017/jfm.2019.1029).
  * Fantuzzi et al. [2018] G. Fantuzzi, A. Pershin, and A. Wynn.  Bounds on heat transfer for Bénard–Marangoni convection at infinite Prandtl number.  _J. Fluid Mech._ , 837:562–596, 2018.  (doi:10.1017/jfm.2017.858).
  * Doering and Constantin [1998] C. R. Doering and P. Constantin.  Bounds for heat transport in a porous layer.  _J. Fluid Mech._ , 376:263–296, 1998.  (doi:10.1017/S002211209800281X).
  * Otero et al. [2004] J. Otero, L. A. Dontcheva, H. Johnston, R. A. Worthing, A. Kurganov, G. Petrova, and C. R. Doering.  High-Rayleigh-number convection in a fluid-saturated porous layer.  _J. Fluid Mech._ , 500:263–281, 2004.  (doi:10.1017/S0022112003007298).
  * Wen et al. [2012] B. Wen, N. Dianati, E. Lunasin, G. P. Chini, and C. R. Doering.  New upper bounds and reduced dynamical modeling for Rayleigh–Bénard convection in a fluid saturated porous layer.  _Commun. Nonlinear Sci. Numer. Simul._ , 17(5):2191–2199, 2012.  (doi:10.1016/j.cnsns.2011.06.039).
  * Wen et al. [2013] B. Wen, G. P. Chini, N. Dianati, and C. R. Doering.  Computational approaches to aspect-ratio-dependent upper bounds and heat flux in porous medium convection.  _Phys. Lett. A_ , 377(41):2931–2938, 2013.  (doi:10.1016/j.physleta.2013.09.009).
  * Wen and Chini [2018] B. Wen and G. P. Chini.  Inclined porous medium convection at large Rayleigh number.  _J. Fluid Mech._ , 837:670–702, 2018.  (doi:10.1017/jfm.2017.863).
  * Lu et al. [2004] L. Lu, C. R. Doering, and F. H. Busse.  Bounds on convection driven by internal heating.  _J. Math. Phys._ , 45(7):2967–2986, 2004.  (doi:10.1063/1.1760845).
  * Whitehead and Doering [2011b] J. P. Whitehead and C. R. Doering.  Internal heating driven convection at infinite Prandtl number.  _J. Math. Phys._ , 52(9):093101, 2011b.  (doi:10.1063/1.3637032).
  * Goluskin [2015] D. Goluskin.  Internally heated convection beneath a poor conductor.  _J. Fluid Mech._ , 771:36–56, 2015.  (doi:10.1017/jfm.2015.140).
  * Arslan et al. [2021a] A. Arslan, G. Fantuzzi, J. Craske, and A. Wynn.  Bounds on heat transport for convection driven by internal heating.  _J. Fluid Mech._ , 919:A15, 2021a.  (doi:10.1017/jfm.2021.360).
  * Arslan et al. [2021b] A. Arslan, G. Fantuzzi, J. Craske, and A. Wynn.  Bounds for internally heated convection with fixed boundary heat flux.  _J. Fluid Mech._ , 992:R1, 2021b.  (doi:10.1017/jfm.2021.527).
  * Kumar et al. [2021] A. Kumar, A. Arslan, G. Fantuzzi, J. Craske, and A. Wynn.  Analytical bounds on the heat transport in internally heated convection.  arXiv:2110.10344 [physics.flu-dyn], 2021.
  * Balmforth et al. [2006] N. J. Balmforth, S. A. Ghadge, A. Kettapun, and S. D. Mandre.  Bounds on double-diffusive convection.  _J. Fluid Mech._ , 569:29–50, 2006.  (doi:10.1017/S0022112006002230).
  * Siggers et al. [2004] J. H. Siggers, R. R. Kerswell, and N. J. Balmforth.  Bounds on horizontal convection.  _J. Fluid Mech._ , 517:55–70, 2004.  (doi:10.1017/S0022112004000497).
  * Marchioro [1994] C. Marchioro.  Remark on the energy dissipation in shear driven turbulence.  _Phys. D_ , 74(3–4):395–398, 1994.  (doi:10.1016/0167-2789(94)90203-8).
  * Nicodemus et al. [1997a] R. Nicodemus, S. Grossmann, and M. Holthaus.  Improved variational principle for bounds on energy dissipation in turbulent shear flow.  _Phys. D_ , 101(1-2):178–190, 1997a.  (doi:10.1016/S0167-2789(96)00210-2).
  * Nicodemus et al. [1998a] R. Nicodemus, S. Grossmann, and M. Holthaus.  The background flow method. Part 2. Asymptotic theory of dissipation bounds.  _J. Fluid Mech._ , 363:301–323, 1998a.  (doi:10.1017/S0022112098001177).
  * Doering et al. [2000] C. R. Doering, E. A. Spiegel, and R. A. Worthing.  Energy dissipation in a shear layer with suction.  _Phys. Fluids_ , 12(8):1955–1968, 2000.  (doi:10.1063/1.870443).
  * Kerswell [2002] R. R. Kerswell.  Upper bounds on general dissipation functionals in turbulent shear flows: revisiting the ‘efficiency’ functional.  _J. Fluid Mech._ , 461:239–275, 2002.  (doi:10.1017/S0022112002008303).
  * Hagstrom and Doering [2014] G. I. Hagstrom and C. R. Doering.  Bounds on surface stress-driven shear flow.  _J. Nonlinear Sci._ , 24(1):185–199, 2014.  (doi:10.1007/s00332-013-9183-4).
  * Nicodemus et al. [1997b] R. Nicodemus, S. Grossmann, and M. Holthaus.  Variational bound on energy dissipation in plane Couette flow.  _Phys. Rev. E_ , 56(6):6774–6786, 1997b.  (doi:10.1103/PhysRevE.56.6774).
  * Nicodemus et al. [1998b] R. Nicodemus, S. Grossmann, and M. Holthaus.  The background flow method. Part 1. Constructive approach to bounds on energy dissipation.  _J. Fluid Mech._ , 363:281–300, 1998b.  (doi:10.1017/S0022112098001165).
  * Tang et al. [2004] W. Tang, C.-C. P. Caulfield, and W. R. Young.  Bounds on dissipation in stress-driven flow.  _J. Fluid Mech._ , 510:333–352, 2004.  (doi:10.1017/S0022112004009589).
  * Fantuzzi and Wynn [2016a] G. Fantuzzi and A. Wynn.  Optimal bounds with semidefinite programming: An application to stress-driven shear flows.  _Phys. Rev. E_ , 93(4):043308, 2016a.  (doi:10.1103/PhysRevE.93.043308).
  * Lee et al. [2019] H. Lee, B. Wen, and C. R. Doering.  Improved upper bounds on the energy dissipation rate for shear flow with injection and suction.  _Phys. Fluids_ , 31(8):085102, 2019.  (doi:10.1063/1.5109059).
  * Constantin [1994] P. Constantin.  Geometric statistics in turbulence.  _SIAM Rev._ , 36(1):73–98, 1994.  (doi:10.1137/1036004).
  * Gallet et al. [2010] B. Gallet, C. R. Doering, and E. A. Spiegel.  Destabilizing Taylor-Couette flow with suction.  _Phys. Fluids_ , 22(3):034105, 2010.  (doi:10.1063/1.3319824).
  * Ding and Marensi [2019] Z. Ding and E. Marensi.  Upper bound on angular momentum transport in Taylor–Couette flow.  _Phys. Rev. E_ , 100:063109, 2019.  (doi:10.1103/PhysRevE.100.063109).
  * Kumar [2020] A. Kumar.  Pressure-driven flows in helical pipes: Bounds on flow rate and friction factor.  _J. Fluid Mech._ , 904:A5, 2020.  (doi:10.1017/jfm.2020.655).
  * Kerswell [1996] R. R. Kerswell.  Upper bounds on the energy dissipation in turbulent precession.  _J. Fluid Mech._ , 321:335–370, 1996.  (doi:10.1017/S0022112096007756).
  * Kumar and Garaud [2020] A. Kumar and P. Garaud.  Bound on the drag coefficient for a flat plate in a uniform flow.  _J. Fluid Mech._ , 900:A6(1–24), 2020.  (doi:10.1017/jfm.2020.477).
  * Tilgner [2021] A. Tilgner.  A rigorous bound on the scaling of dissipation with velocity amplitude in flow past a sphere.  _J. Fluid Mech._ , 916:1–10, 2021.  (doi:10.1017/jfm.2021.239).
  * Chernyshenko et al. [2014] S. I. Chernyshenko, P. J. Goulart, D. Huang, and A. Papachristodoulou.  Polynomial sum of squares in fluid dynamics: a review with a look ahead.  _Philos. Trans. Roy. Soc. A_ , 372(2020):20130350, 2014.  (doi:10.1098/rsta.2013.0350).
  * Fantuzzi et al. [2016] G. Fantuzzi, D. Goluskin, D. Huang, and S. I. Chernyshenko.  Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization.  _SIAM J. Appl. Dyn. Syst._ , 15(4):1962–1988, 2016.  (doi:10.1137/15M1053347).
  * Goluskin and Fantuzzi [2019] D. Goluskin and G. Fantuzzi.  Bounds on mean energy in the Kuramoto–Sivashinsky equation computed using semidefinite programming.  _Nonlinearity_ , 32(5):1705–1730, 2019.  (doi:10.1088/1361-6544/ab018b).
  * Tobasco et al. [2018] I. Tobasco, D. Goluskin, and C. R. Doering.  Optimal bounds and extremal trajectories for time averages in nonlinear dynamical systems.  _Phys. Lett. A_ , 382(6):382–386, 2018.  (doi:10.1016/j.physleta.2017.12.023).
  * Rosa and Temam [2020] R. M. S. Rosa and R. M. Temam.  Optimal minimax bounds for time and ensemble averages of dissipative infinite-dimensional systems with applications to the incompressible Navier–Stokes equations.  arXiv:2010.06730 [math.AP], 2020\.
  * Chernyshenko [2017] S. I. Chernyshenko.  Relationship between the methods of bounding time averages.  arXiv:1704.02475 [phys.flu-dyn], 2017.
  * Fantuzzi and Wynn [2015] G. Fantuzzi and A. Wynn.  Construction of an optimal background profile for the Kuramoto–Sivashinsky equation using semidefinite programming.  _Phys. Lett. A_ , 379(1-2):23–32, 2015.  (doi:10.1016/j.physleta.2014.10.039).
  * Constantin and Doering [1995b] P. Constantin and C. R. Doering.  Variational bounds in dissipative systems.  _Phys. D_ , 82(3):221–228, 1995b.  (doi:10.1016/0167-2789(94)00237-K).
  * Goluskin [2016] D. Goluskin.  _Internally heated convection and Rayleigh–Bénard convection_.  Springer Briefs in Applied Sciences and Technologies.  (doi:10.1007/978-3-319-23941-5).
  * Foias et al. [1987] C. Foias, O. Manley, and R. Temam.  Attractors for the Bénard problem: existence and physical bounds on their fractal dimension.  _Nonlinear Anal._ , 11(8):939–967, 1987.  (doi:10.1016/0362-546X(87)90061-7).
  * Fantuzzi and Wynn [2016b] G. Fantuzzi and A. Wynn.  Semidefinite relaxation of a class of quadratic integral inequalities.  In: IEEE, _Proc. 55th IEEE Conference on Decision and Control_ , vol. 2, pp. 6192–6197, Las Vegas, USA, 2016b.  (doi:10.1109/CDC.2016.7799221).
  * Fantuzzi et al. [2017] G. Fantuzzi, A. Wynn, P. J. Goulart, and A. Papachristodoulou.  Optimization with affine homogeneous quadratic integral inequality constraints.  _IEEE Trans. Automat. Control_ , 62(12):6221–6236, 2017.  (doi:10.1109/TAC.2017.2703927).
  * Nesterov and Nemirovskii [1994] Y. Nesterov and A. Nemirovskii.  _Interior-Point Polynomial Algorithms in Convex Programming_.  (doi:10.1137/1.9781611970791).
  * Vandenberghe and Boyd [1996] L. Vandenberghe and S. Boyd.  Semidefinite Programming.  _SIAM Rev._ , 38(1):49–95, 1996.  (doi:10.1137/1038003).
  * Nemirovski [2006] A. Nemirovski.  Advances in convex optimization: Conic programming.  In: _International Congress of Mathematicians_ , vol. 1, pp. 413–444, 2006.  (doi:10.4171/022-1/17).
  * O’Donoghue et al. [2016] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd.  Conic optimization via operator splitting and homogeneous self-dual embedding.  _J. Optim. Theory Appl._ , 169(3):1–27, 2016\.  (doi:10.1007/s10957-016-0892-3).
  * Wen et al. [2010] Z. Wen, D. Goldfarb, and W. Yin.  Alternating direction augmented Lagrangian methods for semidefinite programming.  _Math. Program. Comput._ , 2(3-4):203–230, 2010\.  (doi:10.1007/s12532-010-0017-1).
  * Gambill [2006] T. N. Gambill.  _Application of uncertainty inequalities to bound the radius of the attractor for the Kuramoto–Sivashinsky equation_.  PhD thesis, University of Illinois at Urbana-Champaign, 2006.  URL: https://www.ideals.illinois.edu/handle/2142/86875.
  * Nesterov [2003] Y. Nesterov.  _Introductory Lectures on Convex Programming_.  (doi:10.1007/978-1-4419-8853-9).
  * Fukuda et al. [2000] M. Fukuda, M. Kojima, K. Murota, and K. Nakata.  Exploiting sparsity in semidefinite programming via matrix completion I: General Framework.  _SIAM J. Optim._ , 11(3):647–674, 2000.  (doi:10.1137/S1052623400366218).
  * Nakata et al. [2003] K. Nakata, K. Fujisawa, M. Fukuda, M. Kojima, and K. Murota.  Exploiting sparsity in semidefinite programming via matrix completion II: implementation and numerical results.  _Math. Program. B_ , 95(2):303–327, 2003.  (doi:10.1007/s10107-002-0351-9).
  * Vandenberghe and Andersen [2015] L. Vandenberghe and M. S. Andersen.  Chordal graphs and semidefinite optimization.  _Found. Trends Optim._ , 1(4):241–433, 2015.  (doi:10.1561/2400000006).
  * Zheng et al. [2020] Y. Zheng, G. Fantuzzi, A. Papachristodoulou, P. J. Goulart, and A. Wynn.  Chordal decomposition in operator-splitting methods for sparse semidefinite programs.  _Math. Program._ , 180:489–532, 2020.  (doi:10.1007/s10107-019-01366-3).
  * Zheng et al. [2021] Y. Zheng, G. Fantuzzi, and A. Papachristodoulou.  Chordal and factor-width decompositions for scalable semidefinite and polynomial optimization.  arXiv:2107.02379 [math.OC], 2021\.
  * Fantuzzi [2018b] G. Fantuzzi.  _Construction of optimal background fields using semidefinite programming_.  Phd thesis, Imperial College London, 2018b.  URL: https://spiral.imperial.ac.uk/handle/10044/1/60642.
  * Boyd and Vandenberghe [2004] S. Boyd and L. Vandenberghe.  _Convex Optimization_.  (doi:10.1017/CBO9780511804441).
  * Burer et al. [2002] S. Burer, R. D. C. Monteiro, and Y. Zhang.  Solving a class of semidefinite programs via nonlinear programming.  _Math. Program. A_ , 93(1):97–122, 2002.  (doi:10.1007/s101070100279).
  * Burer and Monteiro [2003] S. Burer and R. D. C. Monteiro.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  _Math. Program. B_ , 95(2):329–357, 2003.  (doi:10.1007/s10107-002-0352-8).
  * Burer and Monteiro [2005] S. Burer and R. D. C. Monteiro.  Local minima and convergence in low-rank semidefinite programming.  _Math. Program. A_ , 103(3):427–444, 2005.  (doi:10.1007/s10107-004-0564-1).
  * Burer and Choi [2006] S. Burer and C. Choi.  Computational enhancements in low-rank semidefinite programming.  _Optim. Methods Softw._ , 21(3):493–512, 2006\.  (doi:10.1080/10556780500286582).
  * Boumal et al. [2020] N. Boumal, V. Voroninski, and A. S. Bandeira.  Deterministic guarantees for Burer–Monteiro factorizations of smooth semidefinite programs.  _Comm. Pure Appl. Math._ , 73(3):581–608, 2020\.  (doi:10.1002/cpa.21830).
  * Waldspurger and Waters [2020] I. Waldspurger and A. Waters.  Rank optimality for the Burer–Monteiro factorization.  _SIAM J. Optim._ , 30(3):2577–2602, 2020.  (doi:10.1137/19M1255318).
  * Pataki [1998] G. Pataki.  On the rank of extreme matrices in semidefinite programs and the multiplicity of optimal eigenvalues.  _Math. Oper. Res._ , 23(2):339–358, 1998.  (doi:10.1287/moor.23.2.339).