1. Introduction
The fundamental challenge of understanding thermally driven flow in the strongly nonlinear regime has puzzled the fluid dynamics community for more than a century (Doering Reference Doering2019). The typical set-up that is studied theoretically consists of a fluid confined between two horizontal plates that is heated from below and cooled at the top. This is the so-called Rayleigh–Bénard convection (RBC) problem after Rayleigh's proposed model (Rayleigh Reference Rayleigh1916) for Bénard's experiment on buoyancy-driven thermal convection (Bénard Reference Bénard1900, Reference Bénard1927). The classical problem depends on two dimensionless parameters: the Rayleigh number 
$\mathit {Ra}$, which measures the driving effect of buoyancy relative to the stabilising effects of viscosity and thermal diffusivity; and the Prandtl number 
$\mathit {Pr}$, which is the ratio of kinematic viscosity to thermal diffusivity. The global flow properties may be characterised by the Nusselt number 
$\mathit {Nu}$, a further dimensionless parameter that measures the total heat flux relative to the purely conductive heat flux. A way to progress our understanding on RBC is to understand the dependence of the dynamics on 
$\mathit {Ra}$, 
$\mathit {Pr}$ and the boundary conditions.
Bolgiano (Reference Bolgiano1959) and Obukhov (Reference Obukhov1959) proposed the so-called Bolgiano–Obukhov (BO) phenomenology for stably stratified turbulence, in which buoyancy balances inertia in the momentum equation and the potential energy flux is approximately constant for length scales larger than the Bolgiano scale. These assumptions lead to kinetic and potential energy spectra of the form 
$\hat {E}_u(k) \propto k^{-11/5}$ and 
$\hat {E}_\theta (k) \propto k^{-7/5}$, respectively, where 
$k$ is the wavenumber. This scaling law is considered to be universal under certain conditions, irrespective of the specific details of the flow geometry, boundary conditions or forcing mechanisms. However, deviations from universality can occur. Despite originally being proposed for homogeneous stably stratified flows, BO scaling has been partly identified in experimental and numerical studies of inhomogeneous (i.e. with no-slip or free-slip boundary conditions on the top and bottom plates) three-dimensional (3-D) RBC in different geometries (Wu et al. Reference Wu, Kadanoff, Libchaber and Sano1990; Calzavarini, Toschi & Tripiccione Reference Calzavarini, Toschi and Tripiccione2002; Mishra & Verma Reference Mishra and Verma2010; Kaczorowski & Xia Reference Kaczorowski and Xia2013). However, the existence of the BO scaling is still debatable (Lohse & Xia Reference Lohse and Xia2010; Kunnen & Clercx Reference Kunnen and Clercx2014), with some studies of inhomogeneous and homogeneous (i.e. with periodic boundary conditions) 3-D RBC reporting (Kolmogorov Reference Kolmogorov1941) power-law scaling 
$k^{-5/3}$ for the kinetic and potential energy spectra (Borue & Orszag Reference Borue and Orszag1997; Verma, Kumar & Pandey Reference Verma, Kumar and Pandey2017; Verma Reference Verma2018). Similar controversy remains in simulations of two-dimensional (2-D) RBC, with some studies partly identifying BO scaling (Toh & Suzuki Reference Toh and Suzuki1994; Mazzino Reference Mazzino2017; Xie & Huang Reference Xie and Huang2022; Samuel & Verma Reference Samuel and Verma2024), while some argue against the validity of the BO scaling (Biskamp & Schwarz Reference Biskamp and Schwarz1997; Biskamp, Hallatschek & Schwarz Reference Biskamp, Hallatschek and Schwarz2001; Celani, Mazzino & Vergassola Reference Celani, Mazzino and Vergassola2001; Celani et al. Reference Celani, Matsumoto, Mazzino and Vergassola2002).
For the strongly nonlinear regime of thermal convection, which is of paramount importance for geophysical and astrophysical applications, there are two competing theories for the behaviour of 
$\mathit {Nu}$ as 
$\mathit {Ra}$ tends to infinity for arbitrary 
$\mathit {Pr}$. These two proposed asymptotic scaling laws are the ‘classical’ theory 
$\mathit {Nu} \propto \mathit {Pr}^0\mathit {Ra}^{1/3}$ by Malkus (Reference Malkus1954) and the ‘ultimate’ theory 
$\mathit {Nu} \propto \mathit {Pr}^{1/2}\mathit {Ra}^{1/2}$, at least for 
$\mathit {Pr} \leqslant 1$, by Kraichnan (Reference Kraichnan1962) and Spiegel (Reference Spiegel1971). The Rayleigh number at which the transition to the ultimate scaling is presumed to occur is not known, and laboratory experiments and numerical simulations (Niemela et al. Reference Niemela, Skrbek, Sreenivasan and Donnelly2000; Chavanne et al. Reference Chavanne, Chillà, Chabaud, Castaing and Hébral2001; Niemela & Sreenivasan Reference Niemela and Sreenivasan2003; Urban, Musilová & Skrbek Reference Urban, Musilová and Skrbek2011; He et al. Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012; Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018; Doering, Toppaladoddi & Wettlaufer Reference Doering, Toppaladoddi and Wettlaufer2019) have reported different power laws for a wide range of Rayleigh numbers.
In the analogy of RBC with geophysical phenomena, the top and bottom boundaries are often absent, particularly when the focus is on understanding the dynamics of the bulk flow. In this study, we choose the most obvious theoretical approach to bypass boundary layer effects by considering a fully periodic domain (Toh & Suzuki Reference Toh and Suzuki1994; Borue & Orszag Reference Borue and Orszag1997; Biskamp et al. Reference Biskamp, Hallatschek and Schwarz2001; Celani et al. Reference Celani, Matsumoto, Mazzino and Vergassola2002; Ng et al. Reference Ng, Ooi, Lohse and Chung2018) for the Rayleigh–Bénard problem, with an imposed constant vertical temperature gradient. For this homogeneous RBC set-up it has been claimed that 
$Nu \propto Ra^{1/2}$ (Lohse & Toschi Reference Lohse and Toschi2003; Calzavarini et al. Reference Calzavarini, Lohse, Toschi and Tripiccione2005), in line with the ultimate theory, and this scaling is also suggested from simulations of axially periodic RBC in a vertical cylinder (Schmidt et al. Reference Schmidt, Calzavarini, Lohse, Toschi and Verzicco2012).
Homogeneous RBC, however, exhibits exponentially growing solutions in the form of axially uniform vertical jets, called ‘elevator modes’ (Baines & Gill Reference Baines and Gill1969; Batchelor & Nitsche Reference Batchelor and Nitsche1991; Calzavarini et al. Reference Calzavarini, Doering, Gibbon, Lohse, Tanabe and Toschi2006; Liu & Zikanov Reference Liu and Zikanov2015). As soon as these modes grow to a significant amplitude, they experience secondary instabilities, ultimately leading to statistically stationary solutions (Calzavarini et al. Reference Calzavarini, Doering, Gibbon, Lohse, Tanabe and Toschi2006; Schmidt et al. Reference Schmidt, Calzavarini, Lohse, Toschi and Verzicco2012). In recent numerical simulations, the elevator modes were suppressed by introducing an artificial horizontal buoyancy field (Xie & Huang Reference Xie and Huang2022) or large-scale friction (Barral & Dubrulle Reference Barral and Dubrulle2023). The inverse cascade that is observed in 2-D homogeneous RBC (Toh & Suzuki Reference Toh and Suzuki1994; Verma Reference Verma2018; Xie & Huang Reference Xie and Huang2022) is another source of energy for the large-scale modes that grow to extreme values, forming a condensate whose amplitude saturates when the viscous dissipation at the largest scale balances the energy injection (Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007; Boffetta & Ecke Reference Boffetta and Ecke2012). So, in this study, to avoid the unbounded growth of energy we include large-scale dissipation to mimic the effect of friction when there are boundaries and to be able to reach statistically stationary solutions.
The aim of this paper is to bring more insight to the effect of the control parameters and the boundaries on the dynamics of 2-D RBC. Most of the attention on 2-D RBC focuses on the Rayleigh number dependence of the dynamics, with only a few studies considering the effects of the Prandtl number (Calzavarini et al. Reference Calzavarini, Lohse, Toschi and Tripiccione2005; Zhang, Zhou & Sun Reference Zhang, Zhou and Sun2017; He et al. Reference He, Fang, Gao, Huang and Bao2021). So, we extensively study the effects of the Prandtl number using numerical simulations of 2-D RBC driven by a constant temperature gradient in a periodic domain, to effectively remove the boundaries. We also consider high Rayleigh number simulations with normal viscosity and hyperviscosity to permit large scale separation, which is crucial for the analysis of the multi-scale dynamics, that will also shed light on the long standing debate of the power-law scalings of the spectra.
The paper is structured as follows. Section 2 contains the dynamical equations, numerical methods and the definitions of the global and spectral observables under study. The results of our simulations are presented in § 3, where we analyse how the global observables scale with the Prandtl and Rayleigh numbers and then investigate the effects of these dimensionless parameters on the spectral dynamics. In Appendix A we discuss the behaviour of the system at zero and infinite Prandtl number and in § 4 we summarise our conclusions.
2. Problem description
2.1. Governing equations
We consider 2-D RBC of a fluid heated from below in a periodic square cell 
$(x,y)\in [0,L]^2$. The temperature 
$T(x,y,t)$ is decomposed as 
$T=-\Delta T y/L+\theta$, where 
$\Delta T/L$ is the constant imposed temperature gradient and the temperature perturbation 
$\theta (x,y,t)$ satisfies periodic boundary conditions. As usual, for simplicity we employ the Oberbeck–Boussinesq approximation (Oberbeck Reference Oberbeck1879; Boussinesq Reference Boussinesq1903; Tritton Reference Tritton2012), in which the kinematic viscosity 
$\nu$ and the thermal diffusivity 
$\kappa$ are taken to be constant while the temperature dependence of the fluid density 
$\rho$ is neglected, except in the buoyancy term of the momentum equation.
The governing equations of the problem in two dimensions can be written in terms of 
$\theta (x,y,t)$ and the streamfunction 
$\psi (x, y, t)$ as follows:
$$\begin{gather} \partial_t \nabla^2 \psi + \{\psi, \nabla^2 \psi\} = \alpha g \partial_x \theta + \nu {(-\nabla^2)^{n+1}} \psi + \mu \psi, \end{gather}$$
$$\begin{gather}\partial_t \theta + \{\psi, \theta \} = \frac{\Delta T}{L} \partial_x \psi - \kappa {(-\nabla^2)^n} \theta, \end{gather}$$
where 
$\{A, B\} = \partial _x A \partial _y B - \partial _y A \partial _x B$ is the standard Poisson bracket, 
$\alpha$ is the thermal expansion coefficient and 
$g$ is the gravitational acceleration. To prevent the formation of a large-scale condensate in the presence of an inverse cascade, to avoid the elevator modes and to reach a turbulent stationary regime, we supplement our system with a large-scale dissipative term 
$\mu \psi$ that is responsible for saturating the inverse cascade. We consider both normal and hyperviscosity by raising the Laplacian to the powers of 
$n=1$ and 
$n=4$, respectively. The hyperviscous case, albeit not physically realisable, gives a wider inertial range, as diffusive and viscous terms kick in abruptly at much smaller scales compared with the normal viscosity case. In the ideal limit of 
$\nu \to 0$, 
$\kappa \to 0$ and 
$\mu \to 0$, the quantity that is conserved is 
$E_u(t) - ({agL}/{\Delta T})E_\theta (t)$, where the instantaneous kinetic energy 
$E_u(t)$ and the potential energy 
$E_\theta (t)$ are defined by
\begin{equation} E_u(t) =\tfrac{1}{2}\langle{|\nabla\psi|^2}\rangle_\textbf{x}, \quad E_\theta(t) = \tfrac{1}{2}\langle{\theta^2}\rangle_\textbf{x}, \end{equation}
with the angle brackets 
$\langle {\cdot }\rangle _\textbf{x}$ here denoting average in space. Spatio-temporal averages will be denoted as 
$\langle {\cdot }\rangle$.
Equations (2.1) depend on three dimensionless parameters, namely
\begin{equation} \mathit{Pr} = \frac{\nu}{\kappa}, \quad \mathit{Ra} = \frac{\alpha g\Delta T L^{4n-1}}{\nu \kappa}, \quad \mathit{Rh} = \mu \left(\frac{L^5}{\alpha g \Delta T}\right)^{1/2}, \end{equation}which are the Prandtl number, Rayleigh number and friction Reynolds number, respectively.
We perform direct numerical simulations of (2.1) using the pseudospectral method (Orszag & Gottlieb Reference Orszag and Gottlieb1977). We decompose the streamfunction into basis functions with Fourier modes in both the 
$x$ and 
$y$ directions, viz.
\begin{equation} \psi({\boldsymbol{x}},t) = \sum^{N/2}_{{\boldsymbol{k}} ={-}N/2} \hat{\psi}_{\boldsymbol{k}}(t) {\rm e}^{{\rm i} {\boldsymbol{k}} \boldsymbol{\cdot}\boldsymbol{x}}, \end{equation}
where 
$\hat { \psi }_{\boldsymbol {k}}$ is the amplitude of the 
${\boldsymbol {k}} = (k_x, k_y)$ mode of 
$\psi$, and 
$N$ denotes the number of aliased modes in the 
$x$- and 
$y$-directions. We decompose 
$\theta$ in the same way. A third-order Runge–Kutta scheme is used for time advancement and the aliasing errors are removed with the two-thirds dealiasing rule (Gómez, Mininni & Dmitruk Reference Gómez, Mininni and Dmitruk2005). In both the normal and hyperviscous simulations, we find that 
$\mathit {Rh}=(2{\rm \pi} )^{5/2}$ yields a saturated turbulent state that dissipates enough kinetic energy at large scales such that the kinetic energy spectrum peaks at 
$k = 2$ without overdamping the system. So, we fix 
$\mathit {Rh} = (2{\rm \pi} )^{5/2} \simeq 100$ while varying the Rayleigh and Prandtl numbers in the ranges 
$6.2 \times 10^7 \leqslant \mathit {Ra} \leqslant 6.2 \times 10^{11}$ and 
$10^{-3} \leqslant \mathit {Pr} \leqslant 10^2$. To model the large Rayleigh number dynamics, we set 
$\mathit {Ra} = 9.4 \times 10^{49}$ in our hyperviscous simulations. Figure 1 shows the parameter values simulated in the 
$(\mathit {Ra},\mathit {Pr})$-plane as well as the resolution, 
$N$, used in each case. The criterion we used to ensure that our runs are well resolved was to verify that the peaks of the potential and kinetic energy dissipation spectra are at least one order of magnitude above the tail of the spectra. Time-averaged quantities are computed after the system has reached a statistically stationary regime. We have accumulated statistics over long time integrations to ensure that our statistics are well converged.
Figure 1. Parameter values simulated in the 
$(\mathit {Ra},\mathit {Pr})$-plane with the resolution, 
$N$, used at each instant colour coded in the legend. The black dashed line separates runs with normal viscosity (below line) from those with hyperviscosity (above line).
2.2. Global and spectral observables
Next, we briefly outline the global and spectral observables that will be explored in our numerical simulations below. The energy spectra of the velocity field 
$\hat {E}_u(k,t)$ and the temperature field 
$\hat {E}_{\theta }(k,t)$, referred to as the kinetic energy and potential energy spectra, are defined as
$$\begin{gather} \hat{E}_u(k,t) = \frac{1}{2}\sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} k^2\left| \hat{\psi}_{\boldsymbol{k}}(t) \right|^2, \end{gather}$$
$$\begin{gather}\hat{E}_{\theta}(k,t) = \frac{1}{2}\sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} \left| \hat{\theta}_{\boldsymbol{k}}(t) \right|^2, \end{gather}$$
where the sum is performed over the Fourier modes with wavenumber amplitude 
$k = |{\boldsymbol {k}}| = \sqrt {k_x^2 + k_y^2}$ in a shell of width 
$\Delta k = 2{\rm \pi} /L$. Using the Fourier transform, one can derive the evolution equations of kinetic and potential energy spectra from (2.1), namely
$$\begin{gather} \partial_t \hat{E}_u(k,t) ={-} \partial_k \varPi_u(k,t) - D_\nu(k,t) - D_\mu(k,t) + \alpha g F_B(k,t), \end{gather}$$
$$\begin{gather}\partial_t \hat{E}_{\theta}(k,t) ={-}\partial_k \varPi_\theta(k,t) - D_\kappa(k,t) + \frac{\Delta T}{L}F_B(k,t). \end{gather}$$ The energy flux 
$\varPi$ is a measure of the nonlinear cascades in turbulence (Alexakis & Biferale Reference Alexakis and Biferale2018). The energy flux for a circle of radius 
$k$ in the 2-D wavenumber space is the total energy transferred from the modes within the circle to the modes outside the circle. Consequently, we define the flux of kinetic energy 
$\varPi _u(k,t)$ and potential energy 
$\varPi _{\theta }(k,t)$ as
$$\begin{gather} \varPi_u(k,t) = \sum_{k' \leqslant k} T_u(k',t), \end{gather}$$
$$\begin{gather}\varPi_{\theta}(k,t) = \sum_{k' \leqslant k} T_\theta(k',t), \end{gather}$$
where 
$T_u(k,t)$ and 
$T_\theta (k,t)$ are the nonlinear kinetic and potential energy transfer across 
$k$
$$\begin{gather} T_u(k,t) ={-}\sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} \hat{\psi}^*_{\boldsymbol{k}}(t) \widehat{\{\psi,\nabla^2 \psi \}}_{\boldsymbol{k}}(t), \end{gather}$$
$$\begin{gather}T_\theta(k,t) = \sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} \hat{\theta}^*_{\boldsymbol{k}}(t) \widehat{\{\psi,\theta \}}_{\boldsymbol{k}}(t). \end{gather}$$
The notation 
$\widehat {\{ . \}}_{\boldsymbol {k}}$ represents the Fourier mode of the Poisson bracket expanded using (2.4), and the asterisk denotes complex conjugation.
The spectra of the small-scale viscous dissipation 
$D_\nu (k,t)$, the large-scale frictional dissipation 
$D_\mu (k,t)$ and the thermal dissipation 
$D_\kappa (k,t)$ are defined as
$$\begin{gather} D_\nu(k,t) = 2 \nu k^{2n}\hat{E}_u(k,t), \end{gather}$$
$$\begin{gather}D_\mu(k,t) = 2 \mu k^{{-}2} \hat{E}_u(k,t), \end{gather}$$
$$\begin{gather}D_\kappa(k,t) = 2 \kappa k^{2n} \hat{E}_\theta(k,t), \end{gather}$$and the spectrum of the buoyancy term is given by
\begin{equation} F_B(k,t) = \sum_{k \leqslant |{\boldsymbol{k}}| < k+\Delta k} {\rm i} k_x \widehat{\psi^*}_{\boldsymbol{k}}(t) \hat{\theta}_{\boldsymbol{k}}(t). \end{equation}The Nusselt number is a dimensionless measure of the averaged vertical heat flux relative to the conductive heat flux. The form of the conductive heat flux is determined by the thermal diffusive term. Hence, the Nusselt number is defined mathematically by
\begin{equation} \mathit{Nu} = 1 + \frac{\langle{\partial_x \psi \theta}\rangle}{\kappa \Delta T / L^{2n-1}}. \end{equation}Using the above definition, one can derive the following exact relations for the kinetic and potential energy balances in the statistically stationary regime, as in Shraiman & Siggia (Reference Shraiman and Siggia1990):
$$\begin{gather} \epsilon_u = \epsilon_\nu + \epsilon_\mu = \frac{\nu^3}{L^{6n-2}}(\mathit{Nu} - 1)\frac{\mathit{Ra}}{\mathit{Pr}^2}, \end{gather}$$
$$\begin{gather}\epsilon_\theta = \epsilon_\kappa = \frac{\kappa \Delta T^2}{L^{2n}} (\mathit{Nu} - 1), \end{gather}$$
where 
$\epsilon _u = \alpha g \langle {\psi \partial _x \theta }\rangle$ is the injection rate of kinetic energy due to buoyancy, 
$\epsilon _\nu = \nu {\langle \psi \nabla ^{2(n+1)} \psi \rangle }$ is the viscous dissipation rate, 
$\epsilon _\mu = \mu {\langle \psi ^2\rangle }$ is the large-scale dissipation rate, 
$\epsilon _\theta = ({\Delta T}/{L}) \langle {\partial _x \psi \,\theta }\rangle$ is the injection rate of potential energy due to buoyancy and 
$\epsilon _\kappa = \kappa {\langle \theta \nabla ^{2n} \theta \rangle }$ is the thermal dissipation rate.
2.3. Elevator modes
Upon linearising (2.1) about the conductive state (
$\psi =\theta =0$), we find that infinitesimal solutions with 
$\psi (\boldsymbol {x},t)=\exp ({{\rm i} \boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}+\sigma t})$ are possible, provided the normalised linear growth rate 
$\sigma$ satisfies the relation
\begin{equation} \left(\sigma+k^{2n}\right) \left(\sigma+\mathit{Pr} k^{2n}+\frac{\mathit{Rh}\sqrt{\mathit{Ra}\mathit{Pr}}}{k^2}\right) =\frac{\mathit{Ra}\mathit{Pr} (2{\rm \pi} k_x)^2}{k^2}, \end{equation}
which has two real roots for 
$\sigma$. From these two roots, one is positive if and only if
\begin{equation} \mathit{Rh}\sqrt{\frac{\mathit{Ra}}{\mathit{Pr}}}<\frac{\mathit{Ra}(2{\rm \pi} k_x)^2}{k^{2n}}-k^{2n+2}. \end{equation}
One can show that 
$\sigma$ is a monotonic decreasing function of 
$k_y$, so the most dangerous modes are independent of 
$y$, with
\begin{equation} \psi(\boldsymbol{x},t)=\exp({{\rm i}k_x x+\sigma t}),\quad k_x\in\mathbb{Z}_{{>}0}. \end{equation}Indeed, such a unidirectional mode satisfies the nonlinear governing equations (2.1) exactly, because the nonlinear Poisson bracket terms are identically zero.
Although the maximum growth rate does not necessarily occur at the minimum wavenumber 
$k_x=1$, it is straightforward to show that (2.14) is satisfied for some 
$(k_x,k_y)$ if and only if it is satisfied at 
$(k_x,k_y)=(1,0)$. In other words, exact solutions of the problem (2.1) that grow exponentially without bound exist whenever
\begin{equation} \mathit{Rh}\sqrt{\frac{\mathit{Ra}}{\mathit{Pr}}}< \frac{\mathit{Ra}}{(2{\rm \pi})^{2(n-1)}}-(2{\rm \pi})^{2n+2}. \end{equation}
In practice, unless the initial conditions are equal to an elevator mode, we find that elevator modes do not dominate the dynamics at finite Prandtl number if 
$\mathit {Rh}$ is large enough. So, we fix 
$\mathit {Rh}=(2{\rm \pi} )^{5/2}$, which is sufficient to prevent elevator modes from growing and to allow the flow to reach turbulent saturated states without overdamping the system. For zero and infinite Prandtl number convection, elevator modes are seemingly more attractive, even for randomly generated initial conditions (see Appendix A). Thus, our analysis focuses on finite values of Prandtl number.
3. Results
3.1. Global observables
In figures 2(a)–2(c) we show how the global observables in the kinetic and potential energy balances (2.12) vary with 
$\mathit {Pr}$ keeping 
$\mathit {Ra} = 6.2 \times 10^{11}$, while in (d–f) we show how the same quantities vary with 
$\mathit {Ra}$ while keeping 
$\mathit {Pr}=1$. In both cases, we use normal viscosity (
$n=1$) and keep 
$\mathit {Rh}=(2{\rm \pi} )^{5/2}$ constant. Firstly, figures 2(a) and 2(d) show that the kinetic and potential energies remain virtually constant while 
$\mathit {Pr}$ and 
$\mathit {Ra}$ vary by at least four orders of magnitude. Secondly, figures 2(b) and 2(e) show that 
$\epsilon _u \approx \epsilon _\mu$, which indicates that the majority of the kinetic energy, injected by buoyancy in the flow, is dissipated at large scales. This effect is caused by the inverse cascade of kinetic energy, which will be investigated in more detail in § 3.2. Note that 
$\epsilon _u$ and 
$\epsilon _\mu$ are almost independent of both 
$\mathit {Pr}$ and 
$\mathit {Ra}$, while the viscous dissipation scales like
\begin{equation} \epsilon_\nu \propto \mathit{Pr}^{1/2} \mathit{Ra}^{{-}1/4}. \end{equation}
Finally, in figures 2(c) and 2(f), we see that 
$\epsilon _\theta =\epsilon _\kappa$, as required by the potential energy balance (2.12b), and both observables, are also virtually independent of both 
$\mathit {Pr}$ and 
$\mathit {Ra}$.
Figure 2. Spatio-temporal averages of the kinetic and potential energies and the terms in the kinetic and potential energy balances (2.12) as functions of 
$\mathit {Pr}$ with 
$\mathit {Ra} = 6.2 \times 10^{11}$ fixed (a–c) and as functions of 
$\mathit {Ra}$ with 
$\mathit {Pr} = 1$ fixed (d–f). We use normal viscosity (
$n=1$) and keep 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$ fixed in all plots.
To explain the scaling of the viscous dissipation rate (3.1), we now look at how the enstrophy 
$\langle {\omega ^2}\rangle$ varies with the Prandtl and Rayleigh numbers, where 
$\omega = \nabla ^2 \psi$ is the vorticity of the flow. In figure 3(a) we plot enstrophy as a function of 
$\mathit {Pr}$ with 
$\mathit {Ra} = 6.2 \times 10^{11}$ fixed, and in figure 3(b) as a function of 
$\mathit {Ra}$ with 
$\mathit {Pr} = 1$ fixed. We keep 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$ fixed throughout. From figure 3(a) we observe that enstrophy can be considered approximately independent of 
$\mathit {Pr}$, because it varies by less than a factor of two over the five decades of Prandtl numbers considered, while figure 3(b) demonstrates that enstrophy scales like 
$\mathit {Ra}^{1/4}$, i.e.
\begin{equation} \langle{\omega^2}\rangle \propto \mathit{Pr}^0 \mathit{Ra}^{1/4}. \end{equation}
With normal viscosity, by definition we have 
$\epsilon _\nu = \nu \langle {\psi \nabla ^4 \psi }\rangle = \nu \langle {\omega ^2}\rangle$ and, writing 
$\nu \propto (\mathit {Pr}/\mathit {Ra})^{1/2}$, we thus obtain the scaling of (3.1) that is observed in figure 2.
Figure 3. Enstrophy 
$\langle {\omega ^2}\rangle$ as a function of 
$\mathit {Pr}$ with 
$\mathit {Ra} = 6.2 \times 10^{11}$ fixed (a), and as a function of 
$\mathit {Ra}$ with 
$\mathit {Pr} = 1$ fixed (b). We keep 
$n=1$ (normal viscosity) and 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$ fixed in all plots.
According to figure 2, the energy injection rates 
$\epsilon _u$ and 
$\epsilon _\theta$ are both approximately constant in the range of 
$\mathit {Pr}$ and 
$\mathit {Ra}$ we considered. With 
$\mathit {Ra}\mathit {Pr} \gg 1$ and 
$n = 1$, the net energy balances (2.12) thus produce the Nusselt number scaling
\begin{equation} \mathit{Nu} \propto \mathit{Pr}^{1/2}\mathit{Ra}^{1/2}. \end{equation}
This relation is in agreement with the ultimate scaling. In figure 4, we plot the Nusselt number compensated by the classical scaling, 
$\mathit {Nu} \propto \mathit {Pr}^{0}\mathit {Ra}^{1/3}$, and by the ultimate scaling, 
$\mathit {Nu} \propto \mathit {Pr}^{1/2}\mathit {Ra}^{1/2}$. The ultimate scaling provides a much more convincing collapse of the data, with the fit becoming progressively more accurate as 
$\mathit {Ra}$ increases for all values of 
$Pr$ considered. Indeed, the ultimate scaling in terms of Rayleigh number dependence was expected to be exhibited by our simulations, as we have effectively removed the boundary layers by applying periodic boundary conditions, similarly to previous studies (Lohse & Toschi Reference Lohse and Toschi2003; Calzavarini et al. Reference Calzavarini, Lohse, Toschi and Tripiccione2005). The ultimate scaling asserts (3.3) at least for 
$\mathit {Pr} \leqslant 1$ and this is what is expected in 3-D RBC with boundaries. Nevertheless, figure 4 demonstrates that in homogeneous 2-D RBC the Nusselt number depends on all the values of Prandtl number we considered, even for 
$\mathit {Pr} > 1$.
Figure 4. The Nusselt number compensated with the ‘classical’ and the ‘ultimate’ scaling in (a,b), respectively. The Prandtl numbers are colour coded in the legend.
3.2. Spectra and cascades
In this section, we examine the multi-scale dynamics of the flows using spectra. To eliminate finite Rayleigh number effects and to have large enough scale separation, we focus on the hyperviscous simulations (i.e. 
$n=4$ in (2.1)) at 
$\mathit {Ra} = 9.4 \times 10^{49}$ and 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$. Results with normal viscosity (i.e. 
$n=1$ in (2.1)) at 
$\mathit {Ra} = 6.2 \times 10^{11}$ are also presented for comparison. According to the BO scaling (Bolgiano Reference Bolgiano1959; Obukhov Reference Obukhov1959), the ratio of the kinetic to the potential energy spectra scales as 
$\hat {E}_u(k) / \hat {E}_\theta (k) \propto k^{-4/5}$, since 
$\hat {E}_u(k) \propto k^{-11/5}$ and 
$\hat {E}_\theta (k) \propto k^{-7/5}$. In figure 5 we plot 
$\hat {E}_u(k) / \hat {E}_\theta (k)$ compensated by 
$k^{4/5}$ for the hyperviscous simulations and, in the inset, for the runs with normal viscosity. Instead of finding a wavenumber range where this scaling is valid, we observe a 
$k^{-0.3}$ power law for all Prandtl numbers considered, suggesting that BO scaling is not followed in our simulations. For the normal viscosity simulations we find the 
$k^{-0.3}$ power law again to be followed, but within a narrower wavenumber range (see inset of figure 5). Figure 5 also shows that, for 
$\mathit {Pr} \ll 1$, the kinetic energy is much larger than the potential energy at large wavenumbers. This is expected as the small scales are dominated by thermal diffusivity when 
$\mathit {Pr} \ll 1$, and so the potential energy is dissipated much more effectively than the kinetic energy. The opposite is true for 
$\mathit {Pr} \gg 1$.
Figure 5. The ratio of the kinetic energy to the potential energy spectra 
$\hat {E}_u(k)/\hat {E}_\theta (k)$ compensated by 
$k^{4/5}$ for the hyperviscous simulations with 
$\mathit {Ra} = 9.4 \times 10^{49}$, 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend. The inset shows the same ratio but for runs with normal viscosity at 
$\mathit {Ra} = 6.2 \times 10^{11}$.
In figures 6(a) and 6(b), we plot the kinetic and potential energy spectra compensated by the power laws proposed by the BO phenomenology. The hyperviscous runs are shown in the main plots, and normal viscosity runs in the insets. From the hyperviscous runs we observe that the compensated kinetic energy spectrum 
$k^{11/5}\hat {E}_u(k)$ is flat within a range of wavenumbers and hence it is close to the BO scaling, with the best fit being 
$\hat {E}_u(k) \propto k^{-2.3}$. However, the BO scaling is not evident from the inset showing the runs with normal viscosity due finite Rayleigh number effects. On the other hand, it is clear that the compensated potential energy spectrum 
$k^{7/5}\hat {E}_\theta (k)$ does not follow the BO scaling for either the hyperviscous and normal viscosity runs, with the best fit being 
$\hat {E}_\theta (k)\propto k^{-1.2}$. The power-law exponents of the spectra we observe differ from those in 3-D RBC with periodic boundary conditions, where the kinetic and potential energy exhibit 
$k^{-5/3}$ power-law spectra (Borue & Orszag Reference Borue and Orszag1997), similar to those observed in passive scalar turbulence (Warhaft Reference Warhaft2000).
Figure 6. The time-averaged energy spectra compensated by best-fit power laws (a,b), and spectral fluxes normalised by the time-averaged dissipation rates (c,d) for the runs with hyperviscosity at 
$\mathit {Ra} = 9.4 \times 10^{49}$, 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend. The insets show the same quantities but for runs with normal viscosity at 
$\mathit {Ra} = 6.2 \times 10^{11}$.
To understand the dynamics of the turbulent cascades, we plot the associated kinetic and potential energy fluxes normalised by the respective time-averaged injection rates of energy due to buoyancy in figures 6(c) and 6(d). The positive potential energy flux suggests a strong direct cascade, while the weak negative kinetic energy flux suggests an inverse cascade. For 
$\mathit {Pr} = 1$ these types of cascades are in agreement with Xie & Huang (Reference Xie and Huang2022). The negative kinetic energy flux peaks at low wavenumbers, while the potential energy flux peaks at high wavenumbers. The inverse cascade of kinetic energy is not affected by the Prandtl number; however, the direct cascade of potential energy moves to higher wavenumbers along with the peak of 
$\varPi _\theta (k)$ as 
$\mathit {Pr}$ increases. We emphasise the wavenumber dependence of kinetic and potential energy fluxes, with 
$\partial _k \varPi _u(k) > 0$ and 
$\partial _k \varPi _\theta (k) > 0$ in the intermediate range of wavenumbers for all of the Prandtl numbers considered.
In figure 7, we present the time-averaged spectra of the magnitudes of the terms in the kinetic and potential energy balances (2.6) for runs with hyperviscosity and three different values of 
$\mathit {Pr}\in \{10^{-2},1,10^2\}$. The corresponding plots with normal viscosity and 
$\mathit {Ra}=6.2\times 10^{11}$ are shown in figure 8. In both figures, the red dots indicate where the inertial flux terms 
$\partial _k\varPi _u$ and 
$\partial _k\varPi _\theta$ become negative. In the kinetic energy balance, we identify three distinct wavenumber ranges, labelled I to III in the plots. In region I, the kinetic energy injected by buoyancy is dissipated by the large-scale friction. In region II, the inertial term balances buoyancy, which is positive for all wavenumbers in this region, i.e.
\begin{equation} \partial_k \varPi_u(k) \approx \alpha g F_B(k) > 0, \end{equation}
with 
$F_B(k) \propto k^{-1.5}$. Equation (3.4) implies the 
$k$ dependence of the kinetic energy flux we see in figure 6(c), which is expected based on the BO phenomenology. Note that region II is largest for small Prandtl numbers, especially in the runs with normal viscosity, as evident from the results presented in figures 8(a)–8(c). In region III, the balances between terms depend on the Prandtl number. For 
$\mathit {Pr} = 10^{-2}$ buoyancy decays rapidly and so small-scale viscous dissipation is balanced by the inertial term, which is negative in this range of wavenumbers. As 
$\mathit {Pr}$ increases, buoyancy becomes more significant in the balance of region III between the small-scale viscous dissipation and the inertial term. For 
$\mathit {Pr} \gtrsim 10^2$, small-scale viscous dissipation seems to be balanced by buoyancy rather than by the inertial term. This effect is shown more clearly in the runs with normal viscosity shown in figure 8(c), where the small-scale dissipation range is much larger.
Figure 7. Time-averaged spectra of the terms in the kinetic energy balance (2.6a) (a–c); and the potential energy balance (2.6b) (d–f) for the runs with hyperviscosity, 
${\mathit {Ra} = 9.4 \times 10^{49}}$, 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$ with 
$\mathit {Pr} = 10^{-2}$ (a,d), 
$\mathit {Pr} = 1$ (b,e) and 
$\mathit {Pr} = 10^{2}$ (c,f). The terms displayed in the legends are defined in (2.7), (2.9) and (2.10). In (a–c) and (d–f), we observe three and two distinct dominant balances, respectively, annotated I–III and A–B. The red dots indicate where the inertial terms become negative.
Figure 8. Time-averaged spectra of the terms in the kinetic energy balance (2.6a) (a–c); and the potential energy balance (2.6b) (d–f) for the runs with normal viscosity, 
${\mathit {Ra} = 6.2 \times 10^{11}}$, 
$\mathit {Rh} = (2{\rm \pi} )^{5/2}$ with 
$\mathit {Pr} = 10^{-2}$ (a,d), 
$\mathit {Pr} = 1$ (b,e) and 
$\mathit {Pr} = 10^{2}$ (c,f). The terms displayed in the legend are defined in (2.7), (2.9) and (2.10). In (a–c) and (d–f), we observe three and two distinct dominant balances, respectively, annotated I–III and A–B. The red dots indicate where the inertial terms become negative.
In the potential energy balance, we identify two distinct wavenumber ranges, labelled A and B in the plots (see (d–f) in figures 7 and 8). In these plots we observe the inertial term to be balanced by buoyancy in region A and by small-scale thermal dissipation in region B. In other words, the potential energy injected by buoyancy in region A is cascaded to larger wavenumbers, where it is dissipated by thermal diffusivity.
We recall that the red dots in figures 7 and 8 show where the inertial terms of the kinetic and potential energy become negative. This sign change occurs primarily in region III for 
$\partial _k \varPi _u(k)$ and region B for 
$\partial _k \varPi _\theta (k)$, the latter corresponding to the large negative gradient of 
$\varPi _\theta (k)$ observed in figure 6 at large wavenumbers. Near the boundary between regions A and B in figures 7 and 8(d–f), we observe that 
$\partial _k \varPi _\theta (k)$ exhibits fluctuations between positive and negative values over these wavenumbers. However, for the majority of the wavenumbers in region A we find that
\begin{equation} \partial_k \varPi_\theta(k) \approx \frac{\Delta T}{L} F_B(k) > 0. \end{equation}
This relation implies the 
$k$ dependence of the potential energy flux we observe in figure 6(d) for all the Prandtl numbers we studied. Note that in the wavenumber range where region II and A overlap, we have 
$\partial _k \varPi _u(k)/\alpha g \approx ({L}/{\Delta T}) \partial _k \varPi _\theta (k) \approx F_B(k)$, such that 
$\varPi _u(k) - ({\alpha g L}/{\Delta T}) \varPi _\theta (k)$ is approximately constant in the intermediate range of wavenumbers as 
$F_B(k)$ cancels out. This is the inertial range of scales, where the viscous, diffusive and friction effects can be neglected and so the energy flux is constant. This is expected for the ideal invariant 
$E_u(t) - ({agL}/{\Delta T} )E_\theta (t)$, which is conserved in the limit of 
$\nu \to 0$, 
$\kappa \to 0$ and 
$\mu \to 0$. In Appendix B we plot the absolute value of the time-averaged spectra 
$|\hat {E}_u(k) - ({agL}/{\Delta T}) \hat {E}_\theta (k)|$ and the flux 
$\varPi _u(k) - ({\alpha g L}/{\Delta T}) \varPi _\theta (k)$ for the hyperviscous runs at 
$Ra = 9.4 \times 10^{49}$, 
$Rh = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend (see figure 10). The spectra 
$|\hat {E}_u(k) - ({agL}/{\Delta T} )\hat {E}_\theta (k)|$ scale like 
$k^{-1}$ and the flux 
$\varPi _u(k) - ({\alpha g L}/{\Delta T} )\varPi _\theta (k)$ is constant as expected for all values of 
$Pr$ considered.
4. Conclusions
In this paper we study the effects of varying the Prandtl and Rayleigh numbers on the dynamics of 2-D RBC without boundaries, i.e. with periodic boundary conditions. For zero and infinite Prandtl number convection, we find that large-scale elevator modes dominate the dynamics, unless large-scale dissipation is made so strong as to suppress turbulent convection. Such elevator modes have long been known (Baines & Gill Reference Baines and Gill1969). In all parameter values simulated, the inequality (2.16) holds, implying that the system admits exact single-mode solutions that grow exponentially without bound. Whether or not the solution blows up must depend on the initial conditions, in a way that is not currently understood. Instead, we focus on finite Prandtl numbers, where we find that nonlinear interactions between modes allow the system to reach turbulent stationary states.
Examining the Prandtl and Rayleigh number dependences of the terms in the kinetic and potential energy balances, we find that the enstrophy scales as 
$\langle {\omega ^2}\rangle \propto \mathit {Pr}^0 \mathit {Ra}^{1/4}$ and hence the small-scale viscous dissipation scales as 
$\epsilon _\nu = \nu \langle {\omega ^2}\rangle \propto \mathit {Pr}^{1/2}\mathit {Ra}^{-1/4}$. On the other hand, we observe that the injection rate of kinetic energy 
$\epsilon _u$ due to buoyancy is effectively independent of both the Prandtl and the Rayleigh number. Using this observation, we find that 
$\mathit {Nu} \propto \mathit {Pr}^{1/2}\mathit {Ra}^{1/2}$, which agrees with the so-called ultimate scaling. The collapse of the data using the ultimate scaling becomes progressively more accurate as 
$\mathit {Ra}$ increases for all values of 
$\mathit {Pr}$ considered.
Analysing the kinetic and potential energy spectral fluxes, we find an inverse cascade of kinetic energy and a direct cascade of potential energy. The inverse cascade is independent of the Prandtl number, while the peak of the potential energy flux moves to higher wavenumbers as 
$\mathit {Pr}$ increases. The kinetic and potential energy fluxes, 
$\varPi _u(k)$ and 
$\varPi _\theta (k)$, are not constant because both of their derivatives in 
$k$ are balanced by the buoyancy term 
$F_B(k)$, which is predominantly positive in the intermediate range of wavenumbers. These two balances imply a positive slope in 
$k$ for the fluxes. On the other hand, the potential energy flux was recently claimed to be constant in the inhomogeneous set-up of 2-D RBC (Samuel & Verma Reference Samuel and Verma2024). Nevertheless, their observations of the potential energy spectrum did not follow the BO scaling. Although we observe no range of wavenumbers where either 
$\varPi _u(k)$ or 
$\varPi _\theta (k)$ is constant, we find that they are connected by the relation 
$\varPi _u(k) - ({\alpha g L}/{\Delta T}) \varPi _\theta (k) \approx \text {constant}$ in the overlap between regions II and A.
The kinetic energy spectra scale as 
$\hat {E}_u(k) \propto k^{-2.3}$, which is close to 
$k^{-11/5}$ behaviour of the BO phenomenology. However, the potential energy spectra scale as 
$\hat {E}_\theta (k) \propto k^{-1.2}$, which deviates from the 
$k^{-7/5}$ scaling predicted by the BO arguments. The deviation from the BO phenomenology is also clear when we test the scaling 
$\hat {E}_u(k) / \hat {E}_\theta (k) k^{4/5} \approx \text {const.}$, which is clearly not followed by our spectra, which scale as 
$\hat {E}_u(k) / \hat {E}_\theta (k) k^{4/5} \propto k^{-0.3}$. For the hyperviscous simulations, the observed power laws in the inertial range of the kinetic and potential energy spectra do not show any dependence on the Prandtl number. The only dependence we observe is at the dissipative range of wavenumbers, where the viscosity and the thermal diffusivity dominate the dynamics. In the spectra from the normal viscosity simulations, the effects of the Prandtl number are more significant due to the comparatively low scale separation. Hence, the wavenumber range over which a power-law behaviour can be observed is truncated.
This study clearly demonstrates the necessity for large scale separation to be able to make clearer conclusions on the spectral dynamics and the power-law exponents of 2-D RBC. This requirement makes similar studies in three dimensions more challenging. The development of a phenomenology where buoyancy acts as a broadband spectral forcing is required to interpret the current observations. Numerical simulations at Prandtl and Rayleigh numbers outside the ranges we investigated, i.e. 
$\mathit {Pr} < 10^{-3}$, 
$\mathit {Pr} > 10^2$ and 
$\mathit {Ra} > 10^{12}$, are challenging but would be of great interest to see if they agree with the hyperviscous simulations we have performed and to provide a more complete picture of the asymptotic regime of buoyancy-driven turbulent convection.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Zero and infinite Prandtl number convection
To simplify the analysis let us write the non-dimensional form of (2.1) in accordance with
\begin{equation} {\boldsymbol{x}} \sim L, \quad t \sim \frac{L^{2n}}{\kappa}, \quad \psi \sim \frac{\kappa}{L^{2n-2}}, \quad \theta \sim \Delta T, \end{equation}which yields
$$\begin{gather} \partial_t \nabla^2 \psi + \{\psi, \nabla^2 \psi\} = \mathit{Ra}\mathit{Pr} \partial_x \theta + \mathit{Pr} {(-\nabla^2)^{n+1}} \psi + \mathit{Rh} \sqrt{\frac{\mathit{Ra}}{\mathit{Pr}}}\, \psi, \end{gather}$$
$$\begin{gather}\partial_t \theta + \{\psi, \theta \} = \partial_x \psi - {(-\nabla^2)^{n}} \theta. \end{gather}$$
As we take 
$\mathit {Pr} \to \infty$, we ensure that 
$\mathit {Rh} \sqrt {\mathit {Ra}/\mathit {Pr}}$ is finite to maintain the effects of the large-scale dissipation. The system (A2) thus reduces to
$$\begin{gather} \mathit{Ra} \partial_x \theta ={-}\left[{(-\nabla^2)^{n+1}} + \mathit{Rh} \sqrt{\frac{\mathit{Ra}}{\mathit{Pr}}}\right]\psi, \end{gather}$$
$$\begin{gather}\partial_t \theta + \{\psi, \theta \} = \partial_x \psi - {(-\nabla^2)^n} \theta. \end{gather}$$In the complementary limit of zero Prandtl number, we have to rescale the variables according to
\begin{equation} \{\psi,\theta,t \} \mapsto \{\mathit{Pr} \,\psi,\mathit{Pr}\,\theta,t/\mathit{Pr} \}, \end{equation}
before letting 
$\mathit {Pr} \to 0$, which removes the time derivative and advective term from the heat transport equation (A2b). As we take 
$\mathit {Pr} \to 0$, we maintain the effects of the large-scale dissipation by ensuring that 
$\mathit {Rh} \sqrt {\mathit {Ra}/\mathit {Pr}}$ is finite. This process reduces the system (A2) to
$$\begin{gather} \partial_t \nabla^2 \psi + \{\psi, \nabla^2 \psi\} = \mathit{Ra} \partial_x \theta + {(-\nabla^2)^{n+1}} \psi + \mathit{Rh} \sqrt{\frac{\mathit{Ra}}{\mathit{Pr}}} \psi, \end{gather}$$
$$\begin{gather}\partial_x \psi = {(-\nabla^2)^n} \theta. \end{gather}$$ In figure 9 we plot time series of the kinetic energy, 
$E_u$ in both the 
$\mathit {Pr} \to 0$ and 
$\mathit {Pr} \to \infty$ limits. We use normal viscosity (i.e. 
$n=1$), we set 
$\mathit {Ra} =10^7$ for four different values of 
$\mathit {Rh} \sqrt {\mathit {Ra}/\mathit {Pr}}$ and the simulations are initialised with random initial data. The flow converges to a statistically steady state only in the case where 
$\mathit {Pr} \to \infty$ and 
$\mathit {Rh} \sqrt {\mathit {Ra}/\mathit {Pr}} = 9 \times 10^6$, a value only 
$10\,\%$ smaller than the maximum value given by (2.16) at which all convection is suppressed. For all other parameter values attempted, an elevator mode takes over the dynamics and the energy grows without bound.
Figure 9. Time series of the kinetic energy, 
$E_u$, in the 
$\mathit {Pr} \to 0$ and 
$\mathit {Pr} \to \infty$ limits with normal viscosity (
$n=1$), at 
$\mathit {Ra} =10^7$ and four different values of 
$\mathit {Rh} \sqrt {\mathit {Ra}/\mathit {Pr}}$. The time series are initiated with random initial conditions.
Figure 10. Time-averaged spectra (a) 
$|E_u(k) - (\alpha g L/\Delta T) E_\theta (k)|$, where a 
$k^{-1}$ power law has been included as a dashed line to guide the eye, and (b) 
$\varPi _u(k) - (\alpha g L/\Delta T) \varPi _\theta (k)$ for the runs with hyperviscosity at 
$Ra = 9.4 \times 10^{49}$, 
$Rh = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend.
We have explored various initial conditions, different from the elevator modes, to find bounded solutions in the extreme Prandtl number limits. However, we are not able to reliably obtain turbulent saturated states in these limits, unless the large-scale dissipation is made very strong, suppressing turbulent convection.
Appendix B. Spectra and flux of the ideal invariant
Here, we plot the absolute value of the time-averaged spectra of the ideal invariant 
$|\hat {E}_u(k) - ({agL}/{\Delta T}) \hat {E}_\theta (k)|$ and the flux 
$\varPi _u(k) - ({\alpha g L}/{\Delta T}) \varPi _\theta (k)$ for the hyperviscous runs at 
$Ra = 9.4 \times 10^{49}$, 
$Rh = (2{\rm \pi} )^{5/2}$ and different Prandtl numbers as colour coded in the legend.
Open access
