2 Mathematical formulation

2.1 Set-up

We consider the two-dimensional version of the Boussinesq equations (1.1) and (1.2) where $\boldsymbol{u}=u\hat{\boldsymbol{x}}+w\hat{\boldsymbol{z}}$ over a box $(x,z)\in [-{\textstyle \frac{1}{2}}L,{\textstyle \frac{1}{2}}L]\times [0,1]$ together with the following stress-free and isothermal boundary conditions

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}=w=0,\quad T=1,\quad \text{at }z=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}=w=0,\quad T=0,\quad \text{at }z=1, & \displaystyle\end{eqnarray}$$

following Hassanzadeh et al. (Reference Hassanzadeh, Chini and Doering2014). Applying the background method, we decompose the temperature field as

(2.3) $$\begin{eqnarray}T=\unicode[STIX]{x1D70F}(x,z)+\unicode[STIX]{x1D703}(x,z,t),\end{eqnarray}$$

where the (steady) background temperature $\unicode[STIX]{x1D70F}$ carries the boundary conditions of $T$ (i.e. $\unicode[STIX]{x1D70F}|_{z=0}=1$ and $\unicode[STIX]{x1D70F}|_{z=1}=0$ ) so that the perturbation field $\unicode[STIX]{x1D703}$ vanishes at $z=0,1$ . The time-averaged heat transport is characterised by the time-averaged Nusselt number $Nu$

(2.4) $$\begin{eqnarray}Nu:=\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{T}\frac{1}{L}\int _{-L/2}^{L/2}\,\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}\biggl|_{z=1}\,\text{d}x\,\text{d}t=\langle |\unicode[STIX]{x1D735}T|^{2}\rangle =1+\langle wT\rangle ,\end{eqnarray}$$

where the spatial-temporal average defined in (1.5) becomes specifically

$$\begin{eqnarray}\langle \,(\ldots )\,\rangle :=\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{T}\int _{0}^{1}\frac{1}{L}\int _{-L/2}^{L/2}\,(\ldots )\,\text{d}x\,\text{d}z\,\text{d}t.\end{eqnarray}$$

To find the maximum heat transport possible over all solutions to the Boussinesq equations, we construct the Lagrangian

(2.5) $$\begin{eqnarray}\displaystyle \mathscr{L} & = & \displaystyle \langle |\unicode[STIX]{x1D735}T|^{2}\rangle -{\displaystyle \frac{a}{\unicode[STIX]{x1D70E}Ra}}\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\mathcal{N}}}\rangle -b\langle \unicode[STIX]{x1D703}\,{\mathcal{H}}\rangle ,\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & = & \displaystyle \langle |\unicode[STIX]{x1D735}T|^{2}\rangle -\frac{a}{\unicode[STIX]{x1D70E}Ra}\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\mathcal{N}}}\rangle -b\langle T\,{\mathcal{H}}\rangle +b\langle \unicode[STIX]{x1D70F}\,{\mathcal{H}}\rangle ,\end{eqnarray}$$

where $a/\unicode[STIX]{x1D70E}Ra$ is the Lagrange multiplier imposing the global constraint $\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\mathcal{N}}}\rangle =0$ , $b$ is a Lagrange multiplier imposing the global constraint $\langle T\,{\mathcal{H}}\rangle =0$ and $b\unicode[STIX]{x1D70F}(x,z)$ is the Lagrange multiplier field imposing the time-averaged heat equation pointwise in the domain. The inclusion of $b$ is actually redundant given the constraint imposed by $\unicode[STIX]{x1D70F}$ implies $\langle T\,{\mathcal{H}}\rangle =0$ so the value of $b$ is chosen for convenience. Expression (2.6) can be rewritten using integration by parts and the fact that $\langle wT\rangle =\mathscr{L}-1$ (see (2.4)) for solutions of the Boussinesq equations as

(2.7) $$\begin{eqnarray}\mathscr{L}=\frac{1}{1-a}[\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -a]-\frac{1}{1-a}\mathscr{G},\end{eqnarray}$$

where setting $b=2$ makes

(2.8) $$\begin{eqnarray}\mathscr{G}:=\left\langle \frac{a}{Ra}|\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D703}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}\right\rangle ,\end{eqnarray}$$

a purely quadratic form in $\unicode[STIX]{x1D703}$ and $\boldsymbol{u}$ .

2.2 Bounds on $Nu$

The key realisation is that if $\mathscr{G}\geqslant 0$ for all $(\boldsymbol{u},\unicode[STIX]{x1D703})\in \unicode[STIX]{x1D6F1}$ (the set of incompressible velocity and temperature fields which satisfy homogeneous versions of the boundary conditions (2.1) and (2.2)), which is a spectral constraint on $\unicode[STIX]{x1D70F}$ , and $a\in (0,1)$ , we then have the bound

(2.9) $$\begin{eqnarray}Nu\leqslant \frac{1}{1-a}[\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -a].\end{eqnarray}$$

The challenge is then to find the lowest such bound by minimising over all $(\unicode[STIX]{x1D70F},a)$ which satisfy this spectral constraint, i.e.

(2.10) $$\begin{eqnarray}(\unicode[STIX]{x1D70F},a)\in \unicode[STIX]{x1D6FA}:=\{(\unicode[STIX]{x1D70F},a)\,|\,\mathscr{G}(\boldsymbol{u},\unicode[STIX]{x1D703};\unicode[STIX]{x1D70F},a)\geqslant 0\,\forall (\boldsymbol{u},\,\unicode[STIX]{x1D703})\in \unicode[STIX]{x1D6F1}\}.\end{eqnarray}$$

After introducing a streamfunction, $(u,w)=(\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}z,-\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}x)$ , the constraint that

(2.11) $$\begin{eqnarray}\mathscr{G}=\left\langle \frac{a}{Ra}|\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}|^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D703}J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D713})\right\rangle \geqslant 0,\end{eqnarray}$$

where

(2.12) $$\begin{eqnarray}J(A,B):=\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}x}\end{eqnarray}$$

is equivalent to requiring that all of the eigenvalues $\unicode[STIX]{x1D706}$ of the linear problem

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}-J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D713}), & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}=\frac{a}{Ra}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}-J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(with boundary conditions $\unicode[STIX]{x1D713}=\text{d}^{2}\unicode[STIX]{x1D713}/\text{d}z^{2}=\unicode[STIX]{x1D703}=0$ for $z=\{0,1\}$ ) are negative semi-definite.

2.3 Convexity and uniqueness

The Euler–Lagrange equations for making the Lagrangian $\mathscr{L}$ in (2.7) stationary are

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}-J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D713}), & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{a}{Ra}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}-J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70F}-J(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713}), & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -1={\displaystyle \frac{(1-a)}{Ra}}\langle |\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}|^{2}\rangle & \displaystyle\end{eqnarray}$$

and, as a nonlinear set of equations, can have many solutions. However, only solutions with $(\unicode[STIX]{x1D70F},a)\in \unicode[STIX]{x1D6FA}$ yield a bound through the value of $\mathscr{L}$ generated. Due to the convexity of $\unicode[STIX]{x1D6FA}$ (i.e. if $(\unicode[STIX]{x1D70F}_{1},a_{1})$ and $(\unicode[STIX]{x1D70F}_{2},a_{2})$ are in $\unicode[STIX]{x1D6FA}$ then so is $\unicode[STIX]{x1D707}(\unicode[STIX]{x1D70F}_{1},a_{1})+(1-\unicode[STIX]{x1D707})(\unicode[STIX]{x1D70F}_{2},a_{2})$ for $\unicode[STIX]{x1D707}\in (0,1)$ ), and the fact that the objective functional

(2.19) $$\begin{eqnarray}f(\unicode[STIX]{x1D70F},a):=\frac{1}{1-a}[\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -a],\end{eqnarray}$$

to be minimised is a strictly convex functional (the terms second order in $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D6FF}a$ in the difference $f(\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F},a+\unicode[STIX]{x1D6FF}a)-f(\unicode[STIX]{x1D70F},a)$ , specifically

(2.20) $$\begin{eqnarray}\frac{1}{(1-a)^{2}}\langle \unicode[STIX]{x1D735}|(1-a)\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}+(\unicode[STIX]{x1D70F}+z-1)\unicode[STIX]{x1D6FF}a|^{2}\rangle ,\end{eqnarray}$$

are positive definite), there is in fact at most one solution which satisfies the spectral constraint. This solution, hereafter referred to as the optimal solution, is what is sought.

2.4 Numerical approach

Recently, Wen et al. (Reference Wen, Chini, Kerswell and Doering2015) have proved that when $\unicode[STIX]{x1D70F}$ is one-dimensional, i.e. $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}(z)$ , appropriately augmenting the (steady) Euler–Lagrange equations with time derivatives leads to a system where the optimal solution is a unique attracting steady state. This proof carries over to two-dimensional background fields $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}(x,z)$ in the three-dimensional Rayleigh–Bénard problem but not in the two-dimensional problem (see appendix A for details) where the dimensionality of the background field then matches that of the physical fields. This means that any steady attractor which emerges from time stepping using $\unicode[STIX]{x1D70F}(x,z)$ in the 2-D problem is not guaranteed to be the required optimal solution. The time stepping approach can still be used if it is married with a spectral constraint check but then there is always the prospect of rerunning with different initial conditions until the optimal solution is found. Given this, we chose instead to use the branch continuation approach – Newton’s method with parametric continuation – starting from the energy stability bifurcation point as performed in PK03. While very robust, this has the general disadvantage of only being able to continuously trace optimal solutions from the energy stability bifurcation as $Ra$ varies meaning that any new unconnected optimals cannot be found at a given $Ra$ . This is not a problem here as the aforementioned uniqueness of the optimal solution means that no other optimal solution branches exist.

We consider periodic boundary conditions in $x$ and, exactly as in Hassanzadeh et al. (Reference Hassanzadeh, Chini and Doering2014), assume that the streamfunction $\unicode[STIX]{x1D713}$ is odd (or antisymmetric), while $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D70F}$ are even (or symmetric) about $x=0$ by seeking the solution of (2.15)–(2.18) in the following form:

(2.21a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D713}=\mathop{\sum }_{m=1}^{M}\unicode[STIX]{x1D713}_{m}(z)\sin (m\unicode[STIX]{x1D6FC}x),\quad \unicode[STIX]{x1D703}=\mathop{\sum }_{m=0}^{M}\unicode[STIX]{x1D703}_{m}(z)\cos (m\unicode[STIX]{x1D6FC}x),\quad \unicode[STIX]{x1D70F}=\mathop{\sum }_{m=0}^{M}\unicode[STIX]{x1D70F}_{m}(z)\cos (m\unicode[STIX]{x1D6FC}x).\end{eqnarray}$$

We will find that this choice prevents a 1-D background optimal even though this is allowed by the boundary conditions and imposed symmetry. Here $\unicode[STIX]{x1D6FC}:=2\unicode[STIX]{x03C0}/L$ and $\unicode[STIX]{x1D713}_{m}$ , $\unicode[STIX]{x1D703}_{m}$ , $\unicode[STIX]{x1D70F}_{m}$ are expanded in Chebyshev polynomials, $T_{n}$ ,

(2.22) $$\begin{eqnarray}[\unicode[STIX]{x1D713}_{m}(z),\unicode[STIX]{x1D703}_{m}(z),\unicode[STIX]{x1D70F}_{m}(z)]=\mathop{\sum }_{n=0}^{N}[\unicode[STIX]{x1D713}_{mn},\unicode[STIX]{x1D703}_{mn},\unicode[STIX]{x1D70F}_{mn}]T_{n}(z),\end{eqnarray}$$

where $T_{n}(z):=\cos (n\cos ^{-1}(2z-1))$ . Resolution varies from $(N,M)=(30,30)$ to (80, 80) to ensure numerical accuracy as $Ra$ increases and $L$ changes.

Figure 1. (a) The Nusselt number versus the Rayleigh number tracked from the energy stability bifurcation point. The solution to the background upper bounding problem is optimised over the domain length and the data are in excellent agreement with Hassanzadeh et al. (Reference Hassanzadeh, Chini and Doering2014) (their data courtesy of Dr A. Souza). Insets show the flow streamfunctions at $Ra=10^{5}$ and 10⁶. (b) The optimal domain size $L^{\ast }$ versus the Rayleigh number. The bullet is the energy stability bifurcation point.

Figure 2. The first ( $\unicode[STIX]{x1D706}_{1}$ ) and second ( $\unicode[STIX]{x1D706}_{2}$ ) largest eigenvalues of the spectral constraint for $Ra<4468.8$ . At $Ra=4468.8$ , where they cross, the aspect ratio is $L^{\ast }=2.234$ .

3 Connecting to Hassanzadeh et al. (Reference Hassanzadeh, Chini and Doering2014)

The conductive temperature profile $\unicode[STIX]{x1D70F}=1-z$ is a spectrally stable background field until $Ra=27\unicode[STIX]{x03C0}^{4}/4$ in a domain of size $L=2\sqrt{2}$ when energy instability first starts. Ensuring that the marginal fluctuation fields $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})$ – hereafter a mode – stay marginal (called pinning as was done in PK03), the optimal solution was then tracked up to $Ra=10^{7}$ with the domain size $L=L^{\ast }(Ra)$ simultaneously optimised to yield the highest heat flux at a given $Ra$ (see figure 1). The calculated $Nu$ values correspond exactly with those found by Hassanzadeh et al. (Reference Hassanzadeh, Chini and Doering2014) (as do flow fields computed at $Ra=10^{5}$ and 10⁶; see the inset of figure 1). This indicates that Hassanzadeh et al.’s (Reference Hassanzadeh, Chini and Doering2014) wall-to-wall transport approach is equivalent to the background method when a single mode is considered.

Figure 3. At $Ra=4468.8$ , the first mode (a) $\unicode[STIX]{x1D713}_{1}$ and (b) $\unicode[STIX]{x1D703}_{1}$ and the new second mode (c) $\unicode[STIX]{x1D713}_{2}$ and (d) $\unicode[STIX]{x1D703}_{2}$ .

In their wall-to-wall optimal control approach, however, Hassanzadeh et al. (Reference Hassanzadeh, Chini and Doering2014) had no way of identifying whether their local optimal was in fact the global optimal. It should be sufficiently close to the energy stability point but experience in other related problems (e.g. PK03) suggests that further modes in the spectral constraint eventually become marginal as $Ra$ increases. The optimal solution should subsequently modify itself to keep these new modes marginal with concomitant adjustments in the $Nu$ -scaling. Fortunately, in the background approach, the spectral constraint provides a check on whether a given Euler–Lagrange solution is the optimal solution. Solving the eigenvalue problem (2.13)–(2.14) for disturbances which are also periodic over $[0,L^{\ast }(Ra)]$ demonstrates that the eigenvalue ( $\unicode[STIX]{x1D706}_{1}$ in figure 2) of the first mode is pinned at 0, while a second mode becomes marginal at $Ra=4468.8$ for an aspect ratio $L^{\ast }=2.234$ . This suggests that Hassanzadeh et al.’s (Reference Hassanzadeh, Chini and Doering2014) result is either not a bound for $Ra>4468.8$ or that the background bound has started to overestimate the actual maximal flux or, perhaps least likely, both. This divergence in results is apparently not because any further 2-D bifurcations have been missed in the wall-to-wall calculations (G. Chini, private communication 2019) but more a reflection of the ‘duality gap’ suggested by Souza et al. (Reference Souza, Tobasco and Doering2019) being realised. Figure 3(a,b) shows the first mode $(\unicode[STIX]{x1D713}_{1},\unicode[STIX]{x1D703}_{1})$ with wavenumber $\unicode[STIX]{x1D6FC}_{1}:=2\unicode[STIX]{x03C0}/L^{\ast }$ so the flow field contains one pair of convection cells. The second mode ( $\unicode[STIX]{x1D713}_{2},\unicode[STIX]{x1D703}_{2}$ ) with $\unicode[STIX]{x1D6FC}_{2}=2\unicode[STIX]{x1D6FC}_{1}$ illustrated in figure 3(c,d) has two pairs of convection cells. The optimal background field at $Ra=2000$ is shown in figure 4(a) and the now non-optimal 1-mode solution at $Ra=20\,000$ is shown in figure 4(b). In both cases the field is weakly two-dimensional, indicating that the first mode consists of non-monochromatic (i.e. non-single $\unicode[STIX]{x1D6FC}$ ) velocity and temperature fields. The emergence of the second mode at $Ra=4468.8$ indicates that the background profile is now degenerate in a way which has important implications for solving the Euler–Lagrange equations for higher $Ra$ while respecting the spectral constraint. We discuss this issue in the following section.

Figure 4. (a) The optimal background field $\unicode[STIX]{x1D70F}$ plotted at $Ra=2000$ and (b) the non-optimal 1-mode solution at $Ra=20\,000$ .

4 Multi-modal optimals

When a new mode becomes marginal in the spectral constraint as the background field $\unicode[STIX]{x1D70F}$ evolves with $Ra$ , a further ‘pinning’ constraint needs to be added to keep the new mode marginal in the spectral constraint as $Ra$ increases further. This procedure is thoroughly discussed in Doering & Constantin (Reference Doering and Constantin1996) and implemented in PK03 for a background field of lower dimensionality than the fluctuation field. In this situation, an example of which is using a 1-D case $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}(z)$ in the 2-D Rayleigh–Bénard problem, the fluctuation field can be Fourier transformed over the spatial dimension(s) across which $\unicode[STIX]{x1D70F}$ is invariant and then considered parameterised by the Fourier wavenumber $k$ . Different spectrally marginal fluctuation fields have different $k$ and are then naturally orthogonal under averaging over this spatial dimension. This means that the Euler–Lagrange equations (2.15)–(2.18),

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{j}-J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D713}_{j}),\quad j=1,\ldots ,N & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{a}{Ra}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}_{i}-J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703}_{i}),\quad j=1,\ldots ,N & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D70F}_{zz}-\mathop{\sum }_{j=1}^{N}\overline{J(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D713}_{j})}, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \langle |\unicode[STIX]{x1D70F}_{z}|^{2}\rangle -1={\displaystyle \frac{(1-a)}{Ra}}\mathop{\sum }_{j=1}^{N}\langle |\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}_{j}|^{2}\rangle , & \displaystyle\end{eqnarray}$$

(the overbar represents averaging over $x$ ) can simply be extended to include the new marginal mode

(4.5) $$\begin{eqnarray}(\unicode[STIX]{x1D703}_{N+1},\unicode[STIX]{x1D713}_{N+1})(x,z)=(\hat{\unicode[STIX]{x1D703}}_{N+1},\hat{\unicode[STIX]{x1D713}}_{N+1})(z)\text{e}^{\text{i}k_{N+1}x}\end{eqnarray}$$

when it appears. Equivalently, the Lagrangian is just

(4.6) $$\begin{eqnarray}\mathscr{L}=\frac{1}{1-a}[\langle {\unicode[STIX]{x1D70F}_{z}}^{2}\rangle -a]-\frac{1}{1-a}\mathscr{G},\end{eqnarray}$$

where $\mathscr{G}$ naturally partitions into the contributions from the various marginal modes as follows:

(4.7) $$\begin{eqnarray}\mathscr{G}=\left\langle \frac{a}{Ra}|\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}|^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D703}J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D713})\right\rangle =\mathop{\sum }_{j=1}^{N+1}\mathscr{G}_{j},\end{eqnarray}$$

where

(4.8) $$\begin{eqnarray}\mathscr{G}_{j}:=\left\langle \frac{a}{Ra}|\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}_{j}|^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{j}|^{2}+2\unicode[STIX]{x1D703}J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D713}_{j})\right\rangle .\end{eqnarray}$$

The appearance of a new spectrally marginal mode merely extends the set of wavenumbers contributing to the definition of the fluctuation field by one,

(4.9) $$\begin{eqnarray}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})(x,z):=\mathop{\sum }_{j=1}^{N+1}(\unicode[STIX]{x1D713}_{j},\unicode[STIX]{x1D703}_{j})=\mathop{\sum }_{j=1}^{N+1}(\hat{\unicode[STIX]{x1D713}}_{j},\hat{\unicode[STIX]{x1D703}}_{j})(z)\text{e}^{\text{i}k_{j}x}.\end{eqnarray}$$

Importantly, this means it is possible to talk about the unique optimal solution of the variational problem which satisfies the imposed physical constraints as being

(4.10) $$\begin{eqnarray}(\unicode[STIX]{x1D713},T)(x,z)=(0,\unicode[STIX]{x1D70F})(z)+\mathop{\sum }_{j=1}^{N}(\hat{\unicode[STIX]{x1D713}}_{j},\hat{\unicode[STIX]{x1D703}}_{j})(z)\text{e}^{\text{i}k_{j}x},\end{eqnarray}$$

i.e. the spectral constraint is satisfied at a saddle point of $\mathscr{L}$ .

This pleasing situation in which the marginal fluctuation fields have a physical interpretation changes, however, when the dimensionality of the background field equals the dimensionality of the problem (the case here), or, pathologically, there is more than one marginal mode for a given wavenumber (see chap. 3 of Fantuzzi Reference Fantuzzi2018). In these scenarios, the natural orthogonality property of different marginal fluctuation fields disappears with the result that the physical meaning of the fluctuation fields is lost. To see this, the key is to realise that pinning the marginal fluctuation fields is done (Doering & Constantin Reference Doering and Constantin1996) as before by writing the Lagrangian as

(4.11) $$\begin{eqnarray}\mathscr{L}=\frac{1}{1-a}(\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -a)-\frac{1}{1-a}\mathop{\sum }_{j=1}\mathscr{G}_{j}.\end{eqnarray}$$

The constraint that each $\mathscr{G}_{j}$ vanishes pins the $j$ th mode to be marginal (the Lagrange multiplier imposing this is absorbed into the amplitude of the $j$ th marginal fluctuation field) while $\mathscr{G}>0$ for all other fluctuation fields. However, since the modes $(\unicode[STIX]{x1D713}_{j},\unicode[STIX]{x1D703}_{j})$ are not now orthogonal,

(4.12) $$\begin{eqnarray}\mathop{\sum }_{j=1}^{N}\mathscr{G}_{j}\neq \mathscr{G}:=\left\langle \frac{b}{Ra}|\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}|^{2}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D703}J(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D713})\right\rangle ,\end{eqnarray}$$

for $N\geqslant 2$ where

(4.13) $$\begin{eqnarray}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})(x,z)=\mathop{\sum }_{j=1}^{N}(\unicode[STIX]{x1D713}_{j},\unicode[STIX]{x1D703}_{j})(x,z)\end{eqnarray}$$

is taken as the total optimal fluctuation field. In fact, $\mathscr{G}>0$ and so this total optimal field is not a solution of the heat equation. The clear implication is that the spectral constraint is not satisfied for $N\geqslant 2$ at any saddle point of the Lagrangian (2.7) where the steady heat equation is imposed. Consequently, the optimisation procedure is forced to find an optimal away from the saddle points of the Lagrangian (2.7) where the spectral constraint is satisfied to deliver a bound.

The fact that the form of the Lagrangian in (4.11) is different from that in (2.7) warrants further explanation. The Lagrangian in (2.7) is constructed in the usual way from the functional to be maximised (the heat flux) and the constraints to be applied with the highest stationary point being the maximal heat flux under the imposed constraints. The background method adds a further constraint by restricting the set of fields over which the Lagrangian is made stationary. Specifically, the allowed set (2.10) guarantees that any $(\unicode[STIX]{x1D70F},a)\in \unicode[STIX]{x1D6FA}$ will give a value of the Lagrangian greater or equal to that at the highest stationary point: see the discussion in § 2.2. When the background temperature field is of lower dimension than the physical problem – specifically here $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}(z)$ in 2-D Rayleigh–Bénard convection, this extra requirement emerges naturally from the Lagrangian (2.7) since

(4.14) $$\begin{eqnarray}\mathscr{G}=\mathop{\sum }_{j=1}^{N}\mathscr{G}_{j}\end{eqnarray}$$

as the contributions to $\mathscr{G}$ from the various critical wavenumbers naturally separate. Put another way, the Lagrangians (4.11) and (2.7) are one and the same in this situation. There is, however, a further check that all other wavenumbers give $\mathscr{G}>0$ , which gets added onto the process of making the Lagrangian (2.7) stationary. Now, if $\unicode[STIX]{x1D70F}$ has the same dimension as the physical problem, then the spectral problem also has the same dimension as the physical problem and (4.14) no longer holds (unless, trivially, $N=1$ ). In this case the Lagrangian (2.7) has to be modified to that in (4.11) and the optimal solution no longer corresponds with the highest stationary point of (2.7). Instead, the highest stationary point of (4.11) will overestimate the highest stationary point of (2.7) as soon as $N\geqslant 2$ and in particular, the constraints imposed directly in (2.7) will not be satisfied.

From a different perspective, Souza et al. (Reference Souza, Tobasco and Doering2019) have also recently argued that this should happen when exploring the connection between the wall-to-wall approach (a max–min problem) with the associated background method (a min–max problem). A duality gap means that

(4.15) $$\begin{eqnarray}(\text{wall-to-wall})\quad \sup _{\unicode[STIX]{x1D703},\boldsymbol{u}}\,\inf _{\unicode[STIX]{x1D70F}}\mathscr{L}<\inf _{\unicode[STIX]{x1D70F}}\sup _{\unicode[STIX]{x1D703},\boldsymbol{u}}\mathscr{L}\quad (\text{background}),\end{eqnarray}$$

(making the connection $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D70F}-(1-z)$ and $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D703}$ with the variables used by Souza et al. (Reference Souza, Tobasco and Doering2019)) where the optimal solution to the wall-to-wall problem is achieved at a stationary point of $\mathscr{L}$ , thereby implying that to the background method is not. They also supply a simple quadratic polynomial in five variables to illustrate the phenomenon. The calculations described in the next section confirm that this gap starts to exist as soon as $N=2$ .

Figure 5. (a) The Nusselt number $Nu$ versus the Rayleigh number $Ra$ at fixed aspect ratio $L=2\sqrt{2}$ . (b) $Nu$ versus $\unicode[STIX]{x1D6FC}=2\unicode[STIX]{x03C0}/L$ .

5 Extending Hassanzedah et al. with a symmetric 2-D background field $\unicode[STIX]{x1D70F}(x,z)$

To explore multi-modal bounding solutions, a first series of computations was done in the fixed domain $L=2\sqrt{2}$ . In this geometry, the first mode appears at $Ra=27\unicode[STIX]{x03C0}^{4}/4$ (the energy stability threshold), the second mode at $Ra=3,075$ and the third mode at $Ra=24\,650$ . The 1-mode and 2-mode optimal solution branches could be easily continued up to $Ra=10^{5}$ whereas the 3-mode solution branch proved difficult to continue much beyond $Ra>40\,000$ due to numerical issues: see figure 5(a). The 3-mode solution, which provides an upper bound in this geometry over at least the range $24\,650\leqslant Ra\leqslant 40\,000$ , presents only a modest correction to the 2-mode optimal solution which is no longer a bound for these $Ra$ .

A second series of computations were then carried out to investigate the dependence of the $Nu$ -bound on the aspect ratio $L$ . Three different $Ra$ values were chosen to explore the dependence of the bound on $L$ : $Ra=5000$ and 10 000 where the bound is given by a 2-mode solution, and $Ra=25\,000$ where the bound is given by a 3-mode solution. In all three cases, the largest bound is achieved as the aspect ratio $L\rightarrow \infty$ ; see figure 5(b). This is very different from the optimal control results of Hassanzadeh et al. (Reference Hassanzadeh, Chini and Doering2014), where the optimal aspect ratio scales like $Ra^{-1/4}$ and so vanishes as $Ra\rightarrow \infty$ .

Figure 6 shows the structure of the two modes at $Ra=10^{4}$ . The fluctuation fields $\unicode[STIX]{x1D713}_{i}$ and $\unicode[STIX]{x1D703}_{i}$ for both $i=1$ and 2 have a convection roll structure and increasing $L$ just means that more of the rolls fit into the domain. On closer inspection it is clear that the rolls are slightly different near to $x=0$ and $x=\pm \frac{1}{2}L$ where they are forced to have a certain symmetry (symmetry around $x=0$ and periodicity over a length $L$ force symmetry about $x=\pm \frac{1}{2}L$ as well). When the domain is short, e.g. $L=\unicode[STIX]{x03C0}$ , the background field is clearly two-dimensional as seen in figure 7. However, as $L$ increases to $L=8\unicode[STIX]{x03C0}$ , the background field become predominantly one-dimensional away from the imposed lines of symmetry at $x=0$ and $x=\pm \frac{1}{2}L$ (the ends of the domain shown). Plotting the streamfunctions $\unicode[STIX]{x1D713}_{1}$ and $\unicode[STIX]{x1D713}_{2}$ over this long domain – see figure 8 – confirms that the convection cells are similar away from the symmetry lines (‘zone 1’ in figure 8) where $\unicode[STIX]{x1D70F}$ is predominantly one-dimensional but are quite different close to the symmetry lines (‘zone 2’) where $\unicode[STIX]{x1D70F}$ is clearly two-dimensional.

The structure of the optimal fields (both background and fluctuation) and the fact that the bound is maximised as $L\rightarrow \infty$ indicate that the optimal solution is trying to minimise the effect of the imposed symmetry requirements at $x=0$ and $x=\pm \frac{1}{2}L$ . Without this imposed symmetry, the problem becomes translationally invariant and the optimal solution must be one-dimensional by the convexity result in § 2.3 (Doering & Constantin Reference Doering and Constantin1996). There is another simple way to see this. Since the bounding functional $f(\unicode[STIX]{x1D70F},a)$ (see (2.19)) is strictly convex in both $\unicode[STIX]{x1D70F}(x,z)$ and $a$ , any 2-D solution $(\unicode[STIX]{x1D70F}_{2D}(x,z),a)\in \unicode[STIX]{x1D6FA}$ satisfies

(5.1) $$\begin{eqnarray}f\left(\frac{1}{N}\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D70F}_{2D}\left(x+\frac{jL}{N},z\right),a\right)<\frac{1}{N}\mathop{\sum }_{j=1}^{N}f\left(\unicode[STIX]{x1D70F}_{2D}\left(x+\frac{jL}{N},z\right),a\right)\end{eqnarray}$$

by Jensen’s inequality. Taking the limit $N\rightarrow \infty$ on the left-hand side and using translational invariance of the problem on the right-hand side leads to

(5.2) $$\begin{eqnarray}f\left(\unicode[STIX]{x1D70F}_{1\text{D}}(z):=\frac{1}{L}\int _{0}^{L}\unicode[STIX]{x1D70F}_{2D}(x,z)\,\text{d}x,a\right)<f(\unicode[STIX]{x1D70F}_{2D}(x,z),a),\end{eqnarray}$$

so that a 1-D background field always produces a better bound than a 2-D field. The results in figure 7 indicate that this is what the optimal solution to the background problem is trying to achieve.

Figure 6. (a) The first mode ( $\unicode[STIX]{x1D713}_{1},\unicode[STIX]{x1D703}_{1}$ ); (b) the second mode ( $\unicode[STIX]{x1D713}_{2},\unicode[STIX]{x1D703}_{2}$ ) at $Ra=10^{4}$ (only two critical modes are present for this $Ra$ ).

Figure 7. The background field plotted at different aspect ratios at $Ra=10^{4}$ .

Figure 8. The two critical flow fields at $Ra=10^{4}$ and $L=8\unicode[STIX]{x03C0}$ : (a) $\unicode[STIX]{x1D713}_{1}$ and (b) $\unicode[STIX]{x1D713}_{2}$ .

5.1 Lifting the symmetry: 1-D background field

Calculation of the optimal solution assuming from the onset that the background field is one-dimensional simplifies the computation since the fluctuation fields can then be parametrised by their single wavenumber in $x$ (as in (4.5)). In this case, rather than setting a domain periodicity and insisting the fluctuation wavenumbers be consistent with this, the wavenumbers themselves can be optimised over as real continuous variables meaning, in effect, that $L$ is infinite. For example, the Euler–Lagrange equation corresponding to the $m$ th wavenumber $k_{m}$ is

(5.3) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}\mathscr{L}/\unicode[STIX]{x1D6FF}k_{m}:=-2\int _{0}^{1}ak_{m}(u_{m}^{2}+w_{m}^{2})+k_{m}\unicode[STIX]{x1D703}_{m}^{2}\,\text{d}z+\int _{0}^{1}p_{m}u_{m}\,\text{d}z=0.\end{eqnarray}$$

With this formulation, Newton’s method with branch continuation proved much faster than the time stepping approach. It took approximately 4 h CPU time on a 2.8 GHz laptop using Newton’s method to obtain the optimal solution from $Ra=27\unicode[STIX]{x03C0}^{4}/4$ up to $Ra=5\times 10^{8}$ while the time stepping approach took at least a day to generate a single point at $Ra=5\times 10^{8}$ . However, when the domain is fixed, Newton’s method becomes very inefficient as the critical wavenumbers $k_{m}$ are discrete and cannot be tracked using a (continuous) continuation method: in this case, time stepping is the better choice. The numerical solution of the one-dimensional background problem gives the upper bound of $Nu\leqslant 0.055Ra^{1/2}$ , as shown in figure 9(a) with five critical modes present by $Ra=10^{9}$ (see figure 9(b)). This result has the same scaling exponent as the non-slip result $Nu\leqslant 0.026Ra^{1/2}$ of PK03 but with a larger numerical coefficient as should be expected for stress-free boundary conditions. The prior work of Wen et al. (Reference Wen, Chini, Kerswell and Doering2015) indicates that adding a further enstrophy constraint (possible only in stress-free 2-D convection) significantly improves the bound obtained here down to $Nu\leqslant 0.106Ra^{5/12}$ .

Figure 9. (a) The upper bound of $Nu$ versus the Rayleigh number $Ra$ with $Nu\leqslant 0.055Ra^{1/2}$ in the asymptotic regime. (b) The bifurcation diagram of critical wavenumber $k_{m}$ versus the Rayleigh number.

It is worth briefly discussing how the critical wavenumbers which appear in the 1-D and 2-D background field calculations are related. Figure 9(b) indicates that at $Ra=10^{4}$ there is only one critical wavenumber $k_{1}=3.284$ for the 1-D background problem. However, for the 2-D (symmetric) background problem, there are two critical modes, as seen in figure 6 and figure 8, with both having an approximate wavenumber ≈3. 3. Both these modes are forced to be antisymmetric (and so in phase) about $x=0$ and $x=\pm L/2$ (zone 2 in figure 8) but away from these points endeavour to be approximately $\unicode[STIX]{x03C0}/2$ out of phase (zone 1 in figure 8). With this phase difference together with matching amplitudes so

(5.4a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{1}=f(z)\sin (k_{1}x),\quad \unicode[STIX]{x1D703}_{1}=g(z)\cos (k_{1}x), & & \displaystyle\end{eqnarray}$$

(5.5a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{2}=f(z)\sin (k_{1}x+\unicode[STIX]{x03C0}/2),\quad \unicode[STIX]{x1D703}_{2}=g(z)\cos (k_{1}x+\unicode[STIX]{x03C0}/2), & & \displaystyle\end{eqnarray}$$

the nonlinear term in (2.17)

(5.6) $$\begin{eqnarray}\mathop{\sum }_{i=1}^{2}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{i}}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{i}}{\unicode[STIX]{x2202}z}\right]=-k_{1}\left[\frac{\text{d}f}{\text{d}z}g+\frac{\text{d}g}{\text{d}z}f\right]\end{eqnarray}$$

is $x$ independent and therefore can only drive a 1-D background field. From another perspective, the 1-D background field problem really has two modes with $k=3.284$ but only one needs to be tracked as the nonlinear term is horizontally averaged (e.g. see (4.3)) ensuring that the background field stays one-dimensional.

At $Ra=25\,000$ there is even a third mode in the 2-D background problem compared to still only 1 mode in the 1-D background problem (the second wavenumber $k_{2}$ appears at $Ra\approx 26\,450$ ). Figure 10 shows that in fact $\unicode[STIX]{x1D713}_{3}$ is only significant in zone 2 where the imposed symmetries dominate. In zone 1, where the background field is essentially a 1-D profile, $\unicode[STIX]{x1D713}_{3}$ vanishes.

Figure 10. (a–c) The profiles of the three critical flow fields $\unicode[STIX]{x1D713}_{1}$ , $\unicode[STIX]{x1D713}_{2}$ , $\unicode[STIX]{x1D713}_{3}$ at $Ra=25\,000$ , $L=4\unicode[STIX]{x03C0}$ . (d) The profile of the two-dimensional optimal background field.

The conclusion of the computations so far is that using a 2-D background temperature field does not improve the upper bound produced by a 1-D background temperature field, a realisation also reached independently in Doering & Tobasco (Reference Doering and Tobasco2019) (see their lemma 6.1). The next obvious question is whether adding a background velocity field helps either and we now turn our attention to this.

6 Imposing the steady momentum equation: $\unicode[STIX]{x1D753}\neq \mathbf{0}$

In this section, we attempt to improve the bound by using a background temperature field and a background velocity field of the same dimension as the physical problem, which means that the full, albeit steady, momentum and heat equations are imposed as constraints. Importantly, the optimisation problem is no longer convex and so permitting $\unicode[STIX]{x1D753}\neq \mathbf{0}$ could produce a better (reduced) upper bound. However, numerical calculations suggest otherwise with $\unicode[STIX]{x1D753}$ remaining zero after every bifurcation. To attempt to explain this, we use an inductive bifurcation analysis to show that if $\unicode[STIX]{x1D753}=\mathbf{0}$ before a bifurcation then it remains $\mathbf{0}$ after it too, implying that the continuous branch of optimals found by branch tracking out of the energy stability point always has $\unicode[STIX]{x1D753}=\mathbf{0}$ .

The analysis begins by constructing the following Lagrangian:

(6.1) $$\begin{eqnarray}\mathscr{L}=\langle |\unicode[STIX]{x1D735}T|^{2}\rangle -\frac{a}{\unicode[STIX]{x1D70E}Ra}\langle \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\mathcal{N}}}\rangle -2\langle \unicode[STIX]{x1D703}{\mathcal{H}}\rangle ,\end{eqnarray}$$

which, after introducing the extended background decomposition

(6.2a,b ) $$\begin{eqnarray}\boldsymbol{u}=\unicode[STIX]{x1D753}+\boldsymbol{v},\quad T=\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D703}\end{eqnarray}$$

(now both $\unicode[STIX]{x1D753}$ and $\boldsymbol{v}$ vanish on the boundaries $z=0,1$ ), can be rewritten as

(6.3) $$\begin{eqnarray}\mathscr{L}=\frac{\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -a}{1-a}-\frac{a}{1-a}\langle \unicode[STIX]{x1D719}_{z}\unicode[STIX]{x1D70F}\rangle -\frac{1}{1-a}\mathscr{G},\end{eqnarray}$$

( $\unicode[STIX]{x1D719}_{z}:=\unicode[STIX]{x1D753}\boldsymbol{\cdot }\hat{\boldsymbol{z}}$ ), where

$$\begin{eqnarray}\displaystyle \mathscr{G}(\boldsymbol{v},\unicode[STIX]{x1D703}) & := & \displaystyle \langle 2\unicode[STIX]{x1D703}(\unicode[STIX]{x1D753}+\boldsymbol{v})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}\rangle +a\langle \unicode[STIX]{x1D719}_{z}\unicode[STIX]{x1D703}\rangle \nonumber\\ \displaystyle & & \displaystyle +\,\left\langle \frac{a}{\unicode[STIX]{x1D70E}Ra}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+\frac{a}{\unicode[STIX]{x1D70E}Ra}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D753}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+\frac{a}{Ra}|\unicode[STIX]{x1D735}\boldsymbol{v}|^{2}-\frac{a}{Ra}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D753}\right\rangle \nonumber\end{eqnarray}$$

(note $\mathscr{G}$ depends parametrically on $\unicode[STIX]{x1D70F}$ , $\unicode[STIX]{x1D753}$ , $a$ , $\unicode[STIX]{x1D70E}$ and $Ra$ but this is suppressed for clarity). If $\inf _{\boldsymbol{v},\unicode[STIX]{x1D703}}\mathscr{G}$ exists (and necessarily $0<a<1$ ), a bound is then given by

(6.4) $$\begin{eqnarray}Nu\leqslant \frac{\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -a}{1-a}-\frac{a}{1-a}\langle \unicode[STIX]{x1D719}_{z}\unicode[STIX]{x1D70F}\rangle -\frac{1}{1-a}\inf _{\boldsymbol{v},\unicode[STIX]{x1D703}}\mathscr{G}(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D753}).\end{eqnarray}$$

Minimisation of $\mathscr{G}$ with respect to incompressible $\boldsymbol{v}$ and $\unicode[STIX]{x1D703}$ requires

(6.5) $$\begin{eqnarray}\displaystyle & \displaystyle -{\displaystyle \frac{2a}{Ra}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+\frac{a}{\unicode[STIX]{x1D70E}Ra}\boldsymbol{v}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}^{\text{T}})+\frac{a}{\unicode[STIX]{x1D70E}Ra}\unicode[STIX]{x1D753}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}-\frac{a}{\unicode[STIX]{x1D70E}Ra}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D753}+2\unicode[STIX]{x1D703}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D735}p=0,\quad & \displaystyle\end{eqnarray}$$

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

(6.7) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}+a\unicode[STIX]{x1D719}_{z}+2(\boldsymbol{v}+\unicode[STIX]{x1D753})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}=0, & \displaystyle\end{eqnarray}$$

the solution of which is denoted as $(\boldsymbol{v}_{0},\unicode[STIX]{x1D703}_{0})$ . The Lagrangian can then be written as

(6.8) $$\begin{eqnarray}\mathscr{L}=\frac{\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -a}{1-a}-\frac{a}{1-a}\langle \unicode[STIX]{x1D719}_{z}\unicode[STIX]{x1D70F}\rangle -\frac{1}{1-a}\left\{\mathscr{G}(\boldsymbol{v}_{0},\unicode[STIX]{x1D703}_{0};\unicode[STIX]{x1D70F},\unicode[STIX]{x1D753})+\mathop{\sum }_{i=1}^{N}\mathscr{H}(\boldsymbol{v}_{i},\unicode[STIX]{x1D703}_{i})\right\},\end{eqnarray}$$

where $\boldsymbol{v}_{0}$ , $\boldsymbol{v}_{i}$ ( $i=1\ldots N$ ) and $\unicode[STIX]{x1D753}$ are incompressible fields and

(6.9) $$\begin{eqnarray}\mathscr{H}(\boldsymbol{v},\unicode[STIX]{x1D703}):=\left\langle 2\unicode[STIX]{x1D703}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\frac{a}{\unicode[STIX]{x1D70E}Ra}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+\frac{a}{Ra}|\unicode[STIX]{x1D735}\boldsymbol{v}|^{2}\right\rangle\end{eqnarray}$$

is a purely quadratic functional of $(\boldsymbol{v},\unicode[STIX]{x1D703})$ which must be positive semi-definite – the spectral constraint – for $\inf \mathscr{G}$ to exist. The fields $(\boldsymbol{v}_{i},\unicode[STIX]{x1D703}_{i})$ are marginal in that $\mathscr{H}(\boldsymbol{v}_{i},\unicode[STIX]{x1D703}_{i})=0$ and their number $N$ increases with $Ra$ . The aim is to minimise the upper bound over $\unicode[STIX]{x1D70F}$ , $\unicode[STIX]{x1D753}$ and $a$ at fixed $\unicode[STIX]{x1D70E}$ and $Ra$ subject to this spectral constraint. The Euler–Lagrange equations are as follows: the spectral constraint equations for $\boldsymbol{v}_{i}$ and $\unicode[STIX]{x1D703}_{i}$

(6.10) $$\begin{eqnarray}\displaystyle \displaystyle -(1-a)Ra\frac{\unicode[STIX]{x1D6FF}\mathscr{L}}{\unicode[STIX]{x1D6FF}\boldsymbol{v}_{i}} & := & \displaystyle -2a\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}+\frac{a}{\unicode[STIX]{x1D70E}}\boldsymbol{v}_{i}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}^{\text{T}})\nonumber\\ \displaystyle & & \displaystyle +\,2Ra\unicode[STIX]{x1D703}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D735}p_{i}=\mathbf{0},\quad i=1\ldots N,\end{eqnarray}$$

(6.11) $$\begin{eqnarray}\displaystyle \displaystyle -(1-a)\frac{\unicode[STIX]{x1D6FF}\mathscr{L}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}_{i}}:=-2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{i}+2\boldsymbol{v}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}=0,\quad i=1\ldots N; & & \displaystyle\end{eqnarray}$$

the forced field equations for $\boldsymbol{v}_{0}$ and $\unicode[STIX]{x1D703}_{0}$

(6.12) $$\begin{eqnarray}\displaystyle -(1-a)Ra{\displaystyle \frac{\unicode[STIX]{x1D6FF}\mathscr{L}}{\unicode[STIX]{x1D6FF}\boldsymbol{v}_{0}}} & := & \displaystyle -2a\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}+{\displaystyle \frac{a}{\unicode[STIX]{x1D70E}}}\boldsymbol{v}_{0}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}^{\text{T}})\nonumber\\ \displaystyle & & \displaystyle +\,2Ra\unicode[STIX]{x1D703}_{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D735}p_{0}+\frac{a}{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D753}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}-a\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D753}=\mathbf{0};\end{eqnarray}$$

(6.13) $$\begin{eqnarray}\displaystyle \displaystyle -(1-a){\displaystyle \frac{\unicode[STIX]{x1D6FF}\mathscr{L}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}_{0}}}:=-2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{0}+2\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+2\unicode[STIX]{x1D753}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+a\unicode[STIX]{x1D719}_{z}=0; & & \displaystyle\end{eqnarray}$$

the background field equations

(6.14) $$\begin{eqnarray}\displaystyle (1-a)Ra{\displaystyle \frac{\unicode[STIX]{x1D6FF}\mathscr{L}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D753}}} & := & \displaystyle a\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}-aRa(\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D703}_{0})\boldsymbol{e}_{z}-2Ra\unicode[STIX]{x1D703}_{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D735}q+\frac{a}{\unicode[STIX]{x1D70E}}\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{0}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{a}{\unicode[STIX]{x1D70E}}}(\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}^{\text{T}}-\unicode[STIX]{x1D753}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{0})+{\displaystyle \frac{a}{\unicode[STIX]{x1D70E}}}\mathop{\sum }_{i=1}^{N}\boldsymbol{v}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}=\mathbf{0},\end{eqnarray}$$

(6.15) $$\begin{eqnarray}(1-a)\frac{\unicode[STIX]{x1D6FF}\mathscr{L}}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}}:=-2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70F}-a\unicode[STIX]{x1D719}_{z}+2(\boldsymbol{v}_{0}+\unicode[STIX]{x1D753})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{0}+2\mathop{\sum }_{i=1}^{N}\boldsymbol{v}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}=0;\end{eqnarray}$$

and finally the balance parameter equation

(6.16) $$\begin{eqnarray}\displaystyle & & \displaystyle (1-a)\frac{\unicode[STIX]{x1D6FF}\mathscr{L}}{\unicode[STIX]{x1D6FF}a}:=\mathscr{L}-1-\langle \unicode[STIX]{x1D719}_{z}\unicode[STIX]{x1D70F}\rangle -\left\{\mathop{\sum }_{i=1}\left\langle \frac{1}{\unicode[STIX]{x1D70E}Ra}\boldsymbol{v}_{i}\boldsymbol{\cdot }\boldsymbol{v}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+\frac{1}{Ra}|\unicode[STIX]{x1D735}\boldsymbol{v}_{i}|^{2}\right\rangle \right.\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left.\left\langle \unicode[STIX]{x1D719}_{z}\unicode[STIX]{x1D703}_{0}+{\displaystyle \frac{1}{Ra}}|\unicode[STIX]{x1D735}\boldsymbol{v}_{0}|^{2}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70E}Ra}}\boldsymbol{v}_{0}\boldsymbol{\cdot }\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70E}Ra}}\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D753}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}-\frac{1}{Ra}\boldsymbol{v}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D753}\right\rangle \right\}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The pressure-like quantities have been rescaled as follows: $Ra\,p_{i}\rightarrow p_{i}$ , $Ra\,p_{0}\rightarrow p_{0}$ and $(1-a)Ra\,q\rightarrow q$ ) and incompressibility conditions on $\boldsymbol{v}_{0}$ , $\boldsymbol{v}_{i}$ and $\unicode[STIX]{x1D753}$ are left implicit. A key point here is that the forced field pair $(\boldsymbol{v}_{0},\unicode[STIX]{x1D703}_{0})$ is not marginal in the spectral constraint (otherwise equations (6.12) and (6.13) could not be satisfied) and, since $\mathscr{H}$ is positive semi-definite, $\mathscr{H}(\boldsymbol{v}_{0},\unicode[STIX]{x1D703}_{0})>0$ .

6.1 The first bifurcation point

The solution at the first critical point $Ra_{c}=27\unicode[STIX]{x03C0}^{4}/4$ is $\unicode[STIX]{x1D70F}=1-z$ and $\unicode[STIX]{x1D753}=0$ , $(\boldsymbol{v}_{0},\unicode[STIX]{x1D703}_{0})=(\mathbf{0},0)$ and $a=1$ . At $Ra=Ra_{c}$ , the spectral constraint becomes marginal for the first time, i.e. there is a non-trivial solution to the spectral problem

(6.17) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}+2Ra\unicode[STIX]{x1D703}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D735}p_{i}=0, & \displaystyle\end{eqnarray}$$

(6.18) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{i}=0, & \displaystyle\end{eqnarray}$$

(6.19) $$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{i}+\boldsymbol{v}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}=0. & \displaystyle\end{eqnarray}$$

There are two different modes (using symmetries),

(6.20) $$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{v}_{1},\unicode[STIX]{x1D703}_{1})=A_{1}(U(z)\sin (kx)\boldsymbol{e}_{x}+W(z)\cos (kx)\boldsymbol{e}_{z},\unicode[STIX]{x1D6E9}(z)\cos (kx)), & \displaystyle\end{eqnarray}$$

(6.21) $$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{v}_{2},\unicode[STIX]{x1D703}_{2})=A_{2}(U(z)\cos (kx)\boldsymbol{e}_{x}-W(z)\sin (kx)\boldsymbol{e}_{z},-.\unicode[STIX]{x1D6E9}(z)\sin (kx)). & \displaystyle\end{eqnarray}$$

Since $\text{d}\unicode[STIX]{x1D70F}/\text{d}z=-1$ , the structure in $z$ is simple: $U:=\unicode[STIX]{x03C0}\cos (\unicode[STIX]{x03C0}z)$ , $W:=-k\sin (\unicode[STIX]{x03C0}z)$ and $\unicode[STIX]{x1D6E9}:=-\frac{\sqrt{2}}{3\unicode[STIX]{x03C0}}\sin (\unicode[STIX]{x03C0}z)$ where $k=\unicode[STIX]{x03C0}/\sqrt{2}$ . Slightly away from the critical point, $Ra=Ra_{c}+\unicode[STIX]{x1D700}$ , the fields need to be expanded as follows:

$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D700}\unicode[STIX]{x1D70F}_{1}+\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D70F}_{2}+\ldots , & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x1D753}=\unicode[STIX]{x1D700}\unicode[STIX]{x1D753}_{1}+\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D753}_{2}+\unicode[STIX]{x1D700}^{3}\unicode[STIX]{x1D753}_{3}\ldots , & \displaystyle \nonumber\\ \displaystyle & \displaystyle \boldsymbol{v}_{0}=\unicode[STIX]{x1D700}\boldsymbol{v}_{0}^{1}+\unicode[STIX]{x1D700}^{2}\boldsymbol{v}_{0}^{2}+\cdots \,, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x1D703}_{0}=\unicode[STIX]{x1D700}\unicode[STIX]{x1D703}_{0}^{1}+\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D703}_{0}^{2}+\cdots \,, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \boldsymbol{v}_{i}=\unicode[STIX]{x1D700}^{1/2}\boldsymbol{v}_{i}^{0}+\unicode[STIX]{x1D700}^{3/2}\boldsymbol{v}_{i}^{1}+\unicode[STIX]{x1D700}^{5/2}\boldsymbol{v}_{i}^{2}+\cdots \,, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x1D703}_{i}=\unicode[STIX]{x1D700}^{1/2}\unicode[STIX]{x1D703}_{i}^{0}+\unicode[STIX]{x1D700}^{3/2}\unicode[STIX]{x1D703}_{i}^{1}+\unicode[STIX]{x1D700}^{5/2}\unicode[STIX]{x1D703}_{i}^{2}+\cdots \,, & \displaystyle \nonumber\\ \displaystyle & \displaystyle a=a_{0}+\unicode[STIX]{x1D700}a_{1}+\unicode[STIX]{x1D700}^{2}a_{2}+\cdots \,, & \displaystyle \nonumber\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{0}:=1-z$ and $a_{0}:=1$ .

6.1.1 Leading order (e.g. the problem for $v_{i}^{0}$ and $\unicode[STIX]{x1D703}_{i}^{0}$ )

To leading order, the spectral constraint is satisfied by $(\boldsymbol{v}_{i}^{0},\unicode[STIX]{x1D703}_{i}^{0})$ defined in (6.20–6.21) and these fields force the other leading-order equations for the background fields

(6.22) $$\begin{eqnarray}\displaystyle & \displaystyle 2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70F}_{1}+\unicode[STIX]{x1D719}_{1z}=2\mathop{\sum }_{i=1,2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{0}, & \displaystyle\end{eqnarray}$$

(6.23) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}^{1}-Ra_{c}(\unicode[STIX]{x1D70F}_{1}+a_{1}\unicode[STIX]{x1D70F}_{0}-\unicode[STIX]{x1D703}_{0}^{1})\boldsymbol{e}_{z}-\unicode[STIX]{x1D70F}_{0}\boldsymbol{e}_{z}+\unicode[STIX]{x1D735}q=\underbrace{-\frac{1}{\unicode[STIX]{x1D70E}}\mathop{\sum }_{i=1,2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}}_{balanced\,by\,pressure}, & \displaystyle\end{eqnarray}$$

( $\unicode[STIX]{x1D719}_{1z}:=\unicode[STIX]{x1D753}_{1}\boldsymbol{\cdot }\hat{\boldsymbol{z}}$ ), which are coupled with the forced field equations

(6.24) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}^{1}+2Ra_{c}\unicode[STIX]{x1D703}_{0}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D735}p=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D753}_{1}, & \displaystyle\end{eqnarray}$$

(6.25) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{0}^{1}+2\boldsymbol{v}_{0}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}=\unicode[STIX]{x1D719}_{1z}. & \displaystyle\end{eqnarray}$$

The forcing term in (6.22) is

(6.26) $$\begin{eqnarray}\mathop{\sum }_{i=1}^{2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{0}=\frac{\unicode[STIX]{x03C0}(A_{1}^{2}+A_{2}^{2})}{6}\sin (2\unicode[STIX]{x03C0}z),\end{eqnarray}$$

and in (6.23),

(6.27) $$\begin{eqnarray}\mathop{\sum }_{i=1}^{2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}=\frac{k\unicode[STIX]{x03C0}^{2}}{2}(A_{1}^{2}-A_{2}^{2})\sin (2kx)\boldsymbol{e}_{x}+\frac{k^{2}\unicode[STIX]{x03C0}}{2}(A_{1}^{2}+A_{2}^{2})\sin (2\unicode[STIX]{x03C0}z)\boldsymbol{e}_{z}\end{eqnarray}$$

(recall $k=\unicode[STIX]{x03C0}/\sqrt{2}$ ), so that simply

(6.28a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D753}_{1}=\boldsymbol{v}_{0}^{1}=\mathbf{0},\quad \unicode[STIX]{x1D703}_{0}^{1}=0,\quad q=c_{1}(A_{1}^{2}-A_{2}^{2})\cos (2kx)+c_{2}(A_{1}^{2}+A_{2}^{2})\cos (2\unicode[STIX]{x03C0}z),\end{eqnarray}$$

(6.29) $$\begin{eqnarray}\text{and}\quad \unicode[STIX]{x1D70F}_{1}=c_{3}(A_{1}^{2}+A_{2}^{2})\sin (2\unicode[STIX]{x03C0}z),\end{eqnarray}$$

where $c_{1},c_{2}$ and $c_{3}$ are specific constants. Finally, the leading-order balance (which is at $O(\unicode[STIX]{x1D700})$ ) in the balance parameter equation (6.16) is

(6.30) $$\begin{eqnarray}-\frac{1}{a_{1}}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}|^{2}\rangle -\frac{1}{Ra_{c}}\mathop{\sum }_{i=1,2}\langle |\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}|^{2}\rangle =0,\end{eqnarray}$$

which relates $A_{1}^{2}+A_{2}^{2}$ and $a_{1}$ .

6.1.2 Next order (e.g. the problem for $v_{i}^{1}$ and $\unicode[STIX]{x1D703}_{i}^{1}$ )

A further piece of information to identify the leading-order fields comes from a solvability condition on the spectral constraint equations at next order (( $O(\unicode[STIX]{x1D700}^{3/2})$ ) which is

(6.31) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{1}+2Ra_{c}\unicode[STIX]{x1D703}_{i}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D735}p=2a_{1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{0}-2Ra_{c}\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}-2\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}, & \displaystyle\end{eqnarray}$$

(6.32) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{i}^{1}+2\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}=-2\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}. & \displaystyle\end{eqnarray}$$

Formally, this has two solvability conditions,

(6.33) $$\begin{eqnarray}\langle \boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }(2a_{1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{0}-2Ra_{c}\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}-2\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0})-2\unicode[STIX]{x1D703}_{i}^{0}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}\rangle =0,\quad i=1,2,\end{eqnarray}$$

but they are equivalent since $\unicode[STIX]{x1D70F}_{1}$ is one-dimensional (i.e. solely a function of $z$ ) with the resulting condition providing a second equation linking $a_{1}$ and $A_{1}^{2}+A_{2}^{2}$ . These ((6.30) and (6.33)) can then be solved to give

(6.34a,b ) $$\begin{eqnarray}A_{1}^{2}+A_{2}^{2}=24,\quad a_{1}=-{\textstyle \frac{9}{4}}\end{eqnarray}$$

and as a consequence $\unicode[STIX]{x1D70F}_{1}=-\frac{1}{\unicode[STIX]{x03C0}}\sin (2\unicode[STIX]{x03C0}z)$ . The fields $(\boldsymbol{v}_{i}^{1},\unicode[STIX]{x1D703}_{i}^{1})$ depend linearly on $\boldsymbol{v}_{i}^{0}$ , $\unicode[STIX]{x1D703}_{i}^{0}$ and so can be written as

(6.35) $$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{v}_{1}^{1},\unicode[STIX]{x1D703}_{1}^{1})=A_{1}({\mathcal{U}}(z)\sin (kx)\boldsymbol{e}_{x}+{\mathcal{W}}(z)\cos (kx)\boldsymbol{e}_{z},{\mathcal{T}}(z)\cos (kx)), & \displaystyle\end{eqnarray}$$

(6.36) $$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{v}_{2}^{1},\unicode[STIX]{x1D703}_{2}^{1})=A_{2}({\mathcal{U}}(z)\cos (kx)\boldsymbol{e}_{x}-{\mathcal{W}}(z)\sin (kx)\boldsymbol{e}_{z},-{\mathcal{T}}(z)\sin (kx)), & \displaystyle\end{eqnarray}$$

where ${\mathcal{U}}\neq U$ , ${\mathcal{W}}\neq W$ . These fields, along with $(\boldsymbol{v}_{i}^{0},\unicode[STIX]{x1D703}_{i}^{0})$ , drive the higher-order equations governing further corrections to the background fields and the forced fields. These are

(6.37) $$\begin{eqnarray}2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70F}_{2}+\unicode[STIX]{x1D719}_{2z}=\underbrace{2\mathop{\sum }_{i=1,2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{1}+\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{0}}_{driving\,term},\end{eqnarray}$$

(6.38) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}^{2}-Ra_{c}(a_{0}\unicode[STIX]{x1D70F}_{2}+a_{2}\unicode[STIX]{x1D70F}_{0}-\unicode[STIX]{x1D703}_{0}^{2})\boldsymbol{e}_{z}+\unicode[STIX]{x1D735}q\nonumber\\ \displaystyle & & \displaystyle \quad =\underbrace{([a_{1}Ra_{c}+1]\unicode[STIX]{x1D70F}_{1}+a_{1}\unicode[STIX]{x1D70F}_{0})\boldsymbol{e}_{z}-\frac{a_{1}}{\unicode[STIX]{x1D70E}}\mathop{\sum }_{i=1}^{2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}}_{balanced\,by\,pressure}-\underbrace{\frac{1}{\unicode[STIX]{x1D70E}}\mathop{\sum }_{i=1}^{2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{1}+\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}}_{driving\,term}=0,\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(6.39) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}^{2}+2Ra_{c}\unicode[STIX]{x1D703}_{0}^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D735}p_{i}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D753}_{2}, & \displaystyle\end{eqnarray}$$

(6.40) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{0}^{2}+2\boldsymbol{v}_{0}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}=\unicode[STIX]{x1D719}_{2z}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{2z}:=\unicode[STIX]{x1D753}_{2}\boldsymbol{\cdot }\hat{\boldsymbol{z}}$ . The apparent driving term $-a_{1}/\unicode[STIX]{x1D70E}\sum _{i=1}^{2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}$ can be balanced by the pressure term (see (6.27)) as can $([a_{1}Ra_{c}+1]\unicode[STIX]{x1D70F}_{1}+a_{1}\unicode[STIX]{x1D70F}_{0})\boldsymbol{e}_{z}$ . Also, importantly for what follows, $Ra_{c}a_{2}\unicode[STIX]{x1D70F}_{0}\boldsymbol{e}_{z}$ in (6.38) can also be absorbed into the pressure term, which means that $\unicode[STIX]{x1D70F}_{2}$ and $\unicode[STIX]{x1D753}_{2}$ do not depend on $a_{2}$ (this is crucial for the argument surrounding (6.50) below). This leaves the driving term for the 2-D background temperature field

(6.41) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \mathop{\sum }_{i=1,2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{1}+\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{0}=12\left(-kU{\mathcal{T}}-k{\mathcal{U}}\unicode[STIX]{x1D6E9}+W\frac{\text{d}{\mathcal{T}}}{\text{d}z}+{\mathcal{W}}\frac{\text{d}\unicode[STIX]{x1D6E9}}{\text{d}z}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{2}(A_{1}^{2}-A_{2}^{2})\left(kU{\mathcal{T}}+k{\mathcal{U}}\unicode[STIX]{x1D6E9}+W\frac{\text{d}{\mathcal{T}}}{\text{d}z}+{\mathcal{W}}\frac{\text{d}\unicode[STIX]{x1D6E9}}{\text{d}z}\right)\cos (2kx),\end{eqnarray}$$

and the driving term in (6.38) for $\boldsymbol{v}_{0}^{2}$ :

(6.42) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \mathop{\sum }_{i=1}^{2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{1}+\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}=\underbrace{12\left(-kU{\mathcal{W}}-k{\mathcal{U}}W+\frac{\text{d}(W{\mathcal{W}})}{\text{d}z}\right)\boldsymbol{e}_{z}}_{balanced\,by\,pressure}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{2}(A_{1}^{2}-A_{2}^{2})\left(2kU{\mathcal{U}}+W\frac{\text{d}{\mathcal{U}}}{\text{d}z}+{\mathcal{W}}\frac{\text{d}U}{\text{d}z}\right)\sin (2kx)\boldsymbol{e}_{x}.\end{eqnarray}$$

Given the form of these driving terms, $\unicode[STIX]{x1D70F}_{2}$ can be split into two parts: a 1-D part which depends only on $z$ and a 2-D part proportional to $A_{1}^{2}-A_{2}^{2}$ which has both $x$ and $z$ dependences, whereas the remaining corrections $\unicode[STIX]{x1D753}_{2}$ , $\boldsymbol{v}_{0}^{2}$ and $\unicode[STIX]{x1D703}_{0}^{2}$ only have a 2-D part proportional to $A_{1}^{2}-A_{2}^{2}$ , i.e.

(6.43) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{2}=\unicode[STIX]{x1D70F}_{2}^{1D}(z)+\unicode[STIX]{x1D70F}_{2}^{2D}(x,z):=P_{1}(z)+(A_{1}^{2}-A_{2}^{2})P_{2}(z)\cos (2kx), & \displaystyle\end{eqnarray}$$

(6.44) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D753}_{2}=(A_{1}^{2}-A_{2}^{2})[G_{1}(z)\sin (2kx)\boldsymbol{e}_{x}+G_{2}(z)\cos (2kx)\boldsymbol{e}_{z}] & \displaystyle\end{eqnarray}$$

(the expressions for $\boldsymbol{v}_{0}^{2}$ and $\unicode[STIX]{x1D703}_{0}^{2}$ are not needed in what follows and hence suppressed). At this point $\unicode[STIX]{x1D753}_{2}$ is now known as a function of $A_{1}^{2}-A_{2}^{2}$ . Further information about $A_{1}$ and $A_{2}$ comes from solvability conditions at the next order of the spectral constraint ( $O(\unicode[STIX]{x1D700}^{5/2})$ ).

Before pursuing this, we remark that the next order ( $O(\unicode[STIX]{x1D700}^{2})$ ) of the balance (6.16) involves the higher-order unknown $\unicode[STIX]{x1D753}_{3}$ and so at this order $a_{2}$ is unspecified. In fact $a_{2}$ is also set by solvability conditions at $O(\unicode[STIX]{x1D700}^{5/2})$ of the spectral constraint to which we now turn.

6.1.3 Solvability at $O(\unicode[STIX]{x1D700}^{5/2})$ in the problem for $v_{i}^{2}$ and $\unicode[STIX]{x1D703}_{i}^{2}$

The spectral constraint equations at $O(\unicode[STIX]{x1D700}^{5/2})$ are

(6.45) $$\begin{eqnarray}\displaystyle & & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{2}+2Ra_{c}\unicode[STIX]{x1D703}_{i}^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D735}p_{i}=2a_{1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{1}+2a_{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{0}-2Ra_{c}(\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{2}+\unicode[STIX]{x1D703}_{i}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1})\nonumber\\ \displaystyle & & \displaystyle \qquad -\,2(\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}+\unicode[STIX]{x1D703}_{i}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0})-{\displaystyle \frac{1}{\unicode[STIX]{x1D70E}}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{2}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{2}^{T})\boldsymbol{\cdot }\boldsymbol{v}_{i}^{0},\end{eqnarray}$$

(6.46) $$\begin{eqnarray}-2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{i}^{2}+2\boldsymbol{v}_{i}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}=-2\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{2}-2\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}.\end{eqnarray}$$

The operator on the left-hand side of (6.45)–(6.46) is self-adjoint and annihilates the leading-order fields $(\boldsymbol{v}_{j}^{0},\unicode[STIX]{x1D703}_{j}^{0})$ $j=1,2$ . As a result, there are solvability conditions for $(\boldsymbol{v}_{0}^{2},\unicode[STIX]{x1D703}_{0}^{2})$ of the form

(6.47) $$\begin{eqnarray}\displaystyle & & \displaystyle \left\langle \boldsymbol{v}_{j}^{0}\boldsymbol{\cdot }\left[\vphantom{\frac{1}{2}}2a_{1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{1}+2a_{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{0}-2Ra_{c}(\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{2}+\unicode[STIX]{x1D703}_{i}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1})-2(\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}+\unicode[STIX]{x1D703}_{i}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0})\right.\right.\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\left.\left.\frac{1}{\unicode[STIX]{x1D70E}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{2}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{2}^{T})\boldsymbol{\cdot }\boldsymbol{v}_{i}^{0}\right]-2\unicode[STIX]{x1D703}_{j}^{0}(\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{2}+\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1})\right\rangle =0.\end{eqnarray}$$

Taking $i=j$ (the $i\neq j$ conditions vanish trivially), this can be rearranged to

(6.48) $$\begin{eqnarray}\text{Term}_{1}(i)+(A_{1}^{2}-A_{2}^{2})\text{Term}_{2}(i)=0,\quad i=1,2,\end{eqnarray}$$

where

(6.49) $$\begin{eqnarray}(A_{1}^{2}-A_{2}^{2})\text{Term}_{2}(i):=-\frac{1}{A_{i}^{2}}\left\langle \frac{2}{\unicode[STIX]{x1D70E}}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{2}\boldsymbol{\cdot }\boldsymbol{v}_{i}^{0}+2(Ra_{c}+1)\unicode[STIX]{x1D703}_{i}^{0}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{2}^{2D}\right\rangle .\end{eqnarray}$$

Crucially, $\text{Term}_{1}(1)=\text{Term}_{1}(2)$ whereas $\text{Term}_{2}(1)=-\text{Term}_{2}(2)$ so that (6.48) implies that

(6.50) $$\begin{eqnarray}\text{Term}_{1}(1)=(A_{1}^{2}-A_{2}^{2})\text{Term}_{2}(1)=0.\end{eqnarray}$$

The unspecified coefficient $a_{2}$ only enters $\text{Term}_{1}$ and so is set by the condition this vanishes. In contrast, there are no free constants in $\text{Term}_{2}$ which is non-zero in our computations (although we have been unable to prove this is always the case). In this situation, $A_{1}^{2}-A_{2}^{2}=0$ is forced instead, which eliminates at a stroke all two-dimensional fields in the bifurcation analysis. Consequently, the background flow field remains zero and the background temperature field stays one-dimensional after the first bifurcation.

6.2 Subsequent bifurcations

Now, we consider subsequent bifurcations to establish that if $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}(z),\unicode[STIX]{x1D753}=0$ exists before then that situation persists after the bifurcation. The approach is inductive: assume $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D70F}(z),\unicode[STIX]{x1D753}=0$ after $m$ bifurcations and consider the $(m+1)$ th bifurcation at $Ra=Ra_{c}^{(m+1)}$ where two new neutral modes appear so that there are now $2(m+1)$ critical modes in the spectral constraint. Defining $\unicode[STIX]{x1D700}:=Ra-Ra_{c}^{(m+1)}$ we expand

(6.51) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D70F},\unicode[STIX]{x1D753})=(\unicode[STIX]{x1D70F}_{0}(z)+\unicode[STIX]{x1D700}\unicode[STIX]{x1D70F}_{1}(x,z)+\cdots \,,\unicode[STIX]{x1D700}\unicode[STIX]{x1D753}_{0}(x,z)+\cdots \,), & \displaystyle\end{eqnarray}$$

(6.52) $$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{v}_{0},\unicode[STIX]{x1D703}_{0})=(\unicode[STIX]{x1D700}\boldsymbol{v}_{0}^{0}+\cdots \,,\unicode[STIX]{x1D700}\unicode[STIX]{x1D703}_{0}^{0}+\cdots \,), & \displaystyle\end{eqnarray}$$

(6.53) $$\begin{eqnarray}\displaystyle & \displaystyle a=a_{0}+\unicode[STIX]{x1D700}a_{1}+\cdots \,, & \displaystyle\end{eqnarray}$$

(6.54) $$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{v}_{i},\unicode[STIX]{x1D703}_{i})=\left\{\begin{array}{@{}lll@{}}(\!\boldsymbol{v}_{i}^{0}+\unicode[STIX]{x1D700}\boldsymbol{v}_{i}^{1}+\cdots \,, & \unicode[STIX]{x1D703}_{i}^{0}+\unicode[STIX]{x1D700}\unicode[STIX]{x1D703}_{i}^{1}+\cdots \!), & \quad i=1,2,\ldots ,2m,\\ (\!\unicode[STIX]{x1D700}^{1/2}\boldsymbol{v}_{i}^{0}+\unicode[STIX]{x1D700}^{3/2}\boldsymbol{v}_{i}^{1}+\cdots \,, & \unicode[STIX]{x1D700}^{1/2}\unicode[STIX]{x1D703}_{i}^{0}+\unicode[STIX]{x1D700}^{3/2}\unicode[STIX]{x1D703}_{i}^{1}+\cdots \!), & i=2m+1,2m+2,\end{array}\right. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the leading fields $\unicode[STIX]{x1D70F}_{0}(z)$ , $(\boldsymbol{v}_{i}^{0},\unicode[STIX]{x1D703}_{i}^{0})$ ( $i=1,\ldots ,2m+2$ ) and $a_{0}$ are all known. In particular, the $i$ th wavenumber $k_{i}$ ( $i=1,2,\ldots ,m$ ), is associated with two modes which, to leading order, are

(6.55) $$\begin{eqnarray}(\boldsymbol{v}_{2i-1}^{0},\unicode[STIX]{x1D703}_{2i-1}^{0})=A_{i}(U_{i}(z)\sin (k_{i}x)\boldsymbol{e}_{x}+W_{i}(z)\cos (k_{i}x)\boldsymbol{e}_{z},\unicode[STIX]{x1D6E9}_{i}(z)\cos (k_{i}x))\end{eqnarray}$$

and

(6.56) $$\begin{eqnarray}(\boldsymbol{v}_{2i}^{0},\unicode[STIX]{x1D703}_{2i}^{0})=B_{i}(U_{i}(z)\cos (k_{i}x)\boldsymbol{e}_{x}-W_{i}(z)\sin (k_{i}x)\boldsymbol{e}_{z},-\unicode[STIX]{x1D6E9}_{i}(z)\sin (k_{i}x)),\end{eqnarray}$$

where $A_{i}^{2}=B_{i}^{2}$ for $i=1,2,\ldots ,m$ . The two new modes emerging at the $(m+1)$ th bifurcation point can be assumed to have the following general 3-D forms:

(6.57) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}_{2m+1}^{0}=A_{m+1}\left(\begin{array}{@{}c@{}}U_{m+1}(z)\sin (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x})\\ V_{m+1}(z)\sin (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x})\\ W_{m+1}(z)\cos (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x})\end{array}\right), & \displaystyle\end{eqnarray}$$

(6.58) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{2m+1}^{0}=A_{m+1}\unicode[STIX]{x1D6E9}_{m+1}(z)\cos (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x}) & \displaystyle\end{eqnarray}$$

and

(6.59) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}_{2m+2}^{0}=B_{m+1}\left(\begin{array}{@{}c@{}}U_{m+1}(z)\cos (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x})\\ V_{m+1}(z)\cos (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x})\\ -W_{m+1}(z)\sin (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x})\end{array}\right), & \displaystyle\end{eqnarray}$$

(6.60) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{2m+2}^{0}=-B_{m+1}\unicode[STIX]{x1D6E9}_{m+1}(z)\sin (\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x}), & \displaystyle\end{eqnarray}$$

where $\boldsymbol{k}_{m+1}=(k_{x},k_{y},0)$ . The objective in what follows is to show that $A_{m+1}^{2}=B_{m+1}^{2}$ after the bifurcation so that the optimisation problem remains one-dimensional.

At leading order ( $O(\unicode[STIX]{x1D716})$ ) in the forced field and background field equations

(6.61) $$\begin{eqnarray}\displaystyle & \displaystyle -2a_{0}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}^{0}+2Ra_{c}^{m+1}\unicode[STIX]{x1D703}_{0}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D735}p-a_{0}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D753}_{0}=\mathbf{0}, & \displaystyle\end{eqnarray}$$

(6.62) $$\begin{eqnarray}\displaystyle & \displaystyle -2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{0}^{0}+2\boldsymbol{v}_{0}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+2\unicode[STIX]{x1D753}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+a_{0}\unicode[STIX]{x1D719}_{0z}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{0z}:=\unicode[STIX]{x1D753}_{\mathbf{0}}\boldsymbol{\cdot }\hat{\boldsymbol{z}}$ ,

(6.63) $$\begin{eqnarray}\displaystyle & & \displaystyle a_{0}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{0}^{0}-Ra_{c}^{m+1}(a_{0}\unicode[STIX]{x1D70F}_{1}+a_{1}\unicode[STIX]{x1D70F}_{0}+a_{0}\unicode[STIX]{x1D703}_{0}^{0})\boldsymbol{e}_{z}-a_{0}\unicode[STIX]{x1D70F}_{0}\boldsymbol{e}_{z}-2Ra_{c}^{m+1}\unicode[STIX]{x1D703}_{0}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D735}q\nonumber\\ \displaystyle & & \displaystyle +\,\frac{a_{0}}{\unicode[STIX]{x1D70E}}\mathop{\sum }_{i=1}^{2m}\biggl(\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{1}+\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}\biggr)=-\underbrace{\frac{a_{1}}{\unicode[STIX]{x1D70E}}\mathop{\sum }_{i=1}^{2m}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}}_{one\text{-}dimensional}-\underbrace{\frac{a_{0}}{\unicode[STIX]{x1D70E}}\mathop{\sum }_{i=2m+1}^{2m+2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}}_{driving\,term},\quad\end{eqnarray}$$

(6.64) $$\begin{eqnarray}\displaystyle -2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70F}_{1}-a_{0}\unicode[STIX]{x1D719}_{0z}+2\displaystyle \mathop{\sum }_{i=1}^{2m}\biggl(\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{1}+\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{0}\biggr)=-\underbrace{2\mathop{\sum }_{i=2m+1}^{2m+2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{0}}_{driving\,term} & & \displaystyle\end{eqnarray}$$

and again it is implicit that $\boldsymbol{v}_{0}$ , $\boldsymbol{v}_{i}$ and $\unicode[STIX]{x1D753}$ are incompressible fields. The spectral constraint for modes $i=1,2,\ldots ,2m$ at $O(\unicode[STIX]{x1D700})$ and for modes $i=2m+1,2m+2$ at $O(\unicode[STIX]{x1D700}^{3/2})$ is

(6.65) $$\begin{eqnarray}\displaystyle -2a_{0}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{1}+2Ra_{c}^{m+1}\unicode[STIX]{x1D703}_{i}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}+\unicode[STIX]{x1D735}p & = & \displaystyle 2a_{1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{0}-{\displaystyle \frac{a_{0}}{\unicode[STIX]{x1D70E}}}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{0}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{0}^{\text{T}})\nonumber\\ \displaystyle & & \displaystyle -\,2Ra_{c}^{m+1}\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}-2\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0},\end{eqnarray}$$

(6.66) $$\begin{eqnarray}-2\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}_{i}^{1}+2\boldsymbol{v}_{i}^{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}=-2\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}.\end{eqnarray}$$

The system of (6.61)–(6.66) is linear in $\boldsymbol{v}_{0}^{0}$ , $\unicode[STIX]{x1D703}_{0}^{0}$ , $\unicode[STIX]{x1D753}_{0}$ , $\boldsymbol{v}_{i}^{1}$ and $\unicode[STIX]{x1D703}_{i}^{1}$ ( $i=1,2,\ldots ,2m+2$ ). The emergent critical modes at $Ra_{c}^{m+1}$ give rise to the new driving terms in (6.63) and (6.64)

(6.67) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{i=2m+1}^{2m+2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}_{i}^{0}=\frac{1}{2}(A_{m+1}^{2}+B_{m+1}^{2})\unicode[STIX]{x1D735}[W_{m+1}^{2}(z)]\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{2}(A_{m+1}^{2}-B_{m+1}^{2})\left[\begin{array}{@{}c@{}}-U_{m+1}{\displaystyle \frac{\text{d}W_{m+1}}{\text{d}z}}+{\displaystyle \frac{\text{d}U_{m+1}}{\text{d}z}}W_{m+1}\\[12.0pt] -V_{m+1}{\displaystyle \frac{\text{d}W_{m+1}}{\text{d}z}}+{\displaystyle \frac{\text{d}V_{m+1}}{\text{d}z}}W_{m+1}\\ 0\end{array}\right]\sin (2\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x}),\end{eqnarray}$$

(6.68) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{i=2m+1}^{2m+2}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{i}^{0}=\frac{1}{2}(A_{m+1}^{2}+B_{m+1}^{2})\frac{\text{d}(W_{m+1}\unicode[STIX]{x1D6E9}_{m+1})}{\text{d}z}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{2}(A_{m+1}^{2}-B_{m+1}^{2})\left(W_{m+1}\frac{\text{d}\unicode[STIX]{x1D6E9}_{m+1}}{\text{d}z}-\frac{\text{d}W_{m+1}}{\text{d}z}\unicode[STIX]{x1D6E9}_{m+1}\right)\cos (2\boldsymbol{k}_{m+1}\boldsymbol{\cdot }\boldsymbol{x}),\end{eqnarray}$$

which have a 1-D part proportional to $A_{m+1}^{2}+B_{m+1}^{2}$ and a non-1-D part proportional to $A_{m+1}^{2}-B_{m+1}^{2}$ . If $A_{m+1}^{2}=B_{m+1}^{2}$ then $\unicode[STIX]{x1D70F}_{1}=\unicode[STIX]{x1D70F}_{1}(z)$ and $\unicode[STIX]{x1D753}=\mathbf{0}$ using the arguments of § 6.1. Hence, for $A_{m+1}^{2}\neq B_{m+1}^{2}$ and using (6.64), we can assume a solution structure of the form

(6.69)

(6.70)

where $\unicode[STIX]{x1D70F}_{1}^{\ast }(x,z)$ and $\unicode[STIX]{x1D731}^{\ast }(x,z)$ collect all the other wavenumber dependences on $x$ in $\unicode[STIX]{x1D70F}_{1}$ and $\unicode[STIX]{x1D753}_{0}$ respectively (this is unimportant in what follows). The key is then examining the solvability conditions

(6.71) $$\begin{eqnarray}\displaystyle \langle \!\boldsymbol{v}_{j}^{0}\boldsymbol{\cdot }\!\left(2a_{1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}_{i}^{0}-\frac{a_{0}}{\unicode[STIX]{x1D70E}}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{0}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{0}^{\text{T}})-2Ra_{c}^{m+1}\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}-2\unicode[STIX]{x1D703}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{0}\!\right)-2\unicode[STIX]{x1D703}_{i}^{0}\boldsymbol{v}_{j}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}\!\rangle =0 & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

( $i,j,=1,\ldots m+1$ ) on the corrections $\boldsymbol{v}_{i}^{1}$ to all the critical modes $\boldsymbol{v}_{i}^{0}$ of the spectral constraint. These set the amplitudes $A_{i}^{2}\,(=B_{i}^{2})$ ( $i=1,\ldots m$ ) and $(A_{m+1},B_{m+1})$ ( $a_{1}$ is determined by the balance parameter equation at $O(\unicode[STIX]{x1D700})$ ).

To establish that the background fields stay one-dimensional, it is sufficient to focus on the solvability conditions for the new critical modes ( $i=2m+1$ and $2m+2$ ). Here, the solvability condition is explicitly

(6.72) $$\begin{eqnarray}\displaystyle & & \displaystyle -a_{1}\int _{0}^{1}\boldsymbol{k}_{m+1}^{2}(U_{m+1}^{2}+V_{m+1}^{2}+W_{m+1}^{2})+\left({\displaystyle \frac{\text{d}U_{m+1}}{\text{d}z}}\right)^{2}+\left(\frac{\text{d}V_{m+1}}{\text{d}z}\right)^{2}+\left({\displaystyle \frac{\text{d}W_{m+1}}{\text{d}z}}\right)^{2}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\int _{0}^{1}W_{m+1}\unicode[STIX]{x1D6E9}_{m+1}{\displaystyle \frac{\text{d}\unicode[STIX]{x1D70F}_{0}}{\text{d}z}}\,\text{d}z+\mathop{\sum }_{j=1}^{m+1}\left[(A_{j}^{2}+B_{j}^{2})(Ra_{c}^{m+1}+1)\int _{0}^{1}W_{j}\unicode[STIX]{x1D6E9}_{j}{\displaystyle \frac{\text{d}P_{j}}{\text{d}z}}\,\text{d}z\right]\nonumber\\ \displaystyle & & \displaystyle \qquad +\,(A_{m+1}^{2}-B_{m+1}^{2})\text{Term}(i)=0,\end{eqnarray}$$

where

(6.73) $$\begin{eqnarray}\displaystyle (A_{m+1}^{2}-B_{m+1}^{2})\text{Term}(i):=-\frac{2}{c}\left\langle \frac{a_{0}}{\unicode[STIX]{x1D70E}}\boldsymbol{v}_{i}^{0}\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{0}\boldsymbol{\cdot }\boldsymbol{v}_{i}^{0}+(Ra_{c}^{m+1}+1)\unicode[STIX]{x1D703}_{i}^{0}\boldsymbol{v}_{i}^{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70F}_{1}^{2D}\right\rangle , & & \displaystyle\end{eqnarray}$$

with $c:=A_{m+1}^{2}$ for $i=2m+1$ or $B_{m+1}^{2}$ for $i=2m+2$ . Crucially, it is straightforward to show that

(6.74) $$\begin{eqnarray}\text{Term}(2m+1)=-\text{Term}(2m+2),\end{eqnarray}$$

so that, as in § 6.1, we must have

(6.75) $$\begin{eqnarray}\displaystyle (A_{m+1}^{2}-B_{m+1}^{2})\text{Term}(i)=0 & & \displaystyle\end{eqnarray}$$

leaving equation (6.72) to determine the value of $A_{m+1}^{2}+B_{m+1}^{2}$ . Our computations indicate $\text{Term}(2m+1)\neq 0$ (although, as in § 6.1, we have been unable to prove this in general) implying that $A_{m+1}^{2}=B_{m+1}^{2}$ . This forces $\unicode[STIX]{x1D70F}_{1}=\unicode[STIX]{x1D70F}_{1}(z)$ and $\unicode[STIX]{x1D753}_{0}=\mathbf{0}$ so that the optimal solution stays one-dimensional after the $(m+1)$ th ( $m\geqslant 0$ ) bifurcation if it is of this form before.

Taken together with the first bifurcation analysis in § 6.1, the result of this section means that the optimal solution is one-dimensional for all $Ra$ and so, surprisingly, there is no benefit in imposing the full steady momentum and heat balances in the upper bound problem.

7 Adding a background velocity field to Wen et al. (Reference Wen, Chini, Kerswell and Doering2015)

Following the success of adding an enstrophy constraint in 2-D stress-free convection (Wen et al. Reference Wen, Chini, Kerswell and Doering2015), an interesting question is whether adding a 1-D background velocity field by using the decomposition

(7.1a,b ) $$\begin{eqnarray}\boldsymbol{u}(x,z)=\unicode[STIX]{x1D719}(z)\boldsymbol{e}_{x}+\boldsymbol{v}(x,z),\quad T(x,z)=\unicode[STIX]{x1D70F}(z)+\unicode[STIX]{x1D703}(x,z)\end{eqnarray}$$

would improve the bound further, since this imposes additional information from the Navier–Stokes equations (this question is not covered by the conclusion in § 6 since the extra enstrophy constraint was not included there). To maximise the heat flux, the Lagrangian

(7.2) $$\begin{eqnarray}Nu=\langle |\unicode[STIX]{x1D735}T|^{2}\rangle -\frac{a}{\unicode[STIX]{x1D70E}Ra}\langle \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\mathcal{N}}}\rangle -\frac{b}{\unicode[STIX]{x1D70E}Ra}\langle \unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{{\mathcal{N}}}\rangle -2\langle \unicode[STIX]{x1D703}\,{\mathcal{H}}\rangle\end{eqnarray}$$

is constructed where $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D714}(x,z)\boldsymbol{e}_{y}:=(\unicode[STIX]{x1D735}\times \boldsymbol{v})$ . After some integration by parts, judicious use of boundary conditions and building in the fact that $Nu=1+\langle wT\rangle$ , this leads to the expression

(7.3) $$\begin{eqnarray}\displaystyle \displaystyle (1-a)Nu+a & = & \displaystyle \langle |\unicode[STIX]{x1D70F}^{\prime }|^{2}\rangle -\left\langle |\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D703}v_{z}\unicode[STIX]{x1D70F}^{\prime }+\frac{a}{\unicode[STIX]{x1D70E}Ra}v_{x}v_{z}\unicode[STIX]{x1D719}^{\prime }+\frac{b}{\unicode[STIX]{x1D70E}Ra}v_{z}\unicode[STIX]{x1D714}\unicode[STIX]{x1D719}^{\prime \prime }-b\unicode[STIX]{x1D714}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.{\displaystyle \frac{a}{Ra}}(|\unicode[STIX]{x1D6FB}\boldsymbol{v}|^{2}-v_{x}\unicode[STIX]{x1D719}^{\prime \prime })+{\displaystyle \frac{b}{Ra}}(|\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D714}|^{2}-\unicode[STIX]{x1D714}\unicode[STIX]{x1D719}^{\prime \prime \prime })\right\rangle ,\end{eqnarray}$$

where $\boldsymbol{v}=v_{x}\,\hat{\boldsymbol{x}}+v_{z}\hat{\boldsymbol{z}}$ . The two linear terms in the second line of this expression – $v_{x}\unicode[STIX]{x1D719}^{\prime \prime }$ and $\unicode[STIX]{x1D714}\unicode[STIX]{x1D719}^{\prime \prime \prime }$ – mean that optimisation over the fluctuation fields $\boldsymbol{v}$ and $\unicode[STIX]{x1D703}$ will give rise to a non-zero contribution to be added to $\langle |\unicode[STIX]{x1D70F}^{\prime }|^{2}\rangle$ . This complication can be avoided (or rather made more explicit) by defining a shifted variable

(7.4) $$\begin{eqnarray}\hat{\boldsymbol{v}}:=\boldsymbol{v}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D719}(z)\boldsymbol{e}_{x},\end{eqnarray}$$

which is possible if $\boldsymbol{v}$ and $\unicode[STIX]{x1D719}$ are both assumed to satisfy (natural) homogeneous boundary conditions and allows both linear terms to be absorbed into perfect squares. As a result of this, the expression becomes

(7.5) $$\begin{eqnarray}\displaystyle Nu & = & \displaystyle {\displaystyle \frac{1}{1-a}}(\langle |\unicode[STIX]{x1D70F}^{\prime }|^{2}\rangle -a)+{\displaystyle \frac{a}{4(1-a)Ra}}\langle |\unicode[STIX]{x1D719}^{\prime }|^{2}\rangle +{\displaystyle \frac{b}{4(1-a)Ra}}\langle |\unicode[STIX]{x1D719}^{\prime \prime }|^{2}\rangle \nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{1}{1-a}}\mathscr{G}(\hat{\boldsymbol{v}},\unicode[STIX]{x1D714},\unicode[STIX]{x1D703};\unicode[STIX]{x1D70F},\unicode[STIX]{x1D719},a,b,Ra,\unicode[STIX]{x1D70E}),\end{eqnarray}$$

where

$$\begin{eqnarray}\displaystyle \mathscr{G} & := & \displaystyle a\left\langle \frac{1}{Ra}|\unicode[STIX]{x1D6FB}\hat{\boldsymbol{v}}|^{2}+\frac{1}{\unicode[STIX]{x1D70E}Ra}\hat{v}_{x}\hat{v}_{z}\unicode[STIX]{x1D719}^{\prime }\right\rangle +b\left\langle \frac{1}{\unicode[STIX]{x1D70E}Ra}\unicode[STIX]{x1D714}\hat{v}_{z}\unicode[STIX]{x1D719}^{\prime \prime }+\frac{1}{Ra}|\unicode[STIX]{x1D735}\hat{\unicode[STIX]{x1D714}}|^{2}-\hat{\unicode[STIX]{x1D714}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}\right\rangle \nonumber\\ \displaystyle & & \displaystyle +\,\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D703}\hat{v}_{z}\unicode[STIX]{x1D70F}^{\prime }\rangle \nonumber\end{eqnarray}$$

is a purely quadratic functional of $\hat{\boldsymbol{v}}$ , $\hat{\unicode[STIX]{x1D714}}:=\boldsymbol{e}_{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \hat{\boldsymbol{v}}=\unicode[STIX]{x1D714}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D719}^{\prime }$ and $\unicode[STIX]{x1D703}$ . Provided the background fields $\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D719}$ are chosen such that $\mathscr{G}\geqslant 0$ for all permissible $\hat{\boldsymbol{v}}$ , $\hat{\unicode[STIX]{x1D714}}$ and $\unicode[STIX]{x1D703}$ , then a bound follows on $Nu$ . Now, it is clear that: (i) the objective functional is strictly convex in the background fields and (ii) the set of allowable background fields is convex (if $(\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D719}_{1})$ and $(\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D719}_{2})$ ensure $\mathscr{G}\geqslant 0$ so does $(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D719})=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D70F}_{1},\unicode[STIX]{x1D719}_{1})+(1-\unicode[STIX]{x1D707})(\unicode[STIX]{x1D70F}_{2},\unicode[STIX]{x1D719}_{2})$ with $0\leqslant \unicode[STIX]{x1D707}\leqslant 1$ ). This implies that the optimiser is unique and is attained for $(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D719})=(\unicode[STIX]{x1D70F},0)$ i.e. the background velocity field vanishes, indicating that the extra information this folds into the optimisation is, in fact, unimportant. Physically, a bifurcation analysis shows that the fluctuation fields are always such as to produce zero Reynolds stress so that no background flow field is generated.

A numerical solution shown in figure 11 using Newton’s method on the Euler–Lagrange equations in an infinitely long domain confirms that $\unicode[STIX]{x1D719}=0$ , as does a bifurcation analysis developed in the same way as the previous section (not shown). The bound compares well with the earlier results of Wen et al. (Reference Wen, Chini, Kerswell and Doering2015) who considered a fixed domain of $L=2\sqrt{2}$ , indicating further that the bound is not that sensitive to the domain size. The fashion in which the necessarily discretised critical modes found by Wen et al. (Reference Wen, Chini, Kerswell and Doering2015) cluster around the (continuous) optimal wavenumbers in our study confirms this conclusion.

Figure 11. (a) The upper bound of $Nu$ versus the Rayleigh number $Ra$ , and $Nu\leqslant 0.107a^{5/12}$ in the asymptotic ultimate regime. (b) The bifurcation diagram of critical wavenumbers $k_{m}$ versus the Rayleigh number (solid lines found here for $L=\infty$ ; dots from Wen et al. (Reference Wen, Chini, Kerswell and Doering2015) for fixed $L=2\sqrt{2}$ ). The prefactor predicted 0.107 is very slightly higher than Wen et al.’s (Reference Wen, Chini, Kerswell and Doering2015) 0.106 as the domain size $L$ is also optimised over here (data courtesy of B. Wen).