2.1 Variations with respect to $\unicode[STIX]{x1D702}$

We start by taking the variation of ${\mathcal{S}}$ with respect to $\unicode[STIX]{x1D702}$ and setting it equal to zero. This results in the Euler–Lagrange equations

(2.5)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D702}=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}\end{eqnarray}$$

for $\unicode[STIX]{x1D702}$. Substituting this back into ${\mathcal{S}}$ results in the constrained functional

(2.6)$$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D702}}\equiv \left\langle 2u_{3}\unicode[STIX]{x1D709}-|\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E5}^{-1}\left(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}\right)|^{2}-|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}+\unicode[STIX]{x1D707}\left(Pe^{2}-|\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\right)+\unicode[STIX]{x1D735}p\boldsymbol{\cdot }\boldsymbol{u}\right\rangle .\end{eqnarray}$$

Stated differently, constraining the variable $\unicode[STIX]{x1D702}$ to satisfy (2.5) preserves the saddles of (2.4) and yields

(2.7)$$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D702}}=\min _{\unicode[STIX]{x1D702}}{\mathcal{S}}.\end{eqnarray}$$

2.2 Variations with respect to $\unicode[STIX]{x1D709}$ and $\boldsymbol{u}$

Next, we take the variation of ${\mathcal{S}}$ with respect to the variable $\unicode[STIX]{x1D709}$ and set it equal to zero. This yields

(2.8)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D709}=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}-u_{3}\end{eqnarray}$$

and, after substituting it back into ${\mathcal{S}}$, we find

(2.9)$$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D709}}=\left\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E5}^{-1}\left(\boldsymbol{u}\boldsymbol{\cdot }\left(\hat{z}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\right)\right)|^{2}+\unicode[STIX]{x1D707}\left(Pe^{2}-|\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\right)+\unicode[STIX]{x1D735}p\boldsymbol{\cdot }\boldsymbol{u}\right\rangle .\end{eqnarray}$$

As before, the same functional is obtained by maximizing ${\mathcal{S}}$ in $\unicode[STIX]{x1D709}$, i.e.

$$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D709}}=\max _{\unicode[STIX]{x1D709}}{\mathcal{S}}.\end{eqnarray}$$

Finally, we take variations over all incompressible flow fields $\boldsymbol{u}$. This produces the constrained functional

(2.10)$$\begin{eqnarray}{\mathcal{S}}_{\boldsymbol{u}}=\left\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}-|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}+\unicode[STIX]{x1D707}Pe^{2}+\frac{1}{\unicode[STIX]{x1D707}}|\unicode[STIX]{x1D735}S^{-1}\left(\unicode[STIX]{x1D709}\hat{z}-\unicode[STIX]{x1D709}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\right)|^{2}\right\rangle ,\end{eqnarray}$$

where $S^{-1}\left(\unicode[STIX]{x1D709}\hat{z}-\unicode[STIX]{x1D709}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\right)$ denotes the unique flow field $\boldsymbol{u}$ solving

(2.11)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}\unicode[STIX]{x0394}\boldsymbol{u}=\unicode[STIX]{x1D709}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}-\unicode[STIX]{x1D709}\hat{z}-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D735}p, & \displaystyle\end{eqnarray}$$

(2.12)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0. & \displaystyle\end{eqnarray}$$

Note that

$$\begin{eqnarray}\displaystyle {\mathcal{S}}_{\boldsymbol{u}}=\max _{\boldsymbol{u}}{\mathcal{S}}. & & \displaystyle \nonumber\end{eqnarray}$$

These observations serve as a starting point in § 5 for establishing upper bounds on transport and motivate our choice of numerical methods in § 3. But first, let us see how the wall-to-wall problem comes out of these manipulations.

2.3 Finding the wall-to-wall problem

Consider the structure of the $S_{\unicode[STIX]{x1D702}}=\min _{\unicode[STIX]{x1D702}}{\mathcal{S}}$ and $S_{\unicode[STIX]{x1D709}}=\max _{\unicode[STIX]{x1D709}}{\mathcal{S}}$ functionals for a fixed incompressible velocity field $\boldsymbol{u}$ with enstrophy $\langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle =Pe^{2}$. The functional $S_{\unicode[STIX]{x1D702}}$ is concave in $\unicode[STIX]{x1D709}$; likewise $S_{\unicode[STIX]{x1D709}}$ is convex in $\unicode[STIX]{x1D702}$. Thus, finding the maximum of the former with respect to $\unicode[STIX]{x1D709}$, or the minimum of the latter with respect to $\unicode[STIX]{x1D702}$, is equivalent to enforcing their Euler–Lagrange equations.

The Euler–Lagrange equation for $S_{\unicode[STIX]{x1D702}}$ in $\unicode[STIX]{x1D709}$ is

(2.13)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D709}=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E5}^{-1}\left(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}\right)-u_{3}.\end{eqnarray}$$

Similarly, for $S_{\unicode[STIX]{x1D709}}$ in $\unicode[STIX]{x1D702}$ we have

(2.14)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D702}=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E5}^{-1}\left(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}-u_{3}\right).\end{eqnarray}$$

These can be written more concisely as the system

(2.15)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D702}=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$

(2.16)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D709}=\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}-u_{3}, & \displaystyle\end{eqnarray}$$

obtained by setting

(2.17)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=\unicode[STIX]{x1D6E5}^{-1}\left(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}\right), & \displaystyle\end{eqnarray}$$

(2.18)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}=\unicode[STIX]{x1D6E5}^{-1}\left(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}-u_{3}\right) & \displaystyle\end{eqnarray}$$

in (2.13) and (2.14). Hence,

(2.19)$$\begin{eqnarray}Nu\{\boldsymbol{u}\}-1=\max _{\unicode[STIX]{x1D709}}\min _{\unicode[STIX]{x1D702}}{\mathcal{S}}=\min _{\unicode[STIX]{x1D702}}\max _{\unicode[STIX]{x1D709}}{\mathcal{S}}\end{eqnarray}$$

for each fixed velocity field $\boldsymbol{u}$.

The formula (2.19), which first appeared in Tobasco & Doering (Reference Tobasco and Doering2017) and was further analysed in Doering & Tobasco (Reference Doering and Tobasco2019), is an exact variational characterization of the Nusselt number. It allows the optimal wall-to-wall transport problem to be stated succinctly as

(2.20)$$\begin{eqnarray}\max _{\boldsymbol{u}}Nu\{\boldsymbol{u}\}-1=\max _{\boldsymbol{u}}\max _{\unicode[STIX]{x1D709}}\min _{\unicode[STIX]{x1D702}}{\mathcal{S}}=\max _{\boldsymbol{u}}\min _{\unicode[STIX]{x1D702}}\max _{\unicode[STIX]{x1D709}}{\mathcal{S}},\end{eqnarray}$$

where $\boldsymbol{u}$ satisfies the given boundary conditions and intensity constraints. In particular, we note that gradient ascent may be applied with impunity to the constrained functional ${\mathcal{S}}_{\unicode[STIX]{x1D702}}$ to compute local maxima. This functional can also be used to prove lower bounds on the Nusselt number without having to solve the advection–diffusion equation. Indeed, plugging in any $\unicode[STIX]{x1D709}$ admissible in the previous manipulations into ${\mathcal{S}}_{\unicode[STIX]{x1D702}}$ yields the lower bound

(2.21)$$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D702}}\{\boldsymbol{u},\unicode[STIX]{x1D709}\}\leqslant Nu\{\boldsymbol{u}\}-1.\end{eqnarray}$$

In fact, reinterpreting for the moment angle brackets as space and time averages, we note that the bound (2.21) holds for time-independent flows as well, an observation that was exploited in Tobasco & Doering (Reference Tobasco and Doering2017) and Doering & Tobasco (Reference Doering and Tobasco2019) with ‘branching’ trial functions to prove the scaling $\max \,Nu\sim Pe^{2/3}$ up to logarithmic corrections as $Pe\rightarrow \infty$. We return to discuss this asymptotic result in the context of our numerical results much further below. Next, we consider the relationship between the wall-to-wall problem and the background method.

2.4 Finding the background method

Consider now the background method which guarantees the absolute upper bound

(2.22)$$\begin{eqnarray}Nu-1\leqslant \min _{\unicode[STIX]{x1D702}(\boldsymbol{x}),\unicode[STIX]{x1D707}}\left\{\begin{array}{@{}ll@{}}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}\rangle +\unicode[STIX]{x1D707}Pe^{2}\quad & \text{if }\langle {\mathcal{Q}}[\boldsymbol{u},\unicode[STIX]{x1D709};\unicode[STIX]{x1D702}]\rangle \leqslant 0~~\forall \,\boldsymbol{u},\unicode[STIX]{x1D709},\\ \infty \quad & \text{otherwise},\end{array}\right.\end{eqnarray}$$

where

(2.23)$$\begin{eqnarray}{\mathcal{Q}}[\boldsymbol{u},\unicode[STIX]{x1D709};\unicode[STIX]{x1D702}]=2\unicode[STIX]{x1D709}\boldsymbol{u}\boldsymbol{\cdot }\left(\hat{z}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\right)-|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}-\unicode[STIX]{x1D707}|\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}.\end{eqnarray}$$

The reader may not immediately recognize this as the familiar background method as it has been applied to Rayleigh–Bénard convection (Doering & Constantin Reference Doering and Constantin1996). Nevertheless, equation (2.22) does follow from applying the usual argument to the wall-to-wall problem. (The resulting bounds carry over to the time-dependent case.)

Let us recall the argument now. Starting with the advection–diffusion equation

(2.24)$$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=\unicode[STIX]{x0394}T\end{eqnarray}$$

and decomposing $T$ as $T=\unicode[STIX]{x1D709}+1-z+\unicode[STIX]{x1D702}$ yields

(2.25)$$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}=\unicode[STIX]{x0394}\unicode[STIX]{x1D709}+\unicode[STIX]{x0394}\unicode[STIX]{x1D702}+u_{3}.\end{eqnarray}$$

Multiplying through by $\unicode[STIX]{x1D709}$ and integrating by parts yields the balance relation

(2.26)$$\begin{eqnarray}\langle \unicode[STIX]{x1D709}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}-u_{3}\unicode[STIX]{x1D709}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}\rangle =0.\end{eqnarray}$$

Now, utilizing

(2.27)$$\begin{eqnarray}\langle wT\rangle =\langle |\unicode[STIX]{x1D735}T|^{2}\rangle -1=\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}+2\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\rangle ,\end{eqnarray}$$

we subtract twice (2.26) from (2.27) to conclude that

(2.28)$$\begin{eqnarray}Nu-1=\left\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}+2\unicode[STIX]{x1D709}\boldsymbol{u}\boldsymbol{\cdot }\left(\hat{z}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\right)-|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}\right\rangle .\end{eqnarray}$$

Introducing a Lagrange multiplier $\unicode[STIX]{x1D707}/2$ for the enstrophy constraint $\langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle =Pe^{2}$ yields

(2.29)$$\begin{eqnarray}\displaystyle Nu-1 & = & \displaystyle \left\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}+2\unicode[STIX]{x1D709}\boldsymbol{u}\boldsymbol{\cdot }\left(\hat{z}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}\right)-|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}\right\rangle +\unicode[STIX]{x1D707}(Pe^{2}-\langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle ),\end{eqnarray}$$

(2.30)$$\begin{eqnarray}\displaystyle & = & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}\rangle +\unicode[STIX]{x1D707}Pe^{2}+\langle {\mathcal{Q}}[\boldsymbol{u},\unicode[STIX]{x1D709};\unicode[STIX]{x1D702}]\rangle\end{eqnarray}$$

and upon performing the operations $\min _{\unicode[STIX]{x1D702}}\max _{\boldsymbol{u},\unicode[STIX]{x1D709}}$ we deduce (2.22).

2.5 A possible duality gap

We can now discuss the relationship between the background method and the wall-to-wall optimal transport problem. Combining the definition of ${\mathcal{S}}$ from (2.4) and the identity (2.29) we see that the background method bound (2.22) can be alternatively written as

(2.31)$$\begin{eqnarray}\max _{\boldsymbol{u}}Nu-1\leqslant \min _{\unicode[STIX]{x1D702}}\max _{\boldsymbol{u},\unicode[STIX]{x1D709}}{\mathcal{S}}.\end{eqnarray}$$

On the left-hand side appears the wall-to-wall optimal transport problem, while on the right-hand side appears the background method. Note this inequality is consistent with the results of § 2.3 since in any case

(2.32)$$\begin{eqnarray}\max _{\boldsymbol{u},\unicode[STIX]{x1D709}}\min _{\unicode[STIX]{x1D702}}{\mathcal{S}}\leqslant \min _{\unicode[STIX]{x1D702}}\max _{\boldsymbol{u},\unicode[STIX]{x1D709}}{\mathcal{S}}\end{eqnarray}$$

regardless of the definition of ${\mathcal{S}}$. Now, if ${\mathcal{S}}$ were convex in $\unicode[STIX]{x1D702}$ and jointly concave in $(\boldsymbol{u},\unicode[STIX]{x1D709})$ one would be led on general grounds via convex duality to conjecture that equality should hold between the left-hand and right-hand sides above, in which case the wall-to-wall problem and the background method would turn out to be equivalent. We instead propose that the opposite situation is true and that

(2.33)$$\begin{eqnarray}\max _{\boldsymbol{u}}Nu-1\neq \min _{\unicode[STIX]{x1D702}}\max _{\boldsymbol{u},\unicode[STIX]{x1D709}}{\mathcal{S}}.\end{eqnarray}$$

A concurrent study, Ding & Kerswell (Reference Ding and Kerswell2019), has provided further evidence that such an inequality holds. This inequality still allows the possibility for these quantities to achieve the same asymptotic scaling as $Pe\rightarrow \infty$, a situation suggested for three-dimensional wall-to-wall optimal transport by the recent numerical scaling $\max \,Nu\sim Pe^{2/3}$ reported for a finite range of $Pe$ in Motoki et al. (Reference Motoki, Kawahara and Shimizu2018a). It would also be consistent with the dimension-independent logarithmic lower bound $\max \,Nu\geqslant C^{\prime }Pe^{2/3}/(\log Pe)^{4/3}$ proved for all large enough $Pe$ in Tobasco & Doering (Reference Tobasco and Doering2017) and Doering & Tobasco (Reference Doering and Tobasco2019).

Let us illustrate why the inequality in (2.33) may be true by considering how the previous manipulations operate on the polynomial

(2.34)$$\begin{eqnarray}\displaystyle & \displaystyle p(\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D70F}_{3},u,v)=(\unicode[STIX]{x1D70F}_{1})^{2}+(\unicode[STIX]{x1D70F}_{2})^{2}+(\unicode[STIX]{x1D70F}_{3})^{2}+q(u,v,\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D70F}_{3}), & \displaystyle\end{eqnarray}$$

(2.35)$$\begin{eqnarray}\displaystyle & \displaystyle q=\left[\begin{array}{@{}cc@{}}u & v\end{array}\right]\left[\begin{array}{@{}cc@{}}2(1-\unicode[STIX]{x1D70F}_{1}+\unicode[STIX]{x1D70F}_{2}) & \unicode[STIX]{x1D70F}_{3}\\ \unicode[STIX]{x1D70F}_{3} & 2(1-\unicode[STIX]{x1D70F}_{1}-\unicode[STIX]{x1D70F}_{2})\end{array}\right]\left[\begin{array}{@{}c@{}}u\\ v\end{array}\right], & \displaystyle\end{eqnarray}$$

which we see as analogous to (2.4). The variable $\unicode[STIX]{x1D70F}_{1}$ is analogous to the zeroth Fourier mode of $\unicode[STIX]{x1D702}$, while $\unicode[STIX]{x1D70F}_{2}$ and $\unicode[STIX]{x1D70F}_{3}$ are to the non-zero Fourier modes. (This polynomial was not derived as a modal truncation of ${\mathcal{S}}$. That would produce a more complicated example.) The variables $u$ and $v$ are analogous to $\unicode[STIX]{x1D709}$ and $\boldsymbol{u}$. The fact is that

(2.36)$$\begin{eqnarray}\max _{u,v}\min _{\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D70F}_{3}}p=4/5<\min _{\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D70F}_{3}}\max _{u,v}p=1,\end{eqnarray}$$

and that the optimizers for the ‘background method’ $\min$$\max$ problem are not saddle points of $p$. In particular, the critical point equations

(2.37a-c)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{1}=u^{2}+v^{2},\quad \unicode[STIX]{x1D70F}_{2}=v^{2}-u^{2},\quad \text{and}\quad \unicode[STIX]{x1D70F}_{3}=-uv,\end{eqnarray}$$

(2.38)$$\begin{eqnarray}\left[\begin{array}{@{}c@{}}0\\ 0\end{array}\right]=\left[\begin{array}{@{}cc@{}}2(1-\unicode[STIX]{x1D70F}_{1}+\unicode[STIX]{x1D70F}_{2}) & \unicode[STIX]{x1D70F}_{3}\\ \unicode[STIX]{x1D70F}_{3} & 2(1-\unicode[STIX]{x1D70F}_{1}-\unicode[STIX]{x1D70F}_{2})\end{array}\right]\left[\begin{array}{@{}c@{}}u\\ v\end{array}\right],\end{eqnarray}$$

fail to be satisfied by solutions of the $\min$$\max$ problem.

To see why this is the case, consider the background method procedure wherein the maximum occurs first. If any of the eigenvalues of (2.35) are positive then the maximum over $u$ and $v$ yields infinity; thus, we must calculate the eigenvalues of (2.35) to see when this occurs. For fixed $(\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D70F}_{3})$, the eigenvalues of the matrix in (2.35) are $\unicode[STIX]{x1D706}=2(1-\unicode[STIX]{x1D70F}_{1})\pm \sqrt{4(\unicode[STIX]{x1D70F}_{2})^{2}+(\unicode[STIX]{x1D70F}_{3})^{2}}$. The only way these eigenvalues are non-positive is if $2(1-\unicode[STIX]{x1D70F}_{1})+\sqrt{4(\unicode[STIX]{x1D70F}_{2})^{2}+(\unicode[STIX]{x1D70F}_{3})^{2}}\leqslant 0$, or equivalently $2\unicode[STIX]{x1D70F}_{1}\geqslant 2+\sqrt{4(\unicode[STIX]{x1D70F}_{2})^{2}+(\unicode[STIX]{x1D70F}_{3})^{2}}$. Hence,

(2.39)$$\begin{eqnarray}\max _{u,v}p(\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D70F}_{3},u,v)=\left\{\begin{array}{@{}ll@{}}(\unicode[STIX]{x1D70F}_{1})^{2}+(\unicode[STIX]{x1D70F}_{2})^{2}+(\unicode[STIX]{x1D70F}_{3})^{2}\quad & \text{ if }2\unicode[STIX]{x1D70F}_{1}\geqslant 2+\sqrt{4(\unicode[STIX]{x1D70F}_{2})^{2}+(\unicode[STIX]{x1D70F}_{3})^{2}},\\ \infty \quad & \text{otherwise}.\end{array}\right.\end{eqnarray}$$

It follows immediately that $\min \max p=1$, and that the minimizer satisfies $\unicode[STIX]{x1D70F}_{1}=1$ and $\unicode[STIX]{x1D70F}_{2}=\unicode[STIX]{x1D70F}_{3}=0$; however, such $\unicode[STIX]{x1D70F}_{1}$, $\unicode[STIX]{x1D70F}_{2}$ and $\unicode[STIX]{x1D70F}_{3}$ cannot be a saddle point of $p$: if $\unicode[STIX]{x1D70F}_{2}=\unicode[STIX]{x1D70F}_{3}=0$ then from (2.37) we see that $u^{2}=v^{2}$ and $uv=0$ so that $u=v=0$, but then the $\unicode[STIX]{x1D70F}_{1}$ equation cannot be satisfied since $\unicode[STIX]{x1D70F}_{1}=1\neq u^{2}+v^{2}$.

Proceeding in the reverse order we find that

(2.40)$$\begin{eqnarray}\min _{\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D70F}_{3}}p=-(u^{2}+v^{2})^{2}-(u^{2}-v^{2})^{2}-(uv)^{2}+2(u^{2}+v^{2}).\end{eqnarray}$$

The minimizing $\unicode[STIX]{x1D70F}$ satisfy

(2.41a-c)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{1}=u^{2}+v^{2},\quad \unicode[STIX]{x1D70F}_{2}=v^{2}-u^{2},\quad \text{and}\quad \unicode[STIX]{x1D70F}_{3}=-uv.\end{eqnarray}$$

The maximum over $u,v$ is given by $u=\pm \sqrt{2/5}$ and $v=\pm \sqrt{2/5}$. Hence, $\max \min p=4/5$.

While in this example the background method procedure (maximum followed by minimum) yields results that are incompatible with the saddle points of (2.34), the wall-to-wall procedure (minimum followed by maximum) does produce saddle points. Returning to the actual wall-to-wall optimal transport problem, we note that the optimal flow fields reported in § 4 exhibit non-trivial non-zero Fourier modes for the variable $\unicode[STIX]{x1D702}$, whereas in the background method optimizers must satisfy $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}(z)$. Indeed, if $\unicode[STIX]{x1D702}$ satisfies the spectral stability constraint $\langle {\mathcal{Q}}\rangle \leqslant 0$ then so does its periodic average $\overline{\unicode[STIX]{x1D702}}$ in $x$, while by Jensen’s inequality $\langle |(\text{d}/\text{d}z)\overline{\unicode[STIX]{x1D702}}|^{2}\rangle \leqslant \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}\rangle$ with equality if and only if $\unicode[STIX]{x1D702}=\overline{\unicode[STIX]{x1D702}}$. These observations strongly suggest that the conjectured gap (2.33) between the wall-to-wall problem and the background method should hold and, in particular, that the spectral stability constraint should fail to be satisfied by the true saddle points of ${\mathcal{S}}$.

2.6 Comparison with the Howard–Busse–Malkus problem

We would be remiss if we did not additionally state the connection of the previous discussions on the wall-to-wall and background method approach with the classic Howard–Busse–Malkus approach put forth in Howard (Reference Howard1963). To see the connection between these, start with ${\mathcal{S}}$ and restrict attention to incompressible flows $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ with enstrophy $\langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle =Pe^{2}$ and functions $\unicode[STIX]{x1D702}(z)$. At this point, it is useful to introduce notation for the horizontal average of a function,

(2.42)$$\begin{eqnarray}\overline{f}\equiv \frac{1}{\unicode[STIX]{x1D6E4}}\int _{0}^{\unicode[STIX]{x1D6E4}}f(x,z)\,\text{d}x.\end{eqnarray}$$

Computing the minimum of ${\mathcal{S}}$ with respect to $\unicode[STIX]{x1D702}(z)$ yields the optimality condition

(2.43)$$\begin{eqnarray}\frac{\text{d}^{2}}{\text{d}z^{2}}\unicode[STIX]{x1D702}=\overline{\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}}=\frac{\text{d}}{\text{d}z}\overline{u_{3}\unicode[STIX]{x1D709}}\end{eqnarray}$$

whose solution is

(2.44)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}(z)=\int _{0}^{z}\overline{u_{3}\unicode[STIX]{x1D709}}(z^{\prime })\,\text{d}z^{\prime }-z\int _{0}^{1}\overline{u_{3}\unicode[STIX]{x1D709}}(z^{\prime })\,\text{d}z^{\prime }=\int _{0}^{z}\overline{u_{3}\unicode[STIX]{x1D709}}(z^{\prime })\,\text{d}z^{\prime }-z\langle u_{3}\unicode[STIX]{x1D709}\rangle , & \displaystyle\end{eqnarray}$$

(2.45)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}z}\unicode[STIX]{x1D702}=\overline{u_{3}\unicode[STIX]{x1D709}}(z)-\langle u_{3}\unicode[STIX]{x1D709}\rangle . & \displaystyle\end{eqnarray}$$

Employing this relation in ${\mathcal{S}}$ yields

(2.46)$$\begin{eqnarray}\min _{\unicode[STIX]{x1D702}(z)}{\mathcal{S}}=\left\langle 2u_{3}\unicode[STIX]{x1D709}-\left(\overline{u_{3}\unicode[STIX]{x1D709}}-\langle u_{3}\unicode[STIX]{x1D709}\rangle \right)^{2}-|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}\right\rangle .\end{eqnarray}$$

Making the change of variables $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70E}$ for a soon to be determined scalar $\unicode[STIX]{x1D6FC}$ transforms this to

(2.47)$$\begin{eqnarray}\min _{\unicode[STIX]{x1D702}(z)}{\mathcal{S}}=\unicode[STIX]{x1D6FC}\left\langle 2u_{3}\unicode[STIX]{x1D70E}\right\rangle -\unicode[STIX]{x1D6FC}^{2}\left\langle \left(\overline{u_{3}\unicode[STIX]{x1D70E}}+\langle u_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}|^{2}\right\rangle ,\end{eqnarray}$$

a quadratic function of $\unicode[STIX]{x1D6FC}$. Maximizing with respect to $\unicode[STIX]{x1D6FC}$ determines the optimal choice

(2.48)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}^{\ast }=\frac{\langle u_{3}\unicode[STIX]{x1D70E}\rangle }{\left\langle \left(\overline{u_{3}\unicode[STIX]{x1D70E}}-\langle u_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}|^{2}\right\rangle },\end{eqnarray}$$

utilizing $\unicode[STIX]{x1D6FC}^{\ast }$ in the above and taking $\boldsymbol{u}=Pe(\boldsymbol{v}/\sqrt{\langle |\unicode[STIX]{x1D735}\boldsymbol{v}|^{2}\rangle })$ results in the multiplicative form of the functional

(2.49)$$\begin{eqnarray}{\mathcal{M}}[\unicode[STIX]{x1D70E},\boldsymbol{v},Pe]=\frac{\left(\langle v_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}}{Pe^{-2}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}|^{2}\rangle \langle |\unicode[STIX]{x1D735}\boldsymbol{v}|^{2}\rangle +\langle \left(\overline{v_{3}\unicode[STIX]{x1D70E}}-\langle v_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}\rangle },\end{eqnarray}$$

exactly as in Howard (Reference Howard1963) under the appropriate retranscriptions. The functional in (2.49) is homogeneous with respect to $\unicode[STIX]{x1D70E}$ and $\boldsymbol{v}$, i.e. ${\mathcal{M}}[\unicode[STIX]{x1D706}\unicode[STIX]{x1D70E},\unicode[STIX]{x1D706}\boldsymbol{v},Pe]={\mathcal{M}}[\unicode[STIX]{x1D70E},\boldsymbol{v},Pe]$. It has in the past served as a starting point for the analysis of maximal heat transport, in particular leading in Howard (Reference Howard1963) to bounds on transport under the author’s assumptions of homogeneity and statistical similarity. The connection of this functional to the background method has been pointed out before (Kerswell Reference Kerswell1998).

However, we point out now that the resulting bounds on transport can be improved beyond those obtained using (2.49) due to the variational representation (2.19) of the wall-to-wall problem. After a similar series of manipulations utilizing all possible $\unicode[STIX]{x1D702}(x,z)$ instead of functions of $z$ alone, we deduce the improved formula

(2.50)$$\begin{eqnarray}Nu-1=\max _{\unicode[STIX]{x1D70E}}\tilde{{\mathcal{M}}}[\unicode[STIX]{x1D70E},\boldsymbol{v};Pe],\end{eqnarray}$$

where

(2.51)$$\begin{eqnarray}\tilde{{\mathcal{M}}}[\unicode[STIX]{x1D70E},\boldsymbol{v};Pe]=\frac{\left(\langle v_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}}{Pe^{-2}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}|^{2}\rangle \langle |\unicode[STIX]{x1D735}\boldsymbol{v}|^{2}\rangle +\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E5}^{-1}\left(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}\right)\right|^{2}\rangle }.\end{eqnarray}$$

Stated succinctly, the Howard–Busse–Malkus problem maximizes (2.49) and the wall-to-wall problem maximizes (2.51). The crucial difference between (2.49) and (2.51) lies in the fact that

(2.52)$$\begin{eqnarray}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E5}^{-1}(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E})|^{2}\rangle =\langle (\overline{v_{3}\unicode[STIX]{x1D70E}}-\langle v_{3}\unicode[STIX]{x1D70E}\rangle )^{2}\rangle +\text{non-negative terms}\end{eqnarray}$$

since $\langle \left(\overline{v_{3}\unicode[STIX]{x1D70E}}-\langle v_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}\rangle$ is just the zeroth Fourier mode contribution to $\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E5}^{-1}(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E})|^{2}\rangle$. Hence,

(2.53)$$\begin{eqnarray}\tilde{{\mathcal{M}}}[\unicode[STIX]{x1D70E},\boldsymbol{v};Pe]\leqslant {\mathcal{M}}[\unicode[STIX]{x1D70E},\boldsymbol{v};Pe]\end{eqnarray}$$

for all $(\unicode[STIX]{x1D70E},\boldsymbol{v})$ and the inequality is strict in most cases. Bounds on transport obtained using the functional $\tilde{{\mathcal{M}}}$ are, therefore, at least as tight as those that have been obtained using Howard’s functional ${\mathcal{M}}$. Heuristically, this occurs because the wall-to-wall problem is more constrained (and hence provides a tighter bound) since it fully enforces the steady advection–diffusion equation whereas the Howard–Busse–Malkus problem is less constrained since it only enforces the horizontally averaged version. We take (2.52) as additional evidence of the conjectured duality gap (2.33) between the wall-to-wall problem and the background method/Howard–Busse–Malkus approach.

3.1 Gradient ascent methods

Here, we outline various methods for solving the Euler–Lagrange equations of the wall-to-wall problem described above. One method is to evolve

(3.13a-e)$$\begin{eqnarray}0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}},\quad 0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}},\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\boldsymbol{u}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\boldsymbol{u}},\quad 0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}p},\quad \text{and}\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D707}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D707}}\end{eqnarray}$$

forward in time. This strictly enforces the advection–diffusion equation and its adjoint at every time step, utilizing $\unicode[STIX]{x1D6FF}{\mathcal{F}}/\unicode[STIX]{x1D6FF}\boldsymbol{u}$ to compute corrections to the flow field. This approach was taken in Motoki et al. (Reference Motoki, Kawahara and Shimizu2018b), and the rate limiting step is the solution of the advection–diffusion equation and its adjoint. In the present work, we take a different approach and compute numerical solutions to (3.3) and (3.8) utilizing two different algorithms.

The first of our algorithms involves a time-stepping procedure of the form

(3.14a-e)$$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D703}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}},\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D711}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}},\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\boldsymbol{u}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\boldsymbol{u}},\quad 0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}p},\quad \text{and}\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D707}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D707}}.\end{eqnarray}$$

This procedure may be understood as follows. For fixed $\boldsymbol{u}$, the equations for $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D711}$ evolve towards the steady state solutions $\unicode[STIX]{x1D6FF}{\mathcal{F}}/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}=0$ and $\unicode[STIX]{x1D6FF}{\mathcal{F}}/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}=0$, hence, these evolutions are a relaxation of fully solving the Euler–Lagrange equations for $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$. The equation for $\unicode[STIX]{x1D707}$ guarantees that we flow towards a flow field with the desired enstrophy $Pe$. The last condition (for a fixed $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$) evolves to a solution of the optimality condition. We enforce incompressibility at every time step and evolve all fields at once.

There is an alternative description of this algorithm in terms of ${\mathcal{S}}$. Focusing on the $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$ equations, we see that evolving

(3.15a,b)$$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D703}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}\quad \text{and}\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D711}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}}\end{eqnarray}$$

is equivalent to evolving

(3.16a,b)$$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D709}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{S}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D709}}\quad \text{and}\quad -\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D702}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{S}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D702}}\end{eqnarray}$$

after taking sums and differences (making use of $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D709}+\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D711}=\unicode[STIX]{x1D709}-\unicode[STIX]{x1D702}$) and rescaling time. Thus, our time-stepping procedure for $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$ is equivalent to simultaneously applying gradient ascent for the concave variable $\unicode[STIX]{x1D709}$ and gradient descent for the convex variable $\unicode[STIX]{x1D702}$ in the ${\mathcal{S}}$ functional. This is similar to the philosophy of Wen et al. (Reference Wen, Chini, Kerswell and Doering2015). In fact, the only modification required to compute two-dimensional no-slip background fields would be to project the $\unicode[STIX]{x1D702}$ variable to the zeroth Fourier mode at each time step, by taking the $x$ average of the right-hand side of the $\unicode[STIX]{x1D702}$ equation.

The second time-stepping method involves computing the saddles of ${\mathcal{S}}$ via

(3.17a-d)$$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D709}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D709}},\quad 0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D702}},\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\boldsymbol{u}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\boldsymbol{u}},\quad \text{and}\quad 0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}p},\end{eqnarray}$$

while fixing $\unicode[STIX]{x1D707}$. The resulting enstrophy depends implicitly on $\unicode[STIX]{x1D707}$. Enforcing $\unicode[STIX]{x1D6FF}{\mathcal{F}}/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D702}=0$ at each time step yields gradient ascent for the local optima of (2.6), an unconstrained variational problem. The reason why this algorithm is efficient is due to the many existing algorithms for quick inversion of the Laplacian and the Helmholtz operator.

We found numerically that both (3.14) and (3.17) yield the same result. We also implemented various other ascent procedures with different preconditioners. For example, we evolved equations of the form

(3.18a-e)$$\begin{eqnarray}-\unicode[STIX]{x0394}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D703}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}},\quad -\unicode[STIX]{x0394}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D711}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}},\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\boldsymbol{u}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\boldsymbol{u}},\quad 0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}p},\quad \text{and}\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D707}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D707}}\end{eqnarray}$$

forward in time, as well as

(3.19a-e)$$\begin{eqnarray}-\unicode[STIX]{x0394}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D703}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}},\quad -\unicode[STIX]{x0394}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D711}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}},\quad -\unicode[STIX]{x0394}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\boldsymbol{u}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\boldsymbol{u}},\quad 0=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}p},\quad \text{and}\quad \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D707}=\frac{\unicode[STIX]{x1D6FF}{\mathcal{F}}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D707}},\end{eqnarray}$$

evolving different components on different time scales. Amongst all of our results, the ones presented in § 4 maximize $Nu$ for a given $Pe$.

3.2 Temporal discretization

Each of the gradient ascent procedures described in § 3 are of the form

(3.20)$$\begin{eqnarray}{\dot{x}}={\mathcal{L}}x+{\mathcal{N}}(x)+f,\end{eqnarray}$$

where $x$ is the state vector, ${\mathcal{L}}$ is an ‘easily’ invertible linear operator, ${\mathcal{N}}$ is a nonlinear operator and $f$ is a forcing function. We follow Viswanath & Tobasco (Reference Viswanath and Tobasco2013) and consider linear multi-step schemes as follows:

(3.21)$$\begin{eqnarray}\frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\left(\unicode[STIX]{x1D6FE}x^{n+1}+\mathop{\sum }_{j=0}^{s-1}a_{j}x^{n-j}\right)=\mathop{\sum }_{j=0}^{s-1}b_{j}{\mathcal{N}}(x^{n-j})+{\mathcal{L}}x^{n+1}+f^{n+1},\end{eqnarray}$$

where $s$ is the order of the time-stepping scheme, $a_{i}$ and $b_{i}$ are parameters and $\unicode[STIX]{x0394}t$ is the time-step size. The parameters values for orders $s=1,2$ and $3$ are

(3.22a-d)$$\begin{eqnarray}s=1,\quad \unicode[STIX]{x1D6FE}=1,\quad a_{0}=-1,\quad b_{0}=1,\end{eqnarray}$$

(3.23a-f)$$\begin{eqnarray}s=2,\quad \unicode[STIX]{x1D6FE}=3/2,\quad a_{0}=-2,\quad a_{1}=1/2,\quad b_{0}=2,\quad b_{1}=-1,\end{eqnarray}$$

(3.24a-h)$$\begin{eqnarray}s=3,\quad \unicode[STIX]{x1D6FE}=11/6,\quad a_{0}=-3,\quad a_{1}=3/2,\quad a_{2}=-1/3,\quad b_{0}=3,\quad b_{1}=-3,\quad b_{2}=1.\end{eqnarray}$$

For example, with $s=1$ and the advection–diffusion equation,

(3.25)$$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D703}=-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}+u_{3}\end{eqnarray}$$

we use

(3.26)$$\begin{eqnarray}\left(\unicode[STIX]{x1D6E5}-\frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\boldsymbol{I}\right)\unicode[STIX]{x1D703}^{n+1}=\boldsymbol{u}^{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}^{n}-w^{n}-\frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D703}^{n}.\end{eqnarray}$$

Thus for each time step we must solve a modified Poisson’s equation of the form

(3.27)$$\begin{eqnarray}\left(\unicode[STIX]{x1D6E5}-c\boldsymbol{I}\right)\unicode[STIX]{x1D703}=f,\end{eqnarray}$$

where $c\geqslant 0$ and we have made the transcriptions

(3.28a-c)$$\begin{eqnarray}\unicode[STIX]{x1D703}^{n+1}\mapsto \unicode[STIX]{x1D703},\quad \frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\mapsto c,\quad \text{and}\quad \boldsymbol{u}^{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}^{n}-w^{n}-\frac{1}{\unicode[STIX]{x0394}t}\unicode[STIX]{x1D703}^{n}\mapsto f.\end{eqnarray}$$

Analogous discretizations are used for $\unicode[STIX]{x1D711},\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D709}$. For updating the optimality condition with $s=1$ one option is to use

(3.29)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\unicode[STIX]{x1D707}\unicode[STIX]{x1D6E5}-\frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\boldsymbol{I}\right)\boldsymbol{u}^{n+1}=-\left(\unicode[STIX]{x1D711}^{n}+\unicode[STIX]{x1D703}^{n}\right)\hat{e}_{3}+\unicode[STIX]{x1D711}^{n}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}^{n}-\frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\boldsymbol{u}^{n}+\unicode[STIX]{x1D735}p^{n+1}, & \displaystyle\end{eqnarray}$$

(3.30)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{n+1}=0. & \displaystyle\end{eqnarray}$$

Each time step involves solving a modified Stokes equation

(3.31)$$\begin{eqnarray}\displaystyle & \displaystyle \left(\unicode[STIX]{x1D6E5}-c\boldsymbol{I}\right)\boldsymbol{u}=\boldsymbol{f}+\unicode[STIX]{x1D735}p, & \displaystyle\end{eqnarray}$$

(3.32)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $c\geqslant 0$ and we have made the transcriptions

(3.33a-d)$$\begin{eqnarray}\boldsymbol{u}^{n+1}\mapsto \boldsymbol{u},\quad p^{n+1}\mapsto p,\quad \frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\mapsto c,\quad \text{and}\quad \left(\unicode[STIX]{x1D711}^{n}+\unicode[STIX]{x1D703}^{n}\right)\hat{e}_{3}+\unicode[STIX]{x1D711}^{n}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}^{n}-\frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}\boldsymbol{u}^{n}\mapsto \boldsymbol{f}.\end{eqnarray}$$

We solve these boundary value problems by utilizing a pseudo-spectral method, where the wall-bounded direction is represented by Chebyshev polynomials and the periodic directions are represented by Fourier series; however, our use of Chebyshev polynomials utilizes spectral integration in the same way as Viswanath (Reference Viswanath2015). The Stokes equation (3.31) was solved utilizing the Kleiser–Schumann algorithm (Kleiser & Schumann Reference Kleiser, Schumann and Hirschel1980). See appendix A for details regarding implementation. A theoretical discussion of the approach is in Viswanath (Reference Viswanath2015). Let us highlight some of the benefits of spectral integration here:

1. (i) Memory efficiency: the total memory occupied is linear in the number of points required to describe the state variable $\boldsymbol{u}$.

2. (ii) Speed: solving the linear equations involves only tridiagonal matrices which are put into lower–upper form at the beginning of the gradient ascent procedure.

3. (iii) Discretization accuracy: everything is represented using Fourier and Chebyshev modes, allowing for spectral accuracy.

4. (iv) Machine accuracy: utilizing spectral integration in factored form (see appendix A) allows us to avoid taking wall-bounded derivatives and only take derivatives in periodic direction, without passing to the nodal domain.

At this point, all the pieces are in place to compute flow fields maximizing heat transfer. By utilizing gradient ascent or other time-marching schemes of § 3.1, it is possible to adapt existing Rayleigh–Bénard codes to find steady maximizing flow fields and compare their thermal transfer to that of natural flows.

4.1 Singular value decomposition

Given the structure of the solutions depicted in figure 1 one may wonder if the corresponding fields are to first approximation separable. With regards to the streamfunction $\unicode[STIX]{x1D713}$ defined by $(-\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D713},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713})=(u_{1},u_{3})$, this would mean

(4.1)$$\begin{eqnarray}\unicode[STIX]{x1D713}(x,z)\approx \unicode[STIX]{x1D6F9}(z)\unicode[STIX]{x1D719}(x).\end{eqnarray}$$

We investigate this possibility by minimizing the functional

(4.2)$$\begin{eqnarray}A[\unicode[STIX]{x1D6F9},\unicode[STIX]{x1D719}]=\frac{1}{\unicode[STIX]{x1D6E4}}\int _{0}^{\unicode[STIX]{x1D6E4}}\int _{0}^{1}(\unicode[STIX]{x1D713}(x,z)-\unicode[STIX]{x1D6F9}(z)\unicode[STIX]{x1D719}(x))^{2}\,\text{d}x\,\text{d}z\end{eqnarray}$$

subject to appropriate constraints on $\unicode[STIX]{x1D6F9}(z)$ and $\unicode[STIX]{x1D719}(x)$, and similarly for other fields such as $\unicode[STIX]{x1D702}$ or $\unicode[STIX]{x1D709}$. (For example, $\unicode[STIX]{x1D719}(x)$ should be periodic and $\unicode[STIX]{x1D6F9}(z)$ should satisfy no-slip boundary conditions.) If the minimal value is sufficiently small, then we may say that the flow fields are nearly separable.

Upon discretization,

(4.3)$$\begin{eqnarray}\displaystyle A[\unicode[STIX]{x1D6F9},\unicode[STIX]{x1D719}]\approx \mathop{\sum }_{i,j}(\unicode[STIX]{x1D713}_{ij}-\unicode[STIX]{x1D6F9}_{j}\unicode[STIX]{x1D719}_{i})^{2}w_{ij},\end{eqnarray}$$

(4.4a-c)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{ij}=\unicode[STIX]{x1D713}(x_{i},z_{j}),\quad \unicode[STIX]{x1D719}_{i}=\unicode[STIX]{x1D719}(x_{i}),\quad \text{and}\quad \unicode[STIX]{x1D6F9}_{j}=\unicode[STIX]{x1D6F9}(z_{j}),\end{eqnarray}$$

where $x_{i}=\unicode[STIX]{x1D6E4}i/n$ for $i=0,\ldots ,n-1$ and $z_{j}={\textstyle \frac{1}{2}}(1+\cos (\unicode[STIX]{x03C0}j/m))$ for $j=0,\ldots ,m$ are the collocation points of the Fourier and Chebyshev discretizations and $w_{ij}$ are weights for approximating the integral by a discrete sum. Spectral accuracy for the approximation (4.3) was obtained by setting

(4.5a-c)$$\begin{eqnarray}w_{ij}=\unicode[STIX]{x0394}x_{i}\unicode[STIX]{x0394}z_{j}\quad \text{with}\quad \unicode[STIX]{x0394}x_{i}=\frac{1}{n}\quad \text{and}\quad \unicode[STIX]{x0394}z_{j}=\frac{1}{2}\sin (a_{j})\frac{1}{m}\mathop{\sum }_{\ell =1}^{m-1}\sin (\ell a_{j})\frac{(1-\cos (\ell \unicode[STIX]{x03C0}))}{\ell },\end{eqnarray}$$

where $a_{j}=\unicode[STIX]{x03C0}j/m$ for $j=0,\ldots ,m$, see Boyd (Reference Boyd2001). Since we are dealing with the average value of the integral there is no factor of $\unicode[STIX]{x1D6E4}$ in $\unicode[STIX]{x0394}x_{i}$.

In the discrete problem corresponding to (4.3) we minimized the weighted Frobenius norm of the matrix $\unicode[STIX]{x1D713}_{ij}$ with respect to an outer product decomposition, hence, for $w_{ij}=1$, the solution to the discrete problem is to take $\unicode[STIX]{x1D6F9}$ and $\unicode[STIX]{x1D719}$ to be equal to the largest singular vectors of the matrix $\unicode[STIX]{x1D713}_{ij}$ multiplied by the largest singular value. For $w_{ij}\neq 1$, we rescaled the problem taking $\unicode[STIX]{x1D713}_{ij}^{\prime }=\unicode[STIX]{x1D713}_{ij}\sqrt{w_{ij}}$ and found the largest singular vectors of $\unicode[STIX]{x1D713}_{ij}^{\prime }$. The difference between uniform and non-uniform $w_{ij}$ may be interpreted as minimizing with respect to different weighted energy norms and, in our numerical experiments, we found little qualitative difference between utilizing one over the other. All singular value decompositions reported here were computed using $w_{ij}=1$.

Figure 3 shows the first three modes in the singular value decomposition of the streamfunction $\unicode[STIX]{x1D713}$. Figure 5 shows the same for the $\unicode[STIX]{x1D709}$ field and the vertical velocity field $u_{3}$. The modes with smaller singular values may be viewed as a preliminary manifestation of a branching-like pattern, perhaps similar to the ones constructed in Tobasco & Doering (Reference Tobasco and Doering2017) and Doering & Tobasco (Reference Doering and Tobasco2019). The structures corresponding to the largest singular value are the direct analogue for the wall-to-wall problem of the single wavenumber solutions for the Howard–Busse–Malkus problem produced in Howard (Reference Howard1963).

As it turns out, these dominant structures carry ${\approx}99\%$ of the overall heat transport in the computed range of $Pe$. More precisely, writing

(4.6)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}(x,z)\approx \unicode[STIX]{x1D6F9}(z)\unicode[STIX]{x1D719}^{(1)}(x), & \displaystyle\end{eqnarray}$$

(4.7)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}(x,z)\approx \unicode[STIX]{x1D6EF}(z)\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}^{(2)}(x), & \displaystyle\end{eqnarray}$$

(4.8)$$\begin{eqnarray}\displaystyle & \displaystyle N_{1}=\langle (\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713})\unicode[STIX]{x1D709}\rangle , & \displaystyle\end{eqnarray}$$

(4.9)$$\begin{eqnarray}\displaystyle & \displaystyle N_{2}=\langle \unicode[STIX]{x1D6F9}(\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}^{(1)})\unicode[STIX]{x1D6EF}\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}^{(2)}\rangle , & \displaystyle\end{eqnarray}$$

our solutions achieved $(N_{1}-N_{2})/N_{1}\leqslant 0.01$ uniformly over all computed Péclet numbers, see figure 4. (This error was even less for smaller Péclet numbers.) The fact that (4.8) is equal to $Nu-1$ for solutions to the wall-to-wall problem follows from the Euler–Lagrange equation (3.8).

Figure 3. The streamfunction $\unicode[STIX]{x1D713}$ and its first three modes at $Pe\approx 2.4\times 10^{4}$. The modes are ordered left to right by the magnitude of their singular values starting with the largest.

Figure 4. The quantity $|N_{1}-N_{2}|/N_{2}\times 100$ as a function of $Pe$. The value remains below $1\%$ throughout the range of numerically computed values of $Pe$.

Figure 5. The field $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D703}+\unicode[STIX]{x1D711}$ (top row), the vertical velocity $u_{3}$ (bottom row) and their first three modes at $Pe\approx 2.4\times 10^{4}$. The outer products are ordered with respect to their singular values.

Figure 6. The first singular vectors in the singular value decomposition of the $\unicode[STIX]{x1D709}$ field (a,b) and the vertical velocity $u_{3}$ (c,d). The respective products of these singular vectors yields the first approximations in figure 5. The Péclet number ${\approx}2.4\times 10^{4}$.

Figure 6 depicts the first singular vectors for $\unicode[STIX]{x1D709}$ and $u_{3}$. While in the Howard–Busse–Malkus problem the corresponding $x$ dependence is perfectly sinusoidal, for the wall-to-wall problem the $x$ dependences of the first singular vectors resemble Jacobi elliptic functions. The similarity between $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}^{(1)}$ and $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}^{(2)}$ in figure 6 motivates the ‘separable ansatz’

(4.10)$$\begin{eqnarray}\displaystyle & \displaystyle u_{1}(x,z)=-\unicode[STIX]{x1D6F9}^{\prime }(z)\unicode[STIX]{x1D719}(x), & \displaystyle\end{eqnarray}$$

(4.11)$$\begin{eqnarray}\displaystyle & \displaystyle u_{3}(x,z)=\unicode[STIX]{x1D6F9}(z)\unicode[STIX]{x1D719}^{\prime }(x), & \displaystyle\end{eqnarray}$$

(4.12)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}(x,z)=\unicode[STIX]{x1D6EF}(z)\unicode[STIX]{x1D719}^{\prime }(x), & \displaystyle\end{eqnarray}$$

which we consider further in § 5.1 below. There, we derive conditional upper bounds extending the single wavenumber analysis from Howard (Reference Howard1963) to the wall-to-wall problem. We will see that the functional from § 2.6 can be bounded from above by $Pe^{6/11}=Pe^{0.\overline{54}}$ amongst all separable ansatzes, in accord with the numerically computed $Nu\sim Pe^{0.54}$ scaling.

5.1 Upper bounds within a separable ansatz

To obtain upper bounds on the Nusselt number within a separable ansatz, we start by recalling that the Howard–Busse–Malkus problem bounds the wall-to-wall problem from above. Indeed, as was shown in §§ 2.3 and 2.6,

(5.1)$$\begin{eqnarray}\max _{\boldsymbol{u}}Nu-1=\max _{\boldsymbol{u},\unicode[STIX]{x1D709}}\min _{\unicode[STIX]{x1D702}}{\mathcal{S}}=\max _{\boldsymbol{v},\unicode[STIX]{x1D70E}}\tilde{{\mathcal{M}}}\leqslant \max _{\boldsymbol{v},\unicode[STIX]{x1D70E}}{\mathcal{M}}\end{eqnarray}$$

since the minimum is only made smaller upon enlarging the class of admissible functions. Thus, we are led to consider the functional

(5.2)$$\begin{eqnarray}{\mathcal{M}}[\unicode[STIX]{x1D70E},\boldsymbol{v},Pe]=\frac{\left(\langle v_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}}{Pe^{-2}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}|^{2}\rangle \langle |\unicode[STIX]{x1D735}\boldsymbol{v}|^{2}\rangle +\langle \left(\overline{v_{3}\unicode[STIX]{x1D70E}}-\langle v_{3}\unicode[STIX]{x1D70E}\rangle \right)^{2}\rangle }\end{eqnarray}$$

from § 2.6, while also restricting ourselves to the separable ansatz

(5.3)$$\begin{eqnarray}\displaystyle & \displaystyle u_{1}(x,z)=-\unicode[STIX]{x1D6F9}^{\prime }(z)\unicode[STIX]{x1D719}(x), & \displaystyle\end{eqnarray}$$

(5.4)$$\begin{eqnarray}\displaystyle & \displaystyle u_{3}(x,z)=\unicode[STIX]{x1D6F9}(z)\unicode[STIX]{x1D719}^{\prime }(x), & \displaystyle\end{eqnarray}$$

(5.5)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}(x,z)=\unicode[STIX]{x1D6EF}(z)\unicode[STIX]{x1D719}^{\prime }(x), & \displaystyle\end{eqnarray}$$

where $^{\prime }$ denotes an ordinary derivative. Note, while the argument that follows is similar to those of Howard (Reference Howard1963) and Doering & Constantin (Reference Doering and Constantin1996), it is more general in scope since here we allow $\unicode[STIX]{x1D719}(x)$ to be any periodic function rather than some well-chosen Fourier mode. Our choice to study this more general case is motivated by the previous discussion of numerical results, which showed how, in the wall-to-wall problem, non-sinusoidal $\unicode[STIX]{x1D719}(x)$ arise. Our goal now is to bound the functional ${\mathcal{M}}$ from above under the assumptions of separability (5.3)–(5.5). Again, we note that this does not provide the unconditional upper bound on ${\mathcal{M}}$ required to deduce bounds on $\text{max }Nu$ according to (5.1). What we lack is a proof that (5.3)–(5.5) are indeed valid assumptions over the given range of $Pe$. Nevertheless, we proceed.

Note that the velocity field is incompressible, and we normalize $\unicode[STIX]{x1D719}$ by requiring

(5.6)$$\begin{eqnarray}\overline{(\unicode[STIX]{x1D719}^{\prime })^{2}}=1.\end{eqnarray}$$

Given the homogeneity of the functional we may impose that

(5.7)$$\begin{eqnarray}\langle u_{3}\unicode[STIX]{x1D709}\rangle =\int _{0}^{1}\unicode[STIX]{x1D6F9}(z)\unicode[STIX]{x1D6EF}(z)\,\text{d}z=1.\end{eqnarray}$$

Given these, we seek to bound the denominator of (5.2) from below so as to produce the desired upper bound.

Note the identities

(5.8)$$\begin{eqnarray}\displaystyle & \displaystyle \langle u_{3}\unicode[STIX]{x1D709}\rangle =\int _{0}^{1}\unicode[STIX]{x1D6F9}(z)\unicode[STIX]{x1D6EF}(z)\,\text{d}z, & \displaystyle\end{eqnarray}$$

(5.9)$$\begin{eqnarray}\displaystyle & \displaystyle \langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle =\int _{0}^{1}\left(\overline{\unicode[STIX]{x1D719}^{2}}(\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}+2(\unicode[STIX]{x1D6F9}^{\prime })^{2}+\overline{(\unicode[STIX]{x1D719}^{\prime \prime })^{2}}(\unicode[STIX]{x1D6F9})^{2}\right)\,\text{d}z, & \displaystyle\end{eqnarray}$$

(5.10)$$\begin{eqnarray}\displaystyle & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}\rangle =\int _{0}^{1}(\unicode[STIX]{x1D6EF}^{\prime })^{2}+\overline{(\unicode[STIX]{x1D719}^{\prime \prime })^{2}}\unicode[STIX]{x1D6EF}^{2}\,\text{d}z, & \displaystyle\end{eqnarray}$$

(5.11)$$\begin{eqnarray}\displaystyle & \displaystyle \langle (\overline{u_{3}\unicode[STIX]{x1D709}}-\langle u_{3}\unicode[STIX]{x1D709}\rangle )^{2}\rangle =\int _{0}^{1}\left(\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6EF}-1\right)^{2}\,\text{d}z, & \displaystyle\end{eqnarray}$$

where for convenience we have abbreviated $\int _{0}^{1}f(z)\,\text{d}z$ as $\int f$. Using these, we obtain the lower bound

(5.12)$$\begin{eqnarray}\displaystyle \langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}\rangle & {\geqslant} & \displaystyle \overline{\unicode[STIX]{x1D719}^{2}}\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}+\overline{\left(\unicode[STIX]{x1D719}^{\prime \prime }\right)^{2}}\int \unicode[STIX]{x1D6F9}^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}+\overline{\unicode[STIX]{x1D719}^{2}}\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\overline{\left(\unicode[STIX]{x1D719}^{\prime \prime }\right)^{2}}\int \unicode[STIX]{x1D6EF}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\left(\overline{(\unicode[STIX]{x1D719}^{\prime \prime })^{2}}\right)^{2}\int \unicode[STIX]{x1D6F9}^{2}\int \unicode[STIX]{x1D6EF}^{2}\end{eqnarray}$$

(5.13)$$\begin{eqnarray}\displaystyle & {\geqslant} & \displaystyle \overline{\unicode[STIX]{x1D719}^{2}}\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}+2\sqrt{\overline{\left(\unicode[STIX]{x1D719}^{\prime \prime }\right)^{2}}\int \unicode[STIX]{x1D6F9}^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}\overline{\unicode[STIX]{x1D719}^{2}}\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\overline{\left(\unicode[STIX]{x1D719}^{\prime \prime }\right)^{2}}\int \unicode[STIX]{x1D6EF}^{2}}\nonumber\\ \displaystyle & & \displaystyle +\,\left(\overline{(\unicode[STIX]{x1D719}^{\prime \prime })^{2}}\right)^{2}\int \unicode[STIX]{x1D6F9}^{2}\int \unicode[STIX]{x1D6EF}^{2},\end{eqnarray}$$

by an elementary Young’s inequality. The interpolation inequality

(5.14)$$\begin{eqnarray}\left[\overline{(\unicode[STIX]{x1D719}^{\prime \prime })^{2}}\right]\left[\overline{\unicode[STIX]{x1D719}^{2}}\right]\geqslant \left(\overline{\unicode[STIX]{x1D719}\unicode[STIX]{x1D719}^{\prime \prime }}\right)^{2}=\left(\overline{(\unicode[STIX]{x1D719}^{\prime })^{2}}\right)^{2}=1\end{eqnarray}$$

follows by combining Cauchy–Schwarz, integration by parts, and the normalization (5.6). Similarly,

(5.15)$$\begin{eqnarray}\int \unicode[STIX]{x1D6EF}^{2}\int \unicode[STIX]{x1D6F9}^{2}\geqslant 1,\end{eqnarray}$$

by Cauchy–Schwarz and the net flux constraint (5.7).

Proceeding, we deduce from (5.14) and (5.15) that

(5.16)$$\begin{eqnarray}\displaystyle \langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|^{2}\rangle & {\geqslant} & \displaystyle \frac{1}{\overline{(\unicode[STIX]{x1D719}^{\prime \prime })^{2}}}\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\!\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}+2\sqrt{\overline{\left(\unicode[STIX]{x1D719}^{\prime \prime }\right)^{2}}\!\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\!\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}}+(\overline{(\unicode[STIX]{x1D719}^{\prime \prime })^{2}})^{2}\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(5.17)$$\begin{eqnarray}\displaystyle & {\geqslant} & \displaystyle \inf _{a>0}\left\{\frac{1}{a^{2}}\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}+a\sqrt{4\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}}+a^{4}\right\}\end{eqnarray}$$

(5.18)$$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{9}{2^{4/3}}\left(\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}\right)^{2/3}.\end{eqnarray}$$

Utilizing Howard’s lemma (the estimate referred to as such in Doering & Constantin (Reference Doering and Constantin1996), equation (6.26)), we can bound the remaining term by

(5.19)$$\begin{eqnarray}\int \left(\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D6EF}-1\right)^{2}\geqslant \frac{16}{15}\left(\frac{3}{2}\right)^{1/2}\left(\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}\right)^{-1/4}.\end{eqnarray}$$

Taking $\unicode[STIX]{x1D6FF}=\left(\int (\unicode[STIX]{x1D6F9}^{\prime \prime })^{2}\int (\unicode[STIX]{x1D6EF}^{\prime })^{2}\right)^{-1/4}$, we see that the denominator of the functional ${\mathcal{M}}$ given in (5.2) – in the separable ansatz – is bounded below by

(5.20)$$\begin{eqnarray}\frac{16}{15}\left(\frac{3}{2}\right)^{1/2}\unicode[STIX]{x1D6FF}+\frac{9Pe^{-2}}{2^{4/3}}\frac{1}{\unicode[STIX]{x1D6FF}^{8/3}}\geqslant \left(\frac{16}{15}\left(\frac{3}{2}\right)^{1/2}\right)^{8/11}\left(\frac{9Pe^{-2}}{2^{4/3}}\right)^{3/11}\left(\left(\frac{3}{8}\right)^{-3/11}+\left(\frac{3}{8}\right)^{8/11}\right)\end{eqnarray}$$

after minimizing over $\unicode[STIX]{x1D6FF}$.

Thus, ${\mathcal{M}}$ is bounded from above by a constant multiple of $Pe^{6/11}=Pe^{0.\overline{54}}$ in the separable ansatz (5.3)–(5.5), in accord with the numerically computed transport scaling of $Pe^{0.54}$. While this is suggestive of the bound $Nu-1\leqslant 0.324Pe^{6/11}$ (roughly a factor of two larger than the numerical prefactor ${\approx}0.15$), it remains to be shown that the separable ansatz holds up to terms negligible in their transport for the given range of $Pe$.