Bounds on heat transfer for Bénard–Marangoni convection at infinite Prandtl number

Abstract

The vertical heat transfer in Bénard–Marangoni convection of a fluid layer with infinite Prandtl number is studied by means of upper bounds on the Nusselt number $Nu$ as a function of the Marangoni number $Ma$. Using the background method for the temperature field, it has recently been proved by Hagstrom & Doering (Phys. Rev. E, vol. 81, 2010, art. 047301) that $Nu\leqslant 0.838Ma^{2/7}$. In this work we extend previous background method analysis to include balance parameters and derive a variational principle for the bound on $Nu$, expressed in terms of a scaled background field, that yields a better bound than Hagstrom & Doering’s formulation at a given $Ma$. Using a piecewise-linear, monotonically decreasing profile we then show that $Nu\leqslant 0.803Ma^{2/7}$, lowering the previous prefactor by 4.2 %. However, we also demonstrate that optimisation of the balance parameters does not affect the asymptotic scaling of the optimal bound achievable with Hagstrom & Doering’s original formulation. We subsequently utilise convex optimisation to optimise the bound on $Nu$ over all admissible background fields, as well as over two smaller families of profiles constrained by monotonicity and convexity. The results show that $Nu\leqslant O(Ma^{2/7}(\ln Ma)^{-1/2})$ when the background field has a non-monotonic boundary layer near the surface, while a power-law bound with exponent $2/7$ is optimal within the class of monotonic background fields. Further analysis of our upper-bounding principle reveals the role of non-monotonicity, and how it may be exploited in a rigorous mathematical argument.

1 Introduction

When the surface of a layer of fluid experiences sufficiently strong local variations in temperature, surface-tension-induced shear stresses drive bulk convective motion. Bénard–Marangoni convection, as it is commonly known, arises in a variety of industrial processes, including drying of thin polymer films (Yiantsios et al. 2015), fusion welding (DebRoy & David 1995), laser cladding (Kumar & Roy 2009) and the growth of single-crystal semiconductors (Lappa 2010, chapter 3 and references therein). Shear-driven convection is also observed in distillation columns (Zuiderweg & Harmens 1958; Patberg et al. 1983) and in differentially heated fluids in microgravity environments, where buoyancy effects are negligible (Lappa 2010, chapter 2).

Despite its widespread applications, the dynamics and heat transfer properties of Bénard–Marangoni convection have been studied far less than those of buoyancy-driven Rayleigh–Bénard convection. One fundamental question that remains largely unanswered is how the net vertical heat transfer across the layer, described by the Nusselt number $Nu$ , depends on the external forcing, measured by the Marangoni number $Ma$ . A phenomenological boundary layer scaling analysis put forward by Pumir & Blumenfeld (1996) predicts a transition from $Nu=O(Ma^{1/4})$ to $Nu=O(Ma^{1/3})$ as laminar convection rolls are replaced by turbulent convection, with prefactors that depend on the Prandtl number $Pr$ – the ratio of the fluid’s kinematic viscosity and its thermal diffusivity. Two-dimensional direct numerical simulations (DNSs) at low $Pr$ and large $Ma$  (Boeck & Thess 1998; Boeck 2005) confirm the $1/3$ scaling exponent for the turbulent regime when free-slip conditions are imposed on the velocity field, but $Nu=O(Ma^{1/5})$ is observed in the no-slip case. Moreover, further DNSs by Boeck & Thess (2001) indicate that Bénard–Marangoni convection in high Prandtl number fluids may not be turbulent even when $Ma$ is $10^{4}$ times the value at which convection first appears. Under the assumption that the observed stationary convection rolls remain stable as $Ma$ is raised when $Pr$ is infinite, the same authors predict that $Nu=O(Ma^{2/9})$ in this limit.

Unfortunately, available experimental data (see Schatz & Neitzel 2001; Eckert & Thess 2006, and references therein) do not reach the highly nonlinear regime, where these scaling laws are thought to apply. An alternative approach to confirm or disprove them is to try and derive rigorous bounds on $Nu$ as a function of $Ma$ directly from the governing equations. This can be done without recourse to statistical hypotheses or closure models using the background method (Doering & Constantin 1992, 1994, 1996; Constantin & Doering 1995a , b ). The essence of the method is to write the temperature field as the sum of a steady ‘background’ component $\unicode[STIX]{x1D70F}$ and a time-dependent fluctuation, and show that if $\unicode[STIX]{x1D70F}$ satisfies a particular nonlinear stability condition, then $Nu$ is bounded as a function of $\unicode[STIX]{x1D70F}$ only. The problem that results is variational in nature: optimise the bound on $Nu$ over all stable background fields.

The background method has been applied extensively to the Rayleigh–Bénard problem in a variety of configurations (see e.g. Doering & Constantin 1996; Otero 2002; Doering, Otto & Reznikoff 2006; Wittenberg & Gao 2010; Whitehead & Doering 2011, 2012; Goluskin & Doering 2016). On the other hand, the only result for Bénard–Marangoni convection is due to Hagstrom & Doering (2010), who used a monotonically decreasing, piecewise-linear background temperature field to prove $Nu\leqslant 0.841\times Ma^{1/2}$ for finite Prandtl number fluids, while $Nu\leqslant 0.838\times Ma^{2/7}$ in the infinite- $Pr$ limit.

This work investigates whether Hagstrom & Doering’s bound for Bénard–Marangoni convection at infinite Prandtl number can be lowered, reducing the gap with the DNS results and phenomenological predictions of Boeck & Thess (2001). The assumption of infinite $Pr$ significantly simplifies the mathematical treatment of the problem, making it amenable to analysis, and still provides an accurate model for large- $Pr$ fluids (Boeck & Thess 2001), including some silicone oils used in experiments (de Bruyn et al. 1996).

Our primary aim is to determine the best possible upper bound on $Nu$ when the background method is applied to the temperature field. To this end, we revisit Hagstrom & Doering’s background method analysis and derive a new upper-bounding variational principle for the Nusselt number that includes two so-called ‘balance parameters’ (Nicodemus, Grossmann & Holthaus 1997). One of these balance parameters can be optimised analytically, while the remaining one and the background temperature field can be combined to formulate a bound on $Nu$ in terms of a scaled background profile. We then employ convex programming to optimise the scaled background field for Marangoni numbers up to $Ma=10^{9}$ , and observe that the optimal bounds take the form $Nu\leqslant O(Ma^{2/7}(\ln Ma)^{-1/2})$ – a logarithmic improvement on Hagstrom & Doering’s bound.

We also seek to identify which features of the optimal scaled background temperature field are key to lowering the bound on $Nu$ . For instance, non-monotonicity plays an important role in the background method analysis for infinite- $Pr$ Rayleigh–Bénard convection (Plasting & Ierley 2005; Doering et al. 2006), and it is natural to ask if the same is true for the Bénard–Marangoni problem. Another important issue is whether one can expect to improve Hagstrom & Doering’s bound using a relatively simple background field, which is amenable to rigorous mathematical analysis. To answer these questions we utilise convex optimisation once again and minimise the bound on $Nu$ over two families of scaled background fields: those that decrease monotonically, and those constrained by convexity. Our results are supported by analysis of the variational principle for the bound, which also suggests a way to proceed with a rigorous mathematical proof.

Numerical optimisation of the bound on $Nu$ is central to this work, and our computational strategy deserves some remarks. Traditionally, the Euler–Lagrange equations for the optimal background field and balance parameters are derived, discretised and solved (see e.g. Plasting & Kerswell 2003; Wen et al. 2013, 2015). Instead, we discretise the variational problem for the bound to obtain a convex conic programme, i.e., a convex optimisation problem in which the variables are constrained to belong to a convex cone. The procedure is similar to that described in previous works by the authors (Fantuzzi & Wynn 2015, 2016a ), however here we use a different discretisation method. The first advantage of this approach is that very efficient software packages are available to solve conic programmes. The second is that additional linear constraints on the background field, such as monotonicity and convexity, can be included in a straightforward way and without any changes to the numerical optimisation algorithm. Conic programming, therefore, enables one to interrogate the bounding principle in a systematic way, in order to inform rigorous mathematical analysis. This applies not only to infinite- $Pr$ Bénard–Marangoni convection, but to any convex upper-bounding variational problem obtained from the application of the background method.

The outline of this work is the following. Section 2 introduces Pearson’s model (Pearson 1958) for Bénard–Marangoni convection at infinite Prandtl number, which is our starting point. We apply the background method with balance parameters to formulate an upper-bounding variational principle for the Nusselt number in § 3, and compare it to the one derived by Hagstrom & Doering (2010) in § 4. Section 5 is devoted to the numerical optimisation of the background fields, and describes our computational approach in detail. We discuss our results in § 6 with the help of additional analysis of the variational problem for the bound. Section 7 concludes the paper.

Our notation will be mostly standard. Upon non-dimensionalising, we consider a two-dimensional, horizontally periodic layer with domain $[0,2\unicode[STIX]{x03C0}]\times [0,1]$ , with $x$ and $z$ denoting the horizontal and vertical coordinates, respectively. The $L^{2}$ and $L^{\infty }$ norms in the $z$ direction of a generic quantity $q(z,\cdot )$ will be denoted by $\Vert \cdot \Vert _{2}$ and $\Vert \cdot \Vert _{\infty }$ , respectively, i.e.

(1.1a,b ) $$\begin{eqnarray}\displaystyle \Vert q\Vert _{2}:=\left(\int _{0}^{1}|q(z,\cdot )|^{2}\,\text{d}z\right)^{1/2},\quad \Vert q\Vert _{\infty }:=\sup _{z\in [0,1]}|q(z,\cdot )|. & & \displaystyle\end{eqnarray}$$

Overlines denote horizontal and infinite-time averages, while angle brackets indicate volume and infinite-time averages, i.e.

(1.2a,b ) $$\begin{eqnarray}\displaystyle \overline{q}(z):=\lim _{{\mathcal{T}}\rightarrow \infty }\frac{1}{{\mathcal{T}}}\int _{0}^{{\mathcal{T}}}\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}q(x,z,t)\,\text{d}x\,\text{d}t,\quad \langle q\rangle :=\int _{0}^{1}\overline{q}(z)\,\text{d}z. & & \displaystyle\end{eqnarray}$$

Since infinite-time averages need not exist in general, one could be more rigorous and replace $\lim$ with $\limsup$ . Note also that $\langle |q(z)|^{2}\rangle =\Vert q\Vert _{2}^{2}$ for any quantity $q$ that depends only on the vertical coordinate $z$ .

2 Pearson’s model

Consider a two-dimensional layer of incompressible fluid of depth $h$ , density $\unicode[STIX]{x1D70C}$ , kinematic viscosity $\unicode[STIX]{x1D708}$ , thermal diffusivity $\unicode[STIX]{x1D705}$ and thermal conductivity $\unicode[STIX]{x1D706}$  (the model and the results may be generalised to the three-dimensional case as described in Hagstrom & Doering (2010)). The fluid is heated from below at constant temperature, and cooled at the surface with a fixed heat flux $q$ . The problem is made non-dimensional using $h$ as the length unit, $h^{2}/\unicode[STIX]{x1D705}$ as the time unit and $qh/\unicode[STIX]{x1D706}$ as the temperature unit. When the Prandtl number $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ is infinite, Pearson’s equations for the fluid’s motion (Pearson 1958) reduce to (Hagstrom & Doering 2010)

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p=\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}T+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=\unicode[STIX]{x1D6FB}^{2}T, & \displaystyle\end{eqnarray}$$

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

where $\boldsymbol{u}(x,z,t)=u(x,z,t)\boldsymbol{i}+w(x,z,t)\boldsymbol{k}$ is the fluid’s velocity, $p(x,z,t)$ is the pressure and $T(x,z,t)$ is the temperature. All variables are assumed to be periodic in the horizontal direction (i.e. along the $x$ axis) with period $2\unicode[STIX]{x03C0}$ , and satisfy the vertical boundary conditions (BCs)

(2.2a-d ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}|_{z=0}=0,\quad w|_{z=1}=0,\quad T|_{z=0}=0,\quad \unicode[STIX]{x2202}_{z}T|_{z=1}=-1. & & \displaystyle\end{eqnarray}$$

The fluid is driven at the top boundary by surface tension forces due to local temperature gradients, which induce motion in the bulk of the layer through the action of viscosity. Mathematically, the situation is described by the additional BC

(2.3) $$\begin{eqnarray}[\unicode[STIX]{x2202}_{z}u+Ma\unicode[STIX]{x2202}_{x}T]_{z=1}=0.\end{eqnarray}$$

The Marangoni number $Ma=\unicode[STIX]{x1D6FE}qh^{2}/(\unicode[STIX]{x1D706}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$ , where $\unicode[STIX]{x1D6FE}$ is the negative of the derivative of the surface tension with respect to the fluid’s temperature, describes the ratio of surface tension to viscous forces, and is the governing non-dimensional parameter of the flow.

The purely conductive state $\boldsymbol{u}(x,z,t)=0$ , $p=\text{constant}$ , $T(x,z,t)=-z$ is asymptotically stable when $Ma\leqslant 66.84$  (Fantuzzi & Wynn 2017), while for $Ma\geqslant 79.61$ it is subject to linear instabilities (Pearson 1958) and convection sets in (Boeck & Thess 1998, 2001). Taking the divergence of (2.1a ) and using incompressibility shows that $\unicode[STIX]{x1D6FB}^{2}p=0$ , so taking the Laplacian of (2.1a ) gives

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{4}\boldsymbol{u}=0.\end{eqnarray}$$

Thus, each component of the ensuing convective velocity is bi-harmonic, and can be determined as a linear function of the temperature field, which forces (2.4) via the BC (2.3). In particular the horizontal Fourier coefficients ${\hat{w}}_{k}(z)$ , $k\in \mathbb{Z}$ , of the vertical velocity $w$ can be computed as a function of the horizontal Fourier coefficients $\hat{T}_{k}(z)$ of the temperature. One finds (Hagstrom & Doering 2010)

(2.5) $$\begin{eqnarray}{\hat{w}}_{k}(z)=-Maf_{k}(z)\hat{T}(1),\quad k\in \mathbb{Z},\end{eqnarray}$$

where $f_{0}(z)=0$ (so ${\hat{w}}_{0}=0$ and $w$ has zero horizontal mean), and

(2.6) $$\begin{eqnarray}f_{k}(z)=\frac{k\sinh k[kz\cosh (kz)-\sinh (kz)+(1-k\coth k)z\sinh (kz)]}{\sinh (2k)-2k},\quad k\in \mathbb{Z}\setminus \{0\}.\end{eqnarray}$$

Note that the function $f_{k}$ satisfies $f_{k}(z)\leqslant 0$ for $z\in [0,1]$ , $f_{k}(0)=0=f_{k}(1)$ and $f_{k}(z)\rightarrow 0$ pointwise for all $z\in (0,1)$ as $k\rightarrow \infty$ (see figure 1; note that the corresponding figure in Hagstrom & Doering’s original paper is incorrect: they plot the negative of $f_{k}$ ).

Figure 1. The function $f_{k}(z)$ for $k=1$ (dotted line), $k=3$ (dashed line), $k=10$ (dot-dashed line) and $k=100$ (solid line).

Convection enhances the vertical heat transport, and since the BC $\unicode[STIX]{x2202}_{z}T|_{z=1}=-1$ prescribes the heat flux through the top surface, the net effect is a reduction in the temperature drop across the layer. The key non-dimensional parameter to quantify this process is the Nusselt number

(2.7) $$\begin{eqnarray}Nu:=-\frac{1}{\overline{T}(1)}=\frac{1}{\langle |\unicode[STIX]{x1D735}T|^{2}\rangle },\end{eqnarray}$$

where $|\unicode[STIX]{x1D735}T|^{2}=(\unicode[STIX]{x2202}_{x}T)^{2}+(\unicode[STIX]{x2202}_{z}T)^{2}$ . The first equality in (2.7) defines the Nusselt number, while the second one can be proved by taking the volume and infinite-time average of $T\times$ (2.1b ), followed by appropriate integrations by parts using (2.1c ) and the BCs (for more details, see Hagstrom & Doering (2010)).

3 An upper-bounding variational principle for the Nusselt number

3.1 The background method with balance parameters

The background method analysis begins by decomposing the temperature variable as

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

where the steady background field $\unicode[STIX]{x1D70F}(z)$ satisfies the BCs

(3.2a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}(0)=0,\quad \unicode[STIX]{x1D70F}^{\prime }(1)=-1, & & \displaystyle\end{eqnarray}$$

while the time-dependent perturbation $\unicode[STIX]{x1D703}(x,z,t)$ is periodic in the horizontal direction and satisfies

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D703}|_{z=0}=0,\quad \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}|_{z=1}=0.\end{eqnarray}$$

Upon substituting this decomposition into (2.1a ) we obtain an evolution equation for the perturbation $\unicode[STIX]{x1D703}$ ,

(3.4) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D703}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D70F}^{\prime \prime }-w\unicode[STIX]{x1D70F}^{\prime }.\end{eqnarray}$$

Averaging $\unicode[STIX]{x1D703}\times$ (3.4) over the volume and infinite time, followed by appropriate integration by parts using (2.1c ) and the BCs for $\unicode[STIX]{x1D703}$ in (3.3), shows that

(3.5) $$\begin{eqnarray}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D70F}^{\prime }\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}\rangle +\overline{\unicode[STIX]{x1D703}}(1)=0.\end{eqnarray}$$

Moreover, substituting (3.1) into (2.7) gives the two identities

(3.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle Nu^{-1}+\overline{\unicode[STIX]{x1D703}}(1)+\unicode[STIX]{x1D70F}(1)=0, & \displaystyle\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle Nu^{-1}-\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D70F}^{\prime }\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle -\Vert \unicode[STIX]{x1D70F}^{\prime }\Vert _{2}^{2}=0. & \displaystyle\end{eqnarray}$$

Taking the linear combination $\unicode[STIX]{x1D6FC}\times$ (3.5) $-\unicode[STIX]{x1D6FD}\times$ (3.6a ) $+$ (3.6b ) for scalar balance parameters $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}\neq 1$ to be determined, using the fact that $\overline{\unicode[STIX]{x1D703}}(1)=\langle \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\rangle$ by virtue of (3.3), and rearranging yields

(3.7) $$\begin{eqnarray}\frac{1}{Nu}=-\frac{\Vert \unicode[STIX]{x1D70F}^{\prime }\Vert _{2}^{2}+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70F}(1)}{\unicode[STIX]{x1D6FD}-1}+\frac{\unicode[STIX]{x1D6FC}-1}{\unicode[STIX]{x1D6FD}-1}{\mathcal{Q}}\{\unicode[STIX]{x1D703},w\},\end{eqnarray}$$

where

(3.8) $$\begin{eqnarray}{\mathcal{Q}}\{\unicode[STIX]{x1D703},w\}=\left\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6FC}-1}\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}+\left(\frac{\unicode[STIX]{x1D6FC}-2}{\unicode[STIX]{x1D6FC}-1}\unicode[STIX]{x1D70F}^{\prime }+\frac{\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6FC}-1}\right)\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}\right\rangle .\end{eqnarray}$$

If the balance parameters are chosen to satisfy

(3.9) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FC}-1}{\unicode[STIX]{x1D6FD}-1}>0\end{eqnarray}$$

we can bound

(3.10) $$\begin{eqnarray}\frac{1}{Nu}\geqslant -\frac{\Vert \unicode[STIX]{x1D70F}^{\prime }\Vert _{2}^{2}+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70F}(1)}{\unicode[STIX]{x1D6FD}-1}+\frac{\unicode[STIX]{x1D6FC}-1}{\unicode[STIX]{x1D6FD}-1}\inf _{\unicode[STIX]{x1D703},w}{\mathcal{Q}}\{\unicode[STIX]{x1D703},w\},\end{eqnarray}$$

where the infimum is taken over all horizontally periodic fields $\unicode[STIX]{x1D703}$ that satisfy the BCs in (3.3) and over all velocity fields $w$ with horizontal Fourier coefficients given by (2.5). The key simplification is that we do not require $\unicode[STIX]{x1D703}$ to satisfy the nonlinear evolution equation (3.4). As a result, we may without any loss of generality restrict our attention to time-independent perturbations, and interpret $\langle \cdot \rangle$ in (3.8) as a volume average.

To compute the infimum in (3.10) we substitute the Fourier expansions for $\unicode[STIX]{x1D703}$ and $w$ into (3.8). Noticing that $\hat{\unicode[STIX]{x1D703}}_{k}=\hat{T}_{k}$ for $k\neq 0$ by virtue of (3.1), and that $f_{0}(\cdot )=0$ in (2.5), the Fourier coefficients ${\hat{w}}_{k}$ can be expressed in terms of $\hat{\unicode[STIX]{x1D703}}_{k}$ as

(3.11) $$\begin{eqnarray}{\hat{w}}_{k}(z)=-Maf_{k}(z)\unicode[STIX]{x1D703}_{k}(1),\quad k\in \mathbb{Z}.\end{eqnarray}$$

Moreover, $\hat{\unicode[STIX]{x1D703}}_{-k}=\hat{\unicode[STIX]{x1D703}}_{k}^{\ast }$ (where $^{\ast }$ denotes complex conjugation) because the Fourier modes must combine into the real-valued temperature perturbation $\unicode[STIX]{x1D703}$ . Consequently, we may rewrite

(3.12) $$\begin{eqnarray}{\mathcal{Q}}\{\unicode[STIX]{x1D703},w\}={\mathcal{Q}}_{0}\{\hat{\unicode[STIX]{x1D703}}_{0}\}+2\mathop{\sum }_{k\geqslant 1}{\mathcal{Q}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}\end{eqnarray}$$

where

(3.13) $$\begin{eqnarray}{\mathcal{Q}}_{0}\{\hat{\unicode[STIX]{x1D703}}_{0}\}:=\int _{0}^{1}\left[|{\hat{\unicode[STIX]{x1D703}}_{0}}^{\prime }(z)|^{2}+\left(\frac{\unicode[STIX]{x1D6FC}-2}{\unicode[STIX]{x1D6FC}-1}\unicode[STIX]{x1D70F}^{\prime }(z)+\frac{\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6FC}-1}\right){\hat{\unicode[STIX]{x1D703}}_{0}}^{\prime }(z)\right]\,\text{d}z,\end{eqnarray}$$

while for $k\geqslant 1$ the last term in (3.8) vanishes and we have

(3.14) $$\begin{eqnarray}{\mathcal{Q}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}:=\int _{0}^{1}\left\{|{\hat{\unicode[STIX]{x1D703}}_{k}}^{\prime }(z)|^{2}+k^{2}|\hat{\unicode[STIX]{x1D703}}_{k}(z)|^{2}-\frac{\unicode[STIX]{x1D6FC}Ma}{\unicode[STIX]{x1D6FC}-1}\unicode[STIX]{x1D70F}^{\prime }(z)f_{k}(z)Re[\hat{\unicode[STIX]{x1D703}}_{k}(1){\hat{\unicode[STIX]{x1D703}}_{k}(z)}^{\ast }]\right\}\,\text{d}z.\end{eqnarray}$$

Now, the infimum of ${\mathcal{Q}}\{\unicode[STIX]{x1D703},w\}$ must be negative semidefinite since ${\mathcal{Q}}\{0,0\}=0$ . Moreover, each functional ${\mathcal{Q}}_{k}$ , $k\geqslant 0$ must be individually lower bounded because among all perturbations $\unicode[STIX]{x1D703}$ , $w$ are those with only one horizontal wavenumber. In light of (3.3), this lower bound must be sought over all complex-valued functions $\hat{\unicode[STIX]{x1D703}}_{k}(z)$ that satisfy $\hat{\unicode[STIX]{x1D703}}_{k}(0)=0=\hat{\unicode[STIX]{x1D703}}_{k}^{\prime }(1)$ . Since ${\mathcal{Q}}_{0}\{0\}=0$ , the infimum of ${\mathcal{Q}}_{0}$ must be negative semidefinite. When $k\geqslant 1$ , instead, ${\mathcal{Q}}_{k}$ is a homogeneous functional and so if it is lower bounded, its infimum must be exactly zero. Consequently,

(3.15) $$\begin{eqnarray}\inf _{\unicode[STIX]{x1D703},w}{\mathcal{Q}}\{\unicode[STIX]{x1D703},w\}=\left\{\begin{array}{@{}ll@{}}\displaystyle \inf _{\hat{\unicode[STIX]{x1D703}}_{0}}{\mathcal{Q}}_{0}\{\hat{\unicode[STIX]{x1D703}}_{0}\}\quad & \text{if }{\mathcal{Q}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}\geqslant 0,k=1,2,\ldots ,\\ -\infty \quad & \text{otherwise.}\end{array}\right.\end{eqnarray}$$

In appendix A we show that

(3.16) $$\begin{eqnarray}\inf _{\hat{\unicode[STIX]{x1D703}}_{0}}{\mathcal{Q}}_{0}\{\hat{\unicode[STIX]{x1D703}}_{0}\}=-\frac{\Vert (\unicode[STIX]{x1D6FC}-2)\unicode[STIX]{x1D70F}^{\prime }+\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD}\Vert _{2}^{2}}{4(\unicode[STIX]{x1D6FC}-1)^{2}}.\end{eqnarray}$$

Substituting this into (3.10), and using the fact that

(3.17) $$\begin{eqnarray}\unicode[STIX]{x1D70F}(1)=\int _{0}^{1}\unicode[STIX]{x1D70F}^{\prime }(z)\,\text{d}z\end{eqnarray}$$

by virtue of (3.2) to simplify the resulting expression, we obtain

(3.18) $$\begin{eqnarray}\frac{1}{Nu}\geqslant -\frac{4\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FD}-1)\unicode[STIX]{x1D70F}(1)+\Vert \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}^{\prime }+\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD}\Vert _{2}^{2}}{4(\unicode[STIX]{x1D6FC}-1)(\unicode[STIX]{x1D6FD}-1)}.\end{eqnarray}$$

This bound is valid if (3.9) holds, and if the background field $\unicode[STIX]{x1D70F}$ is chosen to make the functional ${\mathcal{Q}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}$ in (3.14) positive semidefinite for all (integer) wavenumbers $k\geqslant 1$ . The latter set of constraints can be combined into the single condition that

(3.19) $$\begin{eqnarray}\left\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6FC}-1}\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}\right\rangle \geqslant 0\end{eqnarray}$$

for all perturbations $\unicode[STIX]{x1D703}$ , $w$ with zero horizontal mean that satisfy (3.3) and (3.11). Using well-established terminology, we refer to such $\unicode[STIX]{x1D703}$ and $w$ as admissible perturbations, and to (3.19) as the spectral constraint.

The best possible bound on $Nu$ is then found upon solving the following optimisation problem:

(3.20) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \sup _{\unicode[STIX]{x1D70F}(z),\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}} & \displaystyle -\frac{4\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FD}-1)\unicode[STIX]{x1D70F}(1)+\Vert \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}^{\prime }+\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD}\Vert _{2}^{2}}{4(\unicode[STIX]{x1D6FC}-1)(\unicode[STIX]{x1D6FD}-1)}\\ \displaystyle \text{subject to} & \displaystyle \left\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6FC}-1}\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}\right\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \frac{\unicode[STIX]{x1D6FC}-1}{\unicode[STIX]{x1D6FD}-1}>0,\\ & \displaystyle \unicode[STIX]{x1D70F}(0)=0,\\ & \displaystyle \unicode[STIX]{x1D70F}^{\prime }(1)=-1.\end{array}\end{eqnarray}$$

Note that we look for the supremum of the objective function (rather than its maximum) because the strict inequality $(\unicode[STIX]{x1D6FC}-1)/(\unicode[STIX]{x1D6FD}-1)>0$ may prevent the existence of a maximiser.

3.2 Optimisation over $\unicode[STIX]{x1D6FD}$

The lower bound (3.18) can be optimised over $\unicode[STIX]{x1D6FD}$ in a relatively straightforward way, because the spectral constraint is independent of $\unicode[STIX]{x1D6FD}$ . Upon setting to zero the first derivative of the right-hand side of (3.18) with respect to $\unicode[STIX]{x1D6FD}$ , and using (3.17) to rearrange, we find two stationary values,

(3.21a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{+}=1+\Vert \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}^{\prime }+\unicode[STIX]{x1D6FC}-1\Vert _{2},\quad \unicode[STIX]{x1D6FD}_{-}=1-\Vert \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}^{\prime }+\unicode[STIX]{x1D6FC}-1\Vert _{2}. & & \displaystyle\end{eqnarray}$$

Inspection of the second derivative of the right-hand side of (3.18) with respect to $\unicode[STIX]{x1D6FD}$ reveals that when $\unicode[STIX]{x1D6FC}$ is constrained by (3.9) both choices $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{+}$ and $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{-}$ correspond to a local maximum. Determining the optimal choice of $\unicode[STIX]{x1D6FD}$ therefore requires comparing the values of such local maxima.

After choosing $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{+}$ and re-parametrising $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D706}/(\unicode[STIX]{x1D706}-1)$ – with $\unicode[STIX]{x1D706}>1$ to satisfy (3.9) – we can use (3.17) to rewrite (3.18) as

(3.22) $$\begin{eqnarray}\frac{1}{Nu}\geqslant \frac{1-\Vert \unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}^{\prime }+1\Vert _{2}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}(1)}{2}.\end{eqnarray}$$

The spectral constraint (3.19) can also be expressed in terms of $\unicode[STIX]{x1D706}$ as

(3.23) $$\begin{eqnarray}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w.\end{eqnarray}$$

Upon introducing the scaled background field $\unicode[STIX]{x1D70C}(z)=\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}(z)=\unicode[STIX]{x1D6FC}/(\unicode[STIX]{x1D6FC}-1)\unicode[STIX]{x1D70F}(z)$ , subject to a suitably scaled version of the BCs in (3.2), the optimal bound on $Nu$ corresponding to the choice $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{+}$ is found by solving the variational problem

(3.24) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \sup _{\unicode[STIX]{x1D70C}(z),\unicode[STIX]{x1D706}} & \displaystyle \frac{1-\Vert \unicode[STIX]{x1D70C}^{\prime }+1\Vert _{2}-\unicode[STIX]{x1D70C}(1)}{2}\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D70C}^{\prime }w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \unicode[STIX]{x1D70C}(0)=0,\\ & \displaystyle \unicode[STIX]{x1D70C}^{\prime }(1)=-\unicode[STIX]{x1D706},\\ & \displaystyle \unicode[STIX]{x1D706}>1.\end{array}\end{eqnarray}$$

Similar steps show that the best possible bound on $Nu$ when setting $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{-}$ in (3.18) is given by the solution of an optimisation problem that differs from (3.24) only in the constraint for $\unicode[STIX]{x1D706}$ ,

(3.25) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \sup _{\unicode[STIX]{x1D70C}(z),\unicode[STIX]{x1D706}} & \displaystyle \frac{1-\Vert \unicode[STIX]{x1D70C}^{\prime }+1\Vert _{2}-\unicode[STIX]{x1D70C}(1)}{2},\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D70C}^{\prime }w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \unicode[STIX]{x1D70C}(0)=0,\\ & \displaystyle \unicode[STIX]{x1D70C}^{\prime }(1)=-\unicode[STIX]{x1D706},\\ & \displaystyle \unicode[STIX]{x1D706}<1.\end{array}\end{eqnarray}$$

The key observation at this stage is that the suprema in (3.24) and (3.25) coincide despite the different constraint on $\unicode[STIX]{x1D706}$ , and furthermore they are equal to the optimal value of the variational problem

(3.26) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D70C}(z)} & \displaystyle \frac{1-\Vert \unicode[STIX]{x1D70C}^{\prime }+1\Vert _{2}-\unicode[STIX]{x1D70C}(1)}{2},\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D70C}^{\prime }w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \unicode[STIX]{x1D70C}(0)=0.\end{array}\end{eqnarray}$$

In fact, for any value of $\unicode[STIX]{x1D706}$ we can construct a feasible $\unicode[STIX]{x1D70C}(z)$ for either (3.24) or (3.25) that approximates the solution of (3.26) arbitrarily accurately: simply let $\unicode[STIX]{x1D70C}_{0}(z)$ be an $\unicode[STIX]{x1D700}$ -suboptimal strictly feasible point for (3.26), and choose $\unicode[STIX]{x1D70C}^{\prime }(z)=\unicode[STIX]{x1D70C}_{0}^{\prime }(z)$ in (3.24) or (3.25) except for an infinitesimally thin layer near $z=1$ , where $\unicode[STIX]{x1D70C}^{\prime }(z)=-\unicode[STIX]{x1D706}$ . A rigorous argument follows steps similar to those used in the energy stability analysis of the conductive state (Fantuzzi & Wynn 2017), and is omitted for brevity. The conclusion is satisfactory: the bound on $Nu$ is independent of whether one sets $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{+}$ or $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{-}$ in (3.18).

3.3 An explicit value for the optimal $\unicode[STIX]{x1D6FD}$

The variational principle (3.26) has been obtained by optimising the balance parameter $\unicode[STIX]{x1D6FD}$ as a function of the other balance parameter, $\unicode[STIX]{x1D6FC}$ , and the background field $\unicode[STIX]{x1D70F}(z)$ . Interestingly, the optimality conditions for the solution $\unicode[STIX]{x1D70C}_{\star }(z)$ of (3.26) allow for the derivation of a precise numerical value for the optimal $\unicode[STIX]{x1D6FD}$ even though the optimal $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}(z)$ are unknown. To show this, we introduce a variable $s$ such that $\Vert \unicode[STIX]{x1D70C}^{\prime }+1\Vert _{2}\leqslant s$ and note that (3.26) is equivalent to

(3.27) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D70C}(z),s} & \displaystyle 1-s-\unicode[STIX]{x1D70C}(1),\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D70C}^{\prime }w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \unicode[STIX]{x1D70C}(0)=0,\\ & \displaystyle \Vert \unicode[STIX]{x1D70C}^{\prime }+1\Vert _{2}\leqslant s.\end{array}\end{eqnarray}$$

The feasible set of this problem is convex, so the linear objective function is maximised on the constraint boundary. Since for any given $\unicode[STIX]{x1D70C}(z)$ we can always choose $s=\Vert \unicode[STIX]{x1D70C}^{\prime }+1\Vert _{2}$ , the optimal bound is attained when $\unicode[STIX]{x1D70C}(z)$ is on the boundary of the feasible set of the spectral constraint, i.e. when

(3.28) $$\begin{eqnarray}\inf _{\unicode[STIX]{x1D703},w\neq 0}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D70C}^{\prime }w\unicode[STIX]{x1D703}\rangle =0.\end{eqnarray}$$

Since the spectral constraint is homogeneous in $\unicode[STIX]{x1D703}$ and $w$ , it suffices to restrict our attention to admissible $\unicode[STIX]{x1D703}$ and $w$ satisfying some normalisation condition ${\mathcal{N}}\{\unicode[STIX]{x1D703},w\}=0$ that excludes the zero fields. The optimal scaled background field $\unicode[STIX]{x1D70C}_{\star }(z)$ and the optimal value $s_{\star }$ are then those that maximise the Lagrangian functional

(3.29) $$\begin{eqnarray}\displaystyle {\mathcal{L}}\{\unicode[STIX]{x1D70C},s,\unicode[STIX]{x1D703},w,\unicode[STIX]{x1D701},\unicode[STIX]{x1D702},\unicode[STIX]{x1D707}\} & := & \displaystyle 1-s-\unicode[STIX]{x1D70C}(1)+\unicode[STIX]{x1D701}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D70C}^{\prime }w\unicode[STIX]{x1D703}\rangle \nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D702}(s^{2}-\Vert \unicode[STIX]{x1D70C}^{\prime }+1\Vert _{2}^{2})+\unicode[STIX]{x1D707}{\mathcal{N}}\{\unicode[STIX]{x1D703},w\},\end{eqnarray}$$

where $\unicode[STIX]{x1D701}$ , $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D707}$ are scalar Lagrange multipliers.

Setting to zero the first variation of ${\mathcal{L}}$ with respect to $\unicode[STIX]{x1D70C}(z)$ shows that the optimal scaled background field $\unicode[STIX]{x1D70C}_{\star }(z)$ must satisfy the ‘natural’ boundary condition

(3.30) $$\begin{eqnarray}1+2\unicode[STIX]{x1D702}+2\unicode[STIX]{x1D702}\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)=0.\end{eqnarray}$$

(Of course, $\unicode[STIX]{x1D70C}_{\star }(z)$ must also satisfy an Euler–Lagrange differential equation, but this will not be important here.) Moreover, setting to zero the derivatives of ${\mathcal{L}}$ with respect to $s$ and $\unicode[STIX]{x1D702}$ , and eliminating $s$ yields

(3.31) $$\begin{eqnarray}2\unicode[STIX]{x1D702}\Vert \unicode[STIX]{x1D70C}_{\star }^{\prime }+1\Vert _{2}-1=0.\end{eqnarray}$$

At this point, note that if $\unicode[STIX]{x1D70C}^{\prime }(z)=-1$ the spectral constraint (3.27) reduces to the condition for global ‘energy’ stability of the conduction solution (see e.g. Fantuzzi & Wynn 2017), which cannot be satisfied in the convective regime. Consequently, $\Vert \unicode[STIX]{x1D70C}_{\star }^{\prime }+1\Vert _{2}\neq 0$ and we may use (3.31) to eliminate $\unicode[STIX]{x1D702}$ from (3.30). The optimal scaled background field must therefore satisfy

(3.32) $$\begin{eqnarray}1+\Vert \unicode[STIX]{x1D70C}_{\star }^{\prime }+1\Vert _{2}+\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)=0.\end{eqnarray}$$

In particular, this implies that $-\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)>1$ , so $\unicode[STIX]{x1D70C}_{\star }(z)$ is also the optimal solution of (3.24) with $\unicode[STIX]{x1D706}=-\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)$ . Recollecting the re-parametrisation $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D706}/(\unicode[STIX]{x1D706}-1)$ we conclude that the optimal value of the balance parameter $\unicode[STIX]{x1D6FC}$ , denoted $\unicode[STIX]{x1D6FC}_{\star }$ , is given by

(3.33) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{\star }=\frac{\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)}{\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)+1}.\end{eqnarray}$$

Finally, recalling that (3.24) was obtained by setting $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{+}$ from (3.21) and that $\unicode[STIX]{x1D6FC}_{\star }\unicode[STIX]{x1D70F}_{\star }^{\prime }(z)/(\unicode[STIX]{x1D6FC}_{\star }-1)=\unicode[STIX]{x1D70C}_{\star }^{\prime }(z)$ according to our rescaling, we can apply (3.33) and (3.32) in succession to conclude that the optimal value of the balance parameter $\unicode[STIX]{x1D6FD}$ is

(3.34) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{\star }=1+(\unicode[STIX]{x1D6FC}_{\star }-1)\Vert \unicode[STIX]{x1D70C}_{\star }^{\prime }+1\Vert _{2}=\frac{\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)+1-\Vert \unicode[STIX]{x1D70C}_{\star }^{\prime }+1\Vert _{2}}{\unicode[STIX]{x1D70C}_{\star }^{\prime }(1)+1}=2.\end{eqnarray}$$

4 Relation to Hagstrom & Doering’s variational problem

The bounding principle formulated by Hagstrom & Doering (2010) can be recovered upon setting $\unicode[STIX]{x1D6FC}=2$ and $\unicode[STIX]{x1D6FD}=2$ in (3.20). These values clearly satisfy (3.9), and we have seen that the choice $\unicode[STIX]{x1D6FD}=2$ is optimal. The variational problem for the optimal background field becomes

(4.1) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D70F}(z)} & \displaystyle -\Vert \unicode[STIX]{x1D70F}^{\prime }\Vert _{2}^{2}-2\unicode[STIX]{x1D70F}(1),\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+2\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \unicode[STIX]{x1D70F}(0)=0.\end{array}\end{eqnarray}$$

Strictly speaking we should also enforce the boundary condition $\unicode[STIX]{x1D70F}^{\prime }(1)=-1$ , but this does not limit the choice of $\unicode[STIX]{x1D70F}$ for the same reasons discussed at the end of § 3.2.

To bring (4.1) in contact with (3.26), we change variables to $\unicode[STIX]{x1D711}=2\unicode[STIX]{x1D70F}$ and use the boundary condition $\unicode[STIX]{x1D711}(0)=0$ to rewrite (4.1) as

(4.2) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D711}(z)} & \displaystyle \frac{1-\Vert \unicode[STIX]{x1D711}^{\prime }+1\Vert _{2}^{2}-2\unicode[STIX]{x1D711}(1)}{4},\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+\unicode[STIX]{x1D711}^{\prime }w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \unicode[STIX]{x1D711}(0)=0.\end{array}\end{eqnarray}$$

It is clear that (3.26) and (4.2) have the same feasible set. It is also not difficult to show that the optimal value of (3.26) is no smaller than that of (4.2); in fact, for any feasible $\unicode[STIX]{x1D711}(z)$

(4.3) $$\begin{eqnarray}\frac{1-\Vert \unicode[STIX]{x1D711}^{\prime }+1\Vert _{2}-\unicode[STIX]{x1D711}(1)}{2}-\frac{1-\Vert \unicode[STIX]{x1D711}^{\prime }+1\Vert _{2}^{2}-2\unicode[STIX]{x1D711}(1)}{4}=\left(\frac{1-\Vert \unicode[STIX]{x1D711}^{\prime }+1\Vert _{2}}{2}\right)^{2}\geqslant 0.\end{eqnarray}$$

In particular, using (3.26) it is almost immediate to obtain a 4.2 % improvement for the prefactor of Hagstrom & Doering’s bound, $Nu\leqslant 0.838Ma^{2/7}$ , at least in the limit of infinite Marangoni number: in appendix B we show that

(4.4) $$\begin{eqnarray}Nu\leqslant 0.803\times Ma^{2/7}\quad \text{as }Ma\rightarrow \infty .\end{eqnarray}$$

What is not immediately apparent when comparing (4.2) to (3.26) is that fixing the balance parameters a priori does not change the asymptotic behaviour of the optimal bounds as $Ma\rightarrow \infty$ . To show that this is true, recall from § 3.3 that the choice $\unicode[STIX]{x1D6FD}=2$ is optimal. After fixing $\unicode[STIX]{x1D6FD}=2$ and re-parametrising $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D706}/(\unicode[STIX]{x1D706}-1)$ as in § 3 – with $\unicode[STIX]{x1D706}>1$ to satisfy (3.9) – the bound in (3.18) becomes

(4.5) $$\begin{eqnarray}\frac{1}{Nu}\geqslant 1-\frac{\unicode[STIX]{x1D706}^{2}}{4(\unicode[STIX]{x1D706}-1)}\Vert \unicode[STIX]{x1D70F}^{\prime }+1\Vert _{2}^{2}.\end{eqnarray}$$

From this point onwards, the analysis is analogous to that of the infinite- $Pr$ Rayleigh–Bénard problem (Plasting 2004, chapter 6). First, let $w=Ma\tilde{w}$ and define the scaled Marangoni number $M=\unicode[STIX]{x1D706}Ma$ to rewrite the spectral constraint (3.23) as

(4.6) $$\begin{eqnarray}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+M\unicode[STIX]{x1D70F}^{\prime }\tilde{w}\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},\tilde{w}.\end{eqnarray}$$

Upon rescaling $w=Ma\tilde{w}$ the Marangoni number drops out of equation (3.11), so $\tilde{w}$ is a (linear) function of $\unicode[STIX]{x1D703}$ only and the admissible test functions in (4.6) are independent of $M$ . Then, consider the family of background fields $\unicode[STIX]{x1D70F}_{M}$ , parametrised by the scaled Marangoni number $M$ , that maximises the right-hand side of (4.5) for a fixed value of $\unicode[STIX]{x1D706}$ . In other words, assume that $\unicode[STIX]{x1D70F}_{M}$ solves the variational problem

(4.7) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \min _{\unicode[STIX]{x1D70F}(z)} & \Vert \unicode[STIX]{x1D70F}^{\prime }+1\Vert _{2}^{2}\\ \text{subject to} & \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+M\unicode[STIX]{x1D70F}^{\prime }\tilde{w}\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},\tilde{w}.\end{array}\end{eqnarray}$$

Moreover, suppose $\unicode[STIX]{x1D70E}(M)$ is such that

(4.8) $$\begin{eqnarray}\Vert \unicode[STIX]{x1D70F}_{M}^{\prime }+1\Vert _{2}^{2}=1-\unicode[STIX]{x1D70E}(M).\end{eqnarray}$$

Note that it is reasonable to assume that $\unicode[STIX]{x1D70E}(M)\rightarrow 0$ as $M\rightarrow \infty$ , because we expect that $\unicode[STIX]{x1D70F}^{\prime }(z)\approx 0$ except for thin boundary layers if the scaled spectral constraint (4.6) is to be satisfied. The optimal bound for a given value of $\unicode[STIX]{x1D706}$ is then given by

(4.9) $$\begin{eqnarray}\frac{1}{Nu}\geqslant \frac{\unicode[STIX]{x1D706}^{2}\unicode[STIX]{x1D70E}(M)-(\unicode[STIX]{x1D706}-2)^{2}}{4(\unicode[STIX]{x1D706}-1)}.\end{eqnarray}$$

Using the fact that $\text{d}M/\text{d}\unicode[STIX]{x1D706}=\text{d}(\unicode[STIX]{x1D706}Ma)/\text{d}\unicode[STIX]{x1D706}=M/\unicode[STIX]{x1D706}$ , it is straightforward to show that the right-hand side of the last expression is maximised with respect to $\unicode[STIX]{x1D706}$ when

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{2-2\unicode[STIX]{x1D70E}(M)-M\unicode[STIX]{x1D70E}^{\prime }(M)}{1-\unicode[STIX]{x1D70E}(M)-M\unicode[STIX]{x1D70E}^{\prime }(M)},\end{eqnarray}$$

and that the corresponding bound on the Nusselt number is

(4.11) $$\begin{eqnarray}Nu\leqslant \frac{4-4\unicode[STIX]{x1D70E}(M)-4M\unicode[STIX]{x1D70E}^{\prime }(M)}{4\unicode[STIX]{x1D70E}(M)-[2\unicode[STIX]{x1D70E}(M)+M\unicode[STIX]{x1D70E}^{\prime }(M)]^{2}}.\end{eqnarray}$$

Now, recall that $Nu\geqslant 1$ since convection enhances the purely conductive vertical heat transport. Using the fact that $\unicode[STIX]{x1D706}>1$ and the assumption that $\unicode[STIX]{x1D70E}(M)\rightarrow 0$ as $M\rightarrow 0$ it is then not difficult to see that since $Nu\geqslant 1$ the quantity $M\unicode[STIX]{x1D70E}^{\prime }(M)$ must be uniformly bounded as the scaled Marangoni number $M$ tends to infinity. Consequently, the solution $\unicode[STIX]{x1D706}_{\star }=\unicode[STIX]{x1D706}_{\star }(M)$ of (4.10) satisfies

(4.12) $$\begin{eqnarray}\lim _{M\rightarrow \infty }\unicode[STIX]{x1D706}_{\star }(M)=O(1).\end{eqnarray}$$

This implies that $Ma=O(M)$ as $M\rightarrow \infty$ , meaning that the optimisation over the balance parameter does not influence the asymptotic scaling of the bound on $Nu$ with the Marangoni number.

5 Optimal bounds

We now turn our attention to the numerical solution of the variational problem (3.26). To implement our computational strategy, described in § 5.1 below, it is convenient to change variables once more and let

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(z):=\int _{0}^{z}[\unicode[STIX]{x1D719}(\unicode[STIX]{x1D709})-1]\,\text{d}\unicode[STIX]{x1D709},\end{eqnarray}$$

so the boundary condition $\unicode[STIX]{x1D70C}(0)=0$ is satisfied. Since $\unicode[STIX]{x1D70C}^{\prime }(z)=\unicode[STIX]{x1D719}(z)-1$ , (3.26) can be rewritten as

(5.2) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D719}(z)} & \displaystyle 1-\frac{1}{2}\Vert \unicode[STIX]{x1D719}\Vert _{2}-\frac{1}{2}\int _{0}^{1}\unicode[STIX]{x1D719}(z)\,\text{d}z,\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+(\unicode[STIX]{x1D719}-1)w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w.\end{array}\end{eqnarray}$$

Moreover, we introduce a non-negative variable $s$ such that $\Vert \unicode[STIX]{x1D719}\Vert _{2}\leqslant s$ . After dropping the constant $1$ as well as a factor of $1/2$ from the objective function, it is not difficult to see that the optimal solution of (5.2) is the same as that of the convex problem

(5.3) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D719}(z),s} & \displaystyle -s-\int _{0}^{1}\unicode[STIX]{x1D719}(z)\,\text{d}z,\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+(\unicode[STIX]{x1D719}-1)w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \Vert \unicode[STIX]{x1D719}\Vert _{2}\leqslant s.\end{array}\end{eqnarray}$$

As anticipated in § 1, we are also interested in optimising the bound on $Nu$ over the restricted classes of monotonically decreasing and convex scaled background fields, i.e. such that $\unicode[STIX]{x1D70C}^{\prime }(z)\leqslant 0$ and $\unicode[STIX]{x1D70C}^{\prime \prime }(z)\geqslant 0$ . This is achieved by solving the convex problems

(5.4) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D719}(z),s} & \displaystyle -s-\int _{0}^{1}\unicode[STIX]{x1D719}(z)\,\text{d}z,\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+(\unicode[STIX]{x1D719}-1)w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \Vert \unicode[STIX]{x1D719}\Vert _{2}\leqslant s,\\ & \displaystyle \unicode[STIX]{x1D719}(z)\leqslant 1,\end{array}\end{eqnarray}$$

and

(5.5) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{\unicode[STIX]{x1D719}(z),s} & \displaystyle -s-\int _{0}^{1}\unicode[STIX]{x1D719}(z)\,\text{d}z,\\ \displaystyle \text{subject to} & \displaystyle \langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}+(\unicode[STIX]{x1D719}-1)w\unicode[STIX]{x1D703}\rangle \geqslant 0\quad \forall \text{admissible }\unicode[STIX]{x1D703},w,\\ & \displaystyle \Vert \unicode[STIX]{x1D719}\Vert _{2}\leqslant s,\\ & \displaystyle \unicode[STIX]{x1D719}^{\prime }(z)\geqslant 0.\end{array}\end{eqnarray}$$

5.1 Computational methodology

Our computational methodology is based on the observation that the constraints in (5.3)–(5.5) are the infinite-dimensional equivalent of well-known types of finite-dimensional convex constraints. As already pointed out in previous work (Fantuzzi & Wynn 2016a ) the spectral constraint is the infinite-dimensional equivalent of a linear matrix inequality (LMI), the condition that a symmetric matrix $\unicode[STIX]{x1D64E}$ whose entries are affine with respect to a set of optimisation variables is positive semidefinite (denoted by $\unicode[STIX]{x1D64E}\succcurlyeq 0$ ). The norm constraint $\Vert \unicode[STIX]{x1D719}\Vert _{2}\leqslant s$ , instead, is the infinite-dimensional version of a second-order cone constraint (SOCC), i.e. the requirement that a vector $\boldsymbol{y}\in \mathbb{R}^{n+1}$ and a scalar $s$ satisfy $\Vert \boldsymbol{y}\Vert \leqslant s$ , where $\Vert \cdot \Vert$ denotes the usual Euclidean norm of a vector. Finally, the pointwise constraints $\unicode[STIX]{x1D719}(z)\leqslant 1$ and $\unicode[STIX]{x1D719}^{\prime }(z)\geqslant 0$ are the infinite-dimensional equivalent of element-wise inequalities for a vector $\boldsymbol{y}\in \mathbb{R}^{n+1}$ of the form $\boldsymbol{A}\boldsymbol{y}\leqslant \boldsymbol{b}$ , with $\boldsymbol{A}\in \mathbb{R}^{m\times (n+1)}$ and $\boldsymbol{b}\in \mathbb{R}^{m}$ given. For more details on LMIs and SOCCs we refer the reader to the works by Boyd et al. (1994) and Boyd & Vandenberghe (2004). Optimisation problems with LMIs, SOCCs and element-wise vector inequalities are well-known instances of so-called conic programmes, and can be solved to high accuracy in polynomial time (Vandenberghe & Boyd 1996; Boyd & Vandenberghe 2004). Consequently, problems (5.3)–(5.5) can be solved numerically if we discretise them to obtain conic programmes.

In order to reduce the norm constraint $\Vert \unicode[STIX]{x1D719}\Vert _{2}\leqslant s$ to a SOCC, we introduce a piecewise-linear ansatz for $\unicode[STIX]{x1D719}$ . Given a set of $n+1$ collocation points $0=z_{0}<z_{1}<\cdots <z_{n-1}<z_{n}=1$ , we denote $\unicode[STIX]{x1D719}_{i}=\unicode[STIX]{x1D719}(z_{i})$ for all $i=1,\ldots ,n$ and consider

(5.6) $$\begin{eqnarray}\unicode[STIX]{x1D719}(z)=\mathop{\sum }_{i=0}^{n}\unicode[STIX]{x1D719}_{i}\unicode[STIX]{x1D713}_{i}(z),\end{eqnarray}$$

where $\unicode[STIX]{x1D713}_{i}(z)$ is the unique piecewise-linear function satisfying $\unicode[STIX]{x1D713}_{i}(z_{i})=1$ and vanishing at all other nodes (cf. figure 2). After defining the column vector of nodal values

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D731}:=[\unicode[STIX]{x1D719}_{0},\ldots ,\unicode[STIX]{x1D719}_{n}]^{\text{T}}\in \mathbb{R}^{n+1},\end{eqnarray}$$

it is clear that there exists a positive definite matrix $\unicode[STIX]{x1D64B}=\unicode[STIX]{x1D64D}^{\text{T}}\unicode[STIX]{x1D64D}$ such that

(5.8) $$\begin{eqnarray}\Vert \unicode[STIX]{x1D719}\Vert _{2}=\left(\int _{0}^{1}\mathop{\sum }_{i,j=0}^{n}\unicode[STIX]{x1D719}_{i}\unicode[STIX]{x1D719}_{j}\unicode[STIX]{x1D713}_{i}(z)\unicode[STIX]{x1D713}_{j}(z)\,\text{d}z\right)^{1/2}=(\unicode[STIX]{x1D731}^{\text{T}}\unicode[STIX]{x1D64B}\unicode[STIX]{x1D731})^{1/2}=\Vert \unicode[STIX]{x1D64D}\unicode[STIX]{x1D731}\Vert .\end{eqnarray}$$

The norm constraint $\Vert \unicode[STIX]{x1D719}\Vert _{2}\leqslant s$ then becomes the SOCC $\Vert \unicode[STIX]{x1D64D}\unicode[STIX]{x1D731}\Vert \leqslant s$ .

Figure 2. Sketch of the piecewise-linear function $\unicode[STIX]{x1D713}_{i}(z)$ .

The spectral constraint can be reduced to a set of LMIs in a similar way. Recall from § 3.1 that the spectral constraint is equivalent to the functional ${\mathcal{Q}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}$ in (3.14) being positive semidefinite for all wavenumbers $k\geqslant 1$ and all complex-valued functions $\hat{\unicode[STIX]{x1D703}}_{k}(z)$ satisfying $\hat{\unicode[STIX]{x1D703}}_{k}(0)=0=\hat{\unicode[STIX]{x1D703}}_{k}^{\prime }(1)$ . Recognising that the real and imaginary parts of $\hat{\unicode[STIX]{x1D703}}_{k}$ give identical and independent contributions to ${\mathcal{Q}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}$ , it suffices to restrict our attention to real-valued functions $\hat{\unicode[STIX]{x1D703}}_{k}(z)$ , so we define the space of test functions

(5.9) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}:=\left\{v(z):[0,1]\rightarrow \mathbb{R},\int _{0}^{1}\left(|v^{\prime }(z)|^{2}+|v(z)|^{2}\right)\,\text{d}z<\infty ,v(0)=0,v^{\prime }(1)=0\right\}.\end{eqnarray}$$

Recalling that we have changed variables according to

(5.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6FC}-1}\unicode[STIX]{x1D70F}^{\prime }(z)=\unicode[STIX]{x1D70C}^{\prime }(z)=\unicode[STIX]{x1D719}(z)-1,\end{eqnarray}$$

we can therefore replace the spectral constraint in (5.3)–(5.5) with the infinite set of Fourier-transformed spectral constraints

(5.11) $$\begin{eqnarray}\displaystyle {\mathcal{Q}}_{k}\{v\} & = & \displaystyle \int _{0}^{1}\left\{|v^{\prime }(z)|^{2}+k^{2}|v(z)|^{2}-Ma[\unicode[STIX]{x1D719}(z)-1]f_{k}(z)v(1)v(z)\right\}\,\text{d}z\geqslant 0\nonumber\\ \displaystyle & & \displaystyle \forall v\in \unicode[STIX]{x1D6E4},k=1,2,\ldots .\end{eqnarray}$$

In appendix C we show that ${\mathcal{Q}}_{k}\{v\}$ is positive semidefinite for a candidate $\unicode[STIX]{x1D719}(z)$ whenever

(5.12) $$\begin{eqnarray}k>k_{c}:=\left\lfloor \left(\frac{3\sqrt{3}}{128}\right)^{1/4}Ma^{1/2}\Vert \unicode[STIX]{x1D719}-1\Vert _{\infty }^{1/2}\right\rfloor ,\end{eqnarray}$$

where $\lfloor \cdot \rfloor$ denotes the integer part of a number. The ‘cutoff’ wavenumber $k_{c}$ represents an upper bound on the largest critical wavenumber, i.e. the largest values of $k$ for which the infimum of the functional ${\mathcal{Q}}_{k}$ in (5.11) over non-zero test functions is zero. When $k\leqslant k_{c}$ , instead, we approximate the test function $v$ using the same piecewise-linear ansatz used for $\unicode[STIX]{x1D719}$ , i.e.

(5.13) $$\begin{eqnarray}v(z)=\mathop{\sum }_{i=0}^{n}v_{i}\unicode[STIX]{x1D713}_{i}(z).\end{eqnarray}$$

We also set $v_{0}=0$ and $v_{n}=v_{n-1}$ in order to enforce the boundary conditions $v(0)=0$ and $v^{\prime }(1)=0$ , but we do not do this explicitly in (5.13) to simplify the following discussion. Substituting (5.13) and (5.6) into ${\mathcal{Q}}_{k}\{v\}$ from (5.11) yields

(5.14) $$\begin{eqnarray}\displaystyle & & \displaystyle {\mathcal{Q}}_{k}\{v\}=\mathop{\sum }_{i,j=0}^{n}v_{i}v_{j}\int _{0}^{1}[\unicode[STIX]{x1D713}_{i}(z)^{\prime }\unicode[STIX]{x1D713}_{j}(z)^{\prime }+k^{2}\unicode[STIX]{x1D713}_{i}(z)\unicode[STIX]{x1D713}_{j}(z)]\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle \quad +\,Ma\mathop{\sum }_{i=0}^{n}v_{n}v_{i}\int _{0}^{1}\unicode[STIX]{x1D713}_{i}(z)f_{k}(z)\,\text{d}z-Ma\mathop{\sum }_{i,j=0}^{n}\unicode[STIX]{x1D719}_{i}v_{n}v_{j}\int _{0}^{1}\unicode[STIX]{x1D713}_{i}(z)\unicode[STIX]{x1D713}_{j}(z)f_{k}(z)\,\text{d}z.\end{eqnarray}$$

Recollecting that we have set $v_{0}=0$ and $v_{n}=v_{n-1}$ , the right-hand side of (5.14) is a quadratic form of the vector of nodal values $\boldsymbol{v}:=[v_{1},\ldots ,v_{n-1}]^{\text{T}}$ , and there exists an $(n-1)\times (n-1)$ symmetric matrix $\unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})$ , affine with respect to $\unicode[STIX]{x1D731}$ , such that ${\mathcal{Q}}_{k}\{v\}=\boldsymbol{v}^{\text{T}}\unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})\boldsymbol{v}$ . Consequently, for each wavenumber $k\leqslant k_{c}$ the Fourier-transformed spectral constraint ${\mathcal{Q}}_{k}\{v\}\geqslant 0$ can be approximated by the LMI $\unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})\succcurlyeq 0$ .

Finally, it is easy to see that the piecewise-linear approximation (5.6) turns the pointwise inequality $\unicode[STIX]{x1D719}(z)\leqslant 1$ into the $n+1$ constraints $\unicode[STIX]{x1D719}_{i}\leqslant 1$ , $i=0,\ldots ,n$ , which can be written succinctly as the element-wise vector inequality $\unicode[STIX]{x1D731}\leqslant \mathbf{1}$ . Similarly, the condition $\unicode[STIX]{x1D719}^{\prime }(z)\geqslant 0$ becomes a set of $n$ inequalities $\unicode[STIX]{x1D719}_{i-1}-\unicode[STIX]{x1D719}_{i}\leqslant 0$ , $i=1,\ldots ,n$ , which can be written in the vector form

(5.15) $$\begin{eqnarray}\boldsymbol{A}\unicode[STIX]{x1D731}\leqslant \mathbf{0},\quad \boldsymbol{A}:=\left[\begin{array}{@{}cccc@{}}1 & -1\\ & \ddots & \ddots \\ & & 1 & -1\end{array}\right]\in \mathbb{R}^{n\times (n+1)}.\end{eqnarray}$$

After substituting (5.6) into the objective function of (5.3) and defining

(5.16) $$\begin{eqnarray}\boldsymbol{c}:=\left[\int _{0}^{1}\unicode[STIX]{x1D713}_{0}(z)\,\text{d}z,\ldots ,\int _{0}^{1}\unicode[STIX]{x1D713}_{n}(z)\,\text{d}z\right]^{\text{T}},\end{eqnarray}$$

we can therefore approximate the infinite-dimensional variational problem (5.3) with the finite-dimensional conic programme

(5.17) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{s,\unicode[STIX]{x1D731}} & -s-\boldsymbol{c}^{\text{T}}\unicode[STIX]{x1D731}\\ \text{subject to} & \unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})\succcurlyeq 0,\quad k=1,\ldots ,k_{c},\\ & \Vert \unicode[STIX]{x1D64D}\unicode[STIX]{x1D731}\Vert \leqslant s.\end{array}\end{eqnarray}$$

Similarly, (5.4) can be approximated as

(5.18) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{s,\unicode[STIX]{x1D731}} & -s-\boldsymbol{c}^{\text{T}}\unicode[STIX]{x1D731}\\ \text{subject to} & \unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})\succcurlyeq 0,\quad k=1,\ldots ,k_{c},\\ & \Vert \unicode[STIX]{x1D64D}\unicode[STIX]{x1D731}\Vert \leqslant s,\\ & \unicode[STIX]{x1D731}\leqslant \mathbf{1},\end{array}\end{eqnarray}$$

while (5.5) becomes

(5.19) $$\begin{eqnarray}\begin{array}{@{}rl@{}}\displaystyle \max _{s,\unicode[STIX]{x1D731}} & -s-\boldsymbol{c}^{\text{T}}\unicode[STIX]{x1D731}\\ \text{subject to} & \unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})\succcurlyeq 0,\quad k=1,\ldots ,k_{c},\\ & \Vert \unicode[STIX]{x1D64D}\unicode[STIX]{x1D731}\Vert \leqslant s,\\ & \boldsymbol{A}\unicode[STIX]{x1D731}\leqslant \mathbf{0}.\end{array}\end{eqnarray}$$

Before describing our numerical implementation of (5.17)–(5.19) in more detail, let us note some important aspects of our finite-dimensional approximations.

The first observation is that introducing the piecewise-linear ansatz (5.6) for $\unicode[STIX]{x1D719}(z)$ means that only lower bounds on the optimal values of (5.3)–(5.5) can be computed, because the ‘true’ optimal $\unicode[STIX]{x1D719}(z)$ is unlikely to be piecewise linear. Moreover, assuming (5.13) enforces the Fourier-transformed spectral constraint only over a particular subset of the test function space $\unicode[STIX]{x1D6E4}$ . This enlarges the set of feasible functions $\unicode[STIX]{x1D719}(z)$ , so (5.17)–(5.19) yield upper limits for lower bounds of the true optimal values of (5.3)–(5.5), respectively. However, one expects the solutions of each conic programme (5.17)–(5.19) to converge to that of the corresponding maximisation problem (5.3)–(5.5) as the number of collocation points in the spatial discretisation increases.

One could also estimate the error between functions in $\unicode[STIX]{x1D6E4}$ and their finite-dimensional approximation, in order to formulate conic programmes whose optimal solutions bound the optimal value of (5.3)–(5.5) rigorously from below. This is possible if global polynomial approximation is utilised when all but the test functions in the spectral constraint are polynomials (Fantuzzi et al. 2017). One could follow a similar line of reasoning in each sub-interval of our piecewise-linear approximation, with the additional complication that the function $f_{k}$ appearing in the spectral constraint is not polynomial. However, we do not do so here because we do not aim to compute bounds on the Nusselt number to the standard of a computer-assisted proof.

Finally, as already pointed out in § 1, a major advantage of our computational methodology is that the monotonicity and convexity constraints can be implemented in a very straightforward way. On the contrary, optimising the bound on $Nu$ over all monotonic or convex background fields seems considerably more challenging if one follows the classical Euler–Lagrange variational approach, because one has to solve a set of differential equations coupled to an inequality (in fact, a differential inequality in the convex case).

5.2 Implementation details

The conic programmes (5.17)–(5.19) were implemented using the MATLAB toolbox YALMIP (Löfberg 2004) and solved with the conic solver SDPT3 (Toh, Todd & Tütüncü 1999; Tütüncü, Toh & Todd 2003). Sparsity was exploited using chordal decomposition methods (Fukuda et al. 2000; Nakata et al. 2003; Kim et al. 2011). All computations were run on a PC with a 3.40 GHz Intel^® Core™ i7-4770 CPU and 16 Gb of RAM.

As collocation points, we used the Chebyshev nodes $z_{i}=[1-\cos (\unicode[STIX]{x03C0}i/n)]/2$ , $i=0,\ldots ,n$ in the sub-interval $(0.05,0.98)$ , and the finer distribution $z_{i}=[1-\cos (\unicode[STIX]{x03C0}i/4n)]/2$ , $i=0,\ldots ,4n$ in the boundary sub-intervals $[0,0.05]$ and $[0.98,1]$ . After initial experiments we set $n=512$ , giving $873$ collocation points in total; our results, presented in § 5.3, change by less than 0.1 % if a larger $n$ is used.

Chebyshev nodes were chosen as they naturally cluster near the boundaries and help resolve boundary layers near $z=0$ and $z=1$ in the optimal $\unicode[STIX]{x1D719}(z)$ . These are expected even if no boundary conditions are imposed because to maximise the objective function in (5.3) one would like to choose $\unicode[STIX]{x1D719}(z)<0$ , but setting $\unicode[STIX]{x1D719}(z)\approx 1$ in the bulk of the domain is necessary to be able to satisfy the spectral constraint. However, it is possible to have $\unicode[STIX]{x1D719}(z)<0$ in thin layers near the walls because the functions $f_{k}$ , which act as a weight on $\unicode[STIX]{x1D719}$ in the Fourier-transformed spectral constraint (5.11), are small there for all $k$ values (cf. figure 1). These observations are confirmed by the numerical results presented in § 5.3.

While boundary layers can in principle be resolved with a sufficiently fine distribution of Chebyshev points, we preferred to refine the discretisation only near the boundaries using a secondary set of Chebyshev nodes to limit the cost of our computations. A precise assessment of the computational burden of our conic programmes relies on technical details of the sparsity-exploiting methods we used (Fukuda et al. 2000; Nakata et al. 2003; Kim et al. 2011) and is beyond the scope of this work. Here, we simply note that it must grow at least linearly with the number of collocation points. Approximately speaking, in fact, sparsity allows replacing each LMI $\unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})\succcurlyeq 0$ with a set of LMIs on certain $3\times 3$ submatrices of $\unicode[STIX]{x1D64C}_{k}$ . Since the number of rows/columns of $\unicode[STIX]{x1D64C}_{k}$ grows linearly with the number of discretisation points, so does the number of such submatrices. Even assuming optimistically that the computational cost of one such $3\times 3$ LMI is fixed and that handling a much larger number of LMIs has negligible overhead, the overall computational cost can grow no slower than linearly with the number of collocation nodes.

One complication to the implementation of (5.17) is that the cutoff wavenumber $k_{c}$ is not known a priori, but it depends on $\unicode[STIX]{x1D731}$ according to (5.12). We therefore employ the following iterative procedure: find the optimal $\unicode[STIX]{x1D731}$ using an initial guess $k_{0}$ for $k_{c}$ , update the value of $k_{c}$ using (5.12), check if $\unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})$ is positive semidefinite for all $k\leqslant k_{c}$ , and repeat the optimisation with the updated guess for $k_{c}$ if any of these checks fail.

A second hurdle is that solving (5.17) with this iterative procedure becomes expensive when the Marangoni number is large because the cutoff wavenumber $k_{c}$ , and therefore the number of LMI constraints, grows proportionally to $Ma^{1/2}$ . For example, at $Ma=2.5\times 10^{6}$ we find that the optimal $\unicode[STIX]{x1D719}$ satisfies $\Vert \unicode[STIX]{x1D719}-1\Vert _{\infty }=2$ , so (5.12) gives $k_{c}=1003$ ; when all $1003$ LMIs are considered in (5.17), SDPT3 takes more than 4 h to converge on our machine. In an effort to reduce the CPU time requirements, we implemented a trial-and-error procedure, inspired by the numerical continuation method employed by Plasting & Kerswell (2003), in which only a subset of wavenumbers are considered in (5.17). More precisely, we progressively increased the Marangoni number according to the update rule $Ma_{i+1}=Ma_{i}\times 10^{p}$ , which gives $p+1$ logarithmically spaced points between successive powers of 10. Given the critical wavenumbers $k_{1},\ldots ,k_{m}$ at one Marangoni number, we solved the optimisation problem for the next $Ma$ considering only wavenumbers in a window of width $2r$ around each $k_{i}$ , $i=1,\ldots ,m$ , i.e. values of $k$ such that

(5.20) $$\begin{eqnarray}k\in \mathop{\bigcup }_{i=1}^{m}[k_{i}-r,k_{i}+r].\end{eqnarray}$$

We then checked if the optimal solution satisfied $\unicode[STIX]{x1D64C}_{k}(\unicode[STIX]{x1D731})\succcurlyeq 0$ for all remaining wavenumbers up to the cutoff value $k_{c}$ . If any of these checks failed, we added the wavenumber with the largest constraint violation (i.e. corresponding to the matrix $\unicode[STIX]{x1D64C}_{k}$ with the most negative eigenvalue) to the list of critical values and repeated the optimisation.

5.3 Results

The conic programmes (5.17)–(5.19) were successfully solved for Marangoni numbers up to $Ma=10^{9}$ using the procedure described in § 5.2 with $p=19$ and $r=10$ . In all cases, at each value of $Ma$ the optimal $\unicode[STIX]{x1D719}(z)$ was used to recover the optimal scaled background field $\unicode[STIX]{x1D70C}(z)$ and the corresponding bound on the Nusselt number.

The most important results of our computations are the bounds on $Nu$ , which are plotted in figure 3. Also shown for comparison are: the analytical bound $Nu\leqslant 0.803Ma^{2/7}$ from § 4; the DNS results obtained by Boeck & Thess (2001); the conductive value $Nu=1$ , which bounds the Nusselt number from below. The results are plotted in two ways: compensated by a factor of $Ma^{-2/7}$ to aid the visual comparison with the asymptotic scaling of the analytical bound, and compensated by $Ma^{-2/7}(\ln Ma)^{1/2}$ .

The main observation is that while a gap with the DNS data remains, the fully optimal bounds and those computed after enforcing convexity grow more slowly than the analytical bound by $(\ln Ma)^{1/2}$ . In particular, the fully optimal bounds exhibit the asymptotic behaviour

(5.21) $$\begin{eqnarray}Nu\leqslant 1.285Ma^{2/7}(\ln Ma)^{-1/2}.\end{eqnarray}$$

In contrast, when the background field is constrained to decrease monotonically the bound on $Nu$ asymptotes to $0.535Ma^{2/7}$ . This suggests that the analytical bound of § 4 attains the optimal asymptotic scaling available when $\unicode[STIX]{x1D70C}(z)$ is monotonic, but it may be lowered by a logarithm upon construction of a non-monotonic background field.

Figure 3. Comparison between: the fully optimal bounds on the Nusselt number, computed using the solution of (5.3) (solid line); the optimal monotonic bounds, computed using the solution of (5.4) (dot-dashed line); the optimal convex bounds, computed using the solution of (5.5) (thick dotted line). Also shown are the conductive Nusselt number $Nu=1$ (dashed line), the analytical bound $Nu\leqslant 0.803Ma^{2/7}$ proved in § 4 (dotted line), and the DNS data (Boeck & Thess 2001, crosses). In panel (a), the data are compensated by $Ma^{-2/7}$ to facilitate the visual comparison with the asymptotic scaling of the analytical bound. In panel (b), the data are compensated by $Ma^{-2/7}(\ln Ma)^{1/2}$ .

Figure 4. Normalised derivatives of the optimal background fields, $\unicode[STIX]{x1D70C}^{\prime }(z)/|\unicode[STIX]{x1D70C}^{\prime }(0)|$ , obtained with (5.17) (solid line), with (5.18) (dot-dashed lines) and with (5.19) (dotted line) for (a) $Ma=100$ , (b) $Ma=186.12$ , (c) $Ma=10^{3}$ , (d) $Ma=10^{4}$ , (e) $Ma=10^{5}$ , and (f) $Ma=10^{6}$ . Insets in (e,f) show a detailed view of the boundary layers near $z=1$ .

Figure 5. (a) Boundary value $\unicode[STIX]{x1D70C}^{\prime }(0)$ for the fully optimal (solid line), monotonic (dot-dashed line) and convex (dotted line) background fields. All curves almost coincide. (b) Plot of the convergence measure $\unicode[STIX]{x1D70C}^{\prime }(0)+2$ , scaled by $Ma^{2/7}(\ln Ma)^{-1/2}$ (solid line) and by $Ma^{2/7}$ (dashed line), for the fully optimal background fields.

Figure 4 shows the derivative of the optimal scaled background field, computed with each of the conic programmes (5.17)–(5.19), for selected values of $Ma$ . We plot $\unicode[STIX]{x1D70C}^{\prime }(z)$ instead of $\unicode[STIX]{x1D70C}(z)$ because by virtue of (5.1) problems (5.3)–(5.5) can be rewritten in terms of $\unicode[STIX]{x1D70C}^{\prime }(z)$ alone. Since $\unicode[STIX]{x1D70C}(z)$ can be recovered by integration using the boundary condition $\unicode[STIX]{x1D70C}(0)=0$ , the derivative $\unicode[STIX]{x1D70C}^{\prime }(z)$ is the actual decision variable in (5.3)–(5.5). Moreover, to ease the comparison the profiles have been normalised by the magnitude of the boundary value $\unicode[STIX]{x1D70C}^{\prime }(0)$ , which converges to $-2$ as $Ma$ grows, as illustrated in figure 5(a). Figure 5(b) demonstrates that in the fully optimal case the convergence is logarithmic; this was also observed when convexity was imposed, while power-law convergence was observed for the monotonic profiles (these results are not shown for brevity). Such evidence corroborates our conjecture that (5.21) is the correct functional form of the optimal bound on $Nu$ .

As illustrated by figure 4, the optimal $\unicode[STIX]{x1D70C}^{\prime }(z)$ is negative for $Ma\leqslant 186.12$ , meaning that the corresponding scaled background field decreases monotonically for sufficiently small Marangoni numbers. When $Ma$ is raised, all profiles are characterised by boundary layers separated by a bulk region where $\unicode[STIX]{x1D70C}^{\prime }(z)\approx 0$ . Note that the transition to the bulk region is not smooth when monotonicity or convexity are enforced, and this is the main reason for preferring the piecewise-linear approximations of § 5.1 to the global polynomial approximation used in previous works (Fantuzzi & Wynn 2016a , b ; Fantuzzi et al. 2017).

In the fully optimal case, $\unicode[STIX]{x1D70C}^{\prime }(z)$ changes sign inside both boundary layers to reach positive local maxima, so the corresponding scaled background field is characterised by non-monotonic boundary layers. Enforcing monotonicity removes these local maxima and makes the boundary layers thinner, while convexity prevents the local maximum near $z=0$ and makes $\unicode[STIX]{x1D70C}^{\prime }(z)$ constant across the boundary layer near $z=1$ .

Figure 6. Details of the boundary layer structure of the fully optimal scaled background field derivative $\unicode[STIX]{x1D70C}^{\prime }(z)$ for $Ma\geqslant 10^{4}$ . The dot-dashed, dashed and solid lines in (b) indicate the approximate scaling laws (5.23a )–(5.23c ), respectively.

Further details of the boundary layer structure of the fully optimal profiles for $Ma\geqslant 10^{4}$ are given in figure 6 (very similar results for the optimal convex profiles are not shown for brevity). Letting $z_{bot}$ and $z_{top}$ denote the coordinates of the positive local maxima of $\unicode[STIX]{x1D70C}^{\prime }(z)$ near $z=0$ and $z=1$ , respectively, we take $\unicode[STIX]{x1D6FF}:=z_{bot}$ and $\unicode[STIX]{x1D700}:=1-z_{top}$ as a measure of the thickness of each boundary layer. The boundary layer near the bottom of the domain ( $z=0$ ) becomes approximately self-similar at large Marangoni numbers, and least-squares power-law fits to the data in figure 6(a,c) for $Ma\geqslant 10^{7}$ return

(5.22a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\approx 3.8Ma^{-0.26},\quad \unicode[STIX]{x1D70C}^{\prime }(z_{bot})\approx 0.07. & & \displaystyle\end{eqnarray}$$

Note that the scaling exponent of $\unicode[STIX]{x1D6FF}$ is not far from $-2/7\approx -0.286$ , suggesting that the width of the boundary layer near $z=0$ is one of the leading factors determining the scaling of the bound on $Nu$ . In fact, we conjecture that asymptotically $\unicode[STIX]{x1D6FF}=O(Ma^{-2/7}(\ln Ma)^{1/2})$ , such that $Nu=O(\unicode[STIX]{x1D6FF}^{-1})$ , but unfortunately the finite precision of our data does not permit us to clearly identify logarithmic corrections. To obtain more precise values we should solve the conic programmes (5.17)–(5.18) to a level of accuracy beyond the capabilities of SDPT3, as well as study larger Marangoni numbers (this issue will be discussed further in § 6).

The situation is more complicated for the boundary layer near $z=1$ . In figure 6(b) we can identify three distinct regions characterised by different scaling laws for $\unicode[STIX]{x1D700}$ :

(5.23a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}\approx 1.65Ma^{-0.36}\quad \text{for }Ma\lessapprox 5\times 10^{4}, & \displaystyle\end{eqnarray}$$

(5.23b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}\approx 11.8Ma^{-0.54}\quad \text{for }5\times 10^{4}\lessapprox Ma\lessapprox 3\times 10^{6}, & \displaystyle\end{eqnarray}$$

(5.23c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}\approx 24.3Ma^{-0.58}\quad \text{for }Ma\gtrapprox 3\times 10^{6}. & \displaystyle\end{eqnarray}$$

In the first and third regions we could also determine approximate scaling laws for the peak value $\unicode[STIX]{x1D70C}^{\prime }(z_{top})$ :

(5.24a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{\prime }(z_{top})\approx 0.34Ma^{-0.01}\quad \text{for }Ma\lessapprox 5\times 10^{4}, & \displaystyle\end{eqnarray}$$

(5.24b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{\prime }(z_{top})\approx 0.04Ma^{0.16}\quad \text{for }Ma\gtrapprox 3\times 10^{6}. & \displaystyle\end{eqnarray}$$

Once again, these scaling laws are only tentative due to the finite precision to which the conic programmes for the optimal bounds could be solved. However, we remark that the large scatter in the data points in figure 6(b) is simply due to plotting $\unicode[STIX]{x1D700}$ after rescaling by $Ma^{0.54}$ , which at large $Ma$ amplifies small inaccuracies in our numerical data.

Figure 7. Bifurcation diagrams for the critical wavenumbers in: (a) the conic programme (5.17) for the fully optimal background fields; (b) the conic programme (5.18) for the optimal monotonic background fields; (c) the conic programme (5.19) for the optimal convex background fields.

Changes in the scaling of the boundary layer near $z=1$ correspond to bifurcations in the critical wavenumbers for the conic programme (5.17). As illustrated in figure 7(a), new critical wavenumbers appear at large values of $k$ for $Ma\approx 5\times 10^{4}$ and $Ma\approx 3\times 10^{6}$ . Another intermediate branch of critical wavenumbers appears for $Ma\approx 10^{8}$ , but this does not seem to influence the scaling of the boundary layer. Such bifurcations can be explained in terms of the interactions in the Fourier-transformed spectral constraint (5.11) between the boundary layer of $\unicode[STIX]{x1D70C}^{\prime }(z)=\unicode[STIX]{x1D719}(z)-1$ and the function $f_{k}(z)$ , which is almost entirely supported near $z=1$ at large $k$ . As shown in figure 7(b,c), similar bifurcations were observed when solving (5.19) but not when solving (5.18), probably because the boundary layer near $z=1$ of the optimal monotonic background fields is too thin to allow interesting interactions for wavenumbers below the cutoff value $k_{c}$ .

Figure 8. (a) Convergence of the optimal balance parameter $\unicode[STIX]{x1D6FC}_{\star }$ , computed using the optimal solution of (5.3) and (3.33), to the asymptotic value $2$ . (b) Plot of the difference $\unicode[STIX]{x1D6FC}_{\star }-2$ , scaled by $Ma^{2/7}(\ln Ma)^{-1/2}$ (solid line) and by $Ma^{2/7}$ (dashed line).

To conclude this section, we plot in figure 8 the variation with $Ma$ of the optimal balance parameter $\unicode[STIX]{x1D6FC}_{\star }$ , computed using (3.33) and the fully optimal background field. The results are interesting for two reasons. First, the convergence of $\unicode[STIX]{x1D6FC}_{\star }$ to $2$ as $Ma$ is raised is logarithmic, giving further evidence in support of (5.21). Second, the results suggest that the choice $\unicode[STIX]{x1D6FC}=2$ in the original work by Hagstrom & Doering (2010) – presumably motivated only by the convenience of eliminating the linear terms when combining (3.5), (3.6a ) and (3.6b ) in the background method analysis – is optimal, at least as far as the asymptotic behaviour of the bound as $Ma\rightarrow \infty$ is concerned. Contrary to what has been observed in previous works (see e.g. Plasting & Kerswell 2003; Wen et al. 2015), this means that not only does the optimisation of the balance parameters have no influence on the asymptotic scaling of the bound (cf. § 4), but it also does not improve the optimal prefactor available to Hagstrom & Doering’s original upper-bounding principle.

6 Discussion

6.1 Towards an improved bound

The results presented in § 5.3 suggest that Hagstrom & Doering’s bound $Nu\leqslant O(Ma^{2/7})$ may be improved by the logarithmic factor $(\ln Ma)^{-1/2}$ . Despite the strong numerical evidence, however, whether the optimal bound scales logarithmically when $Ma\rightarrow \infty$ remains uncertain due to the limited range of Marangoni numbers spanned in the present investigation (see § 6.2 for more on this issue). In particular, we cannot rule out the occurrence of further bifurcations in the critical wavenumbers that may cause a transition to a pure power-law behaviour with scaling exponent of $2/7$ .

Uncertainty about the true asymptotic scaling notwithstanding, our numerical results demonstrate that if the current analytical bound $Nu\leqslant O(Ma^{2/7})$ can be improved, to do so requires a background temperature profile with non-monotonic boundary layers. More precisely, the optimal convex background fields and the corresponding bounds on $Nu$ are evidence that what is needed is a relatively simple non-monotonic boundary layer near $z=1$ , while non-monotonicity near $z=0$ only lowers the prefactor.

Taking advantage of these observations to improve the bound on $Nu$ analytically, however, is likely to require a careful analysis of the sign-indefinite term in each Fourier-transformed spectral constraint, which we restate here in terms of the variable $\unicode[STIX]{x1D70C}(z)$ in the slightly rearranged form

(6.1) $$\begin{eqnarray}{\mathcal{Q}}_{k}\{v\}=\Vert v^{\prime }\Vert _{2}^{2}+k^{2}\Vert v\Vert _{2}^{2}-Mav(1)\int _{0}^{1}\unicode[STIX]{x1D70C}^{\prime }(z)f_{k}(z)v(z)\,\text{d}z\geqslant 0\quad \forall v\in \unicode[STIX]{x1D6E4}.\end{eqnarray}$$

For example, simply estimating

(6.2) $$\begin{eqnarray}\left|Mav(1)\int _{0}^{1}\unicode[STIX]{x1D70C}^{\prime }(z)f_{k}(z)v(z)\,\text{d}z\right|\leqslant Ma|v(1)|\int _{0}^{1}|\unicode[STIX]{x1D70C}^{\prime }(z)||\,f_{k}(z)||v(z)|\,\text{d}z\end{eqnarray}$$

and requiring (as we have done in appendix B) that

(6.3) $$\begin{eqnarray}\Vert v^{\prime }\Vert _{2}^{2}+k^{2}\Vert v\Vert _{2}^{2}-Ma|v(1)|\int _{0}^{1}|\unicode[STIX]{x1D70C}^{\prime }(z)||\,f_{k}(z)||v(z)|\,\text{d}z\geqslant 0\quad \forall v\in \unicode[STIX]{x1D6E4}\end{eqnarray}$$

forces the optimal $\unicode[STIX]{x1D70C}$ to decrease monotonically. In fact, if $\unicode[STIX]{x1D70C}$ satisfies (6.3) and $\unicode[STIX]{x1D70C}^{\prime }(z)\geqslant 0$ for $z\in {\mathcal{U}}\subset [0,1]$ , the profile

(6.4) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D70C}}^{\prime }(z):=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}^{\prime }(z),\quad & z\in [0,1]\smallsetminus {\mathcal{U}},\\ 0,\quad & z\in {\mathcal{U}},\end{array}\right.\end{eqnarray}$$

also satisfies (6.3), but decreases monotonically and gives a larger objective value in (3.26). In light of the numerical results presented in § 5.3, we expect that any bound obtained using the estimate (6.2) will not be better than $Nu\leqslant O(Ma^{2/7})$ .

A better approach is to reformulate the Fourier-transformed spectral constraint (6.1) before applying any estimates. Without any loss of generality, let $\unicode[STIX]{x1D6FF}\in (0,1)$ and write

(6.5) $$\begin{eqnarray}\unicode[STIX]{x1D70C}^{\prime }(z)=\left\{\begin{array}{@{}ll@{}}g(z),\quad & 0\leqslant z\leqslant \unicode[STIX]{x1D6FF},\\ h(z),\quad & \unicode[STIX]{x1D6FF}\leqslant z\leqslant 1.\end{array}\right.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FF}$ represents the thickness of the boundary layer of the optimal background field near $z=0$ . With this choice, the Fourier-transformed spectral constraint (6.1) becomes

(6.6) $$\begin{eqnarray}\displaystyle {\mathcal{Q}}_{k}\{v\} & = & \displaystyle \Vert v^{\prime }\Vert _{2}^{2}+k^{2}\Vert v\Vert _{2}^{2}-Mav(1)\int _{0}^{\unicode[STIX]{x1D6FF}}g(z)f_{k}(z)v(z)\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle -\,Mav(1)\int _{\unicode[STIX]{x1D6FF}}^{1}h(z)f_{k}(z)v(z)\,\text{d}z\geqslant 0\quad \forall v\in \unicode[STIX]{x1D6E4}.\end{eqnarray}$$

Since this inequality is homogeneous in $v$ and holds when $v(1)=0$ , we may restrict attention to test functions normalised such that $v(1)=1$ . Upon adding and subtracting $Ma\int _{\unicode[STIX]{x1D6FF}}^{1}h(z)f_{k}(z)\,\text{d}z$ we then need to check that

(6.7) $$\begin{eqnarray}\displaystyle & & \displaystyle \Vert v^{\prime }\Vert _{2}^{2}+k^{2}\Vert v\Vert _{2}^{2}-Ma\int _{0}^{\unicode[STIX]{x1D6FF}}g(z)f_{k}(z)v(z)\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle \quad +\,Ma\int _{\unicode[STIX]{x1D6FF}}^{1}h(z)f_{k}(z)[1-v(z)]\,\text{d}z-Ma\int _{\unicode[STIX]{x1D6FF}}^{1}h(z)f_{k}(z)\,\text{d}z\geqslant 0.\end{eqnarray}$$

If $\int _{\unicode[STIX]{x1D6FF}}^{1}h(z)f_{k}(z)\,\text{d}z<0$ , the last term in (6.7) gives a net positive contribution to the spectral constraint, and can be used to control the sign-indefinite terms. Recalling from figure 1 that $f_{k}(z)\leqslant 0$ , this requires $h(z)>0$ over a sufficient portion of the interval $(\unicode[STIX]{x1D6FF},1)$ , so the background field $\unicode[STIX]{x1D70C}$ does not decrease monotonically. Moreover, $h(z)$ should be supported in a boundary layer near $z=1$ to be able to control the fourth term in (6.7). Consequently, a non-monotonic boundary layer near $z=1$ helps to enforce the spectral constraint. The situation is similar in infinite- $Pr$ Rayleigh–Bénard convection (Doering et al. 2006; Otto & Seis 2011), so this observation is perhaps not surprising.

In addition to casting light on the role of the surface boundary layer, identity (6.7) may also offer a starting point to improve the bound $Nu\leqslant O(Ma^{2/7})$ analytically. Recalling the boundary condition $v(0)=0$ and that $v(1)=1$ by virtue of our choice of normalisation for $v$ , one possible approach is to use the fundamental theorem of calculus and the Cauchy–Schwarz inequality to bound

(6.8) $$\begin{eqnarray}\displaystyle \left|\int _{0}^{\unicode[STIX]{x1D6FF}}g(z)f_{k}(z)v(z)\,\text{d}z\right| & {\leqslant} & \displaystyle \int _{0}^{\unicode[STIX]{x1D6FF}}|g(z)f_{k}(z)|\left|\int _{0}^{z}v^{\prime }(t)\,\text{d}t\right|\,\text{d}z\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \Vert v^{\prime }\Vert _{2}\int _{0}^{\unicode[STIX]{x1D6FF}}|g(z)f_{k}(z)|\sqrt{z}\,\text{d}z\end{eqnarray}$$

and

(6.9) $$\begin{eqnarray}\displaystyle \left|\int _{\unicode[STIX]{x1D6FF}}^{1}h(z)f_{k}(z)[1-v(z)]\,\text{d}z\right| & {\leqslant} & \displaystyle \int _{\unicode[STIX]{x1D6FF}}^{1}|h(z)f_{k}(z)|\left|\int _{z}^{1}v^{\prime }(t)\,\text{d}t\right|\,\text{d}z\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \Vert v^{\prime }\Vert _{2}\int _{\unicode[STIX]{x1D6FF}}^{1}|h(z)f_{k}(z)|\sqrt{1-z}\,\text{d}z.\end{eqnarray}$$

Defining

(6.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle a_{k}:=\int _{0}^{\unicode[STIX]{x1D6FF}}|g(z)f_{k}(z)|\sqrt{z}\,\text{d}z, & \displaystyle\end{eqnarray}$$

(6.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle b_{k}:=\int _{\unicode[STIX]{x1D6FF}}^{1}|h(z)f_{k}(z)|\sqrt{1-z}\,\text{d}z, & \displaystyle\end{eqnarray}$$

(6.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle c_{k}:=-\int _{\unicode[STIX]{x1D6FF}}^{1}h(z)f_{k}(z)\,\text{d}z & \displaystyle\end{eqnarray}$$

to ease the notation, a sufficient condition for (6.7) is that

(6.11) $$\begin{eqnarray}\Vert v^{\prime }\Vert _{2}^{2}-Ma(a_{k}+b_{k})\Vert v^{\prime }\Vert _{2}+Ma\,c_{k}\geqslant 0,\end{eqnarray}$$

which in turn is satisfied if

(6.12) $$\begin{eqnarray}a_{k}+b_{k}\leqslant 2\sqrt{\frac{c_{k}}{Ma}}.\end{eqnarray}$$

Given a candidate background field, condition (6.12) can be checked for all wavenumbers up to the ‘cutoff’ wavenumber $k_{c}$ in (5.12).

Improving the bound $Nu\leqslant O(Ma^{2/7})$ , however, may not be straightforward. To illustrate one of the difficulties, let us consider a simple background field. Motivated by figure 5(b) and the shape of the derivatives of the optimal convex background fields in figure 4, we fix

(6.13) $$\begin{eqnarray}\displaystyle g(z)=-2,\quad h(z)=\left\{\begin{array}{@{}ll@{}}0,\quad & \unicode[STIX]{x1D6FF}\leqslant z<1-\unicode[STIX]{x1D700},\\ \unicode[STIX]{x1D6FE},\quad & 1-\unicode[STIX]{x1D700}\leqslant z\leqslant 1,\end{array}\right. & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FE}>0$ a constant (independent of the Marangoni number) and $\unicode[STIX]{x1D700}\ll 1$ but such that $1/\unicode[STIX]{x1D700}\leqslant k_{c}=O(Ma^{1/2})$ . When $k\leqslant 1/\unicode[STIX]{x1D700}$ we can use the Taylor expansions $f_{k}(z)=O(z^{2})$ near $z=0$ and $f_{k}(z)=O(k(z-1))$ near $z=1$ to estimate

(6.14a-c ) $$\begin{eqnarray}\displaystyle a_{k}=O(\unicode[STIX]{x1D6FF}^{7/2}),\quad b_{k}=O(\unicode[STIX]{x1D6FE}k\unicode[STIX]{x1D700}^{5/2}),\quad c_{k}=O(\unicode[STIX]{x1D6FE}k\unicode[STIX]{x1D700}^{2}). & & \displaystyle\end{eqnarray}$$

Using these estimates, (6.12) can be rearranged as

(6.15) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}^{7/2}\leqslant O\left(2\unicode[STIX]{x1D700}\sqrt{\frac{\unicode[STIX]{x1D6FE}k}{Ma}}(1-\sqrt{\unicode[STIX]{x1D6FE}kMa\unicode[STIX]{x1D700}^{3}})\right).\end{eqnarray}$$

When $k=O(1)$ the two sides of (6.15) could be balanced by taking $\unicode[STIX]{x1D700}=O(Ma^{-1/3})$ and $\unicode[STIX]{x1D6FF}=O(Ma^{-5/21})$ , and upon computing the bound on $Nu$ we find

(6.16) $$\begin{eqnarray}Nu\leqslant \frac{2}{2\unicode[STIX]{x1D6FF}-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D700}-\sqrt{\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D6FE}+2)\unicode[STIX]{x1D700}}}=O\left(\frac{1}{\unicode[STIX]{x1D6FF}}\right)=O(Ma^{5/21}).\end{eqnarray}$$

Interestingly, the exponent $5/21\approx 0.238$ is extremely close to that of the best power-law fit $Nu=O(Ma^{0.24})$ to the DNS data by Boeck & Thess (2001, see equation (4) in their paper). In these simulations convection takes the form of stationary rolls with energy only at low wavenumbers, and the deviation from the theoretical asymptotic scaling exponent $2/9\approx 0.222$ can be attributed to the contribution to the heat transfer of the thermal boundary layer near the surface (see the discussion after equation (13) in Boeck & Thess (2001)). Although this contribution is expected to vanish as $Ma\rightarrow \infty$ , the background method could yield a bound that agrees well with observations at least over a finite range of Marangoni numbers if the stability of the rolls were deduced rigorously from the governing equations. Given the lack of such information, however, (6.12) must be satisfied for all wavenumbers up to the cutoff value $k_{c}=O(Ma^{1/2})$ . In particular, setting $k=1/\unicode[STIX]{x1D700}$ (which is no larger than $k_{c}$ by assumption) shows that we must choose $\unicode[STIX]{x1D700}\leqslant O(Ma^{-1/2})$ and $\unicode[STIX]{x1D6FF}\leqslant O(Ma^{-2/7})$ , so the eventual bound on $Nu$ cannot grow more slowly than $O(Ma^{2/7})$ . The issue remains when we let $\unicode[STIX]{x1D6FE}$ increase with $Ma$ to mimic the behaviour of the numerically optimal profiles (cf. figure 6 d), because the apparent gain in (6.15) is exactly outbalanced by the need of testing wavenumbers up to $k_{c}=O(\unicode[STIX]{x1D6FE}Ma^{1/2})$ . We therefore expect that to improve Hagstrom & Doering’s scaling using (6.12) will require careful estimates of $a_{k}$ , $b_{k}$ and $c_{k}$ at large wavenumbers, perhaps in conjunction with a more sophisticated choice of background field.

6.2 Reaching the asymptotic regime: current challenges for conic optimisation

As mentioned at the beginning of § 6.1, the true asymptotic nature of our numerical bound remains uncertain due to the limited range of Marangoni numbers that could be studied. Clearly, this kind of uncertainty is inherent to any kind of numerical investigation irrespective of which computational tools are employed. Nonetheless, the challenges faced by conic programming in reaching the asymptotic regime deserve further discussion.

The main limitation to extending the results presented in § 5.3 to larger values of $Ma$ is computational cost: proceeding from $Ma=10^{8}$ to $Ma=10^{9}$ took more than 48 h on our machine, and to achieve significant further progress would require computational resources beyond those available to the present investigation. One difficulty is that at large Marangoni numbers, checking whether a candidate background field satisfies the Fourier-transformed spectral constraints up to the cutoff wavenumber $k_{c}$ becomes a burden. For example, $k_{c}=20073$ at $Ma=10^{9}$ , meaning that $20073$ eigenvalue decompositions must be computed after each iteration of the wavenumber-tracking procedure described in § 5.2. The situation is worsened by the occurrence of bifurcations in critical wavenumbers, because more iterations are needed to correctly track all critical branches. Performance could be not improved by taking smaller steps in $Ma$ , because doing so slows progress towards higher Marangoni numbers. Increasing the parameter $r$ in (5.20) also does not help much, because the cost of adding more LMIs to our conic programmes at each iteration offsets the reduction in number of iterations required to identify the critical wavenumbers.

A possible solution to the critical wavenumber identification problem could be to apply the time-marching algorithm of Wen et al. (2013, 2015) to the optimality conditions for our conic programmes (5.17)–(5.18). This method has been reported to locate the correct critical wavenumbers efficiently, although convergence to the optimal background field can be slow (Wen et al. 2015). Fast but less accurate solvers for conic programmes (such as SCS by O’Donoghue et al. (2016)) may also have similar benefits and drawbacks, with the additional advantage that finely tuned open-source implementations are readily available. Irrespective of which method is utilised, once the critical wavenumbers have been identified, the optimal solution can be computed using accurate conic programming packages such as SDPT3 (used in this work).

Our numerical method and the possible improvements discussed above can of course be applied beyond Bénard–Marangoni convection. However, to study the asymptotic regime of more complex background method problems will require overcoming some additional obstacles. Spectral constraints with multiple test functions, such as those encountered in shear flows (Plasting & Kerswell 2003; Fantuzzi & Wynn 2016a ) or finite Prandtl number convection (Doering & Constantin 1996; Otero 2002), yield conic programmes with larger LMIs. While current state-of-the-art algorithms for conic programming can handle many small LMIs very efficiently, the computational cost of a single LMI grows as a nonlinear function of its size. This problem is exacerbated for problems with two- and higher-dimensional background fields that cannot be Fourier-transformed in the horizontal directions, because after discretisation one obtains a single LMI instead of a set of smaller, independent LMIs corresponding to each wavevector.

On the other hand, the (current) unfavourable scalability of algorithms for conic programming can be mitigated by taking advantage of special properties of the particular background field problem at hand. For instance, the spectral constraint often presents symmetries that can be exploited to reduce the number of degrees of freedom needed to discretise the background field or the test functions (however, this is not the case for Bénard–Marangoni convection). In addition, the very choice of discretisation method plays an important role because it directly impacts the sparsity of the eventual LMI. The piecewise-linear approximation method considered in this work is particularly attractive in this respect because it results in a chordal sparsity pattern, meaning that the non-zero entries of the LMI approximation of the spectral constraint can be represented by a chordal graph (a thorough discussion of these concepts is beyond the scope of this work, and we refer the interested reader to Fukuda et al. 2000, § 2). The same is true when one uses multidimensional piecewise-polynomial representations in the spirit of finite-element methods. Chordal sparsity enables one to decompose a large LMI into multiple smaller ones, at the expense of introducing extra optimisation variables (Fukuda et al. 2000; Nakata et al. 2003; Kim et al. 2011). This procedure can be automated, for instance using the MATLAB toolbox SparseCoLO (Fujisawa et al. 2009). As mentioned above, current algorithms for conic programming can handle multiple small LMIs much more efficiently than a single large one. Decomposition techniques based on chordal sparsity proved extremely effective in our study of Bénard–Marangoni convection and we expect the same to be true for other background method problems.

Finally, the development of efficient algorithms for large-scale conic programmes and implementations that take advantage of modern parallel computer architectures is a very active area of research (Sun, Andersen & Vandenberghe 2014; Madani, Kalbat & Lavaei 2015; Pakazad et al. 2015; O’Donoghue et al. 2016; Zheng et al. 2017a , b ). While it remains imperative to exploit all available symmetries and sparsity, advances at the algorithmic level promise to extend the ability of conic programming to reach the asymptotic regime of background method problems more complex than the one considered in this work.

7 Conclusion

This work studied the vertical heat transfer in Bénard–Marangoni convection of a fluid layer with infinite Prandtl number by means of rigorous upper bounds on the Nusselt number. First, the background method analysis by Hagstrom & Doering (2010) was extended to include balance parameters and formulate a new variational principle for the bound. Using this we proved that $Nu\leqslant 0.803\times Ma^{2/7}$ , reducing the prefactor of the previous best bound by approximately 4.2 %, but we also showed that optimising the balance parameters does not affect the asymptotic scaling of the optimal bounds compared to Hagstrom & Doering’s original formulation. We then employed conic programming to optimise the bound on $Nu$ over all background fields, as well as over two smaller families constrained by either a monotonicity or a convexity constraint. The main result of our numerical investigation was the observation that the fully optimal bounds have the form $Nu\leqslant O(Ma^{2/7}(\ln Ma)^{-1/2})$ for large Marangoni numbers. We also demonstrated that to achieve a logarithmic bound requires a background field with a non-monotonic boundary layer near the surface of the fluid.

Whether the logarithmic scaling observed numerically can be proved analytically remains an open question, and is the subject of ongoing research. The analysis presented in § 6 suggests a way forward by replacing the spectral constraint on the background field with the sufficient condition (6.12). Using (6.12) is an attractive option because it is easier to check than the spectral constraint for a candidate background field, and the role of non-monotonicity is apparent. Moreover, the fact that enforcing (6.12) at large wavenumbers seems to constrain the bound on $Nu$ is reminiscent of the bifurcations in critical wavenumbers observed in our numerical investigation (cf. figure 7). In summary, condition (6.12) seems to capture the essential features of the spectral constraint.

Should (6.12) prove too strong, the analysis of the energy stability problem (Fantuzzi & Wynn 2017) may be adapted to derive an inequality that exactly enforces each Fourier-transformed spectral constraint. The disadvantage is that such an inequality may not be analytically tractable except for very simple choices of the background field. On the other hand, it may be possible to check this condition numerically and confirm that a candidate background field can indeed achieve a logarithmic bound, leaving ‘only’ the task of constructing the correct estimates to prove so rigorously. Alternatively, one may consider the Lagrangian dual of the variational problem obtained with background method. This amounts to constructing the temperature and velocity fields that maximise the heat transfer subject to the linearised momentum equation, the boundary conditions, and suitably averaged versions of the advection–diffusion equation for the temperature (for a detailed discussion of the duality between these two approaches in the context of Rayleigh–Bénard convection, we refer the reader to Plasting & Ierley (2005)). However, only the fields achieving the maximal heat transfer yield a fully rigorous bound, so the maximisation must be solved exactly. Moreover, compared to the Rayleigh–Bénard problem the construction of a suitable hierarchy of nested boundary layers using Busse’s ‘multi- $\unicode[STIX]{x1D6FC}$ ’ solution method (see for example Busse 1979) is complicated by the lack of vertical symmetry and the Neumann conditions at the surface of the fluid.

Irrespective of how the variational problem for the upper bound on $Nu$ is analysed, however, it is evident from the present numerical investigation that the background method for the temperature field (or its dual formulation) cannot close the gap with the phenomenological prediction $Nu=O(Ma^{2/9})$ by Boeck & Thess (2001). It is possible that Boeck & Thess’s assumption that steady convection rolls remain stable as $Ma\rightarrow \infty$ is incorrect, making a scaling exponent of 2/9 unattainable with any bounding method. To prove so rigorously requires a lower bound on $Nu$ that grows faster than $Ma^{2/9}$ , which can also not be achieved with the background method because the unstable conduction solution saturates the constant lower bound $Nu\geqslant 1$ . Consequently, further numerical simulations in the high- $Ma$ seem essential to investigate the issue. The observation of steady convection rolls would provide further supporting evidence for Boeck & Thess’s phenomenological prediction. Determining the stability of the steady rolls is of interest also to reveal if the bifurcations in critical wavenumbers observed in our computations correspond to yet unobserved physical instabilities.

If Boeck & Thess’s phenomenological prediction is corroborated by further DNSs, to confirm it through rigorous bounds on $Nu$ will necessarily require bounding techniques beyond the background method. Unfortunately, the formulation of a wall-to-wall optimal transport problem in the spirit of Hassanzadeh, Chini & Doering (2014) and Tobasco & Doering (2017) does not appear suited to the study Bénard–Marangoni convection at infinite- $Pr$ . In fact, the optimal transport approach treats the temperature as a passively advected and diffusing scalar, and one looks for the (generally time-dependent) incompressible velocity field that maximises the passive vertical transport of heat subject to a maximum power budget. However, in infinite- $Pr$ Bénard–Marangoni convection the flow velocity is a linear function of the temperature field, which is effectively the only dynamical variable. This coupling is crucial in the background method analysis, so improving our bound on $Nu$ without taking it into account seems unlikely.

It would then be tempting to formulate the ‘ultimate’ optimal wall-to-wall transport problem using the temperature as the decision variable, and let the flow velocity be specified as a function of it. However, this corresponds to searching for the exact solution of the equations of motion (2.1a–c ) with maximal heat transfer, so any significant progress does not appear possible. Difficulties remain when one drops the time dependence: maximising the heat transfer among the steady solutions is not much easier, and in any case the eventual bound would rely on the unproven assumption that unsteady flows cannot transport more heat than steady ones. Nonetheless, the construction of $Ma$ -dependent exact solutions remains of interest because knowledge of a (possibly unstable) flow with Nusselt number $Nu_{ss}$ places a strict limit on what can be achieved by upper-bounding theory. In particular, any bounds that apply equally to all solutions of (2.1a–c ) cannot be better than $Nu_{ss}$ . Moreover, any flow with heat transfer $Nu_{ss}\gg O(Ma^{2/9})$ would demonstrate that Boeck & Thess’s phenomenological scaling applies at most to a particular subset of all possible convective flows.

While improving the rigorous upper bound on $Nu$ using the ‘ultimate’ wall-to-wall optimal transport approach described above appears challenging, it may be possible to consider successively weaker, tractable relaxations of it. The idea stems from the aforementioned realisation that the background method analysis is dual to the problem of maximising the heat transfer over all temperature (and associated velocity) fields that satisfy a set of constraints obtained by averaging the heat equation (Plasting & Ierley 2005). The upper bound on $Nu$ may therefore be improved by including additional constraints implied by the heat equation, but not the heat equation itself. A simple way to do so is through a general bounding framework that encompasses the background method (Chernyshenko et al. 2014; Chernyshenko 2017). The essence of this approach is to construct a functional ${\mathcal{V}}$ of the flow variables subject to a positivity condition akin to the spectral constraint in the background method. Each term in this functional can be interpreted as enforcing a particular constraint implied by the governing equations. Taking ${\mathcal{V}}$ to be the volume average of a quadratic polynomial of the flow variables gives the same bound as the background method (Chernyshenko 2017), but experience with finite-dimensional systems (Fantuzzi et al. 2016; Goluskin 2016) indicates that considering more general functionals – for instance, volume averages of higher-than-quadratic polynomials of the flow variables – could yield significant improvements. Although the construction of suitable functionals may be beyond the reach of purely analytical work, progress can be assisted by computations that utilise conic programming techniques similar to those applied in this paper. Whether the numerical bounds can reach the asymptotic regime is of course highly dependent on the availability of efficient algorithmic tools for conic programming. Promising recent developments in this field (see for example O’Donoghue et al. 2016; Zheng et al. 2017a , b ), however, give us hope that Bénard–Marangoni convection and other turbulent hydrodynamic systems may be studied successfully in the near future.