Bounds on heat transfer for B\'enard-Marangoni convection at infinite Prandtl number

The vertical heat transfer in B\'enard-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 proven by Hagstrom&Doering that $ Nu\leq 0.838\,Ma^{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 \leq 0.803\,Ma^{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 \leq 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.

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 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 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 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.

G. Fantuzzi, A. Pershin and A. Wynn
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 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. 2013Wen et al. , 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π] × [0, 1], with x and z denoting the horizontal and vertical coordinates, respectively. The L 2 and L ∞ norms Bounds on heat transfer for infinite-Pr Bénard-Marangoni convection 565 in the z direction of a generic quantity q(z, ·) will be denoted by · 2 and · ∞ , respectively, i.e. Since infinite-time averages need not exist in general, one could be more rigorous and replace lim with lim sup. Note also that |q(z)| 2 = q 2 2 for any quantity q that depends only on the vertical coordinate z.

Pearson's model
Consider a two-dimensional layer of incompressible fluid of depth h, density ρ, kinematic viscosity ν, thermal diffusivity κ and thermal conductivity λ (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 /κ as the time unit and qh/λ as the temperature unit. When the Prandtl number Pr = ν/κ is infinite, Pearson's equations for the fluid's motion (Pearson 1958) reduce to (Hagstrom & Doering 2010) where u(x, z, t) = u(x, z, t)i + w(x, z, t)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π, and satisfy the vertical boundary conditions (BCs) 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 The Marangoni number Ma = γ qh 2 /(λρνκ), where γ 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 u(x, z, t) = 0, p = constant, T(x, z, t) = −z is asymptotically stable when Ma 66.84 , while for Ma 79.61 it is subject to linear instabilities (Pearson 1958) and convection 566 G. Fantuzzi, A. Pershin and A. Wynn sets in (Boeck & Thess 1998. Taking the divergence of (2.1a) and using incompressibility shows that ∇ 2 p = 0, so taking the Laplacian of (2.1a) gives 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ŵ k (z), k ∈ Z, of the vertical velocity w can be computed as a function of the horizontal Fourier coefficientsT k (z) of the temperature. One finds (Hagstrom & Doering 2010) where f 0 (z) = 0 (soŵ 0 = 0 and w has zero horizontal mean), and Note that the function f k satisfies f k (z) 0 for z ∈ [0, 1], f k (0) = 0 = f k (1) and f k (z) → 0 pointwise for all z ∈ (0, 1) as k → ∞ (see figure 1; note that the corresponding figure in Hagstrom & Doering's original paper is incorrect: they plot the negative of f k ). Convection enhances the vertical heat transport, and since the BC ∂ 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 where |∇T| 2 = (∂ x T) 2 + (∂ 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×(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 Bounds on heat transfer for infinite-Pr Bénard-Marangoni convection 567 where the steady background field τ (z) satisfies the BCs while the time-dependent perturbation θ(x, z, t) is periodic in the horizontal direction and satisfies Upon substituting this decomposition into (2.1a) we obtain an evolution equation for the perturbation θ , (3.4) Averaging θ×(3.4) over the volume and infinite time, followed by appropriate integration by parts using (2.1c) and the BCs for θ in (3.3), shows that Moreover, substituting (3.1) into (2.7) gives the two identities Taking the linear combination α×(3.5)−β×(3.6a)+(3.6b) for scalar balance parameters α, β = 1 to be determined, using the fact that θ (1) = ∂ z θ by virtue of (3.3), and rearranging yields where If the balance parameters are chosen to satisfy where the infimum is taken over all horizontally periodic fields θ 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 θ 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 · in (3.8) as a volume average.

G. Fantuzzi, A. Pershin and A. Wynn
To compute the infimum in (3.10) we substitute the Fourier expansions for θ and w into (3.8). Noticing thatθ k =T k for k = 0 by virtue of (3.1), and that f 0 (·) = 0 in (2.5), the Fourier coefficientsŵ k can be expressed in terms ofθ k aŝ (3.11) Moreover,θ −k =θ * k (where * denotes complex conjugation) because the Fourier modes must combine into the real-valued temperature perturbation θ. Consequently, we may rewrite while for k 1 the last term in (3.8) vanishes and we have (3.14) Now, the infimum of Q{θ, w} must be negative semidefinite since Q{0, 0} = 0. Moreover, each functional Q k , k 0 must be individually lower bounded because among all perturbations θ, w are those with only one horizontal wavenumber. In light of (3.3), this lower bound must be sought over all complex-valued functionsθ k (z) that satisfyθ k (0) = 0 =θ k (1). Since Q 0 {0} = 0, the infimum of Q 0 must be negative semidefinite. When k 1, instead, Q k is a homogeneous functional and so if it is lower bounded, its infimum must be exactly zero. Consequently, (3.15) In appendix A we show that (3.16) Substituting this into (3.10), and using the fact that by virtue of (3.2) to simplify the resulting expression, we obtain This bound is valid if (3.9) holds, and if the background field τ is chosen to make the functional Q k {θ k } in (3.14) positive semidefinite for all (integer) wavenumbers k 1. The latter set of constraints can be combined into the single condition that Bounds on heat transfer for infinite-Pr Bénard-Marangoni convection 569 for all perturbations θ, w with zero horizontal mean that satisfy (3.3) and (3.11). Using well-established terminology, we refer to such θ 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) Note that we look for the supremum of the objective function (rather than its maximum) because the strict inequality (α − 1)/(β − 1) > 0 may prevent the existence of a maximiser.
3.2. Optimisation over β The lower bound (3.18) can be optimised over β in a relatively straightforward way, because the spectral constraint is independent of β. Upon setting to zero the first derivative of the right-hand side of (3.18) with respect to β, and using (3.17) to rearrange, we find two stationary values, Inspection of the second derivative of the right-hand side of (3.18) with respect to β reveals that when α is constrained by (3.9) both choices β = β + and β = β − correspond to a local maximum. Determining the optimal choice of β therefore requires comparing the values of such local maxima. After choosing β = β + and re-parametrising α = λ/(λ − 1) -with λ > 1 to satisfy (3.9) -we can use (3.17) to rewrite (3.18) as 1 Nu The spectral constraint (3.19) can also be expressed in terms of λ as Upon introducing the scaled background field ρ(z) = λτ (z) = α/(α − 1)τ (z), subject to a suitably scaled version of the BCs in (3.2), the optimal bound on Nu corresponding to the choice β = β + is found by solving the variational problem sup ρ(z),λ 1 − ρ + 1 2 − ρ(1) 2 subject to |∇θ| 2 + ρ wθ 0 ∀ admissible θ , w, ρ(0) = 0, ρ (1) = −λ, λ > 1.  Similar steps show that the best possible bound on Nu when setting β = β − in (3.18) is given by the solution of an optimisation problem that differs from (3.24) only in the constraint for λ, (3.25) The key observation at this stage is that the suprema in (3.24) and (3.25) coincide despite the different constraint on λ, and furthermore they are equal to the optimal value of the variational problem max ρ(z) (3.26) In fact, for any value of λ we can construct a feasible ρ(z) for either (3.24) or (3.25) that approximates the solution of (3.26) arbitrarily accurately: simply let ρ 0 (z) be an ε-suboptimal strictly feasible point for (3.26), and choose ρ (z) = ρ 0 (z) in (3.24) or (3.25) except for an infinitesimally thin layer near z = 1, where ρ (z) = −λ. A rigorous argument follows steps similar to those used in the energy stability analysis of the conductive state , and is omitted for brevity. The conclusion is satisfactory: the bound on Nu is independent of whether one sets β = β + or β = β − in (3.18).
3.3. An explicit value for the optimal β The variational principle (3.26) has been obtained by optimising the balance parameter β as a function of the other balance parameter, α, and the background field τ (z). Interestingly, the optimality conditions for the solution ρ (z) of (3.26) allow for the derivation of a precise numerical value for the optimal β even though the optimal α and τ (z) are unknown. To show this, we introduce a variable s such that ρ + 1 2 s and note that (3.26) is equivalent to The feasible set of this problem is convex, so the linear objective function is maximised on the constraint boundary. Since for any given ρ(z) we can always choose s = ρ + 1 2 , the optimal bound is attained when ρ(z) is on the boundary of the feasible set of the spectral constraint, i.e. when inf θ ,w =0 |∇θ| 2 + ρ wθ = 0. (3.28) Since the spectral constraint is homogeneous in θ and w, it suffices to restrict our attention to admissible θ and w satisfying some normalisation condition N {θ , w} = 0 that excludes the zero fields. The optimal scaled background field ρ (z) and the optimal value s are then those that maximise the Lagrangian functional where ζ , η and µ are scalar Lagrange multipliers. Setting to zero the first variation of L with respect to ρ(z) shows that the optimal scaled background field ρ (z) must satisfy the 'natural' boundary condition (3.30) (Of course, ρ (z) must also satisfy an Euler-Lagrange differential equation, but this will not be important here.) Moreover, setting to zero the derivatives of L with respect to s and η, and eliminating s yields At this point, note that if ρ (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, ρ + 1 2 = 0 and we may use (3.31) to eliminate η from (3.30). The optimal scaled background field must therefore satisfy 1 + ρ + 1 2 + ρ (1) = 0. (3.32) In particular, this implies that −ρ (1) > 1, so ρ (z) is also the optimal solution of (3.24) with λ = −ρ (1). Recollecting the re-parametrisation α = λ/(λ − 1) we conclude that the optimal value of the balance parameter α, denoted α , is given by Finally, recalling that (3.24) was obtained by setting β = β + from (3.21) and that α τ (z)/(α − 1) = ρ (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 β is

Relation to Hagstrom & Doering's variational problem
The bounding principle formulated by Hagstrom & Doering (2010) can be recovered upon setting α = 2 and β = 2 in (3.20). These values clearly satisfy (3.9), and we have seen that the choice β = 2 is optimal. The variational problem for the optimal background field becomes

G. Fantuzzi, A. Pershin and A. Wynn
Strictly speaking we should also enforce the boundary condition τ (1) = −1, but this does not limit the choice of τ for the same reasons discussed at the end of § 3.2.
To bring (4.1) in contact with (3.26), we change variables to ϕ = 2τ and use the boundary condition ϕ(0) = 0 to rewrite (4.1) as 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 ϕ(z) In particular, using (3.26) it is almost immediate to obtain a 4.2 % improvement for the prefactor of Hagstrom & Doering's bound, Nu 0.838Ma 2/7 , at least in the limit of infinite Marangoni number: in appendix B we show that 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 → ∞. To show that this is true, recall from § 3.3 that the choice β = 2 is optimal. After fixing β = 2 and re-parametrising α = λ/(λ − 1) as in § 3 -with λ > 1 to satisfy (3.9) -the bound in (3.18) becomes 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 = Maw and define the scaled Marangoni number M = λMa to rewrite the spectral constraint (3.23) as Upon rescaling w = Maw the Marangoni number drops out of equation (3.11), sow is a (linear) function of θ only and the admissible test functions in (4.6) are independent of M. Then, consider the family of background fields τ M , parametrised by the scaled Marangoni number M, that maximises the right-hand side of (4.5) for a fixed value of λ. In other words, assume that τ M solves the variational problem min τ (z) τ + 1 2 2 subject to |∇θ | 2 + Mτ wθ 0 ∀ admissible θ ,w. (4.7) Moreover, suppose σ (M) is such that Note that it is reasonable to assume that σ (M) → 0 as M → ∞, because we expect that τ (z) ≈ 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 λ is then given by 1 Nu Using the fact that dM/dλ = d(λMa)/dλ = M/λ, it is straightforward to show that the right-hand side of the last expression is maximised with respect to λ when and that the corresponding bound on the Nusselt number is (4.11) Now, recall that Nu 1 since convection enhances the purely conductive vertical heat transport. Using the fact that λ > 1 and the assumption that σ (M) → 0 as M → 0 it is then not difficult to see that since Nu 1 the quantity Mσ (M) must be uniformly bounded as the scaled Marangoni number M tends to infinity. Consequently, the solution λ = λ (M) of (4.10) satisfies This implies that Ma = O(M) as M → ∞, meaning that the optimisation over the balance parameter does not influence the asymptotic scaling of the bound on Nu with the Marangoni number.

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 Moreover, we introduce a non-negative variable s such that φ 2 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 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 ρ (z) 0 and ρ (z) 0. This is achieved by solving the convex problems (5.5)

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 finitedimensional 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 S whose entries are affine with respect to a set of optimisation variables is positive semidefinite (denoted by S 0). The norm constraint φ 2 s, instead, is the infinite-dimensional version of a second-order cone constraint (SOCC), i.e. the requirement that a vector y ∈ R n+1 and a scalar s satisfy y s, where · denotes the usual Euclidean norm of a vector. Finally, the pointwise constraints φ(z) 1 and φ (z) 0 are the infinite-dimensional equivalent of element-wise inequalities for a vector y ∈ R n+1 of the form Ay b, with A ∈ R m×(n+1) and b ∈ 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 φ 2 s to a SOCC, we introduce a piecewise-linear ansatz for φ. Given a set of n + 1 collocation points 0 = z 0 < z 1 < · · · < z n−1 < z n = 1, we denote φ i = φ(z i ) for all i = 1, . . . , n and consider where ψ i (z) is the unique piecewise-linear function satisfying ψ i (z i ) = 1 and vanishing at all other nodes (cf. figure 2). After defining the column vector of nodal values Bounds on heat transfer for infinite-Pr Bénard-Marangoni convection it is clear that there exists a positive definite matrix P = R T R such that The norm constraint φ 2 s then becomes the SOCC RΦ s. 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 Q k {θ k } in (3.14) being positive semidefinite for all wavenumbers k 1 and all complex-valued functionsθ k (z) satisfyingθ k (0) = 0 =θ k (1). Recognising that the real and imaginary parts ofθ k give identical and independent contributions to Q k {θ k }, it suffices to restrict our attention to real-valued functionsθ k (z), so we define the space of test functions Recalling that we have changed variables according to we can therefore replace the spectral constraint in (5.3)-(5.5) with the infinite set of Fourier-transformed spectral constraints where · 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 Q k in (5.11) over non-zero test functions is zero. When k k c , instead, we approximate the test function v using the same piecewiselinear ansatz used for φ, i.e.

G. Fantuzzi, A. Pershin and A. Wynn
We also set v 0 = 0 and v n = v n−1 in order to enforce the boundary conditions v(0) = 0 and v (1) = 0, but we do not do this explicitly in (5.13) to simplify the following discussion. Substituting (5.13) and (5.6) into Q k {v} from (5.11) yields 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 v := [v 1 , . . . , v n−1 ] T , and there exists an Consequently, for each wavenumber k k c the Fourier-transformed spectral constraint Q k {v} 0 can be approximated by the LMI Q k (Φ) 0. Finally, it is easy to see that the piecewise-linear approximation (5.6) turns the pointwise inequality φ(z) 1 into the n + 1 constraints φ i 1, i = 0, . . . , n, which can be written succinctly as the element-wise vector inequality Φ 1. Similarly, the condition φ (z) 0 becomes a set of n inequalities φ i−1 − φ i 0, i = 1, . . . , n, which can be written in the vector form   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 φ(z) means that only lower bounds on the optimal values of (5.3)-(5.5) can be computed, because the 'true' optimal φ(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 Γ . This enlarges the set of feasible functions φ(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 Γ and their finitedimensional 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 . 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).
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 φ(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 φ(z) < 0, but setting φ(z) ≈ 1 in the bulk of the domain is necessary to be able to satisfy the spectral constraint. However, it is possible to have φ(z) < 0 in thin layers near the walls because the functions f k , which act as a weight on φ in the Fourier-transformed spectral constraint (5.11), are 578 G. Fantuzzi, A. Pershin and A. Wynn 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 Q k (Φ) 0 with a set of LMIs on certain 3 × 3 submatrices of Q k . Since the number of rows/columns of Q 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 × 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 Φ according to (5.12). We therefore employ the following iterative procedure: find the optimal Φ using an initial guess k 0 for k c , update the value of k c using (5.12), check if Q k (Φ) is positive semidefinite for all k 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 × 10 6 we find that the optimal φ satisfies φ − 1 ∞ = 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 × 10 p , which gives p + 1 logarithmically spaced points between successive powers of 10. Given the critical wavenumbers k 1 , . . . , 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, . . . , m, i.e. values of k such that (5.20) We then checked if the optimal solution satisfied Q k (Φ) 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 Q k with the most negative eigenvalue) to the list of critical values and repeated the optimisation.

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 Downloaded from https://www.cambridge.org/core. 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 Ma 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 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 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 . cases, at each value of Ma the optimal φ(z) was used to recover the optimal scaled background field ρ(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 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 Nu 1.285Ma 2/7 (ln Ma) −1/2 . (5.21) 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 ρ(z) is monotonic, but it may be lowered by a logarithm upon construction of a non-monotonic background field. 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 ρ (z) instead of ρ(z) because by virtue of (5.1) problems (5.3)-(5.5) can be rewritten in terms of ρ (z) alone. Since ρ(z) can be recovered by integration using the boundary condition ρ(0) = 0, the derivative ρ (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 ρ (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 580 G. Fantuzzi, A. Pershin and A. Wynn brevity). Such evidence corroborates our conjecture that (5.21) is the correct functional form of the optimal bound on Nu.
boundary layers separated by a bulk region where ρ (z) ≈ 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;).
In the fully optimal case, ρ (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 ρ (z) constant across the boundary layer near z = 1.
Further details of the boundary layer structure of the fully optimal profiles for Ma 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 ρ (z) near z = 0 and z = 1, respectively, we take δ := z bot and ε := 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 10 7 return δ ≈ 3.8Ma −0.26 , ρ (z bot ) ≈ 0.07. (5.22a,b) Note that the scaling exponent of δ is not far from −2/7 ≈ −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 δ = O(Ma −2/7 (ln Ma) 1/2 ), such that Nu = O(δ −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 ε: ε ≈ 1.65Ma −0.36 for Ma 5 × 10 4 , (5.23a) ε ≈ 11.8Ma −0.54 for 5 × 10 4 Ma 3 × 10 6 , (5.23b) ε ≈ 24.3Ma −0.58 for Ma 3 × 10 6 . (5.23c) In the first and third regions we could also determine approximate scaling laws for the peak value ρ (z top ): ρ (z top ) ≈ 0.34Ma −0.01 for Ma 5 × 10 4 , (5.24a) ρ (z top ) ≈ 0.04Ma 0.16 for Ma 3 × 10 6 . (5.24b) 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 ε after rescaling by Ma 0.54 , which at large Ma amplifies small inaccuracies in our numerical data. 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 ≈ 5 × 10 4 and Ma ≈ 3 × 10 6 . Another intermediate branch of critical wavenumbers appears for Ma ≈ 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 ρ (z) = φ(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 .
To conclude this section, we plot in figure 8 the variation with Ma of the optimal balance parameter α , computed using (3.33) and the fully optimal background field. The results are interesting for two reasons. First, the convergence of α to 2 as Ma is raised is logarithmic, giving further evidence in support of (5.21). Second, the results suggest that the choice α = 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 → ∞ 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.

Towards an improved bound
The results presented in § 5.3 suggest that Hagstrom & Doering's bound Nu 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 → ∞ 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 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 ρ(z) in the slightly rearranged form For example, simply estimating G. Fantuzzi, A. Pershin and A. Wynn and requiring (as we have done in appendix B) that forces the optimal ρ to decrease monotonically. In fact, if ρ satisfies (6.3) and ρ (z) 0 for z ∈ U ⊂ [0, 1], the profilẽ 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 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 δ ∈ (0, 1) and write Here, δ 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 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 1 δ h(z)f k (z) dz we then need to check that If 1 δ h(z)f k (z) dz < 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) 0, this requires h(z) > 0 over a sufficient portion of the interval (δ, 1), so the background field ρ 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 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 to ease the notation, a sufficient condition for (6.7) is that v 2 2 − Ma(a k + b k ) v 2 + Ma c k 0, (6.11) which in turn is satisfied if 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 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) with γ > 0 a constant (independent of the Marangoni number) and ε 1 but such that 1/ε k c = O(Ma 1/2 ). When k 1/ε 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

G. Fantuzzi, A. Pershin and A. Wynn
Using these estimates, (6.12) can be rearranged as When k = O(1) the two sides of (6.15) could be balanced by taking ε = O(Ma −1/3 ) and δ = O(Ma −5/21 ), and upon computing the bound on Nu we find Interestingly, the exponent 5/21 ≈ 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 ≈ 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 → ∞, 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/ε (which is no larger than k c by assumption) shows that we must choose ε O(Ma −1/2 ) and δ 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 γ increase with Ma to mimic the behaviour of the numerically optimal profiles (cf. figure 6d), because the apparent gain in (6.15) is exactly outbalanced by the need of testing wavenumbers up to k c = O(γ 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 Bounds on heat transfer for infinite-Pr Bénard-Marangoni convection 587 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. (2013Wen et al. ( , 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.

G. Fantuzzi, A. Pershin and A. Wynn
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.

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 0.803 × 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 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 ) 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) 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-α' 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 → ∞ 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 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 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 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 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; 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.
for δ > 0 sufficiently small. A rigorous proof is omitted for brevity, but a similar argument can be found in a previous work by the authors (Fantuzzi & Wynn 2017, appendix C).
The minimiser of the Dirichlet problem (A 2) satisfies the Euler-Lagrange equation subject to the BCs v(0) = 0 and v(1) = A, and is given by The corresponding minimum is An expression for the minimum over A is readily found, and it can be rearranged in the form (3.16) after noticing that τ (1) = 1 0 τ (z) dz by virtue of (3.2).

(B 1)
The boundary layer slope R > 0 and thickness δ > 0 should be chosen to satisfy the spectral constraint (3.19) whilst optimising the bound on the Nusselt number, 1 Nu Recall from § 3 that the spectral constraint is equivalent to the quadratic form Q k {θ k } in (3.14) being positive semidefinite for all wavenumbers k 1, and recall that we have changed variables such that α/(α − 1)τ (z) = ρ (z). Although the test functionθ k is complex valued, the contributions of its real and imaginary parts to Q k {θ k } are identical and independent, so it suffices to consider real-valued test functions. We conclude that R and δ must be chosen such that, for all k 1,