1. Introduction
Convection driven by spatially varying heating is attracting attention as a generalisation of internally heated convection (IHC) due to its relevance in studying different regimes of heat transport by turbulence. Recent work demonstrates that varying the heating location can enhance heat transport in bounded domains (Lepot, Aumaître & Gallet Reference Lepot, Aumaître and Gallet2018; Bouillaut et al. Reference Bouillaut, Flesselles, Miquel, Aumaître and Gallet2022; Kazemi, Ostilla-Mónico & Goluskin Reference Kazemi, Ostilla-Mónico and Goluskin2022) and that bounds on the heat transport depend on the supply of potential energy due to variable heating (Song, Fantuzzi & Tobasco Reference Song, Fantuzzi and Tobasco2022). Such studies of non-uniformly heated convection have implications for geophysical fluid dynamics, from mixing in lakes, where solar radiation acts as a spatially varying heat source, to convection in the Earth's mantle and liquid outer core driven by energy released irregularly by radioactive isotopes and secular cooling (Schubert Reference Schubert2015). Furthermore, non-uniform heating is theoretically significant because it induces a much larger class of flows than those by boundary-forced examples such as Rayleigh–Bénard convection (RBC).
This work considers non-uniform IHC between parallel plates with fixed thermal fluxes. Specifically, we assume that the lower plate is a thermal insulator and that the upper plate is a poor conductor in comparison with the fluid's ability to transport heat (see for example Goluskin Reference Goluskin2015). The heat flux through the upper plate is therefore assumed fixed and, to ensure that a statistically stationary state is realisable, is specified to match the total heat input due to internal heating (figure 1a). The Rayleigh number 
$R$ for this system, defined precisely in §2, quantifies the destabilising effect of internal heating relative to the stabilising effect of diffusion. An emergent physical property of the flow on which we will focus is the mean temperature difference, 
$\smash {{\langle {\delta T} \rangle _h}}$, between the lower and upper plates (we use 
$\left \langle {\cdot } \right \rangle$ to represent an infinite time and volume average, with subscript 
$h$ denoting a horizontal average).
Figure 1. (a) Internally heated and (b) RBC with fixed-flux boundaries, where 
$H(z)$ is a positive non-uniform heating profile. In both panels, the red dashed line denotes the conductive temperature profiles and the red solid line the mean temperature profiles in the turbulent regime. In (a), the temperature profiles are for 
$H(z)=1$.
This paper communicates bounds on 
$\smash {{\langle {\delta T} \rangle _h}}$ that depend on the Rayleigh number 
$R$ of the flow and the spatial distribution of the heat sources. More specifically, we prove that, for any solution of the governing Boussinesq equations,
\begin{equation} \smash{{\langle {\delta T} \rangle_h}} := \langle T|_{z=0} \rangle_h - \langle T|_{z=1} \rangle_h \geq \sigma R^{-\gamma}-\mu, \end{equation}
where 
$\sigma,\mu >0$ are constants depending only on the spatial distribution of the internal heating, and 
$\gamma$ is a scaling rate. A corollary of (1.1) is that downward conduction, here defined as the emergent mean temperature of the upper plate being above that of the lower plate, i.e. 
$\langle T|_{z=1} \rangle _h > \langle T|_{z=0} \rangle _h$, is impossible for 
$R\leq R_0= (\sigma /\mu )^{1/\gamma }$, where 
$R_0$ is a critical Rayleigh number that depends on the spatial distribution of the internal heating.
To give context to this result, observe that a flow for which a rigorous characterisation of 
$\smash {{\langle {\delta T} \rangle _h}}$ is known is fixed-flux RBC, shown schematically in figure 1(b). In that case, symmetric heating and cooling rates are prescribed at the upper and lower boundaries (as opposed to internally), and the temperature drop satisfies the rigorous bound 
$\smash {{\langle {\delta T} \rangle _h}} \geq cR^{-1/3}$ for a constant 
$c>0$ (Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002). The physical interpretation of this result is that the upper plate, through which heat leaves the domain, must, on average, remain colder than the lower one. Further, the difference in temperature between the plates cannot decrease at a rate faster than 
$R^{-1/3}$, where 
$R \rightarrow \infty$ corresponds to an increasingly turbulent and homogenised flow. Noting that 
$R$ is a flux-based Rayleigh number (rather than being defined in terms of an imposed temperature difference), the provable decay rate of the lower bound on 
$\smash {{\langle {\delta T} \rangle _h}}$ agrees with the predictions of some phenomenological theories, namely the so-called ultimate regime (Kraichnan Reference Kraichnan1962; Spiegel Reference Spiegel1963). See Doering (Reference Doering2020) and Lohse & Shishkina (Reference Lohse and Shishkina2023) for recent reviews on the ultimate regime in RBC.
There is a striking and crucial distinction between what can currently be proven about the behaviour of the temperature drop 
$\langle \delta T \rangle _h$ in IHC compared with fixed-flux RBC. In particular, while for fixed-flux RBC downward conduction can be ruled out for all Rayleigh numbers (i.e. that 
$\smash {{\langle {\delta T} \rangle _h}} >0$ must hold for all 
$R$), it is not yet known whether there are any heating profiles 
$H(z)$ that lead to the equivalent property in IHC. For example, in the case where internal heating is applied uniformly (
$H \equiv 1$), it was shown by Arslan et al. (Reference Arslan, Fantuzzi, Craske and Wynn2021a) that 
$\smash {{\langle {\delta T} \rangle _h}} \geq 1.6552R^{-1/3}-0.03868$. This only rules out downward conduction up to the critical Rayleigh number 
$R_0 = 78\,389 = 54.437 R_L$, where 
$R_L=1440$ denotes the Rayleigh number at which the flow becomes linearly unstable (Goluskin Reference Goluskin2015).
The bound (1.1), which is the main result of this paper, gives valuable additional information on the gap between the known behaviour of 
$\smash {{\langle {\delta T} \rangle _h}}$ for IHC with an insulating bottom plate and fixed-flux RBC. Specifically, it allows us to understand how the gap depends on the spatial distribution 
$H(z)$ of the internally applied heating via the constants 
$\mu = \mu (H)$ and 
$\sigma = \sigma (H)$ in (1.1). The dependence of this gap on 
$H$ is interesting since it is not unreasonable to expect that downward conduction may be more likely to occur if heating is concentrated close to the upper plate, as in figure 2(a), rather than close to the lower plate, as in figure 2(b). The specific expressions for 
$\mu$ and 
$\sigma$ derived in this paper do not allow downward conduction to be ruled out (i.e. 
$R_0 <\infty$) for any fixed heating profile 
$H$. However, we show that it is possible to construct a sequence of heating profiles 
$H_r$, indexed by an integer 
$r$, in which heating is asymptotically concentrated towards the lower plate as 
$r \rightarrow \infty$, for which
\begin{equation} R_0(H_r) = \left( \frac{\sigma(H_r)}{\mu(H_r)} \right)^{{1}/{\gamma(H_r)}} \rightarrow \infty \quad \text{as} \ r\rightarrow \infty. \end{equation}
In this sense, we can partially bridge the gap between what is known about 
$\smash {{\langle {\delta T} \rangle _h}}$ for fixed-flux RBC and what appears, at least from the perspective of rigorous mathematical analysis, to be the much more intricate behaviour of internally heated convective flows.
Figure 2. Examples of non-uniform heating profiles in a non-dimensionalised domain with illustrative sketches of the conductive temperature profiles (black solid line) for fixed flux boundary conditions. In (a), the heating is localised near the upper boundary and in (b) near the lower boundary while being zero elsewhere, shown with red and white spaces, respectively. For (c) the heating is sinusoidal with a maximum at 
$z=0.25$ (red solid line) and a minimum of zero at 
$z=0.75$ (yellow solid line).
In summary, this paper proves a Rayleigh-dependent lower bound on 
$\smash {{\langle {\delta T} \rangle _h}}$ that depends on the geometric properties of the internally applied heating. As a corollary, for any heating profile, we can find an expression for the critical Rayleigh number 
$R_0$ below which downward conduction is impossible for any solution of the governing equations. To do this, we develop the mathematical approach of Arslan et al. (Reference Arslan, Fantuzzi, Craske and Wynn2021a) for flows driven by arbitrary (non-uniform) heating. The paper is structured as follows: § 2 describes the problem set-up, in § 3 we prove the lower bound (1.1), § 4 highlights implications of the bound on the possibility of downward conduction in geophysical flows and in § 5 we conclude.
2. Set-up
We consider a layer of fluid between two horizontal plates separated by a distance 
$d$ and periodic in the horizontal (
$x$ and 
$y$) directions with periods 
$L_x d$ and 
$L_y d$. The fluid has kinematic viscosity 
$\nu$, thermal diffusivity 
$\kappa$, density 
$\rho$, specific heat capacity 
$c_p$ and thermal expansion coefficient 
$\alpha$. Gravity acts in the negative vertical direction, and the fluid is heated internally at a position-dependent non-negative volumetric rate 
$\tilde {H}\geq 0$. We assume that 
$\tilde {H}$ is a non-negative integrable function on the domain that depends only on the vertical coordinate and satisfies 
$\langle \tilde {H}\rangle >0$ strictly, and define a dimensionless heating as 
$H(z):=\tilde {H}/\langle \tilde {H}\rangle$. Note that 
$\langle H\rangle =\|H\|_{1}=1$ because 
$H(z)\geq 0$. Here and throughout the paper, 
$\lVert\kern1.4pt f \rVert _p$ is the standard 
$L^p$ norm of a function 
$f:[0,1]\rightarrow \mathbb {R}$.
To non-dimensionalise the problem, we use 
$d$ as the characteristic length scale, 
$d^2/\kappa$ as the time scale and 
$d^2\langle \tilde {H}\rangle /\kappa \rho c_p$ as the temperature scale. The velocity of the fluid 
$\boldsymbol {u}(\boldsymbol {x},t)=(u(\boldsymbol {x},t),v(\boldsymbol {x},t),w(\boldsymbol {x},t))$ and temperature 
$T(\boldsymbol {x},t)$ in the non-dimensional domain 
$\varOmega = [0,L_x]\times [0,L_y]\times [0,1]$ are then governed by the Boussinesq equations
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\partial_t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} p = Pr\, \boldsymbol{\nabla}^2 \boldsymbol{u} + Pr\, R\, T \boldsymbol{e}_3 , \end{gather}$$
$$\begin{gather}\partial_t T + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} T = \boldsymbol{\nabla}^2 T + H(z). \end{gather}$$The non-dimensional numbers are the Prandtl and Rayleigh numbers, defined as
\begin{equation} Pr = \frac{\nu}{\kappa}\quad\text{and}\quad R = \frac{g \alpha \langle \tilde{H}\rangle d^5 }{\rho c_p \nu \kappa^2 }. \end{equation}The boundary conditions are no slip for the velocity and fixed flux for the temperature:
$$\begin{gather} \boldsymbol{u}|_{z=\{0,1\}} = 0, \end{gather}$$
$$\begin{gather}\partial_z T|_{z=0} = 0, \quad \partial_z T|_{z=1}={-}1. \end{gather}$$A diagram of the system is shown in figure 1(a).
It proves useful to define a function 
$\eta \in L^2(0,1)$ such that 
$\eta (z)$ measures the total (dimensionless) heat added to those parts of the domain below a height 
$z$:
\begin{equation} \eta(z):= \int^{z}_0 H(z)\, \textrm{d}z. \end{equation}
The value 
$\eta (1)=1$ corresponds to the total heat added to the domain and, therefore, to the negative heat flux applied at the top boundary in (2.3b). We can obtain an identity for the mean vertical heat transport of the system by multiplying (2.1c) by 
$z$ and taking an infinite-time and volume average. Use of (2.4) and a standard application of integration by parts with the boundary conditions (2.3a) and (2.3b) gives the identity
\begin{equation} \langle wT\rangle + \smash{{\langle {\delta T} \rangle_h}} ={-}\int^{1}_{0}(z-1)H\,\textrm{d}z=\int^{1}_{0}\eta(z)\,\textrm{d}z. \end{equation}
The left-hand side of (2.5) is the sum of the mean vertical convective heat flux, 
$\langle wT\rangle$, and the mean conductive heat flux, 
$\smash {{\langle {\delta T} \rangle _h}}$, and balances the potential energy added to the system. To see why, note that 
$-zH$ is the pointwise supply of (dimensionless) potential energy (the negative sign is there because potential energy is created by positive heating from below or negative heating from above). The middle term in (2.5) therefore accounts for the potential energy created by internal heating and cooling from the top boundary (alternatively, the shifted coordinate 
$z\mapsto z-1$, centred on the top boundary, can be seen as a means of removing the cooling at the top boundary from the calculation).
To remove the inhomogeneous boundary conditions on 
$T$, it is convenient to rewrite the temperature field in terms of perturbations 
$\vartheta$ from the conductive profile 
$T_c$,
\begin{equation} T(\boldsymbol{x},t) = \vartheta(\boldsymbol{x},t) + T_c(z). \end{equation}
The steady conductive temperature profile can be found after taking 
$\boldsymbol {u}=0$ and 
$T(\boldsymbol {x},t)=T(\boldsymbol {x})$ in (2.1c). Given the insulating lower boundary condition, 
$T'_c(z) = - \eta (z)$, where primes denote derivatives with respect to 
$z$. Then, (2.1c) in terms of 
$\vartheta$ becomes
$$\begin{gather} \partial_t \vartheta + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta = \boldsymbol{\nabla}^2 \vartheta + w\, \eta(z) , \end{gather}$$
$$\begin{gather}\partial_z \vartheta|_{z=\{0,1\}} = 0 . \end{gather}$$3. Bounding heat transport
To bound 
$\smash {{\langle {\delta T} \rangle _h}}$ using methods that have been successfully applied to uniform IHC (Arslan et al. Reference Arslan, Fantuzzi, Craske and Wynn2021a), we search for an upper bound on 
$\langle wT\rangle$ and use (2.5) to bound 
$\smash {{\langle {\delta T} \rangle _h}}$. We will prove upper bounds on 
${\left \langle {w\vartheta } \right \rangle }$, which is equal to 
$\langle wT\rangle$ because the incompressibility of the velocity field and the boundary conditions (2.3a) imply that any function 
$f:[0,1]\rightarrow \mathbb {C}$ satisfies 
$\left \langle w f \right \rangle = \int ^{1}_0 \langle w \rangle _h f(z)\, \textrm {d}z = 0$. The derivation in § 3.1 is analogous to that in Arslan et al. (Reference Arslan, Fantuzzi, Craske and Wynn2021a), so we only outline the main ideas in this paper.
3.1. The auxiliary functional method
Rigorous bounds on the mean quantities of turbulent flows can be found with the background method (Doering & Constantin Reference Doering and Constantin1994, Reference Doering and Constantin1996; Constantin & Doering Reference Constantin and Doering1995), which we formulate here in the more general framework of the auxiliary functional method (Chernyshenko et al. Reference Chernyshenko, Goulart, Huang and Papachristodoulou2014; Fantuzzi, Arslan & Wynn Reference Fantuzzi, Arslan and Wynn2022). We consider the quadratic auxiliary functional
\begin{equation} \mathcal{V}\{\boldsymbol{u},\vartheta\} = \unicode{x2A0F}_\varOmega \frac{a}{2\,Pr\,R} |\boldsymbol{u}|^2 + \frac{b}{2}\left|\vartheta - \frac{\varphi(z)}{b}\right|^2 \textrm{d}\kern0.07em\boldsymbol{x} , \end{equation}
where 
$a$ and 
$b$ are non-negative scalars and 
$\varphi (z)$ a function, all to be optimised. As demonstrated by Chernyshenko (Reference Chernyshenko2022), the use of quadratic auxiliary functionals is equivalent to the background method, where 
$a$ and 
$b$ are referred to as balance parameters and 
$\varphi (z)/b$ is the background temperature field satisfying the boundary conditions on 
$\vartheta$ in (2.8). With a slight abuse of terminology, we refer to 
$\varphi (z)$ as the background field.
Given that 
$\mathcal {V}\{\boldsymbol {u},\vartheta \}$ remains bounded in time along solutions of (2.1a), (2.1b) and (2.7) for any given initial 
$\boldsymbol {u}$ and 
$\vartheta$, the infinite-time average of the time derivative of 
$\mathcal {V}$ is zero. Using this property, we can write
\begin{equation} {\left\langle {w\vartheta} \right\rangle} = U - \left( U - \left\langle w\vartheta \right\rangle - \limsup_{\tau\rightarrow\infty} \frac{1}{\tau}\int^{\tau}_0\frac{\textrm{d}}{\textrm{d}t} \mathcal{V}\{\boldsymbol{u},\vartheta\} \,\textrm{d}t \right), \end{equation}
and deduce that 
${\left \langle {w \vartheta } \right \rangle } \leq U$, if after rearranging, the functional in the brackets in (3.2) is pointwise in time non-negative. Following computations analogous to Arslan et al. (Reference Arslan, Fantuzzi, Craske and Wynn2021a), the terms in the brackets of (3.2) becomes
\begin{align} \mathcal{S}\{\boldsymbol{u},\vartheta\}:=\unicode{x2A0F}_{\varOmega} \frac{a}{R} |\boldsymbol{\nabla} \boldsymbol{u}|^2 + \beta |\boldsymbol{\nabla} \vartheta|^2 - (a+1+b\eta(z) - \varphi'(z))w\vartheta - \varphi'(z)\partial_z\vartheta + U \, \textrm{d}\kern0.07em\boldsymbol{x}. \end{align}The positivity of (3.3) is demonstrated by first exploiting horizontal periodicity and using the Fourier series
\begin{equation} \begin{bmatrix} \vartheta(x,y,z)\\ \boldsymbol{u}(x,y,z) \end{bmatrix}= \sum_{\boldsymbol{k}} \begin{bmatrix} \hat{\vartheta}_{\boldsymbol{k}}(z)\\ \hat{\boldsymbol{u}}_{\boldsymbol{k}}(z) \end{bmatrix}\textrm{e}^{{\rm i}(k_x x + k_y y)}, \end{equation}
where the sum is over suitable wavenumbers 
$\boldsymbol {k} = (k_x,k_y)$ with magnitude 
$\boldsymbol {k} = \sqrt {k_x^2+ k_y^2}$ and where the complex-valued Fourier amplitudes satisfy the complex conjugate relations 
$\hat {\boldsymbol {u}}_{-\boldsymbol {k}} =\hat {\boldsymbol {u}}^*_{\boldsymbol {k}}$ and 
$\hat {T}_{-\boldsymbol {k}} =\hat {T}^*_{\boldsymbol {k}}$. Then, 
$\mathcal {S}\{\boldsymbol {u},\vartheta \}$ can be lower bounded by
\begin{equation} \mathcal{S}\{\boldsymbol{u},\vartheta\} \geq \mathcal{S}_0\{\hat{\vartheta}_{0}\} + \sum_{k} \mathcal{S}_{\boldsymbol{k}}\{\hat{w}_0, \hat{\vartheta}_0\}, \end{equation}where
\begin{equation} \mathcal{S}_0\{\hat{\vartheta}_0\} = U+ \int^{1}_0 b |\hat{\vartheta}'_0|^2 - \varphi'\hat{\vartheta}'_0 \,\textrm{d}z \end{equation}and
\begin{align} \mathcal{S}_{\boldsymbol{k}}\{\hat{w}_{\boldsymbol{k}},\hat{\vartheta}_{\boldsymbol{k}}\}&= \frac{a}{R\,k^2} \lVert \hat{w}''_{\boldsymbol{k}} \rVert_2^2 + \frac{2a}{R} \lVert \hat{w}'_{\boldsymbol{k}} \rVert_2^2 + \frac{ak^2}{R}\lVert \hat{w}_{\boldsymbol{k}} \rVert_2^2 + b \lVert \hat{\vartheta}'_{\boldsymbol{k}}\rVert_2^2 + b k^2\lVert \hat{\vartheta}_{\boldsymbol{k}} \rVert_2^2 \nonumber\\ &\quad - \int^{1}_0 (a + 1 + b \eta(z) - \varphi'(z) )\,\textrm{Re}\{\hat{w}_{\boldsymbol{k}} \hat{\vartheta}^*_{\boldsymbol{k}} \}\, \textrm{d}z \geq 0. \end{align}
Standard arguments (Arslan et al. Reference Arslan, Fantuzzi, Craske and Wynn2021a,Reference Arslan, Fantuzzi, Craske and Wynnb) imply that one can assume that 
$\hat {w}_{\boldsymbol {k}}$ and 
$\hat {\vartheta }_{\boldsymbol {k}}$ to be real-valued functions.
Ensuring that 
$\mathcal {S}_{\boldsymbol {k}} \geq 0$ for any functions 
$\hat {w}_{\boldsymbol {k}}(z), \hat {\vartheta }_{\boldsymbol {k}}(z)$ satisfying the boundary conditions
$$\begin{gather} \hat{w}_{\boldsymbol{k}}(0) = \hat{w}_{\boldsymbol{k}}(1)= \hat{w}_{\boldsymbol{k}}'(0)= \hat{w}_{\boldsymbol{k}}'(1)=0, \end{gather}$$
$$\begin{gather}\hat{\vartheta}_{\boldsymbol{k}}'(0) = \hat{\vartheta}_{\boldsymbol{k}}'(1) =0, \end{gather}$$
is commonly referred to as a spectral constraint at wavenumber 
$\boldsymbol {k}$. The right-hand side of (3.5) is non-negative if the spectral constraint is satisfied for each wavenumber 
$\boldsymbol {k} \in \mathbb {N}_{0}$. This, in turn, implies the desired property that 
$\mathcal {S}(\boldsymbol {u},\vartheta ) \geq 0$ for any solution to the governing equations. It then follows from (3.2) that 
$\langle wT \rangle \leq U$ for any solution of the governing equations, where 
$U$ is the scalar appearing in (3.6) for the spectral constraint at wavenumber 
$\boldsymbol {k}=0$.
The simple structure of the spectral constraint at wavenumber 
$\boldsymbol {k}=0$ can give useful information on the upper bound 
$U$ that can be proven using the auxiliary function method. In particular, if 
$a,b,\varphi (z)$ can be chosen so that the spectral constraints 
$\mathcal {S}_{\boldsymbol {k}}\{\hat {w}_{\boldsymbol {k}},\hat {\vartheta }_{\boldsymbol {k}}\} \geq 0$ holds for all non-zero wavenumbers, then it is possible to prove a bound and the numerical value of the best provable bound 
$U$ can be found by minimising (3.6) over functions 
$\hat {\vartheta }_{\boldsymbol {0}}$ that satisfy the boundary conditions (3.8b) at 
$\boldsymbol {k}=0$. This gives a best provable bound (after fixing 
$\varphi (z),a$ and 
$b$) of
\begin{equation} U = \int^{1}_{0}\frac{\varphi'(z)^2}{4b}\,\textrm{d}z. \end{equation}
Although the parameter 
$a$ does not appear in (3.9), there is an implicit interplay between 
$a$ and 
$\varphi (z)$ since, given (3.7), both must be chosen appropriately if the spectral constraint at non-zero wavenumbers is to be satisfied. This is clarified in § 3.2.
Finding a good bound within the auxiliary function method then corresponds to identifying 
$a,b$ and 
$\varphi (z)$, which satisfy the spectral constraints for positive wavenumbers while keeping the value of 
$U$ given by (3.9) as small as possible.
3.2. Enforcing the spectral constraint for non-zero wavenumbers
Recall that the cumulative heating function 
$\eta (z) \in L^2(0,1)$ is given by (2.4) and that 
$\eta (z)$ is an increasing function satisfying 
$\eta (0)=0$ and 
$\eta (1)=1$. The aim now is to find parameters 
$a,b, \varphi (z)$ such that the spectral constraints 
$\mathcal {S}_{\boldsymbol {k}}\{\hat {w}_{\boldsymbol {k}},\hat {\vartheta }_{\boldsymbol {k}}\} \geq 0$ hold for every non-zero wavenumber 
$\boldsymbol {k}$. To this end, we choose the background field as
\begin{equation} \varphi'(z) = \frac{(a+1)}{\lVert \eta \rVert_2} \begin{cases} 0, & 0\leq z \leq \delta, \\ \eta(z) + \lVert \eta \rVert_2, & \delta \leq z \leq 1-\delta,\\ 0, & 1-\delta \leq z \leq 1, \end{cases} \end{equation}
where 
$\delta \leq \frac 12$ is the width of the boundary layers of 
$\varphi '(z)$. We also set
\begin{equation} b = \frac{a+1}{\lVert \eta \rVert_2}. \end{equation} The motivation for these parameter choices is as follows. The background field 
$\varphi (z)$ is chosen such that in the ‘bulk region’ 
$\delta \leq z \leq 1-\delta$, the term 
$(a+1+b\eta (z) - \varphi '(z))$ appearing in the final, sign-indefinite, integral of (3.7) is identically equal to zero. Consequently, to verify the spectral constraint 
$\mathcal {S}_{\boldsymbol {k}}\{\hat {w}_{\boldsymbol {k}},\hat {\vartheta }_{\boldsymbol {k}}\} \geq 0$ we only need to estimate the final integral term in (3.7) for integrals over only the ‘boundary layers’ 
$0 \leq z \leq \delta$ and 
$1-\delta \leq z \leq 1$. As we show below this can be achieved by using standard functional estimates that exploit the boundary conditions (3.8a) and (3.8b).
In the boundary layers, for simplicity, 
$\varphi '(z)$ is chosen to be zero to satisfy the assumed boundary conditions on 
$\varphi (z)$ and give no detrimental contribution to the bound 
$U$ in (3.9). Once these choices have been made, the key parameter is the boundary layer width 
$\delta$, and there is a natural tension in that a large value of 
$\delta$ improves (reduces) the bound 
$U$ but makes verification of the spectral constraint increasingly difficult (see (3.16)).
Finally, the motivation for the choice of 
$b$ comes from the fact that for the case of uniform heating (where 
$\eta (z)=z$ and 
$\| \eta \|_2 = 1/\sqrt {3}$) previous work (Arslan et al. Reference Arslan, Fantuzzi, Craske and Wynn2021a) reveals that as 
$R$ increases 
$a$ decreases to zero whereas 
$b\rightarrow \sqrt {3}$. The choice (3.11) respects this relation and maintains algebraic simplicity in the following derivations.
We begin by formulating a simpler sufficient condition for the spectral constraint (3.7) to hold. First, we drop the positive terms in 
$\hat {w}_{\boldsymbol {k}}'(z)$, 
$\hat {w}_{\boldsymbol {k}}(z)$ and 
$\hat {\vartheta }'_{\boldsymbol {k}}(z)$. Then, we plug in our choices (3.10) for 
$\varphi '(z)$ and (3.11) for 
$b$, and use the inequality 
$\eta (z)\leq 1$ to conclude that (3.7) holds if
\begin{equation} \frac{a}{R k^2} \lVert \hat{w}''_{\boldsymbol{k}} \rVert_2^2 + \frac{a+1}{\lVert \eta \rVert_2} k^2 \lVert \hat{\vartheta}_{\boldsymbol{k}}\rVert_2^2 - \frac{a+1}{\lVert \eta \rVert_2} (1+\lVert \eta \rVert_2) \int_{(0,\delta)\cup(1-\delta,1)} |\hat{w}_{\boldsymbol{k}} \hat{\vartheta}_{\boldsymbol{k}}| \,\textrm{d}z \geq 0. \end{equation}
Note that the boundary layers 
$(0,\delta )$ and 
$(1-\delta,1)$ are disjoint because 
$\delta \leq \tfrac 12$ by assumption.
Next, given the boundary conditions on 
$w_{\boldsymbol {k}}(z)$ from (3.8a), the use of the fundamental theorem of calculus and the Cauchy–Schwarz inequality gives
\begin{equation} \hat{w}_{\boldsymbol{k}}(z) = \int^{z}_0 \int^{\xi}_0 \hat{w}_{\boldsymbol{k}}''(\sigma) \,\textrm{d}\sigma \,\textrm{d}\xi \leq \int^{z}_0 \sqrt{\xi} \, \lVert \hat{w}_{\boldsymbol{k}}'' \rVert_{2} \,\textrm{d}\xi = \frac23 z^{3/2}\, \lVert \hat{w}_{\boldsymbol{k}}'' \rVert_2. \end{equation}Using (3.13) and the Cauchy–Schwarz inequality, we can estimate
\begin{align} \int_0^\delta | \hat{w}_{\boldsymbol{k}} \hat{\vartheta}_{\boldsymbol{k}}| \,\textrm{d}z &\leq \int_0^\delta |\hat{w}_{\boldsymbol{k}}| |\hat{\vartheta}_{\boldsymbol{k}}| \,\textrm{d}z \leq \frac23 \lVert \hat{w}_{\boldsymbol{k}}'' \rVert_2\int_0^\delta z^{3/2} |\hat{\vartheta}_{\boldsymbol{k}}| \,\textrm{d}z\nonumber\\ &\leq \frac23 \lVert \hat{w}_{\boldsymbol{k}}'' \rVert_2 \lVert \hat{\vartheta}_{\boldsymbol{k}} \rVert_2\left( \int^{\delta}_0 z^3 \,\textrm{d}z \right)^{1/2} \nonumber\\ &=\frac{\delta^2}{3} \lVert \hat{w}_{\boldsymbol{k}}'' \rVert_2 \lVert \hat{\vartheta}_{\boldsymbol{k}} \rVert_2 . \end{align}
The same estimate applies to 
$\int _{1-\delta }^1 |\hat {w}_{\boldsymbol {k}} \hat {\vartheta }_{\boldsymbol {k}}| \,\textrm {d}z$, so we conclude that (3.12) holds if
\begin{equation} \frac{a}{R k^2} \lVert \hat{w}''_{\boldsymbol{k}} \rVert_2^2 + \frac{a+1}{\lVert \eta \rVert_2} k^2 \lVert \hat{\vartheta}_{\boldsymbol{k}}\rVert_2^2 - \frac{2(a+1)}{3 \lVert \eta \rVert_2} (1+\lVert \eta \rVert_2 ) \delta^2 \lVert \hat{w}''_{\boldsymbol{k}} \rVert_2\lVert \hat{\vartheta}_{\boldsymbol{k}} \rVert_2 \geq 0. \end{equation}
The left-hand side of this inequality is a homogeneous quadratic form in 
$\lVert \hat {w}_{\boldsymbol {k}}'' \rVert _2$ and 
$\lVert \hat {\vartheta }_{\boldsymbol {k}}\rVert _2$ and it is non-negative if it has a non-positive discriminant. Therefore, the spectral constraint is satisfied if
\begin{equation} \delta^4 \leq \frac{ 9 a \lVert \eta \rVert_2}{ (a+1) (1+\lVert \eta \rVert_2)^2 R }. \end{equation}3.3. Estimating the bound
Next, we estimate the bound 
$U$ in (3.9) given our ansätze (3.11) and (3.10). Making use of the fact that 
$\eta (z)\geq 0$ gives, after rearranging,
\begin{align} U &=\frac{ a+1 }{4 \lVert \eta \rVert_2 } \int^{1-\delta}_{\delta} (\eta(z) + \lVert \eta \rVert_2)^2 \,\textrm{d}z \nonumber\\ &= \frac{ a+1 }{4 \lVert \eta \rVert_2 } \int^{1}_{0}(\eta(z) + \lVert \eta \rVert_2)^2 \,\textrm{d}z - \frac{ a+1 }{4 \lVert \eta \rVert_2 } \int_{(0,\delta) \cup (1-\delta,1)} (\eta(z) + \lVert \eta \rVert_2)^2 \,\textrm{d}z \nonumber\\ &\leq \frac12 (a+1) ( \lVert \eta \rVert_2 + \lVert \eta \rVert_1 - \lVert \eta \rVert_2 \delta ). \end{align}
By definition, 
$a$ is real and positive, which means that 
$(a+1)^m\geq 1$ for all 
$m>0$. Therefore, taking 
$\delta$ as large as possible in (3.16) and substituting into (3.17) gives
\begin{align} U &\leq \frac12 (a+1) ( \lVert \eta \rVert_2 + \lVert \eta \rVert_1) - \frac{\sqrt{3}}{2} \frac{\lVert\eta\rVert_2^{5/4}}{(1+\lVert\eta\rVert_2)^{1/2}}(a+1)^{3/4} a^{1/4} R^{{-}1/4}\nonumber\\ &\leq \frac12 (a+1) ( \lVert \eta \rVert_2 + \lVert \eta \rVert_1) - \frac{\sqrt{3}}{2} \frac{\lVert\eta\rVert_2^{5/4}}{(1+\lVert\eta\rVert_2)^{1/2}} a^{1/4} R^{{-}1/4} . \end{align}
We now minimise the last expression in (3.18) over 
$a$ to obtain the best possible estimate on 
$U$. The optimal choice of 
$a$ is
\begin{equation} a^* = \frac{a_0 \lVert \eta \rVert_2^{5/3}}{ (\lVert \eta \rVert_2 + \lVert \eta \rVert_1)^{4/3} (1+\lVert \eta \rVert_2 )^{2/3}} R^{{-}1/3}, \end{equation}
where 
$a_0 = (\frac {9}{256})^{1/3}$. Finally, substituting for 
$a$ in (3.18), rearranging and since 
${\left \langle {w\vartheta } \right \rangle } \leq U$, we get
\begin{equation} {\left\langle {w \vartheta} \right\rangle} \leq \frac12 ( \lVert \eta \rVert_2 + \lVert \eta \rVert_1) - \frac{a_1 \lVert \eta \rVert_2^{5/3}}{(\lVert \eta \rVert_2 + \lVert \eta \rVert_1)^{1/3}(1+\lVert \eta \rVert_2)^{2/3} } R^{{-}1/3}, \end{equation}
where 
$a_1 = \frac {3\sqrt {3}}{8} a_0^{1/4}$. Then, by (2.5),
\begin{equation} \smash{{\langle {\delta T} \rangle_h}} \geq{-}\underbrace{\frac12 ( \lVert \eta \rVert_2 - \lVert \eta \rVert_1)}_{\mu} + \underbrace{\frac{a_1 \lVert \eta \rVert_2^{5/3}}{(\lVert \eta \rVert_2 + \lVert \eta \rVert_1)^{1/3}(1+\lVert \eta \rVert_2)^{2/3} }}_{\sigma} R^{{-}1/3}, \end{equation}
completing the proof of (1.1) for the 
$\eta$-dependent values of 
$\mu$ and 
$\sigma$ indicated in (3.21).
Finally, we must also check that requirement that the choice of 
$\delta$ in the above argument via (3.16) satisfies the constraint 
$\delta \leq \frac 12$ to ensure that the two boundary layers do not interact. This gives a minimum Rayleigh number 
$R=R_m$ above which the final bound (3.21) holds. With 
$\delta$ defined by (3.16), it follows that
\begin{equation} R_m = \frac{18 \lVert \eta \rVert_2^{2} }{(\lVert \eta \rVert_1 + \lVert \eta \rVert_2)(1+ \lVert \eta \rVert_2)^2}. \end{equation} In practice, this does not place a strong restriction on our results. For example, consider the case of uniform heating where 
$H(z)=1$. Here, one has 
$\eta (z)=z$ and obtains from (3.21) the lower bound 
$\smash {{\langle {\delta T} \rangle _h}} \geq -0.03868 + 0.1850 R^{-1/3}$, valid for all 
$R\geq R_m =2.2384$ which is far below the critical Rayleigh number for linear instability of 
$R=1440$. The bound matches the scaling and asymptote of Arslan et al. (Reference Arslan, Fantuzzi, Craske and Wynn2021a). The difference in the constant multiplying the 
$O(R^{-1/3})$ term is due to the choice taken in Arslan et al. (Reference Arslan, Fantuzzi, Craske and Wynn2021a) of a 
$\varphi '(z)$ with linear (rather than constant) boundary layers. To obtain as simple a bound as possible in (3.21), the constant multiplying the 
$R$ scaling was not optimised in this work.
4. Implications of the bound
Inequality (3.21) states that 
$\smash {{\langle {\delta T} \rangle _h}} < 0$ (i.e. downward conduction) is impossible for Rayleigh numbers 
$R$ satisfying
\begin{equation} R_m \leq R \leq R_0 := \frac{8a_1^3 \lVert \eta \rVert_2^5 }{(1+\lVert \eta \rVert_2)^2(\lVert \eta \rVert_2 + \lVert \eta \rVert_1)(\lVert \eta \rVert_2 - \lVert \eta \rVert_1)^3} . \end{equation}
The above expression for 
$R_0$ is obtained simply by setting the right-hand side of (3.21) to zero. The magnitude of 
$R_0$ is plotted as a function of 
$\lVert \eta \rVert _1$ and 
$\lVert \eta \rVert _2$ in figure 3. The domain of the function shown in figure 3 respects the fact that
\begin{equation} \|\eta\|_1 \leq \|\eta\|_2 \leq \sqrt{ \|\eta\|_1} , \end{equation}
for any increasing function 
$\eta (z)$ satisfying 
$\eta (0)=0$ and 
$\eta (1)=1$. The lower bound follows from the Cauchy–Schwarz inequality, while the upper bound is a simple consequence of 
$\eta (z) \leq 1$ for any 
$0 \leq z \leq 1$.
Figure 3. (a) Contour plot of 
$R_0$ as given in (4.1) with the inset (b) highlighting a region (
$[0.925,1]\times [0.95,1]$) where 
$\lVert \eta \rVert _1 \sim \lVert \eta \rVert _2$. In (a) the green cross highlights uniform heating, drawn in (c), the blue cross highlights heating near the lower boundary drawn in (d), the yellow cross highlights heating near the upper boundary (e) and the black cross (
$\times$) highlights sinusoidal heating, ( f).
An important implication of these inequalities is that 
$\mu = -\frac 12 (\|\eta \|_2 - \|\eta \|_1)$ is strictly negative unless 
$\|\eta \|_2 = \|\eta \|_1$. When 
$\mu$ is strictly negative, the corresponding value of 
$R_0$ is finite, meaning that there is a maximum Rayleigh number (indicated in figure 3) up to which downward conduction can be ruled out. Given this observation, it is natural to ask whether there are fixed heating profiles for which 
$\|\eta \|_2=\|\eta \|_1$, then one could provably eliminate the possibility of downward conduction for arbitrarily large Rayleigh numbers. Interestingly, given the restrictions that 
$\eta (0)=0$ and 
$\eta (1)=1$, this is only possible in the limiting case that 
$H(z)$ is a single Dirac measure of unit heating at the lower boundary 
$z=0$.
The inequality in (3.21) applies to all heating profiles that can be specified as a function of 
$z$ and integrated to unity. Figure 3(a) is a contour plot of the critical 
$R_0$ over the open region satisfying the inequalities (4.2). Examples of non-uniform heating profiles from figure 2 appear in figure 3 for cases in which the heating is (d) near the bottom boundary, (e) near the upper boundary and ( f) sinusoidal, with the corresponding value of 
$R_0$ highlighted with crosses in figure 3(a). The profiles (c)–( f) and their corresponding 
$R_0$ can be compared with physical examples of IHC, but before doing so, we note that the precise values of 
$R_0$ computed in figure 3 are not as important as the relative magnitudes of 
$R_0$.
The bounds we obtain can help justify physical observations. In the case of mantle convection, it is difficult to accurately determine the distribution of the radioactive heating sources is unknown, and the Rayleigh number (but is believed to be of the order 
$10^6$ to 
$10^8$ (Schubert Reference Schubert2015)). However, a higher concentration of the heat in the upper mantle, in comparison with a uniform distribution, will give rise to smaller 
$R_0$ and a smaller magnitude of 
$\smash {{\langle {\delta T} \rangle _h}}$ that would align with the observation of the slower-than-expected cooling rate of the Earth's interior (Jaupart et al. Reference Jaupart, Labrosse, Lucazeau and Mareschal2015). Alternatively, convection in the Sun resembles heating similar to (d), of distributed heating near the lower regions of the convective zone. Heating concentrated at the bottom of the domain, as seen in figure 3, results in a value of 
$R_0$ much larger than the case of uniform heating. Moreover, for estimated Rayleigh numbers of 
$10^{22}$ (Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020), downwards conduction can be ruled out for heating concentrated sufficiently closely to the lower boundary.
Finally, we demonstrate an application of the bounds obtained in this work. Instead of starting with a fixed 
$H(z)$ and asking for the corresponding 
$R_0$, we can instead take a given Rayleigh number, motivated by a geophysical or astrophysical flow, and for a family of 
$H(z)$ (figure 3d to f are different one-parameter families of heating profiles) asses which heating profiles in that family guarantee that downwards conduction is ruled out. For example, taking figure 3(c) where heating is uniform in a region 
$(0,\epsilon )$ where 
$\epsilon \leq 1$, by our result in (4.1) we can obtain a value for 
$R_0$ as 
$\epsilon$ decreases, as plotted in figure 4(b). Then, for example, taking a Rayleigh number of 
$10^{22}$, our bound (3.21) indicates that downward conduction is impossible when the heating profile is such that 
$\epsilon < 1.366 \times 10^{-7}$ figure 4(b). For the Sun, where the convective zone is estimated to be of the order of 
$0.713 R_{\odot }$ (where 
$R_{\odot }$ is the solar radius) (Miesch Reference Miesch2005), this corresponds to heating in a region of 6.8 cm. Note that 
$10^{22}$ is an estimated Rayleigh number that arises from including boundary and internal heating in the definition of 
$R$. For a well-chosen 
$R$, (3.21) has direct implications on physical flows, especially if one can establish the optimal constants in the 
$R$ scaling of (3.21). This exercise highlights the use of bounds to assess cases of IHC in nature where downward conduction is ruled out.
Figure 4. Panel (a) plots heating profiles 
$H(z)$, with mean one, where the heating is uniform in a region 
$(0,\epsilon )$ near the lower boundary and zero elsewhere. In (a) eight cases are plotted of 
$\epsilon \in [0.05,0.5]$, where 
$\epsilon = 0.05$ is in yellow (yellow solid line) and 
$\epsilon =0.5$ in black (black solid line). Panel (b) plots 
$\epsilon$ against 
$R_0$, the value below which downwards conduction is ruled out. The vertical dashed line (- - -) corresponds to 
$\epsilon = 1.366 \times 10^{-7}$ intersecting the red line at 
$R_0=10^{22}$. For a given Rayleigh number (i.e. 
$10^{22}$) (b) can be used to find the value of 
$\epsilon$ for heating profiles in (a) below which downwards conduction is impossible.
5. Conclusion
We have obtained lower bounds on the mean vertical conductive heat transport, 
$\smash {{\langle {\delta T} \rangle _h}}$, for non-uniform IHC between parallel plates with fixed-flux thermal boundary conditions. The bounds are expressed in terms of a Rayleigh number and the mean vertical heat flux 
$\eta (z)$, corresponding to the prescribed internal heating 
$H(z)$ integrated upwards from the bottom boundary to a height 
$z$. In particular, we proved with the auxiliary functional method that 
$\smash {{\langle {\delta T} \rangle _h}} \geq -\mu + \sigma R^{-1/3}$ in (3.21), where 
$\mu, \sigma >0$. We were, therefore, able to obtain an explicit Rayleigh number 
$R_0$, in terms of 
$\eta (z)$, below which it is possible to guarantee that 
$\smash {{\langle {\delta T} \rangle _h}} \geq 0$. The results are consistent with those obtained from the particular case of uniform internal heating (
$\eta =z$) (Arslan et al. Reference Arslan, Fantuzzi, Craske and Wynn2021a). As an addendum, while this paper only considers no-slip boundary conditions, by the arguments in Fantuzzi (Reference Fantuzzi2018), the same result of (3.21), albeit with a different value of 
$a_1$, is obtained for stress-free boundary conditions.
Only in the limit of heat being concentrated entirely at the bottom boundary is it possible to guarantee that 
$\smash {{\langle {\delta T} \rangle _h}} \geq 0$ for all 
$R$ (
$\mu =0$ in the previous paragraph). However, in a mathematical sense, such heating profiles are distributions rather than functions of 
$z$, lying beyond the scope of the internal heating profiles considered in this work and constitute the closure, of the open set of pairs 
$(\|\eta \|_{1}, \|\eta \|_{2})$ that determine 
$\mu$ and 
$\sigma$ (shown in figure 3). In this regard, there is an intriguing connection with RBC, for which 
$\smash {{\langle {\delta T} \rangle _h}} \geq 0$ for all 
$R$ (Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002).
While the auxiliary functional method provides a rigorous bound on 
$\smash {{\langle {\delta T} \rangle _h}}$, it does not guarantee the existence of flows that saturate the bound or provide insight into the predominant physical processes responsible for extremising 
$\smash {{\langle {\delta T} \rangle _h}}$. In other words, the question of whether there are solutions to the Boussinesq equations that produce negative mean vertical conductive heat transport (
$\smash {{\langle {\delta T} \rangle _h}} <0$) when 
$R>R_0$ remains an open and interesting problem to investigate in the future. A possible avenue for exploring this question is the method of optimal wall-to-wall transport (Tobasco & Doering Reference Tobasco and Doering2017), which would yield a velocity field that minimises 
$\smash {{\langle {\delta T} \rangle _h}}$ subject to a prescribed set of constraints, along with upper bounds on 
$\smash {{\langle {\delta T} \rangle _h}}$.
Finally, caution should be applied in interpreting the Rayleigh number used to study non-uniform IHC. Our Rayleigh number follows convention by defining temperature and length scales using the volume-averaged internal heating and domain height. However, in cases with heating concentrated towards the top of the domain, the physical interpretation of the Rayleigh number as quantifying the size of destabilising forces relative to the stabilising effects of diffusion is questionable. Although the choice of Rayleigh number is somewhat superficial, an appropriate choice might shed light on the dependence of the bounds on 
$R$ in terms of the underlying physics. Heating and cooling that is applied non-uniformly over the horizontal as well as vertical directions (Song et al. Reference Song, Fantuzzi and Tobasco2022) poses further challenges in this regard, and an opportunity to extend the results reported in this work.
Funding
A.A. acknowledges funding from the Engineering and Physical Sciences Research Council (EPSRC) Centre for Doctoral Training in Fluid Dynamics across Scales (award no. EP/L016230/1), the European Research Council (agreement no. 833848-UEMHP) under the Horizon 2020 program and the Swiss National Science Foundation (grant no. 219247) under the MINT 2023 call. J.C. acknowledges funding from the EPSRC (grant no. EP/V033883/1).
Declaration of interests
The authors report no conflict of interest.
Open access
