1 Introduction
When a layer of fluid heated from below is subject to temperature gradients along its surface, local variations in the surface tension generate a shear stress. This phenomenon, called the Marangoni effect, can set the fluid in motion when the ratio of surface tension forces to viscosity is sufficiently large. The ensuing flow, known as Bénard–Marangoni convection, can produce beautiful surface patterns, as famously observed by Bénard (Reference Bénard1901), and is a paradigm for pattern formation. It also underpins a number of industrial processes, such as fusion welding (DebRoy & David Reference DebRoy and David1995) and the growth of semiconductors (Lappa Reference Lappa2010). Nevertheless, Bénard–Marangoni convection remains poorly understood, especially when compared to its buoyancy-driven counterpart, Rayleigh–Bénard convection.
A fundamental open problem is to determine the vertical heat transport as a function of the thermal forcing and the material parameters of the fluid. In non-dimensional terms, one is interested in how the Nusselt number 
$Nu$ varies with the Marangoni number 
$Ma$, which measures the relative strength of thermally driven surface tension to viscous forces, and the Prandtl number 
$Pr$, given by the ratio between the kinematic viscosity and the thermal diffusivity of the fluid.
For finite Prandtl numbers, a phenomenological argument by Pumir & Blumenfeld (Reference Pumir and Blumenfeld1996) predicts 
$Nu\sim Ma^{1/3}$ with a Prandtl-dependent prefactor when 
$Ma\gg 1$ and the flow is turbulent. Two-dimensional direct numerical simulations (DNS) with stress-free boundaries at low 
$Pr$ support this scaling (Boeck & Thess Reference Boeck and Thess1998), but no-slip boundaries in either two or three dimensions yield smaller powers of 
$Ma$ (Boeck Reference Boeck2005). Two-dimensional free-slip DNS at both high and infinite 
$Pr$ also suggest a smaller exponent. Assuming steady convection rolls are stable at arbitrarily large 
$Ma$, a boundary-layer scaling analysis predicts 
$Nu\sim Ma^{2/9}$ in the infinite-
$Pr$ limit (Boeck & Thess Reference Boeck and Thess2001).
Rigorous results, derived directly from the governing equations without introducing unproven assumptions, are key to substantiate or rule out any of these heuristic scaling arguments. By expressing the temperature field in terms of fluctuations around a carefully chosen steady ‘background’ temperature field, Hagstrom & Doering (Reference Hagstrom and Doering2010) proved that 
$Nu\lesssim Ma^{1/2}$ uniformly in 
$Pr$ when this is finite, and 
$Nu\lesssim Ma^{2/7}$ for 
$Pr=\infty$. These bounds are consistent with all aforementioned theories, but the question remains of whether they are sharp – meaning there exist convective flows that saturate them – or can be improved.
Recently, numerical optimisation of the background temperature field for 
$Ma\leqslant 10^{9}$ suggested that Hagstrom & Doering’s bound for the infinite-
$Pr$ case can be improved at least by a logarithm (Fantuzzi, Pershin & Wynn Reference Fantuzzi, Pershin and Wynn2018). Precisely, the best bound available to the ‘background method’ for 
$Ma\gg 1$ appears to be 
$Nu\lesssim Ma^{2/7}(\ln Ma)^{-1/2}$, although the power of the logarithm remains uncertain due to the limited range of 
$Ma$ spanned by the numerical data. In this work, we prove analytically that logarithmic improvements to a power-law bound with exponent 
$2/7$ are indeed possible. Specifically, we show that
$$\begin{eqnarray}Nu\lesssim Ma^{2/7}(\ln Ma)^{-1/7},\quad \text{when }Ma\gg 1.\end{eqnarray}$$ We do this by combining the careful construction of an asymmetric background temperature field, inspired by the optimal profiles from Fantuzzi et al. (Reference Fantuzzi, Pershin and Wynn2018), with new estimates for the coupling between temperature and vertical velocity. These differ fundamentally from the estimates that apply to infinite-
$Pr$ Rayleigh–Bénard convection (Doering, Otto & Reznikoff Reference Doering, Otto and Reznikoff2006; Whitehead & Doering Reference Whitehead and Doering2011; Whitehead & Wittenberg Reference Whitehead and Wittenberg2014) due to the different boundary conditions (BC) for the velocity field.
2 The model
We consider a 
$d$-dimensional layer of fluid (
$d=2$ or 3) in a box domain with horizontal coordinates 
$\boldsymbol{x}\in \unicode[STIX]{x1D6F1}_{i=1}^{d-1}[0,L_{i}]$ and vertical coordinate 
$z\in [0,1]$. In the infinite-
$Pr$ limit, Pearson’s equations for Bénard–Marangoni convection (Pearson Reference Pearson1958) become
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}T+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T-\unicode[STIX]{x0394}T=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & -\unicode[STIX]{x0394}\boldsymbol{u}+\unicode[STIX]{x1D735}p=0. & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}(\boldsymbol{x},z,t)=(\boldsymbol{v}(\boldsymbol{x},z,t),w(\boldsymbol{x},z,t))$ is the velocity vector field with horizontal and vertical components 
$\boldsymbol{v}$ and 
$w$, respectively, 
$T(\boldsymbol{x},z,t)$ is the scalar temperature and 
$p(\boldsymbol{x},z,t)$ is the scalar pressure. We assume that all variables are periodic in the horizontal directions, while 
$$\begin{eqnarray}\displaystyle & T|_{z=0}=0,\quad \unicode[STIX]{x2202}_{z}T|_{z=1}=-1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{u}|_{z=0}=0,\quad w|_{z=1}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & [\unicode[STIX]{x2202}_{z}\boldsymbol{v}+Ma\,\unicode[STIX]{x1D735}_{\boldsymbol{x}}T]_{z=1}=0, & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D735}_{\boldsymbol{x}}$ denotes the horizontal gradient. The steady solution 
$\boldsymbol{u}=0$, 
$T=-z$, 
$p=\text{const.}$ corresponds to a purely conductive state; it is globally asymptotically stable for 
$Ma\leqslant 66.84$ (Fantuzzi & Wynn Reference Fantuzzi and Wynn2017) and linearly stable for 
$Ma\leqslant 79.61$ (Pearson Reference Pearson1958). For larger Marangoni numbers convection ensues, and the velocity field can be completely slaved to the temperature. Precisely, let 
${\hat{w}}_{\boldsymbol{k}}$ and 
$\hat{T}_{\boldsymbol{k}}$ be any Fourier modes of the vertical velocity and temperature, respectively, with horizontal wavevector 
$\boldsymbol{k}$ of magnitude 
$k$. (These are unique when 
$d=2$ but not when 
$d=3$.) One finds (Hagstrom & Doering Reference Hagstrom and Doering2010)
$$\begin{eqnarray}{\hat{w}}_{\boldsymbol{k}}(z)=-Ma\,f_{k}(z)\hat{T}_{\boldsymbol{k}}(1),\end{eqnarray}$$ where, setting 
$h(x):=(\sinh x)/x$ for convenience,
$$\begin{eqnarray}f_{k}(z)=\frac{1}{2}k^{2}z(z-1)\left\{\frac{h(k)h(kz)-h[k(1-z)]}{h(2k)-1}\right\}.\end{eqnarray}$$ Key to proving (1.1) are the following new bounds for the temperature–velocity coupling in (2.3). They are proved in appendix A and hold for any fixed 
$0\leqslant \unicode[STIX]{x1D6FD}<1$ and 
$k\geqslant 0$. First, for 
$0\leqslant z\leqslant \unicode[STIX]{x1D6FD}$ we have
$$\begin{eqnarray}\left|\,f_{k}(z)\right|\leqslant \frac{1}{6}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FD},k)z^{2},\quad \unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FD},k):=k^{4}\frac{h(k)h(k\unicode[STIX]{x1D6FD})}{h(2k)-1}.\end{eqnarray}$$ Further, for 
$\unicode[STIX]{x1D6FD}\leqslant z\leqslant 1$ we can bound
$$\begin{eqnarray}\frac{\left|\,f_{k}(\unicode[STIX]{x1D6FD})\right|}{1-\unicode[STIX]{x1D6FD}}(1-z)\leqslant \left|\,f_{k}(z)\right|\leqslant \frac{k}{2}(1-z)\text{e}^{-k(1-z)}.\end{eqnarray}$$3 Bound on the Nusselt number
Denote the horizontal and long-time average of a quantity 
$q(\boldsymbol{x},z,t)$ by
$$\begin{eqnarray}\langle q\rangle (z)=\limsup _{{\mathcal{T}}\rightarrow \infty }\frac{1}{{\mathcal{T}}L_{1}\cdots L_{d-1}}\int _{0}^{{\mathcal{T}}}\int _{0}^{L_{1}}\cdots \int _{0}^{L_{d-1}}q(\boldsymbol{x},z,t)\,\text{d}\boldsymbol{x}\,\text{d}t.\end{eqnarray}$$Our interest is to derive a Marangoni-dependent upper bound on the Nusselt number, i.e., the ratio of the total vertical heat flux to the purely conductive one:
$$\begin{eqnarray}Nu:=\frac{\displaystyle \int _{0}^{1}\langle wT-\unicode[STIX]{x2202}_{z}T\rangle \,\text{d}z}{\displaystyle \int _{0}^{1}\langle -\unicode[STIX]{x2202}_{z}T\rangle \,\text{d}z}.\end{eqnarray}$$ To bound 
$Nu$, we follow Hagstrom & Doering (Reference Hagstrom and Doering2010) and write the temperature field as the sum of a steady background field 
$\unicode[STIX]{x1D70F}(z)$, which satisfies the inhomogeneous BC in (2.2a) but is otherwise arbitrary, and a fluctuation 
$\unicode[STIX]{x1D703}(\boldsymbol{x},z,t)$ satisfying
$$\begin{eqnarray}\displaystyle & T(\boldsymbol{x},z,t)=\unicode[STIX]{x1D70F}(z)+\unicode[STIX]{x1D703}(\boldsymbol{x},z,t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70F}(0)=0,\quad \unicode[STIX]{x1D70F}^{\prime }(1)=-1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D703}|_{z=0}=0,\quad \unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}|_{z=1}=0. & \displaystyle\end{eqnarray}$$
$z$. It is shown by Hagstrom & Doering (Reference Hagstrom and Doering2010) that 
$$\begin{eqnarray}Nu^{-1}=\int _{0}^{1}\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\right|^{2}+2\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}\rangle \,\text{d}z-\Vert \unicode[STIX]{x1D70F}^{\prime }\Vert _{2}^{2}-2\unicode[STIX]{x1D70F}(1),\end{eqnarray}$$ where 
$\Vert \boldsymbol{\cdot }\Vert _{2}$ denotes the usual 
$L^{2}$-norm. At this stage, suppose that 
$\unicode[STIX]{x1D70F}$ is chosen such that
$$\begin{eqnarray}{\mathcal{Q}}^{\unicode[STIX]{x1D70F}}\{\unicode[STIX]{x1D703}\}:=\int _{0}^{1}\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\right|^{2}+2\unicode[STIX]{x1D70F}^{\prime }w\unicode[STIX]{x1D703}\rangle \,\text{d}z\geqslant 0\end{eqnarray}$$ for all time-independent trial fields 
$\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}(\boldsymbol{x},z)$ that are horizontally periodic and satisfy (3.3c), with 
$w=w(\boldsymbol{x},z)$ being a function of 
$\unicode[STIX]{x1D703}$ defined in Fourier space according to (2.3). This can be interpreted as a nonlinear stability condition for 
$\unicode[STIX]{x1D70F}$ as if it were a solution to (2.1a)–(2.1c) (see, for example, Malkus (Reference Malkus1954)). Then,
$$\begin{eqnarray}Nu^{-1}\geqslant -\Vert \unicode[STIX]{x1D70F}^{\prime }\Vert _{2}^{2}-2\unicode[STIX]{x1D70F}(1)=1-\Vert \unicode[STIX]{x1D70F}^{\prime }+1\Vert _{2}^{2},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D70F}(0)=0$ is used to obtain the second equality. If the right-hand side is positive, inverting this lower bound produces a finite upper bound on 
$Nu$. A background field 
$\unicode[STIX]{x1D70F}$ is now constructed which gives (1.1) when 
$Ma\gg 1$.
4 Proof of the main result
The boundary condition 
$\unicode[STIX]{x1D70F}(0)=0$ can be dropped because 
$\unicode[STIX]{x1D70F}$ can always be shifted by a constant without affecting (3.5) and (3.6), which depend only on 
$\unicode[STIX]{x1D70F}^{\prime }$. Moreover, the boundary condition 
$\unicode[STIX]{x1D70F}^{\prime }(1)=-1$ can formally be ignored because it can be enforced at the end by modifying 
$\unicode[STIX]{x1D70F}^{\prime }$ in a infinitesimally thin layer near 
$1$ without affecting our bound on 
$Nu$. Given these observations, and motivated by the numerically optimal profiles computed by Fantuzzi et al. (Reference Fantuzzi, Pershin and Wynn2018, see figure 4), we choose
$$\begin{eqnarray}\unicode[STIX]{x1D70F}^{\prime }(z):=\left\{\begin{array}{@{}ll@{}}\displaystyle -1+\left(\frac{z}{\unicode[STIX]{x1D6FF}}\right)^{1/s}=:\unicode[STIX]{x1D702}(z),\quad & \text{for }0\leqslant z\leqslant \unicode[STIX]{x1D6FF},\\ \unicode[STIX]{x1D709}(z),\quad & \text{for }\unicode[STIX]{x1D6FF}\leqslant z\leqslant 1,\end{array}\right.\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}<{\textstyle \frac{1}{2}}$, 
$s>0$, and 
$\unicode[STIX]{x1D709}(z)$ is a non-negative function to be specified later. With 
$\unicode[STIX]{x1D709}(z)=0$ and 
$s\rightarrow 0$ this choice yields the piecewise constant profiles already studied by Hagstrom & Doering (Reference Hagstrom and Doering2010) and Fantuzzi et al. (Reference Fantuzzi, Pershin and Wynn2018).
By expanding 
$\unicode[STIX]{x1D703}$ and 
$w$ as Fourier series in the horizontal directions, using (2.3), and noting that 
$\left|\,f_{k}(z)\right|=-f_{k}(z)$ for 
$0\leqslant z\leqslant 1$, it can be shown (Hagstrom & Doering Reference Hagstrom and Doering2010) that the marginal stability condition (3.5) holds if and only if the quadratic form
$$\begin{eqnarray}{\mathcal{Q}}_{\boldsymbol{k}}^{\unicode[STIX]{x1D70F}}\{\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{ k}}\}:=\int _{0}^{1}[|\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{ k}}^{\prime }(z)|^{2}+k^{2}|\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{ k}}(z)|^{2}+2Ma\,\unicode[STIX]{x1D70F}^{\prime }(z)\,\left|\,f_{k}(z)\right|\,\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(z)\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)]\,\text{d}z\end{eqnarray}$$ is non-negative for all 
$k>0$ and all real-valued functions 
$\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(z)$ subject to
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(0)=0,\quad \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }(1)=0.\end{eqnarray}$$ Since 
${\mathcal{Q}}_{\boldsymbol{k}}^{\unicode[STIX]{x1D70F}}\{\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}\}$ is homogeneous, we may assume without loss of generality that 
$\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)\geqslant 0$.
Using (4.1) and dropping the non-negative term 
$k^{2}|\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(z)|^{2}$ we obtain
$$\begin{eqnarray}{\mathcal{Q}}_{\boldsymbol{k}}^{\unicode[STIX]{x1D70F}}\{\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{ k}}\}\geqslant \Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}^{2}+2Ma\,\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{ k}}(1)\int _{0}^{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D702}(z)|f_{k}(z)|\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(z)\,\text{d}z+2Ma\,\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)\int _{\unicode[STIX]{x1D6FF}}^{1}\unicode[STIX]{x1D709}(z)|f_{k}(z)|\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(z)\,\text{d}z.\end{eqnarray}$$The fundamental theorem of calculus, the BC (4.3) and the Cauchy–Schwarz inequality imply
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(z)=\int _{0}^{z}\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}\leqslant \Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}\sqrt{z}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(z)=\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)-\int _{z}^{1}\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }(\unicode[STIX]{x1D701})\,d\unicode[STIX]{x1D701}\geqslant \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)-\Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}\sqrt{1-z}. & \displaystyle\end{eqnarray}$$
$\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)$ and the function 
$\unicode[STIX]{x1D709}(z)$ are non-negative by assumption, we can use these inequalities to bound 
$$\begin{eqnarray}\displaystyle {\mathcal{Q}}_{\boldsymbol{k}}^{\unicode[STIX]{x1D70F}}\{\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}\} & {\geqslant} & \displaystyle \Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}^{2}+2Ma\,I_{0}(\unicode[STIX]{x1D709},k)\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)^{2}\nonumber\\ \displaystyle & & \displaystyle -\,2Ma\,\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)\Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}\left[\int _{0}^{\unicode[STIX]{x1D6FF}}|\unicode[STIX]{x1D702}(z)f_{k}(z)|\sqrt{z}\,\text{d}z+I_{1/2}(\unicode[STIX]{x1D709},k)\right],\end{eqnarray}$$where we have introduced the notation
$$\begin{eqnarray}I_{\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D709},k)=\int _{\unicode[STIX]{x1D6FF}}^{1}\unicode[STIX]{x1D709}(z)|f_{k}(z)|(1-z)^{\unicode[STIX]{x1D6FD}}\,\text{d}z.\end{eqnarray}$$ Let us now estimate the terms inside the square brackets in (4.6). For the integral over 
$(0,\unicode[STIX]{x1D6FF})$, we use estimate (2.5) with 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FF}$ and the definition of 
$\unicode[STIX]{x1D702}(z)$ from (4.1) to obtain
$$\begin{eqnarray}\displaystyle \int _{0}^{\unicode[STIX]{x1D6FF}}\left|\unicode[STIX]{x1D702}(z)f_{k}(z)\right|\sqrt{z}\,\text{d}z & {\leqslant} & \displaystyle \int _{0}^{\unicode[STIX]{x1D6FF}}\left|-1+\left(\frac{z}{\unicode[STIX]{x1D6FF}}\right)^{1/s}\right|{\displaystyle \frac{1}{6}}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)z^{2}\sqrt{z}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{6}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)\int _{0}^{\unicode[STIX]{x1D6FF}}(z^{5/2}-\unicode[STIX]{x1D6FF}^{-1/s}z^{(1/s)+(5/2)})\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle \frac{2}{21(2+7s)}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)\unicode[STIX]{x1D6FF}^{7/2}.\end{eqnarray}$$ To bound 
$I_{1/2}(\unicode[STIX]{x1D709},k)$, instead, we use the Cauchy–Schwarz inequality:
$$\begin{eqnarray}\displaystyle I_{1/2}(\unicode[STIX]{x1D709},k) & = & \displaystyle \int _{\unicode[STIX]{x1D6FF}}^{1}\unicode[STIX]{x1D709}(z)|f_{k}(z)|\sqrt{1-z}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle \int _{\unicode[STIX]{x1D6FF}}^{1}\sqrt{\unicode[STIX]{x1D709}(z)|f_{k}(z)|}\,\sqrt{\unicode[STIX]{x1D709}(z)|f_{k}(z)|(1-z)}\,\text{d}z\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \left(\int _{\unicode[STIX]{x1D6FF}}^{1}\unicode[STIX]{x1D709}(z)|f_{k}(z)|\,\text{d}z\right)^{1/2}\left(\int _{\unicode[STIX]{x1D6FF}}^{1}\unicode[STIX]{x1D709}(z)|f_{k}(z)|(1-z)\,\text{d}z\right)^{1/2}\nonumber\\ \displaystyle & = & \displaystyle \sqrt{I_{0}(\unicode[STIX]{x1D709},k)}\,\sqrt{I_{1}(\unicode[STIX]{x1D709},k)}.\end{eqnarray}$$Substituting these two estimates into (4.6) we arrive at
$$\begin{eqnarray}\displaystyle {\mathcal{Q}}_{\boldsymbol{k}}^{\unicode[STIX]{x1D70F}}\{\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}\} & {\geqslant} & \displaystyle \Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}^{2}+2Ma\,I_{0}(\unicode[STIX]{x1D709},k)\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)^{2}\nonumber\\ \displaystyle & & \displaystyle -\,2Ma\,\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)\Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}\left[\frac{2\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)\,\unicode[STIX]{x1D6FF}^{7/2}}{21(2+7s)}+\sqrt{I_{0}(\unicode[STIX]{x1D709},k)}\,\sqrt{I_{1}(\unicode[STIX]{x1D709},k)}\right].\end{eqnarray}$$ The right-hand side of this estimate is a quadratic form of type 
$ax^{2}-2bxy+cy^{2}$ with 
$x=\Vert \hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}^{\prime }\Vert _{2}$ and 
$y=\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}(1)$. Quadratic forms are non-negative when their discriminant is negative, meaning 
$\left|b\right|\leqslant \sqrt{ac}$, so 
${\mathcal{Q}}_{\boldsymbol{k}}^{\unicode[STIX]{x1D70F}}\{\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}\}\geqslant 0$ for all admissible fields 
$\hat{\unicode[STIX]{x1D703}}_{\boldsymbol{k}}$ if
$$\begin{eqnarray}Ma\left(\frac{2\,\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)\,\unicode[STIX]{x1D6FF}^{7/2}}{21(2+7s)}+\sqrt{I_{0}(\unicode[STIX]{x1D709},k)}\,\sqrt{I_{1}(\unicode[STIX]{x1D709},k)}\right)\leqslant \sqrt{2Ma\,I_{0}(\unicode[STIX]{x1D709},k)}.\end{eqnarray}$$For simplicity, we rewrite this condition as
$$\begin{eqnarray}\frac{2\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)\unicode[STIX]{x1D6FF}^{7/2}}{21(2+7s)\sqrt{I_{0}(\unicode[STIX]{x1D709},k)}}+\sqrt{I_{1}(\unicode[STIX]{x1D709},k)}\leqslant \sqrt{\frac{2}{Ma}}.\end{eqnarray}$$ To prove a bound on the Nusselt number 
$Nu$ we require inequality (4.12) to hold for all 
$k>0$. A sufficient condition for this is that 
$\unicode[STIX]{x1D709}(z)$ and 
$\unicode[STIX]{x1D6FF}$ be chosen such that, for some constant 
$c\in (0,1)$,
$$\begin{eqnarray}\displaystyle & \displaystyle \sup _{k>0}\sqrt{I_{1}(\unicode[STIX]{x1D709},k)}\leqslant (1-c)\sqrt{\frac{2}{Ma}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \sup _{k>0}\frac{2\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)\unicode[STIX]{x1D6FF}^{7/2}}{21(2+7s)\sqrt{I_{0}(\unicode[STIX]{x1D709},k)}}\leqslant c\sqrt{\frac{2}{Ma}}. & \displaystyle\end{eqnarray}$$Equivalently, after squaring both sides of each condition and rearranging,
$$\begin{eqnarray}\displaystyle & \displaystyle \sup _{k>0}I_{1}(\unicode[STIX]{x1D709},k)\leqslant \frac{2(1-c)^{2}}{Ma}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}^{7}\leqslant \frac{441(2+7s)^{2}c^{2}}{2Ma}\times \inf _{k>0}\frac{I_{0}(\unicode[STIX]{x1D709},k)}{\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)^{2}}. & \displaystyle\end{eqnarray}$$ We will now show that (4.15a) and (4.15b) can be satisfied by a suitable choice of 
$\unicode[STIX]{x1D709}(z)$. Inspired by the numerically optimal background fields in Fantuzzi et al. (Reference Fantuzzi, Pershin and Wynn2018, figure 4) we consider
$$\begin{eqnarray}\unicode[STIX]{x1D709}(z):=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}}{(1-z)^{2}},\quad & \text{for }1-\unicode[STIX]{x1D6FE}\leqslant z\leqslant 1-\unicode[STIX]{x1D700},\\ 0,\quad & \text{otherwise}.\end{array}\right.\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D700}$, 
$\unicode[STIX]{x1D6FE}$ and 
$\unicode[STIX]{x1D714}$ are strictly positive parameters, to be determined as a function of the Marangoni number 
$Ma$ subject to the constraint 
$\unicode[STIX]{x1D700}<\unicode[STIX]{x1D6FE}\leqslant {\textstyle \frac{1}{2}}$.
Upon combining this choice with the upper bound on 
$\left|\,f_{k}\right|$ in (2.6) and the elementary inequality 
$\text{e}^{-k\unicode[STIX]{x1D700}}-\text{e}^{-k\unicode[STIX]{x1D6FE}}\leqslant 1$ we can estimate
$$\begin{eqnarray}I_{1}(\unicode[STIX]{x1D709},k)\leqslant \frac{k\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}}{2}\int _{1-\unicode[STIX]{x1D6FE}}^{1-\unicode[STIX]{x1D700}}\text{e}^{-k(1-z)}\,\text{d}z=\frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}}{2}(\text{e}^{-k\unicode[STIX]{x1D700}}-\text{e}^{-k\unicode[STIX]{x1D6FE}})\leqslant \frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}}{2}.\end{eqnarray}$$ This estimate holds for all 
$k$, so we can bound the left-hand side of (4.15a) from above as
$$\begin{eqnarray}\sup _{k>0}I_{1}(\unicode[STIX]{x1D709},k)\leqslant \frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}}{2}.\end{eqnarray}$$ To estimate the right-hand side of (4.15b) from below, instead, observe that the lower bound on 
$\left|\,f_{k}\right|$ in (2.6) with 
$\unicode[STIX]{x1D6FD}=1-\unicode[STIX]{x1D6FE}$ implies
$$\begin{eqnarray}I_{0}(\unicode[STIX]{x1D709},k)=\int _{1-\unicode[STIX]{x1D6FE}}^{1-\unicode[STIX]{x1D700}}\left|\,f_{k}(z)\right|\frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}}{(1-z)^{2}}\,\text{d}z\geqslant \int _{1-\unicode[STIX]{x1D6FE}}^{1-\unicode[STIX]{x1D700}}\frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}\left|\,f_{k}(1-\unicode[STIX]{x1D6FE})\right|}{\unicode[STIX]{x1D6FE}(1-z)}\,\text{d}z=\frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}}{\unicode[STIX]{x1D6FE}}\left|\,f_{k}(1-\unicode[STIX]{x1D6FE})\right|\ln \left(\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D700}}\right).\end{eqnarray}$$Thus,
$$\begin{eqnarray}\inf _{k>0}\frac{I_{0}(\unicode[STIX]{x1D709},k)}{\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)^{2}}\geqslant \inf _{k>0}\frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}\left|\,f_{k}(1-\unicode[STIX]{x1D6FE})\right|}{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)^{2}}\ln \left(\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D700}}\right).\end{eqnarray}$$ After substituting the expressions for 
$\left|\,f_{k}(1-\unicode[STIX]{x1D6FE})\right|=-f_{k}(1-\unicode[STIX]{x1D6FE})$ and 
$\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},k)$ from (2.4) and (2.5) into the right-hand side of the last inequality and rearranging, we conclude from (4.18) and (4.20) that conditions (4.15a) and (4.15b) hold, respectively, if
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}\leqslant \frac{4(1-c)^{2}}{Ma}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}^{7}\leqslant \frac{441}{4Ma}(1-\unicode[STIX]{x1D6FE})\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}(2+7s)^{2}c^{2}\ln \left(\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D700}}\right)\unicode[STIX]{x1D711}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF}):=\inf _{k>0}\frac{\{h(k)h[k(1-\unicode[STIX]{x1D6FE})]-h(k\unicode[STIX]{x1D6FE})\}[h(2k)-1]}{k^{6}h(k)^{2}h(k\unicode[STIX]{x1D6FF})^{2}}.\end{eqnarray}$$ Observe that the right-hand side of (4.21b) is strictly positive because the function 
$z\mapsto h(z)$ is increasing, so for all 
$\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF}\in [0,{\textstyle \frac{1}{2}}]$ the quantity 
$\unicode[STIX]{x1D711}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF})$ satisfies
$$\begin{eqnarray}0<\unicode[STIX]{x1D711}\left({\textstyle \frac{1}{2}},{\textstyle \frac{1}{2}}\right)\leqslant \unicode[STIX]{x1D711}\left(\unicode[STIX]{x1D6FE},{\textstyle \frac{1}{2}}\right)\leqslant \unicode[STIX]{x1D711}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF})\leqslant \unicode[STIX]{x1D711}(0,0).\end{eqnarray}$$ The analysis we have just carried out shows that the background temperature field 
$\unicode[STIX]{x1D70F}(z)$ defined through (4.1) satisfies the marginal stability constraint (3.5) when 
$\unicode[STIX]{x1D709}(z)$ is as in (4.16), provided that (4.21a) and (4.21b) hold. Let us now turn the attention to the bound on the Nusselt number produced by 
$\unicode[STIX]{x1D70F}$. Substituting (4.1) and (4.16) into (3.6) gives
$$\begin{eqnarray}Nu^{-1}\geqslant \frac{2}{2+s}\unicode[STIX]{x1D6FF}-\frac{\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D700}}{3}\left(1-\frac{\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D6FE}}\right)\left(\frac{6}{\unicode[STIX]{x1D714}}+1+\frac{\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D6FE}}+\frac{\unicode[STIX]{x1D700}^{2}}{\unicode[STIX]{x1D6FE}^{2}}\right).\end{eqnarray}$$ Maximising the right-hand side of (4.24) over 
$\unicode[STIX]{x1D6FF}$, 
$\unicode[STIX]{x1D700}$, 
$\unicode[STIX]{x1D6FE}$, 
$s$, 
$\unicode[STIX]{x1D714}$ and 
$c$ subject to (4.21a), (4.21b) and the constraints 
$\unicode[STIX]{x1D6FF}<{\textstyle \frac{1}{2}}$, 
$\unicode[STIX]{x1D700}<\unicode[STIX]{x1D6FE}\leqslant {\textstyle \frac{1}{2}}$ and 
$0<c<1$ is hard analytically, but can be done numerically. The results, plotted in figure 1, strongly suggest that the optimal upper bound on 
$Nu$ provable via (4.24), (4.21a) and (4.21b) is proportional to 
$Ma^{2/7}(\ln Ma)^{-1/7}$ as 
$Ma\rightarrow \infty$, even though not all of the parameters 
$\unicode[STIX]{x1D6FF}$, 
$\unicode[STIX]{x1D700}$, 
$\unicode[STIX]{x1D6FE}$, 
$\unicode[STIX]{x1D714}$, 
$s$ and 
$c$ exhibit a simple scaling behaviour. Optimisation of these parameters in the limit of infinite Marangoni number is also not easy and will not be pursued in this work. Instead, we prove that
$$\begin{eqnarray}Nu\lesssim Ma^{2/7}(\ln Ma)^{-1/7},\quad \text{as }Ma\rightarrow \infty\end{eqnarray}$$ if we set either 
$\unicode[STIX]{x1D6FE}={\textstyle \frac{1}{2}}$ or 
$\unicode[STIX]{x1D6FE}=(\ln Ma)^{-1}$ (the latter gives a better prefactor) and
$$\begin{eqnarray}s={\textstyle \frac{2}{5}},\quad c={\textstyle \frac{1}{2}},\quad \unicode[STIX]{x1D700}=Ma^{-1/2},\quad \unicode[STIX]{x1D714}=1,\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}=\left[\left(\frac{126}{5Ma}\right)^{2}(1-\unicode[STIX]{x1D6FE})\ln (\unicode[STIX]{x1D6FE}Ma^{1/2})\unicode[STIX]{x1D711}\left(\unicode[STIX]{x1D6FE},\frac{1}{2}\right)\right]^{1/7}.\end{eqnarray}$$(a) Bounds on 
$Nu$ obtained with (4.24) for optimised 
$\unicode[STIX]{x1D6FF}$, 
$\unicode[STIX]{x1D700}$, 
$\unicode[STIX]{x1D6FE}$, 
$s$, 
$\unicode[STIX]{x1D714}$, 
$c$ (
$\cdots \cdots$, red), and with (4.29) for 
$\unicode[STIX]{x1D700}=Ma^{-1/2}$, 
$s={\textstyle \frac{2}{5}}$, 
$c={\textstyle \frac{1}{2}}$ and either 
$\unicode[STIX]{x1D6FE}={\textstyle \frac{1}{2}}$ (——, red) or 
$\unicode[STIX]{x1D6FE}=(\ln Ma)^{-1}$ (– – –, red). Also plotted are the analytical bound 
$Nu\leqslant 0.838Ma^{2/7}$ by Hagstrom & Doering (Reference Hagstrom and Doering2010) (– – –, black), the numerical bound by Fantuzzi et al. (Reference Fantuzzi, Pershin and Wynn2018) (——, black), and DNS data by Boeck & Thess (Reference Boeck and Thess2001) (
$\times$). (b,c) Optimised boundary layers of 
$\unicode[STIX]{x1D70F}^{\prime }$ for 
$Ma=10^{4}$. (d–i) Values of 
$\unicode[STIX]{x1D6FF}$, 
$\unicode[STIX]{x1D700}$, 
$\unicode[STIX]{x1D6FE}$, 
$\unicode[STIX]{x1D714}$, 
$s$ and 
$c$ that optimise (4.24) subject to (4.15a), (4.15b), 
$\unicode[STIX]{x1D6FF}<{\textstyle \frac{1}{2}}$, 
$\unicode[STIX]{x1D700}<\unicode[STIX]{x1D6FE}\leqslant {\textstyle \frac{1}{2}}$ and 
$0<c<1$, as a function of 
$Ma$.
First, for simplicity we strengthen (4.21b) by estimating 
$\unicode[STIX]{x1D711}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF})\geqslant \unicode[STIX]{x1D711}(\unicode[STIX]{x1D6FE},{\textstyle \frac{1}{2}})$, cf. (4.23). Then, it follows from (4.24) that 
$\unicode[STIX]{x1D6FF}$ should be taken as large as the resulting inequality allows. Upon insisting that
$$\begin{eqnarray}\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}^{2}=\frac{4(1-c)^{2}}{Ma}\end{eqnarray}$$ at all 
$Ma$, which is the case for the optimal parameters obtained numerically, we find
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}=\left[\frac{441}{Ma^{2}}c^{2}(1-c)^{2}(2+7s)^{2}(1-\unicode[STIX]{x1D6FE})\ln \left(\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D700}}\right)\unicode[STIX]{x1D711}\left(\unicode[STIX]{x1D6FE},\frac{1}{2}\right)\right]^{1/7}.\end{eqnarray}$$ Substituting this expression back into (4.24) and using (4.27) to eliminate 
$\unicode[STIX]{x1D714}$ yields
$$\begin{eqnarray}Nu^{-1}\geqslant A(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D700},c,s)-B(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D700},c),\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle & \displaystyle A(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D700},c,s):=\frac{2}{2+s}\left[\frac{441}{Ma^{2}}c^{2}(1\!-\!c)^{2}(2+7s)^{2}(1\!-\!\unicode[STIX]{x1D6FE})\ln \left(\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D700}}\right)\unicode[STIX]{x1D711}\left(\!\unicode[STIX]{x1D6FE},\frac{1}{2}\right)\!\right]^{1/7}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle B(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D700},c):=\frac{16(1-c)^{4}}{3Ma^{2}\unicode[STIX]{x1D700}^{3}}\left(1-\frac{\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D6FE}}\right)\left(\frac{3Ma\,\unicode[STIX]{x1D700}^{2}}{2(1-c)^{2}}+1+\frac{\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D6FE}}+\frac{\unicode[STIX]{x1D700}^{2}}{\unicode[STIX]{x1D6FE}^{2}}\right). & \displaystyle\end{eqnarray}$$ To proceed, we make two suboptimal but simple choices. First, to simplify the dependence of 
$B(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D700},c)$ on 
$Ma$ we set 
$\unicode[STIX]{x1D700}=Ma^{-1/2}$. This gives 
$\unicode[STIX]{x1D714}=4(1-c)^{2}$ by (4.27). Second, motivated by our computational results we assume that 
$\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D700}=\unicode[STIX]{x1D6FE}Ma^{1/2}\rightarrow \infty$ as 
$Ma$ tends to infinity. Then, 
$B(\unicode[STIX]{x1D6FE},Ma^{-1/2},c)$ decays to zero faster than 
$A(\unicode[STIX]{x1D6FE},Ma^{-1/2},c,s)$ as 
$Ma$ is raised and we conclude from (4.29) that, asymptotically, 
$Nu\leqslant 1/A(\unicode[STIX]{x1D6FE},Ma^{-1/2},c,s)$. Minimising this asymptotic bound over 
$s$ and 
$c$ simply requires maximising 
$A(\unicode[STIX]{x1D6FE},Ma^{-1/2},c,s)$. This is straightforward and yields 
$s={\textstyle \frac{2}{5}}$ and 
$c={\textstyle \frac{1}{2}}$, the same values approached by the optimal parameters in figure 1(h,i). With these values, equation (4.28) reduces to the value in (4.26e) and the asymptotic bound on 
$Nu$ becomes
$$\begin{eqnarray}Nu\leqslant \frac{6}{5}\left[\frac{126^{2}}{25}\frac{1}{Ma^{2}}(1-\unicode[STIX]{x1D6FE})\ln (\unicode[STIX]{x1D6FE}Ma^{1/2})\unicode[STIX]{x1D711}\left(\unicode[STIX]{x1D6FE},\frac{1}{2}\right)\right]^{-1/7},\quad \text{as }Ma\rightarrow \infty .\end{eqnarray}$$ Minimising this expression over 
$\unicode[STIX]{x1D6FE}$ is not possible analytically, but is also not necessary in order to prove (4.25). For instance, simply setting 
$\unicode[STIX]{x1D6FE}=\frac{1}{2}$ gives
$$\begin{eqnarray}Nu\leqslant \frac{6}{5}\left[\frac{126^{2}}{100}\unicode[STIX]{x1D711}\left(\frac{1}{2},\frac{1}{2}\right)\right]^{-1/7}\times \frac{Ma^{2/7}}{(\ln Ma)^{1/7}},\quad \text{as }Ma\rightarrow \infty .\end{eqnarray}$$ Moreover, in light of (4.23) the prefactor can be improved by letting 
$\unicode[STIX]{x1D6FE}\rightarrow 0$ as 
$Ma\rightarrow \infty$, which asymptotically optimises the term 
$(1-\unicode[STIX]{x1D6FE})\unicode[STIX]{x1D711}(\unicode[STIX]{x1D6FE},{\textstyle \frac{1}{2}})$ in (4.32). The decay of 
$\unicode[STIX]{x1D6FE}$ must be sufficiently slow to ensure that 
$\unicode[STIX]{x1D6FE}Ma^{1/2}\rightarrow \infty$, as assumed above. With 
$\unicode[STIX]{x1D6FE}=(\ln Ma)^{-1}$, for instance,
$$\begin{eqnarray}Nu\leqslant \frac{6}{5}\left[\frac{126^{2}}{50}\unicode[STIX]{x1D711}\left(0,\frac{1}{2}\right)\right]^{-1/7}\times \frac{Ma^{2/7}}{(\ln Ma)^{1/7}},\quad \text{as }Ma\rightarrow \infty .\end{eqnarray}$$ The exact bounds on 
$Nu$ obtained from (4.29) at finite 
$Ma$ for 
$\unicode[STIX]{x1D700}=Ma^{-1/2}$, 
$c={\textstyle \frac{1}{2}}$, 
$s={\textstyle \frac{2}{5}}$ and either 
$\unicode[STIX]{x1D6FE}={\textstyle \frac{1}{2}}$ or 
$\unicode[STIX]{x1D6FE}=(\ln Ma)^{-1}$ are plotted in figure 1(a).
5 Conclusion
In this paper we have derived a new rigorous bound for the Nusselt number in Pearson’s model of Bénard–Marangoni convection at infinite Prandtl number. Specifically, we have proved that 
$Nu\lesssim Ma^{2/7}(\ln Ma)^{-1/7}$ at asymptotically high 
$Ma$, thereby refining a pure power-law bound with exponent 
$2/7$ by Hagstrom & Doering (Reference Hagstrom and Doering2010). The quantitative improvement on this previous result is not large for realistic values of the Marangoni number, but our logarithmic correction is significant for two reasons.
First, its proof relies on a subtle balance between the width of the bottom boundary layer of our background temperature field, which drives the asymptotic scaling of 
$Nu$, and the stabilising effect – with respect to the marginal stability constraint (3.5) – of a thin layer near the fluid’s surface where the temperature increases. Qualitatively similar layers characterise the mean vertical temperature profiles observed in DNS by Boeck & Thess (Reference Boeck and Thess2001, figure 2) and their coupling underpins the phenomenological scaling theory proposed by those authors. It is therefore tempting to conjecture that the heat transport in physically realised flows indeed depends on a subtle interplay between the thermal boundary layers. In order to test this hypothesis thoroughly, it would be desirable to perform numerical simulations at higher Marangoni numbers than those considered by Boeck & Thess (Reference Boeck and Thess2001). Further DNS would also enable one to check if our rigorous bound is sharp and if the assumptions in Boeck & Thess’ scaling argument (most notably, the stability of simple steady convection rolls) should be revised.
Second, our result is the first upper bound proved with the background method that has a logarithmic correction with negative exponent. This is reminiscent of scaling laws obtained for wall-bounded flows through ‘mixing length’ turbulent theories (see, for example, chapter 3 in Doering & Gibbon (Reference Doering and Gibbon1995)). While we are not aware of any such theories being applied to Bénard–Marangoni convection, they have historically motivated the development of rigorous upper-bounding theory in general, and the background method in particular (Doering & Constantin Reference Doering and Constantin1992). In the future, it would be interesting to see if bounds with logarithmic corrections with negative exponent are provable for other flows, starting with extensions of the basic model considered in this work to more general types of thermal boundary conditions (see, for example, Pearson (Reference Pearson1958) and Fantuzzi & Wynn (Reference Fantuzzi and Wynn2017)).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Estimates on 
$f_{k}(z)$
For the lower bound on 
$f_{k}(z)$ in (2.6), observe that the functions 
$h(kz)$ and 
$h[k(1-z)]$ are, respectively, increasing and decreasing on 
$[\unicode[STIX]{x1D6FD},1]$ for any fixed 
$k>0$. This means that the function 
$|f_{k}(z)|(1-z)^{-1}$ increases on 
$[\unicode[STIX]{x1D6FD},1]$, which yields the lower bound.
For the upper bound in (2.6), instead, rewrite (2.4) as
$$\begin{eqnarray}\left|\,f_{k}(z)\right|=\frac{k(1-z)}{2}\text{e}^{-k(1-z)}g_{k}(z)\end{eqnarray}$$with
$$\begin{eqnarray}g_{k}(z):=\frac{h(k)h(kz)-h(k(1-z))}{h(2k)-1}kz\text{e}^{k(1-z)}.\end{eqnarray}$$Differentiation gives
$$\begin{eqnarray}g_{k}^{\prime }(z)=\frac{\text{e}^{2k(1-z)}\ell _{k}(z)}{2(h(2k)-1)(1-z)^{2}},\end{eqnarray}$$ with 
$\ell _{k}(z):=\text{e}^{-2k}(\text{e}^{2kz}-(1-z)^{2})+z^{2}(1-2k)+2z(k-1)$. Now, 
$\ell _{k}(0)=0=\ell _{k}(1)$ and 
$\ell _{k}^{\prime }(0)>0$. Further, 
$\ell _{k}^{\prime }$ is the sum of a convex and a linear function, meaning that 
$\ell _{k}$ has at most two stationary points. Since 
$\ell _{k}^{\prime }(1)=0$, there is at most one stationary point in 
$0<z<1$. Thus, both 
$\ell _{k}(z)\geqslant 0$ and 
$g_{k}^{\prime }(z)\geqslant 0$ for 
$z\in [0,1]$. From this we conclude that
$$\begin{eqnarray}g_{k}(z)\leqslant g_{k}(1)=k[h(k)^{2}-1][h(2k)-1]^{-1}\leqslant 1,\end{eqnarray}$$which, by (A 1), proves the upper bound in (2.6).
Finally, to show (2.5), use the definition of 
$h$ and the inequalities 
$1\leqslant x\coth x\leqslant 1+(x^{2}/3)$, which are valid for 
$x\geqslant 0$, to obtain
$$\begin{eqnarray}(z-1)\frac{h[k(1-z)]}{h(k)h(kz)}=kz[\coth (k)-\coth (kz)]\leqslant z\left(1+\frac{1}{3}k^{2}\right)-1.\end{eqnarray}$$ Combining this estimate with (2.4) and the identity 
$\left|\,f_{k}\right|=-f_{k}$ gives
$$\begin{eqnarray}\left|\,f_{k}(z)\right|\leqslant \frac{1}{2}k^{2}z\frac{h(k)h(kz)}{h(2k)-1}\left[1-z+z\left(1+\frac{1}{3}k^{2}\right)-1\right]=\frac{1}{6}k^{4}\frac{h(k)h(kz)}{h(2k)-1}z^{2}.\end{eqnarray}$$ Since 
$h$ is increasing, so 
$h(kz)\leqslant h(k\unicode[STIX]{x1D6FD})$ for all 
$z\in [0,\unicode[STIX]{x1D6FD}]$, the upper bound (2.5) follows.
