Rigorous scaling laws for internally heated convection at infinite Prandtl number

New bounds are proven on the mean vertical convective heat transport, $\overline{\langle wT \rangle}$, for uniform internally heated (IH) convection in the limit of infinite Prandtl number. For fluid in a horizontally-periodic layer between isothermal boundaries, we show that $\overline{\langle wT \rangle} \leq \frac12 - c R^{-2}$, where $R$ is a nondimensional `flux' Rayleigh number quantifying the strength of internal heating and $c = 216$. Then, $\overline{\langle wT \rangle} = 0$ corresponds to vertical heat transport by conduction alone, while $\overline{\langle wT \rangle}>0$ represents the enhancement of vertical heat transport upwards due to convective motion. If, instead, the lower boundary is a thermal insulator, then we obtain $\overline{\langle wT \rangle} \leq \frac12 - c R^{-4}$, with $c\approx 0.0107$. This result implies that the Nusselt number $Nu$, defined as the ratio of the total-to-conductive heat transport, satisfies $Nu \lesssim R^{4}$. Both bounds are obtained by combining the background method with a minimum principle for the fluid's temperature and with Hardy--Rellich inequalities to exploit the link between the vertical velocity and temperature. In both cases, power-law dependence on $R$ improves the previously best-known bounds, which, although valid at both infinite and finite Prandtl numbers, approach the uniform bound exponentially with $R$.


Introduction
Convective flows driven by internal sources of heat have attracted renewed interest in recent years [1][2][3][4][5][6]. Such flows are commonly encountered in geophysics, where atmospheric convection [7] and mantle convection [8,9] are typical examples. They also exhibit unique features not seen in boundary-driven Rayleigh-Bénard convection: for instance, it has recently been observed experimentally that internally heated (IH) convection can transport heat more efficiently than Rayleigh-Bénard convection [10,11]. Nevertheless, the former remains much less studied. In particular, it remains a largely open challenge to rigorously predict how key statistical properties such as the mean vertical heat flux depend on the heating strength and on the fluid's Prandtl number Pr , defined as the ratio between the fluid's kinematic viscosity ν and its thermal diffusivity κ.
One source of difficulty for mathematical studies of IH convection is that the mean thermal dissipation, the mean viscous dissipation, and the mean vertical convective heat flux cannot all be related to each other via a priori relationships. This is in contrast with Rayleigh-Bénard convection, where such relationships enable one to rigorously bound the convective heat transfer through variational analysis of the mean thermal dissipation [12]. Applying the same strategy to IH flows yields bounds on the mean temperature of the fluid [13][14][15][16] but not the convective heat flux. Recently, this variational strategy was extended by taking into account a minimum principle for the temperature, leading to bounds on the mean convective heat flux that approach a constant exponentially fast as the heating strength is increased [2,5]. Here, we demonstrate that these bounds can be improved to algebraic powers when the Prandtl number is taken to be infinite.
Using standard non-dimensional variables [15], Here, u = (u, v, w) is the fluid velocity in cartesian components, p is the pressure, T is the temperature, and the unit forcing in (1.1c) represents the non-dimensional internal heating rate. The flow is controlled by a 'flux' Rayleigh number that measures the destabilising effect of the heating compared to the stabilising effects of diffusion, R := gαQd 5 ρc p νκ 2 . (1.2) Here g is the acceleration of gravity, ρ is the density, c p is the specific heat capacity, α is the thermal expansion coefficient and Q is the heating rate per unit volume. We consider two separate configurations that differ in the choice of boundary conditions at the top (z = 1) and bottom (z = 0) of the domain. In the first configuration, referred to as IH1 and sketched in Figure 1(a), the velocity satisfies noslip conditions and the temperature of both vertical boundaries is held at a constant value, which can be taken as zero without loss of generality. Hence, we enforce IH1: u| z∈{0,1} = 0, T | z=0 = 0, T | z=1 = 0. (1.3a) In the second configuration, illustrated by Figure 1(b), the bottom plate is replaced by a perfect thermal insulator, giving IH3: We seek bounds on the mean vertical convective heat transport, wT , a quantity that is directly proportional to the viscous dissipation for both the IH1 and IH3 configurations. Throughout this paper overbars denote infinite time averages, while angled brackets denote volume averages: For either set of boundary conditions considered in this paper, it is known that 0 ≤ wT ≤ 1 2 uniformly in R and Pr [17]. The zero lower bound is saturated by the (possibly unstable) state in which the flow does not move and heat is transported vertically by conduction alone. The upper bound of 1 2 , instead, takes on a different meaning depending on the thermal boundary conditions. As will become apparent from (1.6) below, for the IH1 configuration a flow with wT = 1 2 would see all heat escape the domain through the upper boundary. This, however, cannot be achieved at any finite value of R. Precisely, our first main results reveal that the mean vertical heat flux wT is strictly smaller than 1 2 by an amount that cannot decrease faster than quadratically as R is raised.
The bound on wT can be given a clear physical interpretation as a measure of the asymmetry of the heat transport due to heating. Indeed, upon computing z · (1.1c) one can show that the average nondimensional heat fluxes through the top and bottom boundaries, denoted by F T and F B respectively, can be expressed as where · h denotes a spatial average over the horizontal directions alone. Thus, our upper bound on wT immediately implies bounds on F T and F B .
Corollary 1.1. For all sufficiently large R, The results discussed so far apply to the IH1 configuration, where the two horizontal domain boundaries are isothermal. Similar results hold also for the IH3 case, where the bottom boundary is perfectly insulating. In this case, however, the deviation of wT from 1 2 may decay as fast as R −4 . Since the IH3 boundary conditions imply that the conductive heat flux is positive, the effects of convection on the enhancement of heat transport in the system can be described using a Nusselt number, Nu. This is defined as the ratio of the mean total heat flux to the mean conductive heat flux, and can be expressed in terms of wT as (1.9) Thus, our upper bound on wT can be transformed into an upper bound on Nu.
The proofs of Theorems 1 and 2 rely on two key ingredients. The first is a variational problem giving an upper bound on wT . This variational problem is derived by enforcing a minimum principle for the fluid's temperature within the classical "background method" [12,18,19], which for simplicity we formulate using the language of a more general framework for bounding infinite-time averages [20][21][22][23] (see [21,24] for further discussion of the link between the two approaches). Using this minimum principle is essential to obtain bounds on wT that asymptote to 1 2 from below. This has already been shown for the IH1 configuration at finite Pr : for this case, without the miniumum principle one obtains only wT R 1/5 [2], while with it one can prove that wT ≤ 1 2 − O(R 1/5 exp (−R 3/5 )) uniformly in Pr [5]. A similar (but not identical) exponentially-varying bound of wT ≤ 1 2 − O(R −1/5 exp (−R 3/5 )) uniformly in Pr was also obtained for the IH3 configuration [5].
The second key ingredient in our proofs are estimates of Hardy-Rellich type, obtained by observing that the reduced momentum equation (1.1b) determines the vertical velocity field as a function of the temperature field. Specifically, taking the vertical component of the double curl of (1.1b) gives where ∆ h := ∂ 2 x + ∂ 2 y , is the horizontal Laplacian. Using the no-slip boundary conditions with the incompressibility condition (1.1a), the vertical velocity w satisfies (1.11) Equation (1.10) was exploited in Rayleigh-Bénard convection to improve the scaling of upper bounds on Nu [25]. This was achieved by using (1.10) to derive inequalities of the Hardy-Rellich type (see Lemma 4 below) that help the construction of a background field with a logarithmically-varying stable stratification in the bulk [25,26].
Here, we use the same inequalities to construct (different) background fields suited to IH convection, which will enable us to bound wT in the infinite Pr limit.

Bounds for the IH1 configuration
We first consider the IH1 configuration, where the top and bottom plate are held at zero temperature. In §2.1, we show that wT can be bounded from above by constructing suitably constrained functions of the vertical coordinate z. Section 2.2 describes parametric ansätze for such functions, while §2.3 establishes auxiliary results that simplify the verification of the constraints and the evaluation of the bound. We then prove Theorem 1 in §2.4 by prescribing R-dependent values of the free parameters in our ansätze.
To simplify the notation, we introduce two sets of temperature fields that encode the thermal boundary conditions and the pointwise nonnegativity constraint implied by the minimum principle: (2.1a) where n = 1 or 3, depending on the boundary conditions of IH1 and IH3 respectively. In §2, T belongs to H 1 .

Bounding framework
To bound wT we employ the auxiliary function method [20,24]. The method relies on the observation that the time derivative of any bounded functional V{T (t)} along solutions of the Boussinesq equations (1.1) averages to zero over infinite time, so If V is chosen such that the quantity being averaged on the right-hand side is bounded above pointwise in time, then this pointwise bound is also an upper bound on wT . Following analysis at finite Prandtl number [2], we restrict our attention to quadratic functionals taking the form which are parametrized by a positive constant β ∈ R + and a piecewise-differentiable function τ : [0, 1] → R with square-integrable derivative. We require τ to satisfy so the coefficient multiplying T in (2.3) vanishes at z = 0 and z = 1. This choice enables us to integrate by parts without picking up boundary terms when calculating d dt V{T (t)}. To ensure that the resulting expression wT + d dt V{T } can be bounded from above poinwise in time, we also require that the pair (β, τ ) satisfies a condition called the spectral constraint.
If the spectral constraint is satisfied, then it is possible to bound wT from above in terms of τ , β, and another suitably constrained function λ : [0, 1] → R.
Proposition 1 (Bounding framework, IH1). Suppose that the pair (β, τ ) satisfies the spectral constraint and the boundary conditions in (2.4). Further, let λ ∈ L 2 (0, 1) be a nondecreasing function such that λ = −1. Then, Proof. A standard calculation using integration by parts, the incompressibility condition (1.1a), and the boundary conditions on u and T yields The infinite-time average on the right-hand side is bounded above by the largest value of the argument over the global attractor of the infinite-Prandtl-number Boussinesq equations (1.1). The minimum principle for the temperature field implies that the global attractor is contained in the set H + defined in (2.1b). Consequently, we can estimate This upper bound is finite if and only if unless the pair (β, τ ) satisfies the spectral constraint (cf. definition 1), in which case the supremum over T can be evaluated using a technical convex duality argument detailed in appendix A. The result of the argument is that the bound in (2.7) is equivalent to This inequality clearly implies the upper bound on wT stated in the proposition, which is therefore proven.
Remark 3. The function λ that arises when deriving (2.8) from (2.7) can be viewed as a Lagrange multiplier enforcing the pointwise nonnegativity of temperature fields in the set H + . Further details regarding this interpretation (with slightly different notation) are given in [2, §4.4].
Remark 4. The best upper bound on wT provable with our approach is found upon minimizing the expression U (τ, λ, β) over all choices of τ , λ and β that satisfy the conditions of Proposition 1. This is hard to do analytically, but can be done computationally using a variety of numerical schemes (see [24] and references therein). We leave such computations to future work and focus on proving Theorem 1 by constructing suboptimal τ , λ and β analytically.

Ansätze
To prove the upper bound on wT , we seek β > 0, τ (z), and λ(z) that satisfy the conditions of Proposition 1 and make the quantity U (β, τ, λ) as small as possible. To simplify this task, we restrict τ to take the form and λ to be given by (2.10) These piecewise-defined functions, sketched in Figure 2, are fully specified by the bottom boundary layer width δ ∈ (0, 1 2 ), the top boundary layer width ε ∈ (0, 1 2 ), and the parameter A > 0 that determines the amplitude of τ in the bulk of the layer.
We also fix This choice is motivated by the desire to minimize the right-hand side of the inequality which is used later in Lemma 2 by estimating from above the value of the bound U (β, τ, λ) for our choices of τ and λ. For any choice of the parameters δ, ε, and A, the function τ satisfies the boundary conditions in (2.4), while λ is nondecreasing and satisfies the normalization condition λ = −1. To establish Theorem 1 using Proposition 1, we only need to specify parameter values such that U (β, τ, λ) ≤ 1 2 − O(R −2 ) while ensuring that the pair (β, τ ) satisfies the spectral constraint. For the purposes of simplifying the algebra in what follows, we shall fix (2.13) from the outset. As explained in remark 7 below, this choice arises when insisting that the upper estimate on U (β, τ, λ) derived in Lemma 2 be strictly less than 1 2 for some values of δ and ε, at least when all other constraints on these parameters are ignored.

Preliminary estimates
We now derive a series of auxiliary results that make it simpler to specify the boundary layer widths δ and ε. The first result gives estimates on the value of β in (2.11).
Lemma 1 (Estimates on β). Let τ (z), λ(z) and β be given by (2.9), (2.10), and (2.11) with A specified by (2.13). Suppose that the boundary layer widths δ and ε satisfy Then, is key to prove the auxiliary results of this section. The other two restrictions, instead, are introduced to more easily keep track of constants in our estimates, which is necessary to obtain an explicit prefactor for the O(R −2 ) term. We have not attempted to optimize this prefactor.
Remark 6. Condition (2.14c) implies that δ ≤ ε. We will use this fact often in the proofs of this section.
Our second auxiliary result estimates the upper bound U (β, τ, λ) on wT from Proposition 1 in terms of the bottom boundary layer width δ alone.
Remark 7. The right-hand side of (2.19) can be strictly smaller than 1 2 only if A δ 3/2 . It is this observation that dictates the choice of A in (2.13). For any fixed value of R, one should choose A ∼ δ 3/2 with a (possibly R-dependent) prefactor that optimises the balance between the positive and negative terms, subject to constraints on A, δ and all other parameters that ensure the spectral constraint. To simplify our proof, however, we choose to fix this prefactor a priori irrespective of R.
Our final auxiliary result gives sufficient conditions on δ and ε that ensure the spectral constraint (cf. definition 1) is satisfied.
The proof of this result relies on Hardy-Rellich inequalities established in Ref. [14], which extract a positive term from the a priori indefinite term τ wT .
Lemma 4 (Hardy-Rellich inequalities [14]). Let T, w : Ω → R be horizontally periodic functions such that ∆ 2 w = −R∆ h T subject to velocity boundary conditions (1.11). Then, (2.22b) Given our choice of τ from (2.9), we can rewrite the spectral constraint as Here, w is determined as a function of T by solving (1.10) subject to the boundary conditions in (1.11). We shall prove that F and G are individually non-negative. First, let us consider F. The Hardy-Rellich inequality (2.21a) gives Next, we estimate f (z)wT . Since T vanishes at z = 0 by virtue of the thermal boundary conditions in (1.3a), we can use the fundamental theorem of calculus and the Cauchy-Schwarz inequality to estimate |T (·, z)| ≤ √ z(´1 0 |∇T | 2 dz) 1/2 . Squaring both sides and taking the horizontal average of which gives T 2 h ≤ z |∇T | 2 . Then, use of the Cauchy-Schwarz inequality, substitution for T 2 h and Youngs inequality gives Upon using the lower bound on β from (2.15) to estimate the last term from above we obtain This can be substituted into (2.24) along with our choice of A from (2.13) to find To conclude, we show that the term in parentheses is nonnegative when δ satisfies δ ≤ 1 6 and (2.20a). To do this, we observe that the function f (z) in (2.22a) is nonnegative function and that it is nonzero only if z ≤ δ. We can therefore bound it from above on the interval (0, δ) using the estimates 1 1−z ≤ 1 1−δ ≤ 6 5 and, consequently, obtain Using the assumption that δ ≤ 1 6 to estimate the expression in parentheses by its value at δ = 1 6 , followed by an application of assumption (2.20a) to estimate the remaining δ 1/2 term in terms of R we arrive at the desired inequality .
Analogous arguments show that G{T } is nonnegative. Using the Hardy-Rellich inequality (2.21b) we have To estimate the last term, we use (in order) the inequality T 2 h ≤ (1 − z) |∇T | 2 , the Cauchy-Schwarz inequality, and Young's inequality: Using the lower bound on β from (2.15) gives (2.29) To conclude the argument we show that the term in parentheses is non-negative. To demonstrate this, we first estimate g(z) on the interval (1 − ε, 1) from above using the assumption that ε ≤ 1 3 , so 2 3 ≤ z ≤ 1 and 1 z(1−z) ≤ 3 2(1−z) . Thus, Then, we use the assumptions δ ≤ 1 6 and δ ≤ ε (cf. remark 6) to observe that . Combining these estimates with the upper bound on g derived above gives as desired. This concludes the proof of Lemma 3.

Proof of Theorem 1
It is now straightforward to prove the upper bound on wT by specifying boundary layer widths δ and ε that satisfy the conditions of Lemmas 1, 2 and 3. Since the estimate for the resulting upper bound obtained in Lemma 2 is minimized when δ is as large as possible, we choose the largest value consistent with (2.20a), With this choice of δ, conditions (2.14c) and (2.20b) require ε to satisfy 27 648 which is possible for R > R 0 1891.35 (cf. Figure 3). For R > R 0 , any choice of ε in this range is feasible. The optimal value could be determined at the expense of more complicated algebra either by optimizing the full bound U (β, τ, λ), or by deriving better ε-dependent estimates for it. However, we expect that any ε-dependent terms will contribute only higher-order corrections to our bound on wT .
To conclude the proof of Theorem 1, there remains to verify that our choice of δ is no larger than 1 6 and that any ε satisfying (2.31) is no larger than 1 3 . It is easily checked that both conditions hold when R ≥ R 0 (see Figure 3 for an illustration). For all such values of R, therefore, Lemma 2 and our choice of δ yield the upper bound wT ≤ U (β, τ, λ) ≤ 1 2 − cR −2 with c = 216.

Bounds for the IHconfiguration
We now move on to studying IH convection in the IH3 configuration, where the top boundary is maintained at constant (zero) temperature and the bottom boundary is insulating. First, in §3.1, we derive a bounding framework for wT following steps similar to those used for the IH1 case (cf. §2.1). In §3.2 we present ansätze for τ and λ, with which we obtain crucial estimates in §3.3 , which give the bound in §3.4. Throughout this section, T belongs to H 3 . Observe that this changes the set of temperature fields over which the spectral constraint in definition 1 is imposed. The notation H + still denotes the subset of temperature fields in H 3 that are nonnegative pointwise almost everywhere.

Bounding framework
Upper bounds on wT for the IH3 configuration can be derived using a quadratic auxiliary function V{T } similar to that used for the IH1 case. Precisely, we still take V to be defined as in (2.3), where the positive constant β and the piecewisedifferentiable square-integrable function τ (z) are tunable parameters. However, this time we impose only the boundary condition These changes result in the following family of parametrized upper bounds on wT .
Proposition 2 (Bounding framework, IH3). Suppose that the pair (β, τ ) satisfies the spectral constraint (cf. definition 1) and the boundary condition in (3.1). Further, let λ ∈ L 2 (0, 1) be a nondecreasing function such that λ(0) = −1. Then, Proof. Proceeding as in the proof of Proposition 1 shows that The supremum on the right-hand side can be evaluated using the convex duality argument summarized in appendix B, leading to the equivalent inequality This clearly implies the upper bound stated in the proposition.

Ansätze
The procedure for the proof of an upper bound on wT is the same as that employed for isothermal boundaries. We construct β > 0, τ (z) and λ(z) that satisfy the conditions of Proposition 2, while trying to minimize the corresponding bound U (β, τ, λ).
Due to the Neumann boundary condition on T at z = 0, we can no longer employ the Poincaré estimates used in §2.3 to control the sign-indefinite term in the spectral constraint at the bottom boundary. Instead we modify τ (z) in (δ, 1 2 ) to increase slower than logarithmically in z and use results established in [26]. The function τ (z) is hence chosen to have the form On the other hand, the Lagrange multiplier λ(z) is still chosen to be These piecewise functions, sketched in Figure 4, are fully specified by the bottom and top boundary layer widths δ, ε ∈ (0, 1 2 ), the constant A > 0, and the exponent α ∈ (0, 1) driving the behaviour of τ (z) in the bulk.
For β we take This choice is motivated by minimizing the right hand side of the estimate which is used in Lemma 6 below to estimate the value of the bound U (β, τ, λ) from above when τ and λ are define by (3.4) and (3.6) respectively.
For any choice of the parameters δ, ε, A, and α, the function τ satisfies the boundary conditions in (3.1), while λ is nondecreasing and satisfies the condition λ(0) = −1. Thus, to establish Theorem 2 using Proposition 2 we need only specify parameter values such that U (β, τ, λ) ≤ 1 2 − O(R −4 ) while ensuring that (β, τ ) satisfy the spectral constraint. For the purposes of simplifying the algebra in what follows, we shall fix from the outset. This choice arises when insisting that the upper estimate on U (β, τ, λ) derived in Lemma 6 below should be strictly less than 1 2 for suitable choices of δ and ε, at least when all other constraints on these parameters are ignored.

Preliminary Estimates
We now derive auxiliary results that simplify the choice of the exponent α and of the boundary widths δ and ε. The first gives estimates on the value of β in (3.7).
Proof of Lemma 5. It suffices to estimate |τ − λ| 2 1 2 from above and below. For a lower bound, we can substitute our choices of τ and λ and then estimate Dropping the positive term 1/(1−z) from the integrand and integrating the rest gives For every α ∈ ( 1 2 , 1), the second term inside the parentheses can be estimated upon observing that constraints (3.10a-b) imply Thus, we obtain Taking the square root of (3.13) and using (3.9) gives |τ − λ| 2 1 2 ≥ √ 6 9 δ 2 , which combined with (3.7) proves the lower bound on β stated in (3.11).

Lemma 6 (Estimates on
Proof. Inequality (3.8) and the upper bound on β from (3.11) give The result follows upon observing that τ is non-negative, so Remark 9. The right hand side of (3.17) can be strictly smaller than 1 2 when δ 2 is small only if A δ α+3/2 . This observation dictates the choice of A in (3.9). For any fixed value of R, one should choose A ∼ δ α+3/2 with a (possibly R-dependent) prefactor that optimises the balance between the positive and negative terms, subject to constraints on A, δ and all other parameters that ensure the spectral constraint. To simplify our proof, however, we choose to fix this prefactor a priori irrespective of R.
Our final auxiliary result gives the sufficient conditions on δ and ε that ensure the spectral constraint (cf. definition 1) is satisfied.
Unlike the analogous result obtained in §2.3, Lemma 7 cannot be proven using only the Hardy-Rellich inequalities stated in Lemma 4. The lack of a fixed boundary temperature at z = 0, makes it impossible to gain sufficient control on the contribution of the bottom boundary layer to the quadratic form in (2.5). This difficulty can be overcome using the following result, obtained as a particular case of a more general analysis by Whitehead and Wittenberg [26, Eqs. (59) & (77)], which upon setting (in their notation) ν 1 = α 2 − ν 2 and ν 2 = 1 2 −1 + 2(1 − α 2 ) .
Proof of Lemma 7. Let 1 (a,b) denote the indicator function of the interval (a, b) and define the functions (3.22b) Given our choice of τ , we can rewrite the spectral constraint as Observe that F{T } and G{T } are functionals of the temperature field only because w is determined as a function of T by solving (1.10) subject to the boundary conditions in (1.11). We shall prove that F{T } and G{T } are individually non-negative for all temperatures T from the space H 3 , which is sufficient for the spectral constraint to hold.
To prove that F{T } ≥ 0, we apply Lemma 8 with ϕ(z) = f (z)/β and µ = A/β, where f is given by (3.22a) and A given by (3.9). We therefore need to check that To verify this inequality, we first bound from above the weighted integral on the left-hand side. By assumption we have 0 ≤ z ≤ δ ≤ 1 6 and α ∈ ( 1 2 , 1), from which we obtain 1 1−z ≤ 1 z α . Using this estimate and the definition of A from (3.9) we can therefore estimate Using again that 0 ≤ δ ≤ 1 6 and α ∈ ( 1 2 , 1), the bracketed expression can be bounded from above to obtain Next, we bound from below the right-hand-side of (3.24). Using the lower bound on β from Lemma 5, the definition (3.9) of A, and the fact that α ∈ ( 1 2 , 1), we have which is true because δ satisfies (3.19a) by assumption. This proves that F{T } ≥ 0, as desired. We now prove that G{T } is also nonnegative. This can be done following the same steps used in §2.3. The Hardy-Rellich inequality (2.21b) gives To estimate the last term, as before we use the inequality T 2 h ≤ (1 − z) |∇T | 2 , the Cauchy-Schwarz inequality, and Young's inequality: Using the lower bound on β from (3.11) gives which can be substituted into (3.27) along with the value of A from (3.9) to obtain To conclude the argument we need to show that term in the parentheses is nonnegative. To demonstrate this, we first estimate from above the function g(z) given in (3.22b) on the interval (1 − ε, 1). Our assumption that δ ≤ ε implies that 1 − ε ≤ 1 − δ ≤ 1 and ln( 1 ε ) ≤ ln( 1 δ ). Thus, for all δ ≤ 1 6 and α ∈ ( 1 2 , 1) the first term in g(z) can be bounded as To estimate the other terms in g(z), we observe that the assumptions ε ≤ 1 3 and α ∈ ( 1 2 , 1) imply that 1 z α ≤ 1 2(1−α)(1−z) and 1 1−z ≤ 1 2(1−α)(1−z) . Consequently, we arrive at .

Proof of Theorem 2
To prove Theorem 2, we only need to specify R-dependent values for α and for the boundary layer widths δ and ε that satisfy the conditions of Lemmas 5, 6 and 7. Motivated by the desire to minimize the upper bound on U (β, τ, λ) stated in Lemma 6, we choose δ = h(α * )R −2 (3.32) where α * = (1 + √ 13)/6 is the unique maximizer of h(α) on the interval ( 1 2 , 1) (see Figure 5(b)). With these choices, conditions (3.10c) and (3.19b) require ε to satisfy where c 0 , c 1 and c 2 are non-negative constants independent of R. Using the upper bound on B from (3.29), it suffices to find ε such that .
(3.34) Figure 5 shows that suitable values of ε exist when R ≥ R 0 ≈ 2960.89. One can also check that for all such values of R and any ε in the range given by §3.4 one has δ ≤ 1 6 , ε ≤ 1 3 , and δ ≤ ε. We have therefore verified all conditions of Lemmas 5, 6 and 7.
To conclude the proof of Theorem 2, we simply substitute our choice of δ from (3.32) into Lemma 6 to find the upper bound wT ≤ U (β, τ, λ) Remark 10. The top boundary layer width ε is not uniquely determined in our construction. Its optimal value could be obtained by considering more refined estimates on U (β, τ, λ) than Lemma 6, but we expect such estimates to provide only higherorder corrections to the eventual bound on wT .

Conclusions
We have proven upper bounds on the mean vertical convective heat transport wT for two configurations of infinite-Prandtl-number convection driven by uniform internal heating between no-slip boundaries. In the first case, where both boundaries are held at a constant temperature, we find wT ≤ 1 2 − O(R −2 ) for all sufficiently large R (cf. Theorem 1). This result implies that the outward heat fluxes through the top and bottom are bounded by F T ≤ 1 − O(R −2 ) and F B ≥ O(R −2 ), respectively. In the second configuration, where the top boundary remains isothermal but the bottom one is insulating (no-flux condition), we find wT ≤ 1 2 − O(R −4 ) (cf. Theorem 2). In this case, we conclude from (1.9) that the Nusselt number is bounded above by Nu ≤ O(R 4 ). Explicit suboptimal values for the prefactors in the Rayleigh-dependent terms were also obtained (cf. remarks 1 and 2).
All of these results were derived using the background method, which we formulated as a search over quadratic auxiliary functionals of the form (2.3) and augmented using a minimum principle for the fluid's temperature. Similar to previous works on infinite-Prandtl-number Rayleigh-Bénard convection, the background temperature fields used vary linearly in thin boundary layers, and increase either logarithmically (IH1 configuration) or as a power law (IH3 configuration) in the bulk of the fluid layer. This bulk behaviour enables us to use Hardy-Rellich inequalities from [14] (Lemma 4) and an integral estimate from [26] (Lemma 8) that were originally developed in the context of Rayleigh-Bénard convection. In contrast to the latter, however, our background fields lack symmetry in the vertical direction, which reflects the lack of vertical symmetry of IH convection problems.
In our choice of background fields, allowing the bottom boundary layer width δ to be smaller than the top boundary layer width ε is essential to prove Theorems 1 and 2. For the IH1 configuration, forcing δ = ε worsens the R-dependent correction to 1 2 in Theorem 1 to O(R −2 ln −2 (R)). For the IH3 configuration, instead, no upper bound on wT that asymptotes to 1 2 from below as R increases can be obtained with our method of proof if δ = ε. This boundary layer asymmetry contrasts the construction of background fields for IH convection at finite Pr [5], where taking δ = ε appears to bring no qualitative improvement to the exponentially-varying upper bounds on wT . We also stress that the a priori uniform limits on the allowed values of δ and ε imposed throughout §2 and §3 have been chosen with the only goal of simplifying the algebra in our proofs. Varying these limits affects the prefactors of the R-dependent terms in Theorems 1 and 2, as well as the range of R values for which they hold. Both could be optimized further if desired.
One crucial difference between our constructions for the IH1 and IH3 configurations is the leading-order behaviour of the background temperature fields-or, more precisely, of the function τ (z)-as the bottom boundary layer edge is approached from the bulk region. For the IH1 configuration, it suffices for τ to have the same logarithmic behaviour as the background temperature fields used to study Rayleigh-Bénard convection [25]. For the IH3 configuration, however, this choice does not work due to the loss of control on the temperature of the bottom boundary, and we are instead forced to take τ (z) ∼ z 1−α with α ∈ (0, 1). This modification was already used in the context of Rayleigh-Bénard convection between imperfectly conducting boundaries [26], where the optimal exponent α depended logarithmically on the Rayleigh number. Within our proof, instead, the optimal α is a constant. Whether this difference is due to our choice of estimates or the inherent differences between Rayleigh-Bénard and IH convection remains an open question.
More generally, we do not know whether the upper bounds on wT stated in Theorems 1 and 2 are qualitatively sharp. To check if the O(R −2 ) and O(R −4 ) corrections to the asymptotic value of 1 2 are optimal within our bounding framework, one could employ a variation of the computational approach taken in [2] and optimize the tunable parameters τ , β, and λ in full (see also [24] and references therein for more details on the numerical optimization of bounds). A more interesting but also more challenging problem is to identify which convective flows maximize wT and the corresponding optimal scaling of this quantity with R. Considerable insight in this direction can be gained through (i) direct numerical simulations, which to the best of our knowledge are currently lacking; (ii) the calculation of steady but unstable solution of the Boussinesq equations (1.1) that, as recently observed in the context of Rayleigh-Bénard convection [27,28], may transport heat more efficiently than turbulence; and (iii) the explicit design of optimally-cooling flows [6,29,30]. Finally, it would be interesting to investigate if more sophisticated PDE analysis techniques used for Rayleigh-Bénard convection [31] can be extended to IH convection to interpolate between the algebraic bounds on wT proved in this paper for infinite-Pr fluids with the finite-Pr exponential bounds obtained in [5].
Acknowledgements A.A. acknowledges funding by the EPSRC Centre for Doctoral Training in Fluid Dynamics across Scales (award number EP/L016230/1). G.F. was supported by an Imperial College Research Fellowship and would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme "Mathematical aspects of turbulence: where do we stand?" (EPSRC grant number EP/R014604/1) where work on this paper was undertaken.

A Convex duality for the IH1 configuration
The equivalence between (2.7) and (2.8) follow from a relatively standard convex duality argument. It will be enough to show that To establish this identity, we start by rewriting the maximization on the left-hand side as a minimization problem for the Legendre transform of Φ. Recall that the Legendre transform of a functional Λ : H → R is the functional We claim that The first, second and fourth equalities are immediate consequences of the definitions of inf, sup, Ψ, and Ψ * . The third one, instead, follows from the Fenchel-Rockafellar minmax theorem when H is viewed as a Hilbert space with the inner product ∇T 1 · ∇T 2 and norm |∇T | 2 . To apply this theorem as stated in [32, Theorem 1.12], we need to verify that the functionals Φ and Ψ are convex, and that Φ is continuous (with respect to the norm on H) at some T 0 ∈ dom(Φ) ∩ dom(Ψ) = H + .
The convexity of Ψ is obvious, while that of Φ follows from the assumption that the pair (β, τ ) satisfies the spectral constraint (cf. definition 1). To see this, write Φ{T } = Q{T, T } − (τ − βz)∂ z T where Q is the bilinear form and observe that the spectral constraint ensures Q{T, T } ≥ 0 for all T in the linear space H. Thus, for any λ ∈ [0, 1] we can setλ = 1 − λ and estimate proving that Φ is convex. The continuity of Φ at any T 0 ∈ H + follows because Φ is continuous on the whole space H. Indeed, the terms |∇T | 2 and (τ − βz)∂ z T in the expression for Φ{T } are clearly continuous on H. To see that the remaining term is also continuous, it is enough to establish that T k → T in H implies ϕ k := ∆ −2 ∆ h T k → ∆ −2 ∆ h T =: ϕ in H. This can be shown by combining the Poincaré inequalities ϕ 2 |∆ϕ| 2 and T 2 |∇T | 2 with the estimate which implies |∆ϕ| 2 ≤ T 2 . This concludes the proof of Lemma 9.
Next, we prove that since Φ{T } is invariant under horizontal translations of the temperature field T , the minimization of its Legendre transform Φ * can be restricted to functionals µ ∈ H * + that are translation invariant. Specifically, for any real numbers r, s define the translation map S r,s : H → H and its adjoint S * r,s : H * → H * via S r,s T := T (x + r, y + s, z) S * r,s µ := µ • S r,s .
The functional µ ∈ H * is translation invariant if S * r,s µ = µ for all r and s.
Lemma 10. Suppose Φ : H → R satisfies Φ{S r,s T } = Φ{T } for all r, s ∈ R. Then, Proof. It suffices to show that for every ν ∈ H * there exists a translation-invariant µ ∈ H * such that Φ * {µ} = Φ * {ν}. Since every temperature field T ∈ H is horizontally periodic, such a µ can be constructed simply by averaging the functionals S * r,s ν over horizontal translations r and s, i.e., by letting To obtain the reverse inequality observe that, by definition of Φ * , for any ε > 0 there exists T ε ∈ H such that ν{T ε }−Φ{T ε } ≥ Φ * {ν}−ε. Then, since Φ is translation invariant, Upon averaging this inequality over all horizontal shifts r, s ∈ R and applying Jensen's inequality to Φ * , which is concave because it is the supremum of linear functions, we find Φ * (ν) − ε ≤ r,s Φ * S * r,s ν ≤ Φ * r,s S * r,s ν = Φ * (µ).
To establish identity (A.1) we now need to show that its right-hand side coincides with the infimum of Φ * over translation-invariant functionals µ ∈ H * + . For this, we use a characterization of such µ established in [2, Appendix C].

B Convex duality for the IH3 configuration
The equivalence between the upper bounds (3.2) and (3.3) for the IH3 configuration follows from the identity This identity can be proven using a convex duality argument analogous to that in appendix A. Indeed, Lemmas 9 and 10 apply to the functional Φ considered in this section with no changes to their proofs. Consequently, To calculate the Legendre transform Φ * , however, we must replace Lemma 11 with a different characterization of translation-invariant linear functionals µ ∈ H * + . This is due to the different boundary conditions imposed on the temperature space H.

Lemma 12.
Let H * + be the set of positive linear functionals on the temperature space H defined in (2.1a). If µ ∈ H * + is translation invariant, there exists a nondecreasing function λ ∈ L 2 (0, 1) nonnegative almost everywhere and such that µ{T } = −λ(z) ∂ z T .
To conclude the argument, we need to calculate Φ * {µ} for translation-invariant µ ∈ H * + , which by Lemma 12 is given by Since the pair (β, τ ) was assumed to satisfy the spectral constraint, this maximization problem can be restricted to temperature fields that depend only on the vertical coordinate z (this can be proven by splitting T into its horizontal mean η and a perturbation ξ with zero horizontal mean, as outlined at the end of appendix A). The optimal value can then be shown to be Changing the optimization variable on the right-hand side toλ = λ − 1 yields (B.1).