4.1. Bounds via auxiliary functional

A rigorous upper bound on $\langle wT \rangle$ can be derived using the simple observation that the infinite-time average of the time derivative of any bounded function $\mathcal {V}\{\boldsymbol {u},T\}$ along solutions of the governing equations (2.1a–c) vanishes, so

(4.1)

In particular, if $\mathcal {V}$ can be chosen such that

(4.2)

for all time along solutions of (2.1a–c) for some constant $U$, then $\langle wT \rangle \leqslant U$. The goal, then, is to construct a $\mathcal {V}$ that satisfies (

) with the smallest possible $U$.

While the evolution equations (2.1b) and (2.1c) cannot be solved explicitly for all possible initial conditions, when $\mathcal {V}$ is given they can be used to derive an explicit expression for $\mathcal {S}$ as a function of $\boldsymbol {u}$ and $T$ alone. Then, to ensure that (4.2) holds along solutions (2.1a–c) pointwise in time, it suffices to enforce that $\mathcal {S}\{\boldsymbol {u},T\}$ be non-negative when viewed as a functional on the space of time-independent incompressible velocity fields and temperature fields that satisfy the boundary conditions,

(4.3) \begin{equation} \mathcal{H} := \{ (\boldsymbol{u}, T) : \ ({{\rm 2.2}\textit{a},\textit{b}}),\ \text{horizontal periodicity and } \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 \}. \end{equation}

We therefore search for a function $\mathcal {V}$ and constant $U$ such that $\mathcal {S}\{\boldsymbol {u},T\} \geqslant 0$ for all velocity and temperature fields in $\mathcal {H}$. This key relaxation makes this approach tractable and, remarkably, it may not be overly conservative. In fact, if the governing equations in (2.1) were well posed (which is not currently known) and solutions eventually remained in a compact subset $\mathcal {K}$ of $\mathcal {H}$, then optimising $\mathcal {V}$ over a sufficiently general class of functions whilst imposing $\mathcal {S}\{\boldsymbol {u},T\} \geqslant 0$ for all $\boldsymbol {u}$ and $T$ in $\mathcal {K}$, would yield an upper bound $U$ exactly equal to the largest possible value of $\langle wT \rangle$ (Rosa & Temam Reference Rosa and Temam2020).

Unfortunately, the construction of such an optimal $\mathcal {V}$ is currently beyond the reach of both analytical and computational methods. Nevertheless, progress can be made if we restrict the search to quadratic $\mathcal {V}$ in the form

(4.4)

where the non-negative scalars $a$ and $b$ and the function $\psi (z)$ are to be optimised such that (4.2) holds for the smallest possible $U$. As shown by Chernyshenko (Reference Chernyshenko2017), this choice of $\mathcal {V}$ amounts to using the background method (Constantin & Doering Reference Constantin and Doering1995; Doering & Constantin Reference Doering and Constantin1994, Reference Doering and Constantin1996): the profile ${1}/{b}[\psi (z)+z-1]$ corresponds to the background temperature field, while the scalars $a$ and $b$ are the so-called balance parameters. Note that the $z-1$ term could be removed by redefining $\psi (z)$, but we isolate it to simplify the analysis in what follows. Note also that, due to the periodicity in the horizontal directions, a ‘symmetry reduction’ argument, following ideas in Goluskin & Fantuzzi (Reference Goluskin and Fantuzzi2019, appendix A), proves that there is no loss of generality in taking $\psi$ to depend only on the vertical coordinate $z$. Similarly, one can show that the upper bound on $\langle wT \rangle$ cannot be improved by adding to $\mathcal {V}$ a term $\boldsymbol {\phi } \boldsymbol {\cdot } \boldsymbol {u}$, proportional to the velocity field via a (rescaled) incompressible background velocity field $\boldsymbol {\phi }$.

To find an expression for $\mathcal {S}\{\boldsymbol {u},T\}$, we calculate the time derivative of the quadratic $\mathcal {V}$ in (4.4) using the governing equations (2.1b,c), substitute the resulting expression into (4.2), and integrate various terms by parts using incompressibility and the boundary conditions (2.2a,b) to arrive at

(4.5)

The best bound on $\langle wT \rangle$ that can be proved with a quadratic $\mathcal {V}$ is

(4.6)\begin{equation} \langle wT \rangle \leqslant \inf_{U,\psi(z), a, b} \left\{U :\ \mathcal{S}\{\boldsymbol{u},T\} \geqslant 0 \ \forall (\boldsymbol{u},T) \in \mathcal{H} \right\}. \end{equation}

The right-hand side of (4.6) is a linear optimisation problem because the optimisation variables, $U$, $\psi (z)$, $a$ and $b$, enter the constraint $\mathcal {S}\{\boldsymbol {u},T\} \geqslant 0$ and the cost $U$ linearly.

4.2. Fourier expansion

The analysis and numerical implementation of the minimisation problem in (4.6) can be considerably simplified by expanding the horizontally periodic velocity and temperature as Fourier series,

(4.7)\begin{equation} \begin{bmatrix} T(x,y,z)\\\boldsymbol{u}(x,y,z) \end{bmatrix} = \sum_{\boldsymbol{k}} \begin{bmatrix} \hat{T}_{\boldsymbol{k}}(z)\\ \hat{\boldsymbol{u}}_{\boldsymbol{k}}(z) \end{bmatrix} \exp({\textrm{i}(k_x x + k_y y)}). \end{equation}

The sums are over wavevectors $\boldsymbol {k}=(k_x,k_y)$ compatible with the horizontal periods $L_x$ and $L_y$. We denote the magnitude of each wavevector by $k = \sqrt {k_{\smash {x}}^2 + k_{\smash {y}}^2}$. The complex-valued Fourier amplitudes satisfy the complex-conjugate relations $\hat {\boldsymbol {u}}_{-\boldsymbol {k}}=\hat {\boldsymbol {u}}_{\boldsymbol {k}}^*$ and $\hat {T}_{-\boldsymbol {k}} = \hat {T}_{\boldsymbol {k}}^*$, as well as the no-slip and isothermal boundary conditions. Using incompressibility and writing $\hat {w}_{\boldsymbol {k}}$ for the vertical component of $\hat {\boldsymbol {u}}_{\boldsymbol {k}}$, these can be expressed as

(4.8a)\begin{gather} \hat{w}_{\boldsymbol{k}}(0) = \hat{w}_{\boldsymbol{k}}'(0) = \hat{w}_{\boldsymbol{k}}(1) = \hat{w}_{\boldsymbol{k}}'(1)= 0, \end{gather}

(4.8b)\begin{gather}\hat{T}_{\boldsymbol{k}}(0)=\hat{T}_{\boldsymbol{k}}(1) = 0. \end{gather}

After inserting the Fourier expansions (4.7) into (4.5) and applying standard estimates based on the incompressibility condition and Young's inequality to write the horizontal Fourier amplitudes $\hat {u}_{\boldsymbol {k}}$ and $\hat {v}_{\boldsymbol {k}}$ as a function of $\hat {w}_{\boldsymbol {k}}$, the functional $\mathcal {S}\{\boldsymbol {u},T\}$ can be estimated from below as

(4.9)\begin{equation} \mathcal{S}\{\boldsymbol{u},T\} \geqslant \mathcal{S}_{0}\{\hat{T}_0\} + \sum_{\boldsymbol{k} \neq \boldsymbol{0}} \mathcal{S}_{\boldsymbol{k}} \{\hat{w}_{\boldsymbol{k}},\hat{T}_{\boldsymbol{k}}\}, \end{equation}

where

(4.10) \begin{equation} \mathcal{S}_{0}\{\hat{T}_0\} = \int_{0}^{1} [b \vert \hat{T}_{0}' \vert^{2}+ (bz-\psi')\hat{T}_{0}' + \psi] \textrm{d}z + \psi(1)\hat{T}_{0}'(1)-(\psi(0)-1)\hat{T}_{0}'(0) + U - \tfrac12 \end{equation}

and for all $\boldsymbol{k} \neq \boldsymbol{0}$

(4.11) \begin{align} \mathcal{S}_{\boldsymbol{k}}\{\hat{w}_{\boldsymbol{k}},\hat{T}_{\boldsymbol{k}}\} &= \int_{0}^{1} \left[ \frac{a}{R}\left( \frac{1}{k^2} \left\vert \hat{w}_{\boldsymbol{k}}'' \right\vert^{2} + 2\left\vert \hat{w}_{\boldsymbol{k}}' \right\vert^{2} + k^{2}\left\vert \hat{w}_{\boldsymbol{k}} \right\vert^{2} \right) \right. \nonumber\\ &\quad \left. + b\vert \hat{T}_{\boldsymbol{k}}' \vert^{2} + bk^{2}\vert \hat{T}_{\boldsymbol{k}} \vert^{2} - (a-\psi')\hat{w}_{\boldsymbol{k}}\hat{T}_{\boldsymbol{k}}^{*} \vphantom{\frac{1}{k^2}}\right] \textrm{d}z. \end{align}

Equality holds in (4.9) if $\boldsymbol {u}$ has only one non-zero horizontal component because, in this case, the Fourier-transformed incompressibility condition yields an exact relation between $\hat {w}_{\boldsymbol {k}}$ and either $\hat {u}_{\boldsymbol {k}}$ or $\hat {v}_{\boldsymbol {k}}$, so Young's inequality is not needed to eliminate the latter.

Velocity and temperature fields with a single non-zero Fourier mode are admissible in the optimisation problem (4.6), so the right-hand side of (4.9) is non-negative if and only if each term is non-negative. Moreover, since the real and imaginary parts of the Fourier amplitudes $\hat {w}_{\boldsymbol {k}}$ and $\hat {T}_{\boldsymbol {k}}$ give identical and independent contributions to $\mathcal {S}_{\boldsymbol {k}}$, we may assume them to be real without loss of generality. Thus, we may replace the minimisation problem in (4.6) with

(4.12) \begin{equation} \left. \begin{aligned} \inf_{U,\psi(z),a,b} \quad & U\\ \text{subject to} \quad & \mathcal{S}_0\{\hat{T}_0\} \geqslant 0 & & \forall \hat{T}_0: ({{\rm 4.8}\textit{b}}),\\ & \mathcal{S}_{\boldsymbol{k}}\{\hat{w}_{\boldsymbol{k}},\hat{T}_{\boldsymbol{k}}\} \geqslant 0 & & \forall \hat{w}_{\boldsymbol{k}},\hat{T}_{\boldsymbol{k}}: ({{\rm 4.8}\textit{a},\textit{b}}), \quad \forall \boldsymbol{k}\neq 0.\!\!\!\!\!\!\!\ \end{aligned}\right\} \end{equation}

Any choice of $U$, $\psi (z)$, $a$ and $b$ satisfying the constraints yields a rigorous upper bound on the mean vertical convective heat flux $\langle wT \rangle$. Following established terminology, we refer to the inequalities on $\mathcal {S}_{\boldsymbol {k}}$ as the spectral constraints.

Just like (4.6), the optimisation problem in (4.12) is linear and its optimal solution can be approximated using efficient numerical algorithms after discretisation. For computational convenience, however, we simplify the spectral constraints by dropping all non-negative terms that depend explicitly on $k$. Specifically, we replace the spectral constraints with the stronger, but simpler, single condition

(4.13) \begin{equation} \tilde{\mathcal{S}}\{\hat{w}, \hat{T} \} :=\int_0^1 \frac{2a}{R}\left\vert \hat{w}' \right\vert^{2} + b\vert \hat{T}' \vert^{2} - (a-\psi')\hat{w}\hat{T} \,\textrm{d}z \geqslant 0 \quad \forall \hat{T}, \hat{w}:\ ({{\rm 4.8}\textit{a},\textit{b}}) \end{equation}

and solve

(4.14) \begin{equation} \left. \begin{aligned} \inf_{U,\psi(z),a,b} \quad & U\\ \text{subject to} \quad & \mathcal{S}_0\{\hat{T}_0\} \geqslant 0 & & \forall \hat{T}_0: ({{\rm 4.8}\textit{b}}),\\ & \tilde{\mathcal{S}}\{\hat{w},\hat{T}\} \geqslant 0 & & \forall \hat{w},\hat{T}: ({{\rm 4.8}\textit{a},\textit{b}})\!\!\!\!\!\!\!\ \end{aligned}\right\} \end{equation}

instead of (4.12). This simplification leads to suboptimal bounds on $\langle wT \rangle$ at a fixed R but, as discussed in Appendix B, still captures the qualitative behaviour of the optimal ones. On the other hand, considering the simplified spectral constraint (4.13) allows for significant computational savings when optimising bounds numerically, because it removes the need to consider a large set of wavenumbers and it enables implementation using simple piecewise-linear basis functions (cf. Appendix D). This allows for discretisation of (4.14) and its generalisation (4.25) derived in § 4.4 below on very fine meshes, which is essential to resolve sharp boundary layers in $\psi$ accurately.

4.3. Explicit formulation

To simplify the analysis (but not the numerical implementation) of (4.12) it is convenient to eliminate the explicit appearance of $U$. This can be done upon observing that, given any profile $\psi (z)$ and balance parameters $a$, $b$, the smallest $U$ for which $\mathcal {S}_0\{\hat {T}_0\}$ is non-negative over admissible $\hat {T}_0$ is

(4.15a)\begin{align} U^* &= \frac12 + \sup_{\substack{\hat{T}_{\boldsymbol{k}}(0)=0\\ \hat{T}_{\boldsymbol{k}}(1) = 0}} \left\{ -\int_{0}^{1} [b \vert \hat{T}_{0}' \vert^{2}+ (bz-\psi')\hat{T}_{0}' + \psi ]\, \textrm{d}\kern0.7pt z\right.\nonumber\\ &\qquad\qquad \qquad\ \left.\vphantom{\int_{0}^{1} } -\psi(1)\hat{T}_{0}'(1) +(\psi(0)-1)\hat{T}_{0}'(0) \right\}. \end{align}

By modifying any admissible $\hat {T}_0$ in infinitesimally thin layers near $z=0$ and $z=1$, it is possible to show that this value is finite if and only if

(4.15b,c)\begin{equation} \psi(0)=1 \quad \text{and} \quad \psi(1)= 0. \end{equation}

Then, the optimal temperature field in (4.15a) is given by

(4.16)\begin{equation} \hat{T}_0'(z) = \frac{\psi'(z)+1}{2b} - \frac{z}{2} + \frac14, \end{equation}

and we obtain

(4.17)\begin{equation} U^* = \frac12 + \frac{1}{4b}\left\| bz-\frac{b}{2} - \psi'(z) - 1\right\|_2^2 - \int_{0}^{1} \psi\, \textrm{d}z. \end{equation}

Thus, we may replace the minimisation problem (4.12) with the more explicit version

(4.18)\begin{equation} \left. \begin{aligned} \inf_{\psi(z),a,b} \quad & \frac12 + \frac{1}{4b}\left\| bz-\frac{b}{2} - \psi'(z) - 1\right\|_2^2 - \int_{0}^{1}\psi\, \textrm{d}z\\ \text{subject to} \quad & \psi(0)=1,\\ & \psi(1)=0,\\ & \mathcal{S}_{\boldsymbol{k}}\{\hat{w},\hat{T}\} \geqslant 0 \qquad\forall \hat{w},\hat{T}: ({{\rm 4.8}\textit{a},\textit{b}}). \end{aligned}\right\} \end{equation}

Note that, although this formulation is not suitable for numerical implementation because the cost function is not convex with respect to $b$, it is more convenient when attempting to prove an upper bound on $\langle wT \rangle$ analytically.

4.4. Restriction to non-negative temperature fields

The upper bounding principle derived above can be improved by imposing a minimum principle, which guarantees that temperature fields solving the Boussinesq equations (2.1a–c) with homogeneous boundary conditions (2.2b) are non-negative in the domain $\varOmega$ at large time. More precisely, Appendix A proves that the fluid's temperature is non-negative on $\varOmega$ at all times if it is so initially, and the negative part of the temperature decays exponentially quickly otherwise.

Since $\langle wT \rangle$ is determined by the long-time behaviour of the velocity and temperature fields, the minimum principle enables us to replace the upper bound (4.6) with

(4.19)\begin{equation} \langle wT \rangle \leqslant \inf_{U,\psi(z), a, b} \left\{U :\ \mathcal{S}\{\boldsymbol{u},T\} \geqslant 0 \quad \forall (\boldsymbol{u},T) \in \mathcal{H}_+ \right\}, \end{equation}

where the space $\mathcal {H}$ of admissible velocity and temperature fields has been replaced with its subset

(4.20)\begin{equation} \mathcal{H}_+ := \{ (\boldsymbol{u}, T) \in \mathcal{H}: \ T(\boldsymbol{x}) \geqslant 0 \text{ on } \varOmega \}. \end{equation}

The constraint in the modified optimisation problem (4.19) is clearly weaker than the original one in (4.6), so imposing the minimum principle allows for a better bound on $\langle wT \rangle$ in principle. This is indeed the case at large Rayleigh numbers, as shall be demonstrated by numerical results in §§ 5 and 6.

In order to impose the inequality constraint $\mathcal {S}\{\boldsymbol {u},T\} \geqslant 0$ for non-negative temperatures, but relax it when $T$ is negative on a non-negligible subset of the domain, we effectively use a Lagrange multiplier. Specifically, we search for a positive bounded linear functional $\mathcal {L}$ such that $\mathcal {S}\{\boldsymbol {u},T\} \geqslant \mathcal {L}\{T\}$ for all pairs $(\boldsymbol {u},T) \in \mathcal {H}$, which satisfy only the boundary condition and incompressibility. Indeed, the positivity of $\mathcal {L}$ implies $\mathcal {S}\{\boldsymbol {u},T\} \geqslant \mathcal {L}\{T\} \geqslant 0$ if $T$ is non-negative, as desired. When $T$ is negative on a subset of the domain, instead, $\mathcal {L}\{T\}$ need not be positive and the constraint on $\mathcal {S}\{\boldsymbol {u},T\}$ is relaxed.

Analysis in Appendix C proves that there is no loss of generality in taking

(4.21)

where $q:(0,1)\to \mathbb {R}$ is a non-decreasing function, square integrable but not necessarily continuous, to be optimised alongside the bound $U$, the profile $\psi (z)$, and the balance parameter $a,b$. If $q$ were differentiable, one could integrate by parts to obtain

(4.22)

and identify $q'$ as a standard Lagrange multiplier for the condition $T(\boldsymbol {x}) \geqslant 0$; working with (

) simply removes the differentiability requirement from $q$. Moreover, the value of $\mathcal {L}\{T\}$ in (

) does not change when $q$ is shifted by a constant by virtue of the vertical boundary conditions on $T$, so we may normalise $q$ such that

(4.23)\begin{equation} \int_0^1 q(z) \, \textrm{d}z = \psi(1)-\psi(0). \end{equation}

The upper bound (4.19) can therefore be replaced with

(4.24)

Observe that setting $q(z)=-1$ causes the integral to vanish because $T$ is zero at the top and bottom boundaries, so this choice of $q$ results in the upper bound (4.6) derived without the minimum principle for $T$.

The minimisation problem in (4.24) can be expanded using a Fourier series exactly as explained in § 4.2. After estimating the functionals $\mathcal {S}_{\boldsymbol {k}}$ with $\tilde {\mathcal {S}}$, one concludes that $\langle wT \rangle$ is bounded above by the optimal value of the following linear optimisation problem:

(4.25) \begin{equation} \left. \begin{aligned} \inf_{U,a,b,\psi(z),q(z)} \quad & U\\ \text{subject to} \quad & \mathcal{S}_0\{\hat{T}_0\} + \int_0^1 q(z) \hat{T}_0'(z) \,\textrm{d}z \geqslant 0 \quad\forall \,\hat{T}_0: \hat{T}_0(0)=0=\hat{T}_0(1),\\ & \tilde{\mathcal{S}}\{\hat{w},\hat{T}\} \geqslant 0 \quad\forall \hat{w},\hat{T}: ({{\rm 4.8}\textit{a},\textit{b}}),\\ & q(z) \text{ non-decreasing}. \end{aligned}\right\} \end{equation}

If one does not simplify the spectral constraints, one obtains a very similar problem where the inequality on $\tilde {\mathcal {S}}$ is replaced by the same $\boldsymbol {k}$-dependent inequalities appearing in (4.12). This problem gives a quantitatively better bound on $\langle wT \rangle$ but has a higher computational complexity than (4.25) and, as explained in Appendix B, we do not expect qualitative improvements in the behaviour of the bounds with R.

The constant $U$ can be eliminated from (4.25) by following the same procedure outlined in § 4.3 in order to obtain a minimisation problem that is more suitable for analysis, but less convenient for computations. Indeed, the normalisation condition (4.23) implies that the functional $\mathcal {S}_0\{\hat {T}_0\}$ is minimised when

(4.26)\begin{equation} \hat{T}_0'(z) = \frac{\psi'(z) - q(z)}{2b} - \frac{z}{2} + \frac14 \end{equation}

and the best $U$ for given $a$, $b$, $\psi (z)$ and $q(z)$ is

(4.27)\begin{equation} U^* = \frac12 + \frac{1}{4b}\left\| bz-\frac{b}{2} - \psi'(z) + q(z) \right\|_2^2 - \int_{0}^{1} \psi\, \textrm{d}z. \end{equation}

As one would expect, this expression reduces to (4.17) when $q(z)=-1$. It is also possible to show that, since we are only interested in non-negative $\hat {T}_0$, the conditions $\psi (0)=1$ and $\psi (1)=0$ in § 4.3 may be weakened into inequalities $\psi (0)\leqslant 1$ and $\psi (1)\leqslant 0$. Thus, when the minimum principle for the temperature is imposed, the optimisation problem (4.18) relaxes into

(4.28) \begin{equation} \left. \begin{aligned} \inf_{\psi(z),q(z),a,b} \quad & \frac12 + \frac{1}{4b}\left\| bz-\frac{b}{2} - \psi'(z) + q(z) \right\|_2^2 - \int_{0}^{1}\psi\, \textrm{d}z\\ \text{subject to} \quad & \psi(0) \leqslant 1,\, \psi(1) \leqslant 0,\\ & q(z) \text{ non-decreasing, (4.23)},\\ & \tilde{\mathcal{S}}\{\hat{w},\hat{T}\} \geqslant 0 \qquad\forall \hat{w},\hat{T}: ({{\rm 4.8}\textit{a},\textit{b}}). \end{aligned}\right\} \end{equation}

Finally, observe that setting $\psi (z)=0$, $q(z)=0$, $a=0$ and letting $b\to 0$ in this minimisation problem yields a sequence of feasible solutions with optimal cost approaching $1/2$ from above irrespective of R. Thus, the uniform bound $\langle wT \rangle \leqslant 1/2$ proved by Goluskin & Spiegel (Reference Goluskin and Spiegel2012) can be recovered within our approach by taking the auxiliary functional in (4.4) to be

(4.29)

which corresponds to the flow's potential energy measured with respect to the upper boundary. As shown in § 6, however, more general choices can yield better bounds.

5.1. Numerically optimal bounds

Figure 2 compares the numerically optimal upper bounds $U$ on the mean vertical heat transfer $\langle wT \rangle$ to the experimental data by Kulacki & Goldstein (Reference Kulacki and Goldstein1972) and the direct numerical simulation data by Goluskin & Spiegel (Reference Goluskin and Spiegel2012) and Goluskin & van der Poel (Reference Goluskin and van der Poel2016). The bounds were calculated by solving the minimisation problem (4.12) for a fluid layer with horizontal periods $L_x = L_y = 2$, and with the simplified problem (4.14), which is independent of wavevectors and, therefore, of $L_x$ and $L_y$. As expected, the bounds obtained with (4.12) are zero when the Rayleigh number is smaller than the energy stability limit ${\textit {R}}_E \approx 29\,723$, which differs slightly from the value $26\,927$ reported by Goluskin (Reference Goluskin2016) due to our choice of horizontal periods. They then increase monotonically with ${\textit {R}}$, showing the same qualitative behaviour as the bounds computed with the simplified problem (4.14), which reach the value of $1/2$ at ${\textit {R}} = 259\,032$. Both sets of results exceed the uniform upper bound $\langle wT \rangle \leqslant 1/2$ for sufficiently large Rayleigh numbers. The apparent contradiction is due to the fact that the uniform bound relies on the minimum principle for the temperature, which was disregarded in the formulation (4.12).

Figure 2.

(a) Optimal bounds on $\langle wT \rangle$ obtained by solving the wavenumber-dependent problem (4.12) with horizontal periods $L_x = L_y = 2$ (- - -, black dashed line) and the simplified problem (4.14) (—, black solid line). Also plotted are experiments by Kulacki & Goldstein (Reference Kulacki and Goldstein1972) ($\times$, green), two-dimensional DNSs by Goluskin & Spiegel (Reference Goluskin and Spiegel2012) ($\triangle$, red) and three-dimensional DNSs by Goluskin & van der Poel (Reference Goluskin and van der Poel2016) ($\square$, blue). Circles mark values of R at which optimal profiles $\psi (z)$ are plotted in panel (b). The dashed horizontal line (- - -, black dashed line) is the uniform bound $\langle wT \rangle \leqslant 1/2$. (b) Optimal profiles $\psi (z)$ at ${\textit {R}} = 10^4$ (—, orange solid line), $10^5$ (—, red solid line) and $10^6$ (—, brown solid line). (c) Optimal balance parameters $a$ (—, blue solid line, left axis) and $b$ (- - -, red dashed line, right axis) as a function of R.

Numerically optimal profiles of $\psi (z)$ for the simplified bounding problem (4.14) at selected Rayleigh numbers and the variation of the optimal balance parameters $a$ and $b$ with R are illustrated in figures 2(b) and 2(c) respectively. As expected from the structure of the indefinite term $(a-\psi ')\hat {w}\hat {T}$ of the functional $\tilde {\mathcal {S}}$, the derivative $\psi '$ in the bulk of the domain approaches the value of $a$ as R is raised and leads to the formation of two boundary layers. Moreover, the asymmetry of the boundary layers reflects qualitatively the asymmetry of the IH convection problem we are studying, which is characterised by a stable thermal stratification near the bottom boundary ($z=0$) and an unstable one near the top ($z=1$). However, note that while $\psi$ is related to the background temperature field, it is not a physical quantity and need not behave nor scale like the mean temperature in turbulent convection.

Other insightful observations can be made by considering the critical temperature fields $\hat {T}_{0}$, which minimise the functional $\mathcal {S}_0\{\hat {T}_0\}$ for the optimal choice of $\psi$, $a$, $b$ and $U$. These critical temperatures can be recovered upon integrating (4.16), and are plotted in figure 3 for a selection of Rayleigh numbers. As one might expect, when ${\textit {R}}$ is sufficiently small such that $U=0$, $\hat {T}_{0}=\frac {1}{2}z(1-z)$ coincides with the conductive temperature profile. With the onset of convection and increasing R, boundary layers form at $z=0$ and $z=1$ and the maximum of $|\hat {T}_0|$ decreases. The profiles are also consistent with the uniform rigorous bound $\langle T \rangle \leqslant 1/12$ (Goluskin Reference Goluskin2016). However, for sufficiently high Rayleigh numbers they are evidently not related to the horizontal and infinite-time averages of the physical temperature field, because they become negative near $z=0$. Interestingly, as shown in figure 3(d,e), this unphysical behaviour first occurs away from the boundary at ${\textit {R}} = 256\,269$, while the numerical upper bound on $\langle wT \rangle$ reaches the value of $1/2$ only at $259\,032$, when $\hat {T}_0'(0)=0$. The latter is not surprising because the identity $\bar {T}'(0)=1/2-\langle wT\rangle$, derived from (1.1a,b) upon recognising that $\mathcal {F}_0 = \bar {T}'(0)$, implies that an upper bound of 1/2 on $\langle wT \rangle$ is equivalent to a zero lower bound on $\bar {T}'(0)$, which is obtained when $\hat {T}_0'(0)$ vanishes. It is therefore clear that the bounding problems (4.12) and (4.14) fail to improve the uniform bound $\langle wT \rangle \leqslant 1/2$ proved by Goluskin & Spiegel (Reference Goluskin and Spiegel2012) at large R due to a violation of the minimum principle for the temperature, which was not taken into account when formulating them.

Figure 3.

(a–c) Critical temperature $\hat {T}_0(z)$ recovered using (4.16). Colours indicate the Rayleigh number. (a,c) Show details of the boundary layers. (d–f) Detailed view of $\hat {T}_0$ for ${\textit {R}} = 256\,269$, $257\,500$ and $259\,032$. Dashed lines (- - -, black dashed line) are tangent to $\hat {T}_0$ at $z=0$. In (d), $\hat {T}_0$ is non-negative and has minimum of zero inside the layer. In (e), $\hat {T}_0$ is initially positive but has a negative minimum. In (f), $\hat {T}_0'(0)=0$ and there is no positive initial layer. (g) Upper bounds $U$ on $\langle wT \rangle$. Circles mark the values of ${\textit {R}}$ considered in (d–f) and $U=1/2$ at ${\textit {R}} = 259\,032$.

5.2. A new Rayleigh-dependent analytical bound

The numerical results in the previous section demonstrate that optimising $\psi$, $a$ and $b$ cannot improve the uniform bound $\langle wT \rangle \leqslant \frac 12$ at arbitrarily large R. Nevertheless, it is possible to derive better bounds analytically over a finite range of Rayleigh numbers by considering piecewise-linear profiles $\psi (z)$ with two boundary layers, such as the one sketched in figure 4. Even though the numerically optimal profiles in figure 2 show no symmetry with respect to the vertical midpoint $z=\frac 12$, we impose anti-symmetry and take

(5.1) \begin{equation} \psi(z)= \begin{cases} 1- \left(\frac{a+1}{2\delta} -a \right)z ,& 0\leqslant z \leqslant \delta\\ az +\tfrac{1}{2}(1-a), & \delta \leqslant z \leqslant 1-\delta \\ \left(\frac{a+1}{2\delta} -a \right)(1-z), & 1-\delta \leqslant z \leqslant 1, \end{cases} \end{equation}

where $\delta$ is a boundary layer width to be specified later. This considerably simplifies the algebra, at the cost of a quantitatively (but not qualitatively) worse bound on $\langle wT \rangle$. Our goal is to determine values for $\delta$, $a$ and $b$ such that $\psi$ satisfies the constraints in the reduced optimisation problem (4.18), while trying to minimise its objective function.

Figure 4.

General piecewise-linear $\psi (z)$, parametrised by the boundary layer widths $\delta$ and $\varepsilon$, the bulk slope $a$ and the boundary layer height $\sigma$. In our proof, we set $\varepsilon =\delta$ and $\sigma = \tfrac {1}{2}-a (\tfrac {1}{2}-\delta )$ to obtain a profile that is anti-symmetric with respect to $z=1/2$.

To show that the functional $\mathcal {S}_{\boldsymbol {k}}$ is non-negative, observe that the only sign-indefinite term in (4.11) is

(5.2)\begin{equation} \int_0^1 (a-\psi')\hat{w}\hat{T} \, \textrm{d}z = \frac{a+1}{2\delta}\int_{[0,\delta] \cup [1-\delta,1]} \hat{w}\hat{T} \,\textrm{d}z. \end{equation}

As detailed in Appendix E, this term can be estimated using the fundamental theorem of calculus, the Cauchy–Schwarz inequality, and the boundary conditions (4.8a,b) to conclude that $\mathcal {S}_{\boldsymbol {k}}$ is non-negative if

(5.3)\begin{equation} \delta \leqslant \frac{8 \sqrt{2 a b}}{ (a+1) \sqrt{{\textit{R}}}}. \end{equation}

According to the analysis in § 4.3, any choice of $a$, $b$ and $\delta$ satisfying this inequality produces the upper bound

(5.4)\begin{align} \langle wT \rangle &\leqslant \frac{1}{2} + \frac{1}{4b} \left\| bz -\frac{b}{2} - \psi'(z) - 1 \right\|_2^{2} - \int^{1}_{0} \psi(z)\, \textrm{d}z \nonumber\\ &= \frac{b}{48} + \frac{(a+1)^2}{8b\delta} - \frac{(a+1)^2}{4b}. \end{align}

This bound is clearly minimised when $\delta$ is as large as allowed by (5.3), leading to

(5.5)\begin{equation} \langle wT \rangle \leqslant \frac{b}{48} + \frac{ (a+1)^3 \sqrt{{\textit{R}}}}{64 b \sqrt{2 a b}} - \frac{(a+1)^2}{4b}. \end{equation}

The values of $a$ and $b$ minimising the right-hand side of this inequality satisfy

(5.6a)\begin{gather} {-}\sqrt{b} + \frac{\sqrt{R}}{64\sqrt{2a^3}}(a+1)(5a -1 )=0, \end{gather}

(5.6b)\begin{gather}1 + \frac{12(a+1)^2}{b^2} - \frac{9(a+1)^2\sqrt{R}}{8\sqrt{2ab^5}}= 0, \end{gather}

and can be computed numerically for fixed ${\textit {R}}$ to obtain upper bounds plotted as a blue dot-dashed line in figure 5.

Figure 5.

Comparison between the semi-analytical bounds computed with (5.5) and (5.6a,b) ($\text {-}{\cdot }\text {-}$, blue dot dashed line), the explicit bound (5.8) (—, blue solid line), the numerically optimal bounds obtained in § 5.1 with (4.12) (- - -, black dashed line, $L_x=L_y=2$) and with (4.14) ($\bullet$) and experimental and DNS data (see figure 2 for a key of the symbols). A dashed vertical line (- - -, grey dashed line) indicates the smallest energy stability threshold, ${\textit {R}} \approx 26\,926$, while a dashed horizontal line (- - -, grey dashed line) indicates the uniform upper bound $\langle wT \rangle \leqslant \frac 12$.

Fully analytical bounds can instead be obtained if we drop the last term from (5.5). In this case, the optimal $a$ and $b$ are found explicitly as

(5.7a,b)\begin{equation} a = \frac15 \quad\text{and}\quad b = \frac95 \left(\frac{{\textit{R}}}{2}\right)^{1/5},\end{equation}

such that we obtain

(5.8)\begin{equation} \langle wT \rangle \leqslant 2^{-{21}/{5}} R^{1/5} . \end{equation}

This bound, plotted as a solid line in figure 5, is smaller than the uniform bound of 1/2 up to ${\textit {R}} = 2^{16}=65\,536$, which is approximately 2.43 times larger than the energy stability threshold.

6.1. Numerically optimal bounds

Problem (4.25) was discretised into an SDP and solved with the high-precision solver sdpa-gmp (Yamashita et al. Reference Yamashita, Fujisawa, Fukuda, Kobayashi, Nakata and Nakata2012) for $2.0\times 10^5 \leqslant {\textit {R}} \leqslant 3.4 \times 10^5$. The MATLAB toolbox SparseCoLO (Fujisawa et al. Reference Fujisawa, Kim, Kojima, Okamoto and Yamashita2009) was used to exploit sparsity in the SDPs. At each Rayleigh number, we employed the finite-element discretisation approach described in Appendix D on a Chebyshev mesh with at least 6000 piecewise-linear elements, increasing the resolution until the upper bounds changed by less than 1 %. Achieving this at ${\textit {R}} = 3.4\times 10^5$ required approximately $12\,200$ elements. The numerical challenges associated with setting up the SDPs accurately in double precision using SparseCoLO on even finer meshes prevented us from considering a wider the range of Rayleigh numbers.

Figure 6 compares numerical upper bounds on the convective heat flux obtained with (solid line) and without (dot-dashed line) the minimum principle for the temperature, that is, by optimising $q(z)$ or by setting $q(z)=-1$ in (4.25), respectively. The results for the latter case coincide with those described in § 5.1 and shown in figure 2.

Figure 6.

(a) Optimal bounds $U$ on $\langle wT \rangle$ computed by solving (4.25), which incorporates the minimum principle for temperature (—, black solid line). Also plotted ($\text {-}{\cdot }\text {-}$, black dot dashed line) are the bounds computed in § 5 without the minimum principle ($q(z)=-1$). (b) Detail of the region inside the red dashed box (- - -, red dashed line) in (a). (c) Difference between our optimal bounds and the uniform bound $\langle wT \rangle \leqslant 1/2$, shown in all panels as a dotted line (${\cdots }{\cdots }$, grey dotted). In (b,c), a circle at ${\textit {R}}\approx 256\,269$ marks the point at which $q(z)$ begins to vary from $-1$, while stars ($\ast$) mark the Rayleigh numbers at which $\psi$ are plotted in figure 7.

The choice $q(z)=-1$ is optimal for ${\textit {R}} < 256\,269$. For higher R, the numerical upper bounds on $\langle wT \rangle$ with optimised $q(z)$ are strictly better than those with $q(z)=-1$ and, crucially, appear to approach the uniform bound $\frac 12$ from below as ${\textit {R}}$ is raised. Note that, although the deviation from $\frac 12$ is small at the highest values of R that could be handled, it is much larger than the tolerance ($10^{-25}$) used by the multiple-precision SDP solver sdpa-gmp, giving us confidence that it is not a numerical artefact. Moreover, the deviation from $\frac 12$ of the numerical bounds, shown in figure 6(c), appears to decay as a power law, suggesting that the optimal bound available within our bounding framework may have the functional form (1.2). Although the range of Rayleigh numbers spanned by our computations is too small to accurately predict the exponent and the prefactor, it is clear that the decay is much faster than predicted by the heuristic arguments in § 3.

Optimal profiles $\psi (z)$ for selected Rayleigh numbers are shown in figure 7(a,b) and differ significantly from the corresponding profiles in figure 2(b) obtained when the minimum principle for the temperature is disregarded. When the minimum principle is enforced, $\psi$ appears to approach zero almost everywhere as R is raised, but always satisfies $\psi (0)=1$. This leads to the formation of a very thin boundary layer near $z=0$, which at high R consists of three distinct sublayers identified by two points, $z=\delta _0$ and $z=\delta _1$, at which $\psi$ is not differentiable. These are indicated by grey dashed lines in figure 7(b). In the first sublayer, from $z=0$ to $z=\delta _1$, $\psi$ is observed to vary linearly. The second sublayer, $\delta _0 < z < \delta _1$, is observed only for ${\textit {R}} > 2.6\times 10^5$ and we observe that $\psi \sim z^{-1}$ approximately. The third sublayer, from $z=\delta _1$ to the point $z=\delta _2$ at which $\psi$ attains a local minimum, does not have a simple functional form. In the bulk, $\psi$ increases approximately linearly with slope very close to $a$, as was the case in § 5, and the condition $\psi (1)=0$, which emerges as a result of the optimisation and is not imposed a priori, is attained through a small boundary layer of width $\varepsilon$ near $z=1$. We choose the boundary of this layer as the point $z=1-\varepsilon$ at which $\psi$ has a local maximum.

Figure 7.

(a) Optimal $\psi (z)$ in the entire domain, plotted for $0\leqslant \psi \leqslant 0.006$ for visualisation purposes. At all ${\textit {R}}$ values, $\psi (0)=1$. (b) Logarithmic plot of $\psi (z)$, highlighting the behaviour near $z=0$. Coloured lines correspond to the ${\textit {R}}$ highlighted in figure 6, for ${\textit {R}} = 2.6\times 10^5$ (—, yellow solid line), $2.8\times 10^5$ (—, orange solid line), $3.0\times 10^5$ (—, red solid line), $3.2\times 10^5$ (—, brown solid line). Dashed lines (- - -, grey dashed line) mark the boundaries of the sublayers $(0,\delta _0)$, $(\delta _0,\delta _1)$ and $(\delta _1,\delta _2)$, and the sublayer edges are labelled explicitly for ${\textit {R}} = 3.2\times 10^5$. (c) Variation with ${\textit {R}}$ of the optimal balance parameters $a$ (—, blue solid line) and $b$ (- - -, red dashed line) for (4.25). (d) Ratio of balance parameters, $b/a$, when the Lagrange multiplier is active.

Figure 7(c) shows the variation of the balance parameters with ${\textit {R}}$. Below $256\,269$, both coincide with the values plotted in figure 2(c). At higher ${\textit {R}}$, the minimum principle for the temperature becomes active and both balance parameters start to decay rapidly. It would be tempting to conjecture that $a\sim b \sim {\textit {R}}^{-p}$ for some power $p$ but, given the small variation of $b/a$ evident in figure 7(d), we cannot currently exclude that subtly different scaling exponents or higher-order corrections do not play an important role in obtaining an upper bound on $\langle wT \rangle$ that approaches $\frac 12$ asymptotically from below.

Figure 8 shows the structure of $q$, whose (distributional) derivative represents the Lagrange multiplier enforcing the minimum principle. The multiplier, therefore, is active in regions where $q$ is not constant. The choice $q(z)=-1$ is optimal for ${\textit {R}} = 256\,269$. At higher Rayleigh numbers, the multiplier becomes active between $z=\delta _0$ and $z=\delta _1$, which is what causes the second boundary sublayer in $\psi$. In the immediate vicinity of the bottom boundary ($0\leqslant z \leqslant \delta _0$), $q$ is constant and it appears that $q(z) - \psi '(z) \approx b/2$ (cf. figure 10c). Indeed, inspection of the cost function in (4.28) suggests that $q(z) - \psi '(z) = b(1-z)/2$ should be optimal, but we could not identify the very small $z$-dependent correction in our numerical results. In the second sublayer ($\delta _0\leqslant z \leqslant \delta _1$), where the Lagrange multiplier is active, $q(z)\sim -z^{-2}$. Again, this is consistent with the minimisation of the cost function in (4.28), as one expects $q$ to cancel the very large contribution of $\psi '$ near the bottom boundary.

Figure 8.

Variation of $q(z)$ shown on a logarithmic scale to highlight the non-zero region of the Lagrange multiplier $q'(z)$, from ${\textit {R}}= 2.6 - 3.0\times 10^{5}$. The colour bar here applies for all $\psi$, $q$ and $\hat {T}_0$ presented in § 6. The edges of the boundary sublayers $(0,\delta _0)$ and $(\delta _0,\delta _1)$ for the largest ${\textit {R}}$ value are explicitly labelled.

Further validation of our numerical results comes from inspection of the critical temperatures $\hat {T}_0(z)$, which can be recovered using (4.26) and are shown in figure 9(a,b) for selected values of ${\textit {R}}$. As expected, the critical temperatures are non-negative for all $z$ and vanish identically (up to small numerical tolerances) in the region $\delta _0 \leqslant z \leqslant \delta _1$, where $q$ is active. We note that, for a given Rayleigh number, this region is strictly larger than the range of $z$ values for which the critical temperatures in figure 3(c) are negative, indicating that the minimum principle alters the problem in a more subtle way than simply saturating the constraint $T\geqslant 0$.

Figure 9.

Variation with ${\textit {R}}$ of the critical temperature profiles, $\hat {T}_0$, plotted for $0\leqslant z \leqslant 3.5\times 10^{-3}$. (a) Plots of individual $\hat {T}_0$ from ${\textit {R}}= 2.62 \text {--} 2.69\times 10^{5}$. (b) Contour plot of $\hat {T}_0$ vs ${\textit {R}}$. Solid white lines mark the boundary sublayer $\delta _0 \leqslant z \leqslant \delta _1$, within which $q'>0$ and $\hat {T}_0 = 0$. Values of $\hat {T}_0$ below $10^{-12}$ are assumed to be numerical zeros.

To analyse the results further we define the diagnostic function

(6.1)\begin{equation} \chi(z) = \psi'(z) - q(z) , \end{equation}

and rewrite (4.27) as

(6.2)\begin{equation} \langle wT \rangle \leqslant \frac{1}{2} + b\left[ \int^{1}_{0} \frac{1}{4}\left(z-\frac{1}{2} - \frac{\chi(z)}{b} \right)^{2} - \frac{\psi(z)}{b} \,\textrm{d}z\right] . \end{equation}

Figure 10(a,c,d) suggest that $\chi (0)/b \to -\tfrac {1}{2}$ and $\chi (1)/b \to -\tfrac {3}{2}$ as ${\textit {R}}\rightarrow \infty$. In fact, profiles $\chi (z)/b$ for different Rayleigh numbers collapse almost exactly throughout the layer, but subtle corrections are present; for instance, figure 10(b) demonstrates that, in the bulk, the mean value of $\chi (z)/b$ decreases with ${\textit {R}}$. It is also evident that $\chi (z)$ is not exactly constant throughout the bulk, but increases by approximately $10^{-4}$. Since $q(z)$ is constant in this region, we conclude that $\psi '(z)$ is not constant, but displays subtle and thus non-trivial variation.

Figure 10.

(a) Profiles $\chi (z)/b$ for $2.6\times 10^5 \leqslant {\textit {R}} \leqslant 3.0 \times 10^{5}$. (b) Detailed view of $\chi (z)/b$ in the bulk. (c,d) Variation of $\chi (0)/b$ and $\chi (1)/b$ with ${\textit {R}}$. (e) Contributions to the integral of $\psi /b$ in the regions $(0, \delta _0)$ (—, blue solid line), $(\delta _0, \delta _1)$ (—, red solid line), $(\delta _1, \delta _2)$ (—, green solid line), $(\delta _2, 1-\varepsilon )$ (—, yellow solid line) and $(1-\varepsilon , 1)$ (—, purple solid line), as a function of R.

Figure 10(e) illustrates the variation with ${\textit {R}}$ of contributions to the integral of $\psi (z) / b$ from regions $(0,\delta _0)$, $(\delta _0,\delta _1)$, $(\delta _1,\delta _2)$, $(\delta _2,1-\varepsilon )$ and $(1-\varepsilon ,1)$. The largest contribution comes from the bulk ($\delta _2 \leqslant z \leqslant 1-\varepsilon$), but it slowly decreases with ${\textit {R}}$. The same is true of the contribution of the top boundary layer ($1-\varepsilon \leqslant z \leqslant 1$) and the outermost boundary sublayer near $z=0$ ($\delta _1 \leqslant z \leqslant \delta _2$). Only in the second sublayer near $z=0$ does the value of the integral increase with ${\textit {R}}$, suggesting that the integral of $\psi (z)/b$ near this boundary layer may become the dominant term as ${\textit {R}} \to \infty$. While the range of Rayleigh numbers covered by our computations is too small to confirm or disprove this conjecture, it is certain that the integral of $\psi (z)/b$ must remain large enough to offset the positive term in (6.2) in order to obtain a bound on $\langle wT \rangle$ smaller than $1/2$.

Finally, figures 8 and 10 suggest that, for sufficiently large ${\textit {R}}$,

(6.3) \begin{equation} q(z) = \begin{cases} \psi'(z) + \frac{b}{2}, & 0 \leqslant z \leqslant \delta_0,\\ -q_{0}z^{{-}2}, & \delta_0 \leqslant z \leqslant \delta_1,\\ \psi'(1) + \frac{3b}{2}, & \delta_1 \leqslant z \leqslant 1, \end{cases} \end{equation}

for some positive constant $q_0$, and that $a,\psi ,q_0 \sim b$. The next section investigates whether such an ansatz can lead to upper bounds on $\langle wT \rangle$ that are strictly smaller than $1/2$ at all Rayleigh numbers.

6.2. Towards an analytical bound

The numerical evidence presented in § 6.1 suggests that the upper bound on $\langle wT \rangle$ obtained from the optimisation problem (4.28) approaches $\tfrac {1}{2}$ asymptotically from below as $R \rightarrow \infty$. To confirm this observation with a proof requires, for every $R \geqslant 1$, construction of feasible decision variables $a,b,\psi (z),q(z)$ whose corresponding cost is strictly less than $\tfrac {1}{2}$. This section discusses the challenges presented by this goal. Specifically, we show that no construction is possible if one tries to mimic key properties of the numerically optimal decision variables presented in § 6.1 and, at the same time, enforces the spectral constraint using estimates typically used in successful applications of the background method.

To aid the discussion, the following definition introduces three subsets $\mathcal {A}, \mathcal {B}$ and $\mathcal {C}$ of decision variables that capture some of the properties observed from our numerical study.

The set $\mathcal {A}\{\delta ,\varepsilon \}$ contains profiles $\psi$ which possess an initial (and potentially severe) boundary layer in an interval $[0,\delta ]$, then increase in a bulk region $[\delta ,1-\varepsilon ]$, before approaching $\psi (1)=0$ in an upper boundary layer contained in the interval $[1-\varepsilon ,1]$ and in which $\psi '$ is controlled by the balance parameter $b$, as seen in previous results. Optimal decision variables $(a,b, \psi (z),q(z))$ obtained in § 6.1 appear, with compelling evidence, to belong to a set of the form $\mathcal {A}\{\delta ,\varepsilon \}$ with the exception of the differentiability condition $\psi \in C^1[0,1]$. Indeed, the optimal profiles $\psi$ appear to be piecewise differentiable, losing differentiability at two points corresponding to the boundaries of non-constant behaviour of the multiplier $q(z)$ observed in figure 8. However, since no higher derivatives of $\psi$ appear in the optimisation problem (4.28), adding the constraint $(a,b,\psi (z),q(z)) \in \mathcal {A}$ to (4.28) will not change its optimal cost.

The set $\mathcal {B}\{\gamma \}$ contains decision variables for which $\psi '$ and $q$ are constant in some interval centred at $\tfrac {1}{2}$. Figures 8 and 10 reveal that this is not the case for the numerically optimal $\psi$, so the use of $\mathcal {B}\{\gamma \}$ corresponds to a proof which ignores subtle variations from a purely linear profile away from the boundaries. Without further assumptions (say, restriction to fluids with infinite Prandtl number), it is not clear how such variations can be exploited in analytical constructions.

The set $\mathcal {C}\{\delta ,\varepsilon ,R\}$ relates to a choice of profiles $\psi$ and balance parameters $a,b$ for which the constraint $\mathcal {S}_{\boldsymbol {k}} \geqslant 0$ can be proven to hold for a given $R$. In particular, if $a,b,\psi (z)$ satisfy (E11) and it is the case that $\psi '|_{{[\delta ,1-\varepsilon ]}}=a$, then the argument in Appendix E implies that $\tilde {\mathcal {S}} \geqslant 0$. Specifically, $\mathcal {C}\{\delta ,\varepsilon ,R\}$ provides sufficient control of the severity of the boundary layers of $\psi$ for the spectral constraint to be provably satisfied. While crude, estimates of this form in conjunction with constant $\psi '(z)$ in a bulk region $\delta \leqslant z \leqslant 1-\varepsilon$ are employed for almost all analytical constructions of background fields.

We now return to the original question of attempting to upper bound $\langle wT \rangle$ via an analytical construction of feasible decision variables $(a,b,\psi (z),q(z))$ for (4.28). It is not unreasonable, based upon the above evidence, to propose $R$-dependent balance parameters $a=a_R, b=b_R$, boundary layer widths $\delta = \delta _R, \varepsilon = \varepsilon _R$ and profiles $\psi _R,q_R$ which satisfy

(6.5)\begin{equation} (a_R,b_R,\psi_R(z),q_R(z)) \in \mathcal{A}\{ \delta_R,\varepsilon_R\} \cap \mathcal{B}\{\gamma\} \cap \mathcal{C}\{ \delta_R,\varepsilon_R,R\}, \quad R \geqslant 1, \end{equation}

for some $\gamma >0$. The following result shows that, for such a construction, there is a hard lower bound on the optimal cost achievable using (4.28).

Proposition 6.2 Let $0 < \delta < 1-\varepsilon <1$ and $\gamma >0, R \geqslant 1$. Suppose that

(6.6)\begin{equation} (a,b,\psi(z),q(z)) \in \mathcal{A}\{\delta,\varepsilon\} \cap \mathcal{B}\{\gamma\} \cap \mathcal{C}\{\delta,\varepsilon,R\}. \end{equation}

Then

(6.7)\begin{equation} \frac{1}{4b}\left\| bz-\frac{b}{2} - \psi'(z) + q(z) \right\|_2^2 - \int_{0}^{1}\psi \, \text{d}z \geqslant \frac{b}{6} \left( \gamma^3 - 24 \cdot \frac{1 + 2\sqrt{2}}{R^{1/4}} \right). \end{equation}

From the proof in Appendix F, the consequence of Proposition 6.2 is the following. If one constructs feasible decision variables for the optimisation problem (4.28) which satisfy (6.6), then the best achievable bound $U$ for which $\langle wT \rangle \leqslant U$ must satisfy

(6.8)\begin{equation} U \geqslant \frac12 + \frac{b}{6} \left( \gamma^3 - 24 \cdot \frac{1+2\sqrt{2}}{R^{\frac14}} \right). \end{equation}

Hence, using such a construction with the assumption that $\psi '$ and $q$ are constant in a bulk region $(\frac 12-\gamma ,\frac 12+\gamma )$, it is not possible to prove that $\langle wT \rangle < \frac 12$ at arbitrarily high Rayleigh number.

Consequently, one must ask what conditions should be dropped if a rigorous bound $\langle wT \rangle < \frac 12, R \geqslant 1$, is to be found. The numerical evidence presented in § 6.1 suggests that the optimal decision variables do belong to $\mathcal {A}$. Consequently, either $\mathcal {B}$ or $\mathcal {C}$ must be dropped. Figure 8 indicates that optimal multipliers $q(z)$ are constant outside a lower boundary layer. Hence, dropping $\mathcal {B}$ corresponds to choosing a profile with non-constant $\psi '(z)$ in the bulk; the cost function of (4.28) and figure 10(b) suggests that a quadratic ansatz for $\psi (z)$ may be beneficial. Dropping $\mathcal {C}$ corresponds to requiring more sophisticated analysis of the spectral constraint. Using these insights will be the focus of future research.