1. Introduction
Rayleigh–Bénard convection (RBC) is the buoyancy-driven flow in a fluid layer heated from below and cooled from above in the presence of gravity. The emergent convective flow enhances heat flux from the warm bottom boundary to the cool top boundary beyond the conductive flux from diffusion alone. This dimensionless enhancement factor – the ratio of bulk-averaged vertical heat flux from both conduction and convection to the flux from conduction alone – defines the Nusselt number 
$Nu$. In Rayleigh's mathematical model (Rayleigh Reference Rayleigh1916) 
$Nu$ depends on several dimensionless quantities characterizing the problem at hand: (i) what we now call the Rayleigh number 
$Ra$, which is proportional to the imposed temperature drop across the layer, (ii) the fluid's Prandtl number 
$Pr$, which is the ratio of kinematic viscosity to thermal diffusivity and (iii) details of the spatial domain, often captured by an aspect ratio 
$\varGamma$ that is a ratio of a horizontal length scale to the vertical layer height.
Convection is coherent at 
$Ra$ values not too far above the critical value 
$Ra_c$ beyond which the conductive no-flow state is linearly unstable. By coherent we mean flows with few scales present; spatial scales might include a horizontal period and the vertical thickness of boundary layers, and temporally the flow may be steady or time-periodic. Meanwhile, convection is turbulent at the large 
$Ra$ values pertinent to many engineering and scientific applications. Turbulent flows are complex and contain a range of spatial and temporal scales and, in the present context, have thermal and viscous boundary layers at the top and bottom boundaries from which thermal plumes emerge and mix the bulk. In a given domain it is expected that a scaling of 
$Nu$ with respect to both 
$Pr$ and 
$Ra$ will emerge in the 
$Ra\to \infty$ limit (Kadanoff Reference Kadanoff2001).
After nearly a century of increasingly sophisticated mathematical analysis, increasingly resolved direct numerical simulations (DNS) and increasingly refined laboratory experiments, two quantitatively distinct conjectures remain in contention for the heat transport scaling law at large 
$Ra$ (Chillà & Schumacher Reference Chillà and Schumacher2012; Doering Reference Doering2020). The two conjectures follow from heuristic physical arguments that both seem plausible but give incompatible predictions: the ‘classical’ scaling 
$Nu \sim Pr^0 Ra^{1/3}$ and the ‘ultimate’ scaling 
$Nu \sim Pr^{1/2} Ra^{1/2}$, with the latter sometimes including logarithmic-in-
$Ra$ modifications.
For RBC between flat, no-slip, isothermal boundaries, rigorous analysis of the governing equations has yielded upper bounds of the form 
$Nu\le {O}(Ra^{1/2})$ uniformly in 
$Pr$ and 
$\varGamma$ (Howard Reference Howard1963; Doering & Constantin Reference Doering and Constantin1996), but this still allows for either classical or ultimate scaling. Upper bounds that rule out ultimate scaling by being asymptotically smaller than 
${{O}}(Ra^{1/2})$ have been derived in the limit of infinite 
$Pr$ (Doering, Otto & Reznikoff Reference Doering, Otto and Reznikoff2006; Otto & Seis Reference Otto and Seis2011; Whitehead & Doering Reference Whitehead and Doering2012) and for two-dimensional convection between stress-free boundaries (Whitehead & Doering Reference Whitehead and Doering2011). For the no-slip boundaries relevant to experiments, however, it remains an open question whether an upper bound asymptotically smaller than 
$Ra^{1/2}$ is possible.
In view of the problem's stubbornness, a new strategy is called for to determine – or at least to bound – 
$Nu$ as a function of 
$Ra, Pr$ and 
$\varGamma$. Towards that end we have undertaken an indirect approach consisting of two parts. The first part is to study coherent flows for which one can reasonably hope to determine asymptotic heat transport, and the second part is to investigate how transport by those coherent flows compares with transport by turbulent convection. The simplest coherent flows are steady – i.e. time-independent – solutions of the equations of motion. Many such states exist, although they are generally unstable at large 
$Ra$. We focus on what might be called the simplest type of steady states: two-dimensional convection rolls like the counter-rotating pairs shown in figure 1(
$a$,
$b$). In horizontally periodic or infinite domains in two or three dimensions, such rolls bifurcate supercritically from the conductive state in the linear instability identified by Rayleigh (Reference Rayleigh1916). A roll pair of any width-to-height aspect ratio 
$\varGamma$ admitted by the domain exists for sufficiently large 
$Ra$.
Figure 1. Steady convection rolls at 
$Ra=10^9$ and 
$Pr=1$ with (
$a$) 
$\varGamma = 2$ and (
$b$) the value 
$\varGamma ^*\approx 0.235$ that maximizes 
$Nu$ at these values of 
$Ra$ and 
$Pr$. Colour indicates temperature, and streamlines are shown for anticlockwise (solid) and clockwise (dash-dotted) motions. (
$c$) Dependence of 
$Nu$ on the horizontal wavenumber 
$k=2{\rm \pi} /\varGamma$ found by computing steady rolls of various aspect ratios (
${\triangle }$, green). Highlighted points are the two 
$\varGamma$ shown in (
$a$,
$b$) along with the value 
$\varGamma ^*_{loc}\approx 0.0614$ that locally maximizes 
$Nu(\varGamma )$. Cubic spline interpolation (dashed line) is used to find 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$ precisely.
For steady rolls, the dependence of 
$Nu$ on the parameters 
$(\varGamma,Pr, Ra)$ at asymptotically large 
$Ra$ is accessible to computation. As for whether heat transport by steady rolls can be connected to transport by turbulence, there are several reasons for optimism. Relationships between turbulent attractors and the unstable coherent states embedded therein have been established in models of wall-bounded shear flows (Graham & Floryan Reference Graham and Floryan2021), where particular steady states, travelling waves and time-periodic states have been found that closely reflect turbulent flows in terms of integral quantities as well as particular flow structures. Analogous study of RBC began only recently but indeed suggests that certain steady states capture qualitative aspects of turbulent convection (Sondak, Smith & Waleffe Reference Sondak, Smith and Waleffe2015; Waleffe, Boonkasame & Smith Reference Waleffe, Boonkasame and Smith2015; Kooloth, Sondak & Smith Reference Kooloth, Sondak and Smith2021; Motoki, Kawahara & Shimizu Reference Motoki, Kawahara and Shimizu2021). Our findings add to this evidence. The desire to understand and perhaps strengthen the mathematical bound 
$Nu\le {O}(Ra^{1/2})$ is further motivation for studying unstable states since bounds apply to all solutions of the governing equations regardless of stability. It is an open question whether any solutions can achieve ultimate scaling, let alone turbulent solutions.
Here we report numerical computations of steady convection rolls for a 
$Pr=1$ fluid contained between no-slip isothermal top and bottom boundaries. We reach sufficiently large 
$Ra$ values to convincingly reveal several asymptotic scalings of 
$Nu$, depending on the horizontal periods of the rolls. These are the first clearly asymptotic scalings found for any type of flow – steady, turbulent or otherwise – for RBC in the no-slip case. Notably, the largest heat transport among steady rolls of all horizontal periods displays the classical 
$Nu\sim Ra^{1/3}$ scaling. We further observe that 
$Nu$ for these steady rolls is larger than turbulent 
$Nu$ from all laboratory experiments and two- or three-dimensional (2-D or 3-D) simulations at comparable parameters. This observation supports the conjecture that steady states maximize 
$Nu$ among all stable or unstable flows, as was recently verified for a truncated model of RBC (Olson et al. Reference Olson, Goluskin, Schultz and Doering2021) using methods that are not yet applicable to the full governing equations. If steady-roll transport continues to dominate turbulent transport as 
$Ra\to \infty$, then our finding of classical scaling for steady rolls would rule out ultimate scaling of turbulent convection.
The asymptotic scaling of steady rolls is already known in the case of stress-free velocity conditions at the top and bottom boundaries, which were considered for mathematical convenience in Rayleigh's original work. In that case 
$Nu\sim Ra^{1/3}$ as 
$Ra\to \infty$ at fixed 
$Pr$ and 
$\varGamma$, and the aspect ratio of the roll pair maximizing 
$Nu$ at each 
$Ra$ and 
$Pr$ approaches 
$\varGamma \approx 1.9$ (Chini & Cox Reference Chini and Cox2009; Wen et al. Reference Wen, Goluskin, LeDuc, Chini and Doering2020). Recent computations of steady rolls in the
no-slip case for preasymptotic 
$Ra$ values up to 
$10^9$ revealed significant differences from the stress-free problem (Sondak et al. Reference Sondak, Smith and Waleffe2015; Waleffe et al. Reference Waleffe, Boonkasame and Smith2015). The dependence 
$Nu(\varGamma )$ for no-slip rolls at fixed 
$Ra$ and 
$Pr$ can have multiple local maxima, as shown in figure 1(
$c$), and the aspect ratio 
$\varGamma ^*$ that globally maximizes 
$Nu(\varGamma )$ approaches zero rather than a constant as 
$Ra\to \infty$. Steady rolls of 
$Nu$-maximizing aspect ratios 
$\varGamma ^*$ were reported in Sondak et al. (Reference Sondak, Smith and Waleffe2015) for 
$Ra \in [5\times 10^6, 3\times 10^8]$ at 
$Pr=1$, yielding fits of 
$\varGamma ^* \sim Ra^{-0.217}$ and 
$Nu(\varGamma ^*) \sim Ra^{0.31}$. This heat transport scaling is faster than with 
$\varGamma$ fixed: computations in Waleffe et al. (Reference Waleffe, Boonkasame and Smith2015) for 
$Ra \in [5\times 10^5, 5\times 10^6]$ at 
$Pr=7$ with 
$\varGamma = 2$ fixed yield the fit 
$Nu\sim Ra^{0.28}$. These best-fit scaling exponents are, however, not asymptotic.
Steady convection rolls are dynamically unstable at large 
$Ra$ and cannot be found by standard time integration, so we employed a purpose-written code that iteratively solves the time-independent equations. We computed rolls with 
$\varGamma =2$ fixed for 
$Ra\lesssim 2\times 10^{10}$ and with the parameter-dependent aspect ratios 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$ (cf. figure 1) that globally and locally maximize 
$Nu(\varGamma )$, respectively, for 
$Ra\le 10^{14}$. These 
$Ra$ values are evidently large enough to reach asymptotia: the results reported below strongly suggest that fixed-
$\varGamma$ rolls asymptotically transport heat like 
$Nu\sim Ra^{1/4}$ while the ever-narrowing rolls of aspect ratio 
$\varGamma ^*$ achieve the classical 
$Nu\sim Ra^{1/3}$ scaling.
2. Computation of steady-convection-roll solutions
Following Rayleigh (Reference Rayleigh1916), we model RBC using the Boussinesq approximation to the Navier–Stokes equations with constant kinematic viscosity 
$\nu$, thermal diffusivity 
$\kappa$ and coefficient of thermal expansion 
$\alpha$. We non-dimensionalize lengths by the layer height 
$h$, temperatures by the fixed difference 
$\varDelta$ between the boundaries, velocities by the free-fall scale 
$U_f = \sqrt {g\alpha h\varDelta }$, and time by the free-fall time 
$h/U_f$. Calling the horizontal coordinate 
$x$ and the vertical coordinate 
$z$, the gravitational acceleration of magnitude 
$g$ is in the 
$-\boldsymbol {\hat z}$ direction. The evolution equations governing the dimensionless velocity vector 
$\boldsymbol {u} = (u, w)$, temperature 
$T$ and pressure 
$p$ are then
\begin{gather} \partial_t u + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla} p + ({Pr}/{Ra})^{1/2} \;\nabla^2 \boldsymbol{u} + T \boldsymbol{\hat z}, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather}\partial_t T + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = (Pr Ra)^{{-}1/2}\; \nabla^2 T, \end{gather}where
\begin{equation} Ra = \frac{g\alpha h^3\varDelta}{\kappa \nu} \quad \mbox{and} \quad Pr = \frac{\nu}{\kappa}. \end{equation}
The dimensionless spatial domain is 
$(x,z)\in [0,\varGamma ]\times [0,1]$, and all variables are horizontally periodic. The top and bottom boundaries are isothermal with 
$T=0$ and 
$T=1$, respectively, while no-slip conditions require 
$\boldsymbol {u}$ to vanish on both boundaries. The conductive state 
$(\boldsymbol {u},T)=(\boldsymbol 0,1-z)$ becomes unstable when 
$Ra$ increases past the critical value 
$Ra_c\approx 1708$ (Jeffreys Reference Jeffreys1928), at which a roll pair with horizontal period 
$\varGamma \approx 2.016$ bifurcates supercritically. As 
$Ra\to \infty$ the horizontal period of the narrowest marginally stable roll pair decreases as 
${O}(Ra^{-1/4})$, while the horizontal period of the fastest-growing linearly unstable mode decreases more slowly as 
${O}(Ra^{-1/8})$.
In terms of the dimensionless solutions to (2.1), the Nusselt number is
\begin{equation} Nu = 1 + (Pr Ra)^{1/2} \, \langle wT\rangle, \end{equation}
where 
$\langle \cdot \rangle$ denotes an average over the spatial domain and infinite time. For steady states no time average is needed.
To compute rolls at 
$Ra$ values large enough to reach the asymptotic regime, we developed a numerical scheme by adapting the approach of Wen et al. (Reference Wen, Goluskin, LeDuc, Chini and Doering2020) and Wen & Chini (Reference Wen and Chini2018) to the case of no-slip boundary conditions. In these numerics the temperature is represented using the deviation 
$\theta$ from the conductive profile, meaning 
$T=1-z+\theta$, and the velocity is represented using a stream function 
$\psi$, where 
$\boldsymbol {u} = \partial _z \psi \boldsymbol {\hat x} - \partial _x \psi \boldsymbol {\hat z}$ so that the (negative) scalar vorticity is 
$\omega =\partial _x w- \partial _z u = -\nabla ^2 \psi$. In terms of these variables, steady (
$\partial _t=0$) solutions of (2.1) satisfy
\begin{gather} \partial_z\psi\partial_x\omega - \partial_x\psi\partial_z\omega = ({Pr}/{Ra})^{1/2}\; \nabla^2 \omega + \partial_x \theta, \end{gather}
\begin{gather}\nabla^2 \psi ={-}\omega, \end{gather}
\begin{gather}\partial_z\psi\partial_x\theta - \partial_x\psi\partial_z\theta ={-}\partial_x\psi + (Pr Ra)^{{-}1/2}\; \nabla^2 \theta \end{gather}with fixed-temperature and no-slip boundary conditions
\begin{equation} \theta|_{z=0, 1}=0, \quad \psi|_{z=0,1}=0 \quad \text{and} \quad \partial_z \psi|_{z=0,1}=0. \end{equation}
The equality of 
$\psi$ values between each boundary restricts to solutions whose horizontal velocities average to zero over every vertical section.
To compute solutions of the time-independent equations (2.4) and (2.5a–c) by an iterative method, we do not need to impose all boundary conditions precisely on each iteration – the conditions need to hold only for the converged solution. Thus we do not impose (2.5c) exactly, instead using approximate boundary conditions on 
$\omega$ for equation (2.4a). These are derived by Taylor expanding 
$\psi$ about the top and bottom boundaries to find
\begin{gather} \psi|_{z=1-\delta} = \psi|_{z=1} - \left.{\partial_z \psi}\right|_{z=1} \delta + \left.{\partial^2_{z} \psi}\right|_{z=1} \frac{\delta^2}{2} - \left.{\partial^3_z \psi}\right|_{z=1} \frac{\delta^3}{6} + {O}(\delta^4), \end{gather}
\begin{gather}\psi|_{z=\delta} = \psi|_{z=0} + \left.{\partial_z \psi}\right|_{z=0} \delta + \left.{\partial^2_z \psi}\right|_{z=0} \frac{\delta^2}{2} + \left.{\partial^3_z \psi}\right|_{z=0} \frac{\delta^3}{6} + {O}(\delta^4), \end{gather}
where 
$\delta >0$ is small. Combining (2.4b) with (2.5b,c) and neglecting 
${O}(\delta ^4)$ terms in (2.6) give the approximate boundary conditions
\begin{align} \partial_z\omega|_{z=1} - \frac{3}{\delta}\ \omega|_{z=1} - \frac{6}{\delta^3}\ \psi|_{z=1-\delta} = 0, \quad -\partial_z\omega|_{z=0} - \frac{3}{\delta}\ \omega|_{z=0} - \frac{6}{\delta^3}\ \psi|_{z=\delta} = 0. \end{align}
In computations we set 
$\delta$ to be the distance between the boundary and the first interior mesh point.
The time-independent equations (2.4) are solved numerically subject to boundary conditions (2.5a,b) and (2.7) using a Newton–GMRES (generalized minimal residual) iterative scheme. The spatial discretization is spectral, using a Fourier series in 
$x$ and a Chebyshev collocation method in 
$z$ (Trefethen Reference Trefethen2000). All of our computations had at least 20 collocation points in the viscous and thermal boundary layers. At 
$Ra$ just above the linear instability, iterations starting from the unstable eigenmode converge to the steady rolls we seek. At larger 
$Ra$, already-computed steady rolls from nearby 
$Ra$ and 
$\varGamma$ values were used as the initial iterate. Every two to four Newton iterations, we change the boundary values of the iterate to match the 
$\partial _z\psi =0$ boundary condition exactly. Prior to convergence this makes the boundary values slightly inconsistent with the governing equations, but the converged solutions satisfy the equations and the no-slip boundary conditions to high precision. Newton iterations were carried out until the Lebesgue 
$L^2$-norm of the residual of the governing steady equations had a relative magnitude less than 
$10^{-10}$. To accurately locate 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$, rolls were computed at several nearby 
$\varGamma$, and then 
$Nu(\varGamma )$ was interpolated with cubic splines like those in figure 1(
$c$). Details of computational results, including resolutions used, are included in the supplementary material available at https://doi.org/10.1017/jfm.2021.1042.
3. Results
We computed steady 
$Pr=1$ rolls for aspect ratios 
$\varGamma$ encompassing the three distinguished values indicated by figure 1(
$c$): the fixed value 
$\varGamma =2$ and the 
$Ra$-dependent values 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$ that globally and locally maximize 
$Nu$ over 
$\varGamma$. As previously observed by Sondak et al. (Reference Sondak, Smith and Waleffe2015), the 
$Nu(\varGamma )$ curve has a single maximum when 
$Ra$ is small and develops a second local maximum at smaller 
$\varGamma$ when 
$Ra$ increases past roughly 
$2\times 10^5$. The value of 
$Nu$ at this second local maximum remains less than the value at the first, so the picture remains as in figure 1(
$c$) with 
$\varGamma ^*$ on the left and 
$\varGamma ^*_{loc}$ on the right, in contrast to the 
$Pr=10$ and 100 cases (Sondak et al. Reference Sondak, Smith and Waleffe2015). For most 
$Ra$ values we did not compute rolls over a full sweep through 
$\varGamma$ as in figure 1(
$c$), instead searching over 
$\varGamma$ only as needed to locate 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$. The rest of this section reports Nusselt number and Reynolds number scalings for the computed steady rolls, and the supplementary material provides tabulated data.
3.1. Asymptotic heat transport
Figure 2 shows the dependence of 
$Nu$ on 
$Ra$ for steady rolls with aspect ratios 
$\varGamma =2$, 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$. In figure 2(a) 
$Nu-1$ is compensated by 
$Ra^{1/3}$, so the horizontal line approached by rolls of the 
$Nu$-maximizing aspect ratios 
$\varGamma ^*$ corresponds to classical 
$1/3$ scaling. The downward slopes of the data for aspect ratios 2 and 
$\varGamma ^*_{loc}$ correspond to scaling exponents smaller than 1/3. Values of 
$Nu$ at 
$\varGamma ^*$ computed previously for 
$Ra\le 10^9$ (Sondak et al. Reference Sondak, Smith and Waleffe2015; Waleffe, personal communication 2020) are shown in figure 2 also, and they agree with our computations very precisely – e.g. the 
$Ra=10^9$ data point agrees with our value of 
$Nu$ to within 0.0008 %.
Figure 2. (a) Compensated plot of 
$Nu-1$ vs 
$Ra$ for steady rolls with 
$Pr=1$ and aspect ratios of 
$\varGamma =2$, 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$, where the 
$Ra$-dependent values 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$ are where 
$Nu(\varGamma )$ has global and local maxima, respectively (cf. figure 1). Values of 
$Nu$ at 
$\varGamma ^*$ from Sondak et al. (Reference Sondak, Smith and Waleffe2015) and Waleffe (personal communication 2020) are also shown (
$\star$, green). Scaling fits (dashed lines) over the last decade of each data set yield exponents of 
$0.33$, 
$0.29$ and 
$0.25$. (b) Finite-difference approximations of the local scaling exponent 
$\beta _n = {\textrm {d}(\log Nu)}/{\textrm {d}(\log Ra)}$. Exponents of 
$1/3$, 
$2/7$ and 
$1/4$ are shown to guide the eye (grey dotted lines).
Figure 2(b) shows the 
$Ra$-dependent local scaling exponent 
$\beta _n = {\mathrm {d}( \log Nu)}/{\mathrm {d} (\log Ra)}$ of the 
$Nu$–
$Ra$ relation for 
$\varGamma =2$, 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$. This quantity educes small variations not visible in figure 2(a). In particular, for rolls of aspect ratios 
$\varGamma ^*$, the exponent 
$\beta _n$ exhibits a small but rapid change just below 
$Ra = 10^9$, beyond which it smoothly approaches the classical 
$1/3$ exponent that appears to be the 
$Ra\to \infty$ asymptotic behaviour. This rapid change seems to coincide with the velocity becoming vertically uniform outside the boundary layers, as reflected in the streamlines of figure 1(
$b$); further details of the rolls’ structure will be reported elsewhere. Rolls with 
$\varGamma =2$ fixed undergo a similarly rapid change around 
$Ra \approx 2 \times 10^7$ and then approach 
$Nu \sim Ra^{1/4}$ scaling that appears to be asymptotic. Rolls of aspect ratio 
$\varGamma ^*_{loc}$ show intermediate 
$Nu$ scaling whose best-fit exponent over the last decade of data is 0.29.
Figure 3(a) shows the 
$Ra$-dependence of the wavenumber 
$k=2{\rm \pi} /\varGamma$ for 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$, compensated by 
$Ra^{1/5}$. The compensated wavenumbers for 
$\varGamma ^*$ approach a horizontal line, suggesting that the 
$Nu$-maximizing rolls narrow according to the power law 
$\varGamma ^*\sim Ra^{-1/5}$. This narrowing of 
$\varGamma ^*$ is slow relative to the case of RBC in a porous medium, where 
$\varGamma ^*\sim Ra^{-1/2}$ (Wen, Corson & Chini Reference Wen, Corson and Chini2015).
Figure 3. (a) Compensated plot of the fundamental horizontal wavenumber 
$k=2{\rm \pi} /\varGamma$ vs 
$Ra$ for the aspect ratios 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$ that maximize 
$Nu(\varGamma )$ globally and locally, respectively, at 
$Pr=1$. Scaling fits (dashed lines) to 
$2{\rm \pi} /\varGamma ^*$ over 
$Ra\in [10^{10},10^{14}]$ and 
$2{\rm \pi} /\varGamma ^*_{loc}$ over 
$Ra\in [10^{13},10^{14}]$ yield exponents of 0.20 and 0.25, respectively. (b) Finite-difference approximations of the local exponent 
$\beta _k = {\textrm {d}(\log k)}/{\textrm {d}(\log Ra)}$. The values 
$1/4$ and 
$1/5$ (grey dotted lines) agree with the scaling fit exponents to two digits.
Figure 3(b) shows the 
$Ra$-dependence of the local scaling exponent 
$\beta _k = {\textrm {d}(\log k)}/{\textrm {d}(\log Ra)}$. For 
$k=2{\rm \pi} /\varGamma ^*$ the local scaling exponent remains close to 
$1/5$ after the transition around 
$Ra = 10^9$. For 
$k=2{\rm \pi} /\varGamma ^*_{loc}$ the exponent seems to approach 
$1/4$, suggesting that 
$\varGamma ^*_{loc}$ has the same 
$Ra^{-1/4}$ scaling as the narrowest marginally stable mode. Variations in 
$\beta _k$ beyond 
$Ra=10^{12}$ for 
$\varGamma ^*$ are evident, but these might be due to numerical imprecision: 
$Nu$ depends very weakly on 
$\varGamma$ around the maximum of 
$Nu(\varGamma )$, as seen in figure 1(
$c$), so the value of 
$\varGamma ^*$ cannot be determined nearly as precisely as the value of 
$Nu(\varGamma ^*)$.
3.2. Asymptotic kinetic energy
Another emergent quantity central to RBC is the bulk Reynolds number based on root-mean-squared velocity, which in terms of dimensionless solutions to (2.1) is
\begin{equation} Re = \left( \frac{Ra}{Pr} \right)^{1/2} \langle | {\boldsymbol u}|^2 \rangle^{1/2}. \end{equation} Figure 4 depicts the dependence of 
$Re$ on 
$Ra$ for the steady rolls of aspect ratios 
$\varGamma = 2$, 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$. Figure 4(a) shows 
$Re$ compensated by 
$Ra^{1/2}$ while figure 4(b) shows the local scaling exponent 
$\beta _r = {\textrm {d}( \log Re)}/{\textrm {d}( \log Ra)}$. Rolls with the fixed aspect ratio 
$\varGamma =2$ approach the asymptotic scaling 
$Re\sim Ra^{1/2}$ that corresponds to the root-mean-squared velocity being proportional to the free-fall velocity 
$U_f$. For rolls with 
$Nu$-maximizing aspect ratios 
$\varGamma ^*$, the scaling fit over 
$Ra\in [10^{10},10^{14}]$ is 
$Re\sim Ra^{0.47}$, which is quite close to the 
$Re \sim Ra^{0.46}$ scaling observed in recent 3-D DNS up to 
$Ra=10^{15}$ at 
$Pr=1$ in a slender cylinder with a height 10 times its diameter (Iyer et al. Reference Iyer, Scheel, Schumacher and Sreenivasan2020). For the 
$\varGamma ^*_{loc}$ rolls, the scaling exponent of 
$Re$ is indistinguishable from 
$2/5$. The measured exponents (0.50, 0.47, 0.40) are unchanged if 
$Re$ is defined using the pointwise maximum velocity rather than using the root-mean-squared velocity as in (3.1). All three aspect ratios result in smaller speeds than steady rolls between stress-free boundaries, where 
$Re\sim Ra^{2/3}$ for any fixed 
$Pr$ and 
$\varGamma$ (Wen et al. Reference Wen, Goluskin, LeDuc, Chini and Doering2020).
Figure 4. (a) Compensated plot of 
$Re$ vs 
$Ra$ for steady rolls with 
$Pr=1$ and aspect ratios 
$\varGamma =2$, 
$\varGamma ^*$ and 
$\varGamma ^*_{loc}$. Scaling fits yield 
$Re\sim Ra^{0.50}$ for 
$\varGamma =2$ and 
$Re\sim Ra^{0.40}$ for 
$\varGamma ^*_{loc}$ over the last decade of each data set, and 
$Re\sim Ra^{0.47}$ for 
$\varGamma ^*$ over 
$Ra\in [10^{10},10^{14}]$ (dashed lines). (b) Finite-difference approximations to the local exponent 
$\beta _r = {\textrm {d}( \log Re)}/{\textrm {d} (\log Ra)}$. Exponents of 
$1/2$ and 
$2/5$ are shown to guide the eye (grey dotted lines).
4. Comparison with turbulent convection
To compare heat transport by steady rolls with that by turbulent thermal convection, we compiled Nusselt number data from high-
$Ra$ DNS with 
$Pr=1$ or 0.7 and laboratory experiments where the estimated 
$Pr$ is between 0.7 and 1.3. Figure 5 shows these 
$Nu$ values compensated by 
$Ra^{1/3}$, along with 
$Nu$ values of steady convection rolls at the 
$Nu$-maximizing aspect ratios 
$\varGamma ^*$. Strikingly, heat transport by the 
$Nu$-maximizing 2-D steady rolls is larger than transport by turbulent convection in all cases.
Figure 5. 
$Nu$ compensated by 
$Ra^{1/3}$ for steady rolls of 
$Nu$-maximizing aspect ratios 
$\varGamma ^*$ at 
$Pr=1$, along with 
$Nu$ from turbulent 2-D and 3-D DNS and experiments with estimated 
$Pr \in [0.7,1.3]$. For horizontally periodic domains, 2-D DNS with 
$(\varGamma,Pr) = (2,1)$ were done by Johnston & Doering (Reference Johnston and Doering2009) and Zhu et al. (Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018), and 3-D DNS with 
$\varGamma \ge 8$ and 
$Pr = 1$ were done by Vieweg & Schumacher (personal communication 2020) and Krug, Lohse & Stevens (Reference Krug, Lohse and Stevens2020). For DNS in cylinders of diameter-to-height ratio 
$\varGamma _c$, Iyer et al. (Reference Iyer, Scheel, Schumacher and Sreenivasan2020) used 
$(\varGamma _c,Pr)=(0.1,1)$ and Scheel & Schumacher (Reference Scheel and Schumacher2017) used 
$(\varGamma _c,Pr)=(1,0.7)$. For laboratory experiments in cylinders, where the plotted data are truncated according to 
$Pr\in [0.7,1.3]$, the domains and estimated 
$Pr$ ranges are 
$\varGamma _c = 0.5$ and 
$Pr\in [0.7,1.3]$ for Chavanne et al. (Reference Chavanne, Chilla, Chabaud, Castaing and Hebral2001), 
$\varGamma _c = 4$ and 
$Pr\in [0.7,1.27]$ for Niemela & Sreenivasan (Reference Niemela and Sreenivasan2006), 
$\varGamma _c = 0.5$ and 
$Pr\in [0.79,0.86]$ for He et al. (Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012), and 
$\varGamma _c = 1$ and 
$Pr\in [0.95,1.17]$ for Urban et al. (Reference Urban, Hanzelka, Musilová, Králík, La Mantia, Srnka and Skrbek2014). Experiments used working fluids of low-temperature helium gas (Chavanne et al. Reference Chavanne, Chilla, Chabaud, Castaing and Hebral2001; Niemela & Sreenivasan Reference Niemela and Sreenivasan2006; Urban et al. Reference Urban, Hanzelka, Musilová, Králík, La Mantia, Srnka and Skrbek2014) or sulphur hexafluoride (He et al. Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012).
The turbulent data shown in figure 5, as detailed in the figure caption, include DNS in horizontally periodic 2-D and 3-D domains, wherein 2-D steady rolls solve the equations of motion, as well as 3-D DNS and laboratory experiments in cylinders that do not admit 2-D rolls. Values of 
$Nu$ for steady rolls with 
$\varGamma =2$ fixed are omitted from figure 5 for clarity, but they lie below all turbulent values once 
$Ra$ approximately exceeds 
$2\times 10^{9}$ (cf. figure 2), and this gap would only widen at larger 
$Ra$ if their 
$Nu\sim Ra^{1/4}$ scaling persists. The laboratory data sets in figure 5 have unavoidably varying 
$Pr$ values that can be hard to estimate, as well as non-Oberbeck–Boussinesq effects (Urban, Musilová & Skrbek Reference Urban, Musilová and Skrbek2011; Urban et al. Reference Urban, Hanzelka, Kralik, Musilova, Srnka and Skrbek2012, Reference Urban, Hanzelka, Musilová, Králík, La Mantia, Srnka and Skrbek2014). The figure includes only a narrow range of estimated 
$Pr$ values in order to avoid significant non-Oberbeck–Boussinesq effects. When data over a wider range of estimated 
$Pr$ is included, a few data points from the experiments of Chavanne et al. (Reference Chavanne, Chilla, Chabaud, Castaing and Hebral2001) lie above the 
$Nu(\varGamma ^*)$ values of steady rolls, as shown in the supplementary material.
Our finding that steady rolls of 
$Nu$-maximizing aspect ratios apparently display classical 
$Nu\sim Ra^{1/3}$ asymptotic scaling does not ineluctably imply anything about turbulent convection. Taking a dynamical systems point of view, however, steady solutions admitted by the domain are fixed points of (2.1), so they and their unstable manifolds are part of the global attractor. Turbulent trajectories may linger near these fixed points and so inherit some quantitative features (Kooloth et al. Reference Kooloth, Sondak and Smith2021), as has been found for unstable coherent states in shear flows (Nagata Reference Nagata1990; Waleffe Reference Waleffe1998; Wedin & Kerswell Reference Wedin and Kerswell2004; Gibson, Halcrow & Cvitanović Reference Gibson, Halcrow and Cvitanović2008; Suri et al. Reference Suri, Kageorge, Grigoriev and Schatz2020; Graham & Floryan Reference Graham and Floryan2021). Indeed, figure 5 shows scaling similarities between steady and turbulent convection. Further exploration of the global attractor calls for study of 3-D steady flows. Recently computed ‘multi-scale’ 3-D steady states (Motoki et al. Reference Motoki, Kawahara and Shimizu2021) give larger 
$Nu$ values than all 2-D rolls at moderate 
$Ra$, but their scaling at large 
$Ra$ is unknown. Simpler 3-D steady convection patterns remain to be computed as well. Analytically, it is an open challenge to construct approximations of 2-D or 3-D steady flows that are asymptotically accurate as 
$Ra\to \infty$, as has been done for 2-D rolls between stress-free boundaries (Chini & Cox Reference Chini and Cox2009; Wen et al. Reference Wen, Goluskin, LeDuc, Chini and Doering2020). Such constructions could be used to verify that 
$Nu\sim Ra^{1/3}$ is indeed the exact asymptotic scaling for the 
$Nu$-maximizing rolls we have computed, as well as to determine the precise 
$Re$–
$Ra$ scaling relations for rolls of both 
$Nu$- and 
$Re$-maximizing aspect ratios.
More generally, figure 5 highlights the absence of reproducible evidence for ultimate 
$Nu \sim Ra^{1/2}$ scaling, and it raises the intriguing possibility that steady rolls with 
$Nu\sim Ra^{1/3}$ might transport more heat than turbulent convection as 
$Ra\to \infty$. We know of no counterexamples to this hypothesis, including in the case of stress-free boundaries (Wen et al. Reference Wen, Goluskin, LeDuc, Chini and Doering2020). Heat transport by solutions of (2.1) with no-slip isothermal boundaries has been mathematically proved to be limited by 
$Nu\le {O}(Ra^{1/2})$ (Howard Reference Howard1963; Doering & Constantin Reference Doering and Constantin1996), but it remains unknown whether any solutions attain the ultimate scaling of this upper bound. One avenue for pursuing a stronger mathematical statement is to study two conjectures suggested by our computations: that steady convection maximizes 
$Nu$ among all solutions of (2.1) regardless of their stability or time-dependence, and that steady solutions of (2.1) are subject to an upper bound of the form 
$Nu\le {O}(Ra^{1/3})$. Therefore, although numerically computed flows can never determine 
$Ra\to \infty$ scaling definitively, our results suggest a new mathematical approach that may be able to finally resolve the question of asymptotic 
$Nu$ scaling in turbulent convection.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2021.1042.
Acknowledgements
After this manuscript was written our senior author, C.R. Doering, passed away too soon. Beyond his many contributions to the present study, we are forever indebted to him for his mentorship, to say nothing of his many lasting contributions to the field of fluid dynamics. He will be deeply missed by us and many others. We also want to acknowledge helpful discussions about the present work with L.M. Smith, D. Sondak and F. Waleffe.
Funding
This work was supported by US National Science Foundation awards (DMS-1515161, DMS-1813003), Canadian NSERC Discovery Grants Program awards (RGPIN-2018-04263, RGPAS-2018- 522657, DGECR-2018-00371) and computational resources provided by Advanced Research Computing at the University of Michigan.
Declaration of interests
The authors report no conflict of interest.
Open access
