Steady Rayleigh--B\'enard convection between no-slip boundaries

The central open question about Rayleigh--B\'enard convection -- buoyancy-driven flow in a fluid layer heated from below and cooled from above -- is how vertical heat flux depends on the imposed temperature gradient in the strongly nonlinear regime where the flows are typically turbulent. The quantitative challenge is to determine how the Nusselt number $Nu$ depends on the Rayleigh number $Ra$ in the $Ra\to\infty$ limit for fluids of fixed finite Prandtl number $Pr$ in fixed spatial domains. Laboratory experiments, numerical simulations, and analysis of Rayleigh's mathematical model have yet to rule out either of the proposed `classical' $Nu \sim Ra^{1/3}$ or `ultimate' $Nu \sim Ra^{1/2}$ asymptotic scaling theories. Among the many solutions of the equations of motion at high $Ra$ are steady convection rolls that are dynamically unstable but share features of the turbulent attractor. We have computed these steady solutions for $Ra$ up to $10^{14}$ with $Pr=1$ and various horizontal periods. By choosing the horizontal period of these rolls at each $Ra$ to maximize $Nu$, we find that steady convection rolls achieve classical asymptotic scaling. Moreover, they transport more heat than turbulent convection in experiments or simulations at comparable parameters. If heat transport in turbulent convection continues to be dominated by heat transport in steady rolls as $Ra\to\infty$, it cannot achieve the ultimate scaling.


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 N u. In Rayleigh's mathematical model (Rayleigh 1916) N u 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 P r, which is the ratio of kinematic viscosity to thermal diffusivity, and (iii) details of the spatial domain, often captured by an aspect ratio Γ 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 N u with respect to both P r and Ra will emerge in the Ra → ∞ limit (Kadanoff 2001).
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 2012;Doering 2020). The two conjectures follow from heuristic physical arguments that both seem plausible but give incompatible predictions: the 'classical' scaling N u ∼ P r 0 Ra 1/3 and the 'ultimate' scaling N u ∼ P r 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 N u O(Ra 1/2 ) uniformly in P r and Γ (Howard 1963;Doering & Constantin 1996), 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 et al. 2006;Otto & Seis 2011;Whitehead & Doering 2012) and for two-dimensional convection between stress-free boundaries (Whitehead & Doering 2011). 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-N u as a function of Ra, P r, and Γ . Toward 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-independentsolutions 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 (1916). A roll pair of any width-to-height aspect ratio Γ admitted by the domain exists for sufficiently large Ra.
For steady rolls, the dependence of N u on the parameters (Γ, P r, 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 2021), where particular steady states, traveling 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  Figure 1: Steady convection rolls at Ra = 10 9 and P r = 1 with (a) Γ = 2 and (b) the value Γ * ≈ 0.235 that maximizes N u at these values of Ra and P r. Color indicates temperature, and streamlines are shown for counterclockwise (solid) and clockwise (dashdotted) motions. (c) Dependence of N u on the horizontal wavenumber k = 2π/Γ found by computing steady rolls of various aspect ratios ( ). Highlighted points are the two Γ shown in panels (a) and (b) along with the value Γ * loc ≈ 0.0614 that locally maximizes N u(Γ ). Cubic spline interpolation ( ) is used to find Γ * and Γ * loc precisely.
certain steady states capture qualitative aspects of turbulent convection Sondak et al. 2015;Kooloth et al. 2021;Motoki et al. 2021). Our findings add to this evidence. The desire to understand and perhaps strengthen the mathematical bound N u 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 P r = 1 fluid contained between no-slip isothermal top and bottom boundaries. We reach sufficiently large Ra values to convincingly reveal several asymptotic scalings of N u, 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 N u ∼ Ra 1/3 scaling. We further observe that N u for these steady rolls is larger than turbulent N u from all laboratory experiments and two-or threedimensional (2D or 3D) simulations at comparable parameters. This observation supports the conjecture that steady states maximize N u among all stable or unstable flows, as was recently verified for a truncated model of RBC (Olson et al. 2021) using methods that are not yet applicable to the full governing equations. If steady-roll transport continues to dominate turbulent transport as Ra → ∞, 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 N u ∼ Ra 1/3 as Ra → ∞ at fixed P r and Γ , and the aspect ratio of the roll pair maximizing N u at each Ra and P r approaches Γ ≈ 1.9 (Chini & Cox 2009;Wen et al. 2020). Recent computations of steady rolls in the no-slip case for pre-asymptotic Ra values up to 10 9 revealed significant differences from the stress-free problem Sondak et al. 2015). The dependence N u(Γ ) for no-slip rolls at fixed Ra and P r can have multiple local maxima, as shown in figure 1(c), and the aspect ratio Γ * that globally maximizes N u(Γ ) approaches zero rather than a constant as Ra → ∞. Steady rolls of N u-maximizing aspect ratios Γ * were reported in Sondak et al. (2015) for Ra ∈ [5 × 10 6 , 3 × 10 8 ] at P r = 1, yielding fits of Γ * ∼ Ra −0.217 and N u(Γ * ) ∼ Ra 0.31 . This heat transport scaling is faster than with Γ fixed: computations in Waleffe et al. (2015) for Ra ∈ [5 × 10 5 , 5 × 10 6 ] at P r = 7 with Γ = 2 fixed yield the fit N u ∼ 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 Γ = 2 fixed for Ra 2 × 10 10 and with the parameter-dependent aspect ratios Γ * and Γ * loc (cf. figure 1) that globally and locally maximize N u(Γ ), respectively, for Ra 10 14 . These Ra values are evidently large enough to reach asymptotia: the results reported below strongly suggest that fixed-Γ rolls asymptotically transport heat like N u ∼ Ra 1/4 while the ever-narrowing rolls of aspect ratio Γ * achieve the classical N u ∼ Ra 1/3 scaling.

Computation of steady-convection-roll solutions
Following Rayleigh (1916), we model RBC using the Boussinesq approximation to the Navier-Stokes equations with constant kinematic viscosity ν, thermal diffusivity κ, and coefficient of thermal expansion α. We nondimensionalize lengths by the layer height h, temperatures by the fixed difference ∆ between the boundaries, velocities by the freefall scale U f = √ gαh∆, 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 −ẑ direction. The evolution equations governing the dimensionless velocity vector u = (u, w), temperature T , and pressure p are then The dimensionless spatial domain is (x, z) ∈ [0, Γ ] × [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 u to vanish on both boundaries. The conductive state (u, T ) = (0, 1 − z) becomes unstable when Ra increases past the critical value Ra c ≈ 1708 (Jeffreys 1928), at which a roll pair with horizontal period Γ ≈ 2.016 bifurcates supercritically. As Ra → ∞ 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 where · 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. (2020) and Wen & Chini (2018) to the case of no-slip boundary conditions. In these numerics the temperature is represented using the deviation θ from the conductive profile, meaning T = 1 − z + θ, and the velocity is represented using a stream function ψ, where In terms of these variables, steady (∂ t = 0) solutions of (2.1) satisfy with fixed-temperature and no-slip boundary conditions, To compute solutions of the time-independent equations (2.4) and (2.5) by an iterative method, we do not need to impose all boundary conditions precisely on each iterationthe conditions need to hold only for the converged solution. Thus we do not impose (2.5c) exactly, instead using approximate boundary conditions on ω for equation (2.4a). These are derived by Taylor expanding ψ about the top and bottom boundaries to find where δ > 0 is small. Combining equations (2.4b) with (2.5b,c) and neglecting O(δ 4 ) terms in (2.6) give the approximate boundary conditions, In computations we set δ 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 2000). 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 Γ values were used as the initial iterate. Every 2 to 4 Newton iterations, we change the boundary values of the iterate to match the ∂ z ψ = 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 2norm of the residual of the governing steady equations had a relative magnitude less than 10 −10 . To accurately locate Γ * and Γ * loc , rolls were computed at several nearby Γ , and then N u(Γ ) was interpolated with cubic splines like those in figure 1(c). Details of computational results, including resolutions used, are included in the supplementary material.

Results
We computed steady P r = 1 rolls for aspect ratios Γ encompassing the three distinguished values indicated by figure 1(c): the fixed value Γ = 2 and the Ra-dependent values Γ * and Γ * loc that globally and locally maximize N u over Γ . As previously observed by Sondak et al. (2015), the N u(Γ ) curve has a single maximum when Ra is small and Ra for steady rolls with P r = 1 and aspect ratios of Γ = 2, Γ * , and Γ * loc , where the Ra-dependent values Γ * and Γ * loc are where N u(Γ ) has global and local maxima, respectively (cf. figure 1). Values of N u at Γ * from Sondak et al. (2015) and Waleffe (2020) are also shown ( ). Scaling fits ( ) over the last decade of each data set yield exponents of 0.33, 0.29, and 0.25. Bottom: Finite difference approximations of the local scaling exponent β n = d(log N u) d(log Ra) . Exponents of 1/3, 2/7, and 1/4 are shown to guide the eye ( ).
develops a second local maximum at smaller Γ when Ra increases past roughly 2 × 10 5 . The value of N u at this second local maximum remains less than the value at the first, so the picture remains as in figure 1(c) with Γ * on the left and Γ * loc on the right, in contrast to the P r = 10 and 100 cases (Sondak et al. 2015). For most Ra values we did not compute rolls over a full sweep through Γ as in figure 1(c), instead searching over Γ only as needed to locate Γ * and Γ * 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 N u on Ra for steady rolls with aspect ratios Γ = 2, Γ * , and Γ * loc . In the top panel N u − 1 is compensated by Ra 1/3 , so the horizontal line approached by rolls of the N u-maximizing aspect ratios Γ * corresponds to classical 1/3 scaling. The downward slopes of the data for aspect ratios 2 and Γ * loc correspond to scaling exponents smaller than 1/3. Values of N u at Γ * computed previously for Ra 10 9 (Sondak et al. 2015;Waleffe 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 N u to within 0.0008%.
The bottom panel of figure 2 shows the Ra-dependent local scaling exponent β n = d(log N u) d(log Ra) of the N u-Ra relation for Γ = 2, Γ * , and Γ * loc . This quantity educes small variations not visible in the top panel. In particular, for rolls of aspect ratios Γ * , the exponent β 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 → ∞ asymptotic behavior. 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 Γ = 2 fixed undergo a similarly rapid change around Ra ≈ 2 × 10 7 and then approach N u ∼ Ra 1/4 scaling that appears to be asymptotic. Rolls of aspect ratio Γ * loc show intermediate N u scaling whose best-fit exponent over the last decade of data is 0.29.
The top panel of figure 3 shows the Ra-dependence of the wavenumber k = 2π/Γ for Γ * and Γ * loc , compensated by Ra 1/5 . The compensated wavenumbers for Γ * approach a horizontal line, suggesting that the N u-maximizing rolls narrow according to the power law Γ * ∼ Ra −1/5 . This narrowing of Γ * is slow relative to the case of RBC in a porous medium, where Γ * ∼ Ra −1/2 (Wen et al. 2015).
The bottom panel of figure 3 shows the Ra-dependence of the local scaling exponent β k = d(log k) d(log Ra) . For k = 2π/Γ * the local scaling exponent remains close to 1/5 after the transition around Ra = 10 9 . For k = 2π/Γ * loc the exponent seems to approach 1/4, suggesting that Γ * loc has the same Ra −1/4 scaling as the narrowest marginally stable mode. Variations in β k beyond Ra = 10 12 for Γ * are evident, but these might be due to numerical imprecision: N u depends very weakly on Γ around the maximum of N u(Γ ), as seen in figure 1(c), so the value of Γ * cannot be determined nearly as precisely as the value of N u(Γ * ).

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 Ra for steady rolls with P r = 1 and aspect ratios Γ = 2, Γ * , and Γ * loc . Scaling fits yield Re ∼ Ra 0.50 for Γ = 2 and Re ∼ Ra 0.40 for Γ * loc over the last decade of each data set, and Re ∼ Ra 0.47 for Γ * over Ra ∈ [10 10 , 10 14 ] ( ). Bottom: Finite difference approximations to the local exponent β r = d(log Re) d(log Ra) . Exponents of 1/2 and 2/5 are shown to guide the eye ( ). Figure 4 depicts the dependence of Re on Ra for the steady rolls of aspect ratios Γ = 2, Γ * , and Γ * loc . The top panel shows Re compensated by Ra 1/2 while the bottom panel shows the local scaling exponent β r = d(log Re) d(log Ra) . Rolls with the fixed aspect ratio Γ = 2 approach the asymptotic scaling Re ∼ Ra 1/2 that corresponds to the root-mean-squared velocity being proportional to the free-fall velocity U f . For rolls with N u-maximizing aspect ratios Γ * , the scaling fit over Ra ∈ [10 10 , 10 14 ] is Re ∼ Ra 0.47 , which is quite close to the Re ∼ Ra 0.46 scaling observed in recent 3D direct numerical simulations up to Ra = 10 15 at P r = 1 in a slender cylinder with a height 10 times its diameter (Iyer et al. 2020). For the Γ * 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 ∼ Ra 2/3 for any fixed P r and Γ (Wen et al. 2020).

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 P r = 1 or 0.7 and laboratory experiments where the estimated P r is between 0.7 and 1.3. Figure 5 shows these N u values compensated by Ra 1/3 , along with N u values of steady convection rolls at the N u-maximizing aspect ratios Γ * . Strikingly, heat transport by the N u-maximizing 2D steady rolls is larger than transport by turbulent convection in all cases.
The turbulent data shown in figure 5, as detailed in the figure caption, include DNS in horizontally periodic 2D and 3D domains, wherein 2D steady rolls solve the equations of motion, as well as 3D DNS and laboratory experiments in cylinders that do not admit 2D rolls. Values of N u for steady rolls with Γ = 2 fixed are omitted from figure 5 for clarity, but they lie below all turbulent values once Ra approximately exceeds 2 × 10 9 (cf. figure 2), and this gap would only widen at larger Ra if their N u ∼ Ra 1/4 scaling persists. The laboratory data sets in figure 5 have unavoidably varying P r values that can be hard to estimate, as well as non-Oberbeck-Boussinesq effects (Urban et al. 2011(Urban et al. , 2012(Urban et al. , 2014. The figure includes only a narrow range of estimated P r values in order to avoid significant non-Oberbeck-Boussinesq effects. When data over a wider range of estimated P r is included, a few data points from the experiments of Chavanne et al. (2001) lie above the N u(Γ * ) values of steady rolls, as shown in the supplementary material. Our finding that steady rolls of N u-maximizing aspect ratios apparently display classical N u ∼ 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. 2021), as has been found for unstable coherent states in shear flows (Nagata 1990;Waleffe 1998;Wedin & Kerswell 2004;Gibson et al. 2008;Suri et al. 2020;Graham & Floryan 2021). Indeed, figure 5 shows scaling similarities between steady and turbulent convection. Further exploration of the global attractor calls for study of 3D steady flows. Recently computed 'multi-scale' 3D steady states (Motoki et al. 2021) give larger N u values than all 2D rolls at moderate Ra, but their scaling at large Ra is unknown. Simpler 3D steady convection patterns remain to be computed as well. Analytically, it is an open challenge to construct approximations of 2D or 3D steady flows that are asymptotically accurate as Ra → ∞, as has been done for 2D rolls between stress-free boundaries (Chini & Cox 2009;Wen et al. 2020). Such constructions could be used to verify that N u ∼ Ra 1/3 is indeed the exact asymptotic scaling for the N u-maximizing rolls we have computed, as well as to determine the precise Re-Ra scaling relations for rolls of both N u-and Re-maximizing aspect ratios.
More generally, figure 5 highlights the absence of reproducible evidence for ultimate N u ∼ Ra 1/2 scaling, and it raises the intriguing possibility that steady rolls with N u ∼ Ra 1/3 might transport more heat than turbulent convection as Ra → ∞. We know of no counterexamples to this hypothesis, including in the case of stress-free boundaries (Wen et al. 2020). Heat transport by solutions of (2.1) with no-slip isothermal boundaries has been mathematically proved to be limited by N u O(Ra 1/2 ) (Howard 1963;Doering & Constantin 1996), 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 N u 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 N u O(Ra 1/3 ). Therefore, although numerically computed flows can never determine Ra → ∞ scaling definitively, our results suggest a new mathematical approach that may be able to finally resolve the question of asymptotic N u scaling in turbulent convection.

Numerical solutions
Tables 1S to 3S give the N u, k, and Re values for numerical solutions with P r = 1 and Γ = 2, Γ * , and Γ * loc , respectively.  Comparison with turbulent convection Figure 1S is nearly identical to figure 5 in the main text, comparing heat transport by N u-maximizing steady rolls with transport by turbulent convection, except that more experimental data with Prandtl numbers further from 1 are included. In figure 1S the criterion for inclusion is an estimated Prandtl number of P r ∈ [0.5, 2] rather than the range Pr ∈ [0.7, 1.3] in figure 5 of the main text. (In fact all of the estimated Pr are at least 0.6, so Pr ∈ [0.6, 2] in figure 1S.) The working fluids in the experiments-gaseous helium or sulfur hexaflouride-are used near their critical points, leading to coupling and sensitive variation of material parameters that can be difficult to estimate. Faster variation of P r with Ra is associated with increasing non-Oberbeck-Boussinesq effects as well; see Urban et al. (2011Urban et al. ( , 2012Urban et al. ( , 2014 for a discussion of experimental challenges. Data in figure 5 is truncated using the narrower range P r ∈ [0.7, 1.3] mainly to reduce non-Oberbeck-Boussinesq effects-we expect Pr alone to have a more modest effect, even over the wider range [0.5, 2].