1. Introduction
Thermal convection in stellar and planetary interiors and their atmospheres is complex because it is driven by several factors in combination, but researchers often study the idealized model of convection – Rayleigh–Bénard convection (RBC) – where a fluid layer bounded between two horizontal plates is heated and cooled uniformly from the bottom and the top, respectively (Tritton Reference Tritton1977; Siggia Reference Siggia1994; Kadanoff Reference Kadanoff2001; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Chillà & Schumacher Reference Chillà and Schumacher2012). An important parameter governing RBC is the Prandtl number 
$Pr$, which is the ratio of the kinematic viscosity 
$\nu$ and the thermal diffusivity 
$\kappa$ of the fluid. Stellar and planetary convection is often characterized by very small Prandtl numbers; for example, 
$Pr \approx 10^{-6}$ in the Sun (Brandenburg & Subramanian Reference Brandenburg and Subramanian2005; Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020; Garaud Reference Garaud2021) and 
$Pr \approx 10^{-2}\unicode{x2013}10^{-1}$ in the Earth's outer core (Calkins et al. Reference Calkins, Aurnou, Eldredge and Julien2012; Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015; Guervilly, Cardin & Schaeffer Reference Guervilly, Cardin and Schaeffer2019). The Rayleigh number 
$Ra$, which characterizes the strength of the driving buoyancy force relative to viscous and thermal dissipative forces, is another control parameter, and is very high in most natural flows. The ratio of the horizontal and vertical extents of a convective flow is the aspect ratio 
$\varGamma$; a crucial feature of all natural convective flows is that 
$\varGamma \gg 1$, which allows the formation of turbulent superstructures – coherent flow patterns with characteristic scale larger than the depth of the convection layer (Cattaneo, Lenz & Weiss Reference Cattaneo, Lenz and Weiss2001; Rincon, Lignières & Rieutord Reference Rincon, Lignières and Rieutord2005). A possible example is supergranulation on the Sun's surface (Nordlund, Stein & Asplund Reference Nordlund, Stein and Asplund2009; Rincon & Rieutord Reference Rincon and Rieutord2018). Even though RBC incorporates a number of simplifications, such as the Boussinesq approximation (Tritton Reference Tritton1977; Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020), it appears that a study of the flow in extended layers at low Prandtl numbers is highly worthwhile. This is the primary objective of the current work. We are also motivated by the relevance of low-
$Pr$ convection for industrial applications that use liquid metals.
Despite this relevance, low-
$Pr$ turbulent convection has not been explored extensively in experiments mainly because liquid metals such as mercury, gallium and sodium are difficult to handle and optically opaque (Cioni, Ciliberto & Sommeria Reference Cioni, Ciliberto and Sommeria1997; Zürner et al. Reference Zürner, Schindler, Vogt, Eckert and Schumacher2019). Even when these difficulties are circumvented, the lowest Prandtl number that can be explored is 
$Pr \approx 0.006$ for liquid sodium (Horanyi, Krebs & Müller Reference Horanyi, Krebs and Müller1999), which is still some three orders of magnitude higher than in the solar convection zone. Direct numerical simulations (DNS) offer important tools but they, too, are hindered by demanding resolution requirements due to the highly inertial nature of low-
$Pr$ convection (Breuer et al. Reference Breuer, Wessling, Schmalzl and Hansen2004; Schumacher, Götzfried & Scheel Reference Schumacher, Götzfried and Scheel2015; Scheel & Schumacher Reference Scheel and Schumacher2016, Reference Scheel and Schumacher2017; Pandey, Scheel & Schumacher Reference Pandey, Scheel and Schumacher2018; Zwirner et al. Reference Zwirner2020). The fact that the required computational power increases with increasing aspect ratio as 
$\varGamma ^2$ further limits numerical investigations. Yet, Pandey et al. (Reference Pandey, Scheel and Schumacher2018) were able to perform DNS of RBC in a 
$\varGamma = 25$ domain by achieving 
$Pr$ as low as 0.005 at 
$Ra = 10^5$ to study turbulent superstructures. In the present work, we significantly extend the parameter range from Pandey et al. (Reference Pandey, Scheel and Schumacher2018) and Fonda et al. (Reference Fonda, Pandey, Schumacher and Sreenivasan2019) by further decreasing the Prandtl number fivefold, while also increasing the Rayleigh number by two orders of magnitude. The highest Reynolds number achieved in the present work is nearly 
$5.6 \times 10^4$, requiring massively parallel DNS on computational grids of more than 
$5\times 10^{11}$ points.
This work has three main goals: (i) report trends of heat and momentum transfer with respect to Prandtl number over nearly 4 orders of magnitude; (ii) assess the closeness of small-scale statistical properties in the bulk of the convection layer to the classical Kolmogorov-type behaviour (Kolmogorov Reference Kolmogorov1941b; Frisch Reference Frisch1995). This assessment comprises energy spectra, a test of the 4/5th law and an investigation of the dissipative anomaly (Sreenivasan Reference Sreenivasan1984, Reference Sreenivasan1998) in turbulent convection flow with boundaries; (iii) analyse the large-scale circulation patterns, the turbulent superstructures of convection, particularly their trends with decreasing 
$Pr$. Note that the convective flows in the Earth's and stellar interiors are also associated with (differential) rotation and magnetic fields that can lead to strong departures from local isotropy at larger and intermediate scales (Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015), but we shall here focus on the influence of low Prandtl number. These DNS series will thus provide a unique data base for the parametrization of turbulent transport in mesoscale configurations characterized by a degree of large-scale order, with well-resolved thermal and kinetic energy dissipation rates.
It is becoming increasingly clear that many properties of low-
$Pr$ convective flows differ from those at moderate and high Prandtl numbers. For example, the efficacy of low-
$Pr$ flows in transporting heat is lower, and that in transporting momentum higher, than in high-
$Pr$ flows; the disparity between the two increases as 
$Pr$ is lowered (Scheel & Schumacher Reference Scheel and Schumacher2017). Due to high (low) thermal (momentum) diffusivity, low-
$Pr$ convection exhibits coarser thermal structures but the length scales in the velocity field have a broader distribution. This results in an enhanced separation between the energy injection and energy dissipation scales in low-
$Pr$ convective flows (Schumacher et al. Reference Schumacher, Götzfried and Scheel2015). In this regime, the kinetic energy spectrum has been observed to approximate the classical Kolmogorov scaling (Kolmogorov Reference Kolmogorov1941b) with the 
$-5/3$ power in the inertial range (Lohse & Xia Reference Lohse and Xia2010; Mishra & Verma Reference Mishra and Verma2010; Bhattacharya, Verma & Samtaney Reference Bhattacharya, Verma and Samtaney2021). Here, we analyse the kinetic energy spectra in the bulk region of the flow and show that they indeed exhibit the classical Kolmogorov scaling with an inertial range, showing no tendency towards Bolgiano scaling (Bolgiano Reference Bolgiano1959), according to which, the conversion of the kinetic to potential energy leads to the steeper 
$k^{-11/5}$ scaling with the wavenumber k in the inertial range (Lohse & Xia Reference Lohse and Xia2010; Verma, Kumar & Pandey Reference Verma, Kumar and Pandey2017; Verma Reference Verma2018).
Turbulent superstructures of convection can be characterized by a typical spatial scale 
$\lambda$ and a temporal scale 
$\tau$; finer scales evolve much faster than 
$\tau$. Thus, the scales 
$\lambda$ and 
$\tau$ help distinguish the coarse and gradually evolving large-scale patterns from the finer (and faster) turbulent fluctuations (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Krug, Lohse & Stevens Reference Krug, Lohse and Stevens2020). The characteristic scales of superstructures have been observed to depend on 
$Pr$ and 
$Ra$ (Hartlep, Tilgner & Busse Reference Hartlep, Tilgner and Busse2003, Reference Hartlep, Tilgner and Busse2005; von Hardenberg et al. Reference von Hardenberg, Parodi, Passoni, Provenzale and Spiegel2008; Bailon-Cuba, Emran & Schumacher Reference Bailon-Cuba, Emran and Schumacher2010; Emran & Schumacher Reference Emran and Schumacher2015; Pandey et al. Reference Pandey, Scheel and Schumacher2018; Schneide et al. Reference Schneide, Pandey, Padberg-Gehle and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Fonda et al. Reference Fonda, Pandey, Schumacher and Sreenivasan2019; Green et al. Reference Green, Vlaykov, Mellado and Wilczek2020; Krug et al. Reference Krug, Lohse and Stevens2020; Lenzi, von Hardenberg & Provenzale Reference Lenzi, von Hardenberg and Provenzale2021; Pandey, Schumacher & Sreenivasan Reference Pandey, Schumacher and Sreenivasan2021), as well as on the thermal boundary conditions at the horizontal top and bottom plates (Vieweg, Scheel & Schumacher Reference Vieweg, Scheel and Schumacher2021a). The characteristic length scale of superstructures, which is nearly twice the depth 
$H$ of the convection layer at the onset of convection, increases with increasing 
$Ra$ (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018). This dependence on 
$Pr$ is complex, with 
$\lambda (Pr)$ showing a peak near 
$Pr \approx 7$ and decreasing as 
$Pr$ departs from this value (Pandey et al. Reference Pandey, Scheel and Schumacher2018). This decreasing trend of 
$\lambda (Pr)$ continues to hold up to a 
$Pr \approx 0.005$ below which the scales seem to level off at a wavelength of 
$\lambda \gtrsim 3H$.
The remainder of this article is organized as follows. In § 2, we briefly describe the DNS and note the parameter space explored. In § 3, we discuss the flow structures and the scaling of the global transport of heat and momentum, and study in § 4 the vertical profiles of temperature, convective and diffusive transports, as well as dissipation rates. Statistical properties of the flow such as kinetic energy spectra and third-order structure functions are examined in § 5 for their compatibility with Kolmogorov forms. The characterization of turbulent superstructures is presented in § 6. We conclude the main findings in § 7, and present an outlook. Appendices A and B deal with specific tests of sufficient resolution. Appendix C discusses technical detail of the estimates of length and velocity scales of superstructures.
2. Details of DNS
We perform DNS of RBC in a closed rectangular domain with square cross-section of length 
$L = 25H$, where 
$H$ is the depth of the convection layer. We solve the following non-dimensionalized equations incorporating the Oberbeck–Boussinesq (OB) approximation:
$$\begin{gather} \frac{\partial {\boldsymbol u}}{\partial t} + ({\boldsymbol u} \boldsymbol{\cdot} {\boldsymbol \nabla}) {\boldsymbol u} ={-}{\boldsymbol \nabla}p + T \hat{{\boldsymbol z}} + \sqrt{\frac{Pr}{Ra}} \, \nabla^2 {\boldsymbol u}, \end{gather}$$
$$\begin{gather}\frac{\partial T}{\partial t} + ({\boldsymbol u} \boldsymbol{\cdot} {\boldsymbol \nabla}) T = \frac{1}{\sqrt{PrRa}} \, \nabla^2 T, \end{gather}$$
$$\begin{gather}{\boldsymbol \nabla} \boldsymbol{\cdot} {\boldsymbol u} = 0. \end{gather}$$
Here, 
${\boldsymbol u} \equiv (u_x,u_y,u_z)$, 
$T$ and 
$p$ are the velocity, temperature and pressure fields, respectively. The Prandtl number 
$Pr$ and the Rayleigh number 
$Ra = \alpha g \Delta T H^3/\nu \kappa$, where 
$\alpha$ is the coefficient of thermal expansion of the fluid, 
$g$ is the acceleration due to gravity and 
$\Delta T$ is the imposed temperature difference between the horizontal plates. We have used the layer's depth 
$H$, the free-fall velocity 
$u_f = \sqrt {\alpha g \Delta T H}$, the free-fall time 
$t_f = H/u_f$ and 
$\Delta T$ as the non-dimensionalizing length, velocity, time and temperature scales, respectively. We use the no-slip condition on all boundaries. We employ the isothermal condition on the horizontal plates and the adiabatic condition on the sidewalls. This allows direct comparison with other simulations at higher 
$Pr$ (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Fonda et al. Reference Fonda, Pandey, Schumacher and Sreenivasan2019) and controlled laboratory experiments, such as those of Moller et al. (Reference Moller, Käufer, Pandey, Schumacher and Cierpka2022), in exactly the same setting to enable a consistent analysis across the 
$Ra$–
$Pr$ parameter plane.
We use two different solvers for the simulations. For moderate and large 
$Pr$, we use a spectral element solver Nek5000 (Fischer Reference Fischer1997), where the flow domain is divided into a finite number of elements 
$N_e$. The Lagrangian interpolation polynomials of order 
$N$ are further used to expand the turbulence fields within each element (Scheel, Emran & Schumacher Reference Scheel, Emran and Schumacher2013), thus resulting in a total of 
$N_e N^3$ mesh cells in the entire flow domain. As low-
$Pr$ convection is dominated by inertial forces, the resulting flow acquires increased fine structure whose resolution requires more extensive computational resources (Schumacher et al. Reference Schumacher, Götzfried and Scheel2015). Therefore, we performed those simulations using a second-order finite difference solver that requires significantly less working memory at a given grid size. Here, the flow domain is divided into 
$N_x \times N_y \times N_z$ non-uniform mesh cells (Krasnov, Zikanov & Boeck Reference Krasnov, Zikanov and Boeck2011; Liu, Krasnov & Schumacher Reference Liu, Krasnov and Schumacher2018). We have verified that the results obtained from both the solvers agree well with each other by performing two simulations for 
$Pr = 0.005, Ra = 10^5$ and 
$Pr = 0.7, Ra = 10^7$ using both solvers. We refer to Appendix A for a direct comparison of globally averaged and horizontally averaged convective heat fluxes and dissipation rates from the two solvers. Important parameters of all the simulations are provided in table 1.
Important parameters of the simulations in a rectangular box of 
$\varGamma = 25$ with square cross-section; the numbers of mesh cells, 
$N_x \times N_y \times N_z$ and 
$N_e \times N^3$, are for the finite difference and spectral element solvers, respectively; 
$N_{tot}$ represents the total number of mesh cells in units of a billion; 
$Nu$ and 
$Re$ are the volume and time averaged Nusselt and Reynolds numbers, respectively; 
$u_{rms}$ is the root-mean-square velocity computed over the entire volume; 
$\lambda$ and 
$\tau$ are, respectively, the characteristic length and time scales of turbulent superstructures. Runs with superscript 
$a$ are taken from Pandey et al. (Reference Pandey, Scheel and Schumacher2018), while those with superscript 
$b$ are taken from Fonda et al. (Reference Fonda, Pandey, Schumacher and Sreenivasan2019).
3. Flow morphology and global transport
3.1. Structures of velocity and temperature fields
Because the time scales of heat and momentum diffusion processes are very different in low-
$Pr$ convection, the temperature field shows coarser structures than the velocity field. This is illustrated in figure 1, which displays the instantaneous temperature, vertical velocity and local turbulent kinetic energy fields in the mid-horizontal plane, 
$z = H/2$, for the biggest simulations with 
$Pr = 0.001, Ra = 10^7$. Panels (a,d,g) show the fields in the entire cross-section, and (b,c,e,f,h,i) depict the marked magnifications to highlight small-scale structures. The flow pattern of the temperature and vertical velocity fields show similarities at large scales, but the velocity field also consists of very fine structures compared with the highly diffusive temperature field. The finest scale of the turbulent velocity field, denoted as the Kolmogorov scale 
$\eta$, is estimated as 
$\eta = (\nu ^3/\langle \varepsilon _u\rangle _{V,t})^{1/4}$. Here, 
$\langle \varepsilon _u\rangle _{V,t}$ is the combined volume–time average of the kinetic energy dissipation rate field per unit mass, computed at each point by
\begin{equation} \varepsilon_u({\boldsymbol x},t) = \frac{\nu}{2} \sum_{i,j=1}^3\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)^2, \end{equation}
with 
$u_i$ representing the velocity component in the direction of coordinate 
$x_i$. The finest scale of the temperature field is either the Corrsin scale 
$\eta _C=\eta /Pr^{3/4}$, which marks the end of the inertial–convective range for 
$Pr<1$, or the Batchelor scale 
$\eta _B=\eta /Pr^{1/2}$, which marks the end of the viscous–convective range for 
$Pr>1$; see e.g. Sreenivasan & Schumacher (Reference Sreenivasan and Schumacher2010). It is clear that 
$\eta < \eta _C$ when 
$Pr < 1$, and 
$\eta _B < \eta$ when 
$Pr>1$. Thus, the finest scales in the flow at hand are either 
$\eta$ for 
$Pr<1$ or 
$\eta _B$ for 
$Pr>1$.
Turbulent superstructures of convection in a low-
$Pr$ flow with 
$Pr = 0.001$ and 
$Ra = 10^7$. The panels represent instantaneous temperature fields (a–c), vertical velocity (d–f) and turbulent kinetic energy (g–i) in the midplane. In this low-
$Pr$ flow, the thermal energy is primarily contained in large-scales, whereas the kinetic energy is distributed over a broad range of scales. Panels (a,d,g), (b,e,h) and (c, f,i) represent fields of view that are 
$25H \times 25H$, 
$6.25H \times 6.25H$ and 
$1.56H \times 1.56H$, respectively.
The Corrsin scale is nearly 178 times larger than the Kolmogorov scale for 
$Pr = 0.001$. This large difference is clearly visible in figure 1(c, f), where the temperature and vertical velocity fields are shown in a small cross-section of size 
$1.56H \times 1.56H$. The figures reveal that the smallest length scale of the thermal structures – the length scale over which the temperature variation is significant – is of the order of 
$H$, whereas that for velocity structures is much finer. We also find that the dominant structures in the kinetic energy field resemble those in the vertical velocity because of its dominance in the midplane (Pandey et al. Reference Pandey, Scheel and Schumacher2018). This wide range of length scales present in low-
$Pr$ convection engenders a broad inertial range in kinetic energy spectrum, which will be discussed in § 5.2.
To see the effects of 
$Pr$ on flow structures, we show the temperature and the vertical velocity fields for 
$Pr = 0.001$, 
$0.7$ and 
$7$ at 
$Ra = 10^7$ in figure 2. To accentuate small structures, the fields are shown in a quarter (in linear dimension) of the entire cross-section. With increasing 
$Pr$, increasingly finer thermal structures are generated due to decreasing thermal diffusivity. On the other hand, the velocity variation becomes progressively regular as the viscosity increases (or the Reynolds number decreases) with increasing 
$Pr$.
Temperature (a–c) and vertical velocity (d–f) in midplane for 
$Ra = 10^7$. Panels (a,d) are for 
$Pr = 0.001$, (b,e) for 
$Pr = 0.7$ and (c, f) for 
$Pr = 7$. Fields are displayed in a magnified region of dimensions 
$6.25 H \times 6.25 H$ around the centre. Finer temperature contours can be observed for higher 
$Pr$, panel(c), compared with lower 
$Pr$, panel (a). The velocity field for lower 
$Pr$, panel (d), exhibits much finer structures than those in the other two panels, because the Reynolds number decreases as 
$Pr$ increases.
3.2. Heat and momentum transport laws
Convection at low Prandtl numbers differs from its high-
$Pr$ counterpart by reduced heat transport and enhanced momentum transport (Schumacher et al. Reference Schumacher, Götzfried and Scheel2015; Scheel & Schumacher Reference Scheel and Schumacher2016, Reference Scheel and Schumacher2017; Pandey et al. Reference Pandey, Scheel and Schumacher2018; Zürner et al. Reference Zürner, Schindler, Vogt, Eckert and Schumacher2019; Zwirner et al. Reference Zwirner2020). Heat transport is quantified by the Nusselt number 
$Nu$, defined as the ratio of the total heat transport to that by conduction alone. It is computed as
\begin{equation} Nu = 1 + \sqrt{RaPr} \, \langle u_z T \rangle_{V,t} , \end{equation}
where 
$\langle \,\cdot \, \rangle _{V,t}$ denotes again the average over the entire simulation domain and time. We compute 
$Nu$ for all simulations and plot them, for fixed 
$Ra$, as a function of 
$Pr$ in figure 3(a); 
$Nu$ increases up to 
$Pr = 0.7$ but does not change significantly thereafter. A similar trend has also been reported in literature for convection in 
$\varGamma \approx 1$ domains (Verzicco & Camussi Reference Verzicco and Camussi1999; Schmalzl, Breuer & Hansen Reference Schmalzl, Breuer and Hansen2004; van der Poel, Stevens & Lohse Reference van der Poel, Stevens and Lohse2013) and also for 
$\varGamma = 0.1$ (Pandey & Sreenivasan Reference Pandey and Sreenivasan2021). Figure 3(a) indicates that the molecular diffusion becomes an increasingly dominant mode of heat transport as 
$Pr$ decreases. For 
$Ra = 10^5$ and 
$10^6$, we do best fits to the data for 
$Pr \leq 0.7$. The transport laws 
$Nu(Pr)$ for both Rayleigh numbers are given, including error bars, in table 2. In summary, we find that the Nusselt number is consistent with the power law 
$Nu\sim Pr^{0.19}$ for 
$Ra = 10^5$ and 
$Nu\sim Pr^{0.18}$ for 
$Ra=10^6$. These power-law exponents in the extended convection domain are within the range observed in RBC with 
$\varGamma \lesssim 1$, as shown in Pandey & Sreenivasan (Reference Pandey and Sreenivasan2021), where more discussion of the 
$Nu-Pr$ scaling exponent can be found.
(a) Nusselt number for 
$\varGamma = 25$ increases with increasing 
$Pr$ for low Prandtl numbers but does not change much for 
$Pr \geq 0.7$, consistent with Pandey & Sreenivasan (Reference Pandey and Sreenivasan2021). (b) The value of 
$Nu$ as a function of 
$Ra$ increases approximately as 
$Ra^{0.29}$. (c) The Reynolds number based on the root-mean-square velocity decreases with increasing 
$Pr$, with power law exponents consistent with others in the literature. (d) The value of 
$Re$ as a function of 
$Ra$ increases as 
$Ra^{0.5}$ for all the three cases. The dashed lines in panels (b,d) are best fits, summarized in table 2. Solid lines in panels (a,c) are predictions of 
$Nu$ and 
$Re$ from the Grossmann–Lohse theory. Legends in panels (a,b) also apply to the corresponding panels (c,d).
Summary of scaling relations for global heat and momentum transports as functions of 
$Ra$ and 
$Pr$. Note that the scaling laws with respect to 
$Ra$ have been obtained by fits to three data points only; it is clear that more definitive results require larger number of data points.
Normalized values of global heat transport in RBC increase with increasing thermal driving and the rate of increase depends on the Prandtl number (Scheel & Schumacher Reference Scheel and Schumacher2016, Reference Scheel and Schumacher2017). We plot 
$Nu$ for 
$Pr = 0.001, 0.7$ and 7 against 
$Ra$ in figure 3(b), which shows that 
$Nu$ values for 
$Pr = 0.7$ and 7 are similar. There are only three data points, which would not be adequate to establish a new result. Nevertheless, power-law fits to those three points serve to supplement existing results. We obtain approximately 
$Nu \sim Ra^{0.29}$ (see table 2). The exponents for all the three Prandtl numbers are essentially similar and agree with those observed in convection for 
$\varGamma \sim 1$ (Bailon-Cuba et al. Reference Bailon-Cuba, Emran and Schumacher2010; Stevens, Lohse & Verzicco Reference Stevens, Lohse and Verzicco2011; Scheel, Kim & White Reference Scheel, Kim and White2012; Scheel & Schumacher Reference Scheel and Schumacher2014, Reference Scheel and Schumacher2016). It is interesting that the exponent for 
$Pr = 0.001$ is not lower compared with that for 
$Pr \geq 0.7$. A slightly lower scaling exponent of 
$0.27\pm 0.01$ was reported from simulations in closed cylinders for 
$\varGamma =1$ (Scheel & Schumacher Reference Scheel and Schumacher2017) at 
$Pr=0.021$. Recent experiments in strongly turbulent liquid metal convection by Schindler et al. (Reference Schindler, Eckert, Zürner, Schumacher and Vogt2022), at nearly the same Prandtl number but for Rayleigh numbers up to 
$Ra=5\times 10^9$, in a cylinder with 
$\varGamma =1/2$, reported an even smaller scaling exponent of 0.124. It is possible that the constrained large-scale flow in a closed cylinder affects the scaling exponent at low and moderate Rayleigh numbers, see discussion by Pandey & Sreenivasan (Reference Pandey and Sreenivasan2021).
The Reynolds number 
$Re$ quantifies the momentum transport in RBC. We compute it with 
$H$ and 
$u_{rms}$ as the relevant length and velocity scales, as
\begin{equation} Re =
u_{rms} \, \sqrt{\frac{Ra}{Pr}} \quad \mbox{where}\ u_{rms}
= \sqrt{\langle u_i^2 \rangle_{V,t}}
\end{equation}
is the root-mean-square (r.m.s.) velocity. The Reynolds number as a function of 
$Pr$ is plotted in figure 3(c), which reveals that, for fixed 
$Ra$, the flow loses its effectiveness in transporting momentum as 
$Pr$ increases (Käpylä Reference Käpylä2021). In thermal convection, the power-law exponent of 
$Re-Pr$ scaling depends on the range of 
$Pr$; the Reynolds number decreases with increasing 
$Pr$ even when the flow is dominated by inertia. The detailed power-law fits can be found in table 2. As a summary, we get 
$Re \sim Pr^{-0.62}$ and 
$Re \sim Pr^{-0.65}$ for 
$Ra = 10^5$ and 
$Ra = 10^6$, respectively. As for the Nusselt number, the exponents of the 
$Re-Pr$ scaling agree with those observed for 
$\varGamma \lesssim 1$ (Verzicco & Camussi Reference Verzicco and Camussi1999; Pandey & Sreenivasan Reference Pandey and Sreenivasan2021).
The Reynolds number variation with 
$Ra$, plotted in figure 3(d), shows that 
$Re$ is consistently higher for lower Prandtl numbers, manifesting in the enhanced prefactors of the power-law fits which are summarized in table 2. The best fits yield 
$Re \sim Ra^{0.53}$, 
${Re \sim Ra^{0.49}}$ and 
$Re \sim Ra^{0.54}$ for 
$Pr = 0.001, 0.7$ and 7, respectively. Note that the Reynolds number based on the free-fall velocity scales as 
$Ra^{0.50}$. Thus, these scaling exponents suggest that the free-fall velocity for a fixed 
$Pr$ does not depend strongly on 
$Ra$ (see table 1). The scaling exponents are in the same range as in several other studies in the past; the exponent does not, however, decrease with decreasing Prandtl number as found by Scheel & Schumacher (Reference Scheel and Schumacher2017). This might be the result of differences in the aspect ratio, but we reiterate that the present fits use only three data points.
Attempts have been made to predict the global transports in RBC as a function of the control parameters (Shraiman & Siggia Reference Shraiman and Siggia1990; Grossmann & Lohse Reference Grossmann and Lohse2000; Pandey & Verma Reference Pandey and Verma2016). Grossmann & Lohse (Reference Grossmann and Lohse2000) assumed the existence of a large-scale circulation of the order of the size of the convection cell, and proposed a set of coupled equations relating 
$Nu$ and 
$Re$ as functions of 
$Ra$ and 
$Pr$ (Grossmann & Lohse Reference Grossmann and Lohse2001). The equations also include a set of constant coefficients, whose values depend on the aspect ratio of the domain. Using the coefficients provided in Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013) for 
$\varGamma \approx 1$ RBC, we compute 
$Nu$ as a function of 
$Pr$ from the Grossmann–Lohse model and show them as solid curves in figure 3(a). The Nusselt numbers thus estimated are somewhat higher than those computed from the DNS for 
$Ra = 10^5$. On the low-
$Pr$ end, this may be attributed to the fact that the temperature fields for these parameter pairs are dominated by diffusion and barely mixed in the bulk. However, the agreement is better for 
$Ra \geq 10^6$, which indicates that the heat transport in our extended cell is not much different from that in 
$\varGamma = 1$ cells. We also plot 
$Re(Pr)$ from the Grossmann–Lohse model in figure 3(c), and find that there is fair agreement; see also Verma (Reference Verma2018).
4. Vertical profiles across the convection layer
4.1. Temperature and heat flux fields
In the conductive equilibrium state, the vertical temperature gradient is a constant; inhomogeneities in the horizontal directions arise in the convective state, leading to a modification of the linear temperature profile. We compute the mean temperature profiles 
$\langle T \rangle _{A,t}(z)$ and plot them in figure 4. Here, 
$\langle \,\cdot \, \rangle _{A,t}$ stands for the averaging over the entire horizontal cross-section of 
$A=25H\times 25H$ at a fixed height 
$z$ and the full time interval. In a turbulent convective flow, almost the entire temperature drop occurs within the thermal boundary layers (BLs) on the horizontal plates, while the bulk of the flow outside these BLs remains nearly isothermal (and thus well mixed). Figure 4(a) exhibits this feature. However, the slope of the temperature profile in the midplane increases as 
$Pr$ decreases. We plot the profiles for 
$Pr = 0.001$ for all the Rayleigh numbers in figure 4(b). The profile for 
$Ra = 10^5$ departs only weakly from the linear conduction profile despite a high Reynolds number of the flow. The temperature gradient in the central plane decreases with increasing 
$Ra$, and even a Rayleigh number of 
$10^7$ is not enough to generate a well-mixed temperature field in the bulk region for this very low 
$Pr$.
Horizontal and time averages of temperature as a function of the depth for simulations at (a) 
$Ra = 10^6$ and (b) 
$Pr = 0.001$. A well-mixed isothermal region away from the walls occurs only for 
$Pr \geq 0.7$ in (a), whereas a significant temperature gradient in the central region occurs for lower 
$Pr$. Dashed black line in panel (b) corresponds to the dimensionless conduction temperature profile 
$T_{cond} = 1 - z$.
In OB convection, the temperature averaged over the entire flow domain is 
$\Delta T/2$ but fluctuates at each point in the flow. We decompose the temperature field into its mean and fluctuation as
\begin{equation} T({\boldsymbol x},t) = \langle T \rangle_{A,t}(z) + \theta({\boldsymbol x},t) . \end{equation}
Even though the temperature field becomes increasingly diffusive as 
$Pr$ decreases, the fluctuations increase with decreasing 
$Pr$; see figure 5(a) for 
$Ra = 10^7$. The depth variation is captured by the planar temperature fluctuation computed as
\begin{equation} \theta_{rms}(z) = \sqrt{ \langle [T - \langle T \rangle_{A,t}(z)]^2 \rangle_{A,t} } = \sqrt{ \langle T^2(z) \rangle_{A,t} - \langle T(z) \rangle_{A,t}^2}. \end{equation}(a) The r.m.s. temperature fluctuation profiles averaged over the top and bottom halves varying with the distance from the plate for 
$Ra = 10^7$. The peaks in 
$\theta _{rms}(z)$ occur near the thermal BL edges, which are indicated by dashed vertical lines. (b) Vertical profiles of the convective heat flux for 
$Pr = 0.001$. Dashed horizontal lines indicate the global heat flux (
$Nu$) for each case. The convective flux vanishes at the top and bottom plates and is largest in the central plane.
Figure 5(a) shows that 
$\theta _{rms}(z)$ vanishes at the plate due to the imposed isothermal boundary condition. With increasing distance from the bottom plate, however, the strength of fluctuations increases within the thermal BL region. The maxima in 
$\theta _{rms}(z)$ profiles occur near the edge of the thermal BL (computed as 
$0.5H/Nu$), marked as dashed vertical lines in figure 5(a). This suggests that the thermal plumes retain their temperature, while the temperature of the ambient fluid decreases (increases) with increasing distance from the bottom (top) plate. This leads to an increasing contrast between the two components of the flow and is reflected as an increasing 
$\theta _{rms}(z)$ within the BL region (Pandey Reference Pandey2021). In the bulk region, however, 
$\theta _{rms}$ decreases with distance from the plate because the plumes do not retain their identity and begin to mix with the bulk fluid.
The heat transport occurs due to convective as well as diffusive processes, with their ratio varying with depth. To get the total heat flux in a horizontal plane, we average the temperature equation (2.2) in horizontal directions and in time, which leads to
\begin{equation} Nu(z) = \sqrt{RaPr} \langle u_z T \rangle_{A,t} - \frac{\partial \langle T \rangle_{A,t}}{\partial z}=\text{const.} \end{equation}
It is clear from the temperature profiles in figure 4 that the diffusive contribution 
$- \partial \langle T \rangle _{A,t}/\partial z$ should be small in the well-mixed bulk region – increasing towards the plates and becoming largest at the plates. The variation of the convective heat flux 
$\sqrt {RaPr} \langle u_z T \rangle _{A,t}$ with depth in figure 5(b) for 
$Pr = 0.001$ confirms this expectation. The magnitudes of the globally averaged heat flux are indicated as dashed horizontal lines in figure 5(b), showing that the diffusive flux (the distance between the solid curves and the corresponding dashed horizontal lines) is not negligible even in the central region for 
$Ra \leq 10^6$. The diffusive component dominates the total heat flux in the central region for 
$Ra = 10^5$, which is consistent with the highly inefficient convective heat transport; see table 1. However, for 
$Ra = 10^7$, the diffusive contribution diminishes in the central plane. Thus, as 
$Pr$ becomes smaller, one requires increasing 
$Ra$ before turbulent processes become important.
4.2. Thermal and kinetic energy dissipation rates
While the mean temperature 
$\langle T\rangle _{A,t}(z)$ varies sharply near the horizontal plates and weakly in the central region, the vertical mean profile of the thermal dissipation rate field, which is the rate of loss of thermal variance that is computed pointwise by
\begin{equation} \varepsilon_T({\boldsymbol x},t) = \kappa \left[ \left( \frac{\partial T}{\partial x} \right)^2 + \left( \frac{\partial T}{\partial y} \right)^2 + \left( \frac{\partial T}{\partial z} \right)^2 \right], \end{equation}is higher in the vicinity of the horizontal plates and decreases towards the centre (Scheel & Schumacher Reference Scheel and Schumacher2016). We also compute the thermal dissipation rate field defined as
\begin{equation} \varepsilon_\theta({\boldsymbol x},t) = \kappa \left[ \left( \frac{\partial \theta}{\partial x} \right)^2 + \left( \frac{\partial \theta}{\partial y} \right)^2 + \left( \frac{\partial \theta}{\partial z} \right)^2 \right] \end{equation}
to quantify the spatial variation of the temperature fluctuations. The mean profile of the thermal dissipation rate 
$\langle \varepsilon _\theta \rangle _{A,t}(z)$ is plotted in figure 6(a) for 
$Ra = 10^7$. Note that the vertical mean profiles of 
$\varepsilon _T$ and 
$\varepsilon _\theta$ are related by
\begin{equation} \langle \varepsilon_T \rangle_{A,t}(z) = \langle \varepsilon_\theta \rangle_{A,t}(z) + \varepsilon_{\langle T \rangle}(z), \end{equation}
where 
$\varepsilon _{\langle T \rangle } = \kappa \,({\rm d} \langle T \rangle _{A,t}/{\rm d} z)^2$ is the dissipation rate corresponding to the mean temperature profile (Emran & Schumacher Reference Emran and Schumacher2008). In convective flows with well-developed thermal BLs, 
$\varepsilon _{\langle T \rangle }$ contributes primarily to the BLs and negligibly in the bulk. This rapid decrease of 
$\varepsilon _{\langle T \rangle }$ outside the thermal BL region shows a shallow kink in the profiles of 
$\langle \varepsilon _\theta \rangle _{A,t}(z)$. In figure 6(a), we indicate the thermal BL thicknesses for 
$Pr = 0.7$ and 
$Pr = 7$ as dashed vertical lines, and note that the kinks are observed near the edge of the thermal BL. The kink does not appear for 
$Pr = 0.001$ due to the absence of well-developed thermal BLs.
Variation of the horizontally averaged (a) thermal and (b) kinetic energy dissipation rates in the vertical direction for 
$Ra = 10^7$. The profiles are further averaged over the top and bottom halves of the domain to improve the statistics. The dissipation profiles are largest at the plates and decrease towards the central plane; however, 
$\varepsilon _u(z)$ is nearly uniform in the bulk region. The (indistinguishable) dashed vertical lines in panel (a) indicate the edges of the thermal BLs for 
$Pr = 0.7$ and 
$Pr = 7$.
We find that 
$\langle \varepsilon _\theta \rangle _{A,t}(z)$ increases with decreasing 
$Pr$. This is because the volume-averaged thermal dissipation rate is related to the global heat transport (Shraiman & Siggia Reference Shraiman and Siggia1990) as
\begin{equation} \langle \varepsilon_T \rangle_{V,t} = \frac{Nu}{\sqrt{RaPr}} . \end{equation}
As we observe 
$Nu \sim Pr^{0.2}$, this leads to 
$\langle \varepsilon _T \rangle _{V,t} \sim Pr^{-0.3}$ for a fixed 
$Ra$. Thus, the decrease of the thermal dissipation rate with increasing 
$Pr$ is consistent with the 
$Pr$-dependence of the Nusselt number.
We now plot in figure 6(b) the profiles of the viscous dissipation rate defined in (3.1). Similar to 
$\langle \varepsilon _\theta \rangle _{A,t}(z)$, the largest values of 
$\langle \varepsilon _u\rangle _{A,t}(z)$ are found near the horizontal plate owing to the strongly varying velocity field in the vicinity of the plates. Further, the variation of the profiles 
$\langle \varepsilon _u\rangle _{A,t}(z)$ in the bulk region is almost negligible compared with that in the viscous BL region near the plates. Figure 6(b) shows that 
$\langle \varepsilon _u\rangle _{A,t}(z)$ increases with decreasing 
$Pr$ for all 
$z$. Note that the globally averaged viscous dissipation rate is related to the Nusselt number as
\begin{equation} \langle \varepsilon_u \rangle_{V,t} = \frac{Nu-1}{\sqrt{RaPr}} , \end{equation}
and therefore, 
$\langle \varepsilon _u \rangle _{V,t}$ should decrease with increasing 
$Pr$.
5. Characterization of the turbulence in the bulk
5.1. Isotropy in the midplane
Vorobev et al. (Reference Vorobev, Zikanov, Davidson and Knaepen2005) used the ratios
\begin{equation} G_{ij}=\frac{\langle(\partial u_i/\partial z)^2\rangle (1+\delta_{iz})}{\langle(\partial u_i/\partial x_j)^2\rangle (1+\delta_{ij})} \quad\text{with}\ i,j=x,y,z\, \end{equation}
to determine the degree of anisotropy on the level of second-order derivative moments. Flows with no variation in the vertical direction 
$z$ yield 
$G_{ij}\to 0$ (and are thus anisotropic), while 
$G_{ij} = 1$ for perfectly isotropic flows. The coefficient 
$G_{11}$, relating the in-plane derivative to a transverse derivative with respect to the vertical direction is summarized in three horizontal planes in table 3; 
$G_{11}$ remains nearly unity in the bulk region for 
$0.1 \leq z/H \leq 0.9$, but significant departures are found near the horizontal plates. Similar amplitudes follow for other combinations; see also Nath et al. (Reference Nath, Pandey, Kumar and Verma2016). We thus conclude that a plausible case exists for exploring similarities with Kolmogorov turbulence in the bulk region; see Mishra & Verma (Reference Mishra and Verma2010) and Verma et al. (Reference Verma, Kumar and Pandey2017).
The anisotropy coefficient 
$G_{11}$ in three horizontal planes; for definition see (5.1); 
$G_{11}$ as well as the other coefficients 
$G_{ij}$ remain close to unity in the bulk region between 
$z = 0.1H$ and 
$z = 0.9H$, but depart significantly as the horizontal plate is approached and the shear effects dominate.
5.2. Kinetic energy spectra
In three-dimensional turbulent flows, the kinetic energy injected at large length scales cascades towards smaller scales and eventually gets dissipated at the smallest scales by viscous action. We shall not consider the connection to the Onsager conjecture that it may be related to singularities in weak solutions of the Euler equations. In the inertial range – the range of length scales far from both the injection as well as dissipation scales – the kinetic energy spectrum 
$E(k)$, whose integral over all wavenumbers 
$k$ yields the kinetic energy, follows the standard Kolmogorov scaling
\begin{equation} E(k) = K_{K} \varepsilon_u^{2/3} k^{{-}5/3}, \end{equation}
where 
$K_{K}$ is the Kolmogorov constant and 
$\varepsilon _u$ denotes the volume- and time-averaged kinetic energy dissipation rate.
For our purposes, it would be useful to study the behaviour of two-dimensional (2-D) energy spectrum in a horizontal plane. The 2-D Fourier transform of a field 
$f(x,y,z_0)$ in a horizontal plane at 
$z = z_0$ is defined as
\begin{equation} f(x,y,z_0) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \hat{F}(k_x,k_y) \exp({-{\rm i} (k_x x + k_y y)})\,{\rm d}k_x \,{\rm d}k_y, \end{equation}
where 
$\hat {F}(k_x,k_y)$ is the Fourier mode corresponding to the wavevector 
${\boldsymbol k} \equiv (k_x,k_y)$. Thus, the Fourier modes of the velocity field in midplane 
${\boldsymbol u}(x,y,z=H/2) \equiv {\boldsymbol U}(x,y)$ are denoted as 
$\hat {\boldsymbol U}({\boldsymbol k}) \equiv [\hat {U}_x({\boldsymbol k}), \hat {U}_y({\boldsymbol k}), \hat {U}_z({\boldsymbol k})]$. The kinetic energy in a horizontal plane is equal to the sum of the energies of each Fourier mode, i.e.
\begin{align} \frac{1}{2} \langle {\boldsymbol U}^2 \rangle_{A,t} &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{1}{2} |\hat{\boldsymbol U}(k_x,k_y)|^2 \,{\rm d}k_x \,{\rm d}k_y \nonumber\\ &= \int_{0}^{\infty} \left[ \int_{0}^{2{\rm \pi}} \frac{1}{2} |\hat{\boldsymbol U}(k,\phi_k)|^2 \,{\rm d}\phi_k \right] k \,{\rm d}k, \end{align}with
\begin{equation} k = \sqrt{k_x^2+k_y^2} \quad\mbox{and}\quad \phi_k = \arctan(k_y,k_x). \end{equation}
Using the horizontal isotropy of the fields, the expression in the square brackets could be readily integrated to yield 
${\rm \pi} |\hat {\boldsymbol U}(k)|^2$, where 
$|\hat {\boldsymbol U}(k)|^2/2$ is the average kinetic energy of all the Fourier modes lying in an annular region between radii 
$k$ and 
$k+{\rm d}k$. Thus, the average planar kinetic energy becomes
\begin{equation} \frac{1}{2} \langle {\boldsymbol U}^2 \rangle_{A,t} = \int_{0}^{\infty} {\rm \pi}k |\hat{\boldsymbol U}(k)|^2 \,{\rm d}k = \int_{0}^{\infty} E(k) \,{\rm d}k, \end{equation}
where 
$E(k) = {\rm \pi}k |\hat {\boldsymbol U}(k)|^2$ is the 1-D kinetic energy spectrum in a horizontal plane (Peltier et al. Reference Peltier, Wyngaard, Khanna and Brasseur1996).
We compute the energy spectrum in the midplane of the low-
$Pr$ flows for each instantaneous snapshot and then average the instantaneous spectra over all the available snapshots to obtain the mean kinetic energy spectrum. The energy spectra for flows with small-scale universality at different Reynolds numbers should collapse at sufficiently high Reynolds number if they are plotted against 
$k \eta$. Equation (5.2) in terms of the normalized wavenumber 
$k \eta$ reads then (Monin & Yaglom Reference Monin and Yaglom2007) as
\begin{equation} E(k\eta) = K_{K} (\varepsilon_u \nu^5)^{1/4} (k\eta)^{{-}5/3}. \end{equation}
We now plot the normalized energy spectra 
$E(k\eta ) (\varepsilon _u \nu ^5)^{-1/4}$ as a function of 
$k\eta$ in figure 7. The spectra for 
$Pr = 0.001$ are shown in figure 7(a) for all three Rayleigh numbers, whereas the spectra for a fixed 
$Ra = 10^6$ and 
$Pr = 0.001, 0.005$ and 0.021 are displayed in figure 7(b). The collapse is excellent beyond the wavenumber corresponding to the maximum of 
$E(k\eta )$. We show the same spectra in the normalized form 
$E(k) \varepsilon _u^{-2/3} k^{5/3}$ in the insets of figure 7, which confirm that they collapse onto each other for all the low-
$Pr$ flows in low to moderately large wavenumber range. The plateau in the dimensionless wavenumber range 
$k \in [k_0, k_1]$ increases with increasing 
$Ra$. The inset of figure 7(a) shows that the inertial range 
$[k_0, k_1]$ can be estimated to be 
$[2, 60]$, 
$[4, 150]$ and 
$[6, 300]$ for 
$Ra = 10^5, 10^6, 10^7$ and 
$Pr = 0.001$, respectively. The inertial range increases with increasing 
$Ra$ and decreasing 
$Pr$, as expected from the increased Reynolds numbers.
Normalized kinetic energy spectra 
$E(k\eta )(\varepsilon _u \nu ^5)^{-1/4}$ in the midplane as a function of the normalized wavenumber 
$k\eta$ for 
$Pr = 0.001$ (a) and 
$Ra = 10^6$ (b) collapse well in the inertial as well as viscous ranges. Here, 
$\varepsilon _u = \langle \varepsilon _u \rangle _{V,t}$. The short dotted vertical lines indicate the normalized wavenumber 
$k_T \eta$, where 
$k_T = 2{\rm \pi} /\delta _T$ is the wavenumber corresponding to the mean thermal BL thickness 
$\delta _T=H/(2 Nu)$. Insets show that the corresponding spectra normalized with 
$\varepsilon _u^{2/3} k^{-5/3}$ exhibit plateau in the inertial range, which becomes wider with increasing 
$Ra$ or decreasing 
$Pr$. The dashed horizontal line in the insets yields the Kolmogorov constant 
$K_{K} \approx 1.6$ in the planar energy spectra, consistent with the value for isotropic turbulence (Sreenivasan Reference Sreenivasan1995).
The energy spectra normalized in the same way for 
$Pr = 0.001, 0.005$ and 
$0.021$ for a fixed 
$Ra = 10^6$ are shown in the inset of figure 7(b). Again, the normalized spectra collapse quite well. A plateau can be detected for the two lower 
$Pr$ with the inertial range corresponding to 
$[4, 150]$ and 
$[4, 90]$ for 
$Pr = 0.001$ and 
$Pr = 0.005$, respectively. The plateau for all cases corresponds to 
$K_{K} \approx 1.6$, consistent with the experimental and numerical value in isotropic turbulence (Sreenivasan Reference Sreenivasan1995; Yeung & Zhou Reference Yeung and Zhou1997; Ishihara, Gotoh & Kaneda Reference Ishihara, Gotoh and Kaneda2009).
5.3. Third-order structure function
To further explore whether the velocity fluctuations in the bulk of low-
$Pr$ convection are close to Kolmogorov turbulence, we compute the third-order longitudinal structure function defined as
\begin{equation} S_3(r) = \langle (\delta_L u(r))^3 \rangle, \end{equation}
where 
$\langle \,\cdot \, \rangle$ denotes an appropriate averaging, and the longitudinal velocity increment is
\begin{equation} \delta_L u(r) = [{\boldsymbol u}({\boldsymbol x} + {\boldsymbol r}) - {\boldsymbol u}({\boldsymbol x})] \cdot \hat{\boldsymbol r}, \end{equation}
with 
$\hat {\boldsymbol r}={\boldsymbol r}/r$. Kolmogorov (Reference Kolmogorov1941a) showed that 
$S_3(r)$ in high-
$Re$ homogeneous and isotropic turbulent flow is a universal function of the separation 
$r$ and varies as
\begin{equation} S_3(r) ={-}\frac{4}{5} \varepsilon_u r \end{equation}
in the inertial range; here, and for the remainder of the work, 
$\varepsilon _u:=\langle \varepsilon _u \rangle _{A,t}(z=1/2)$. The longitudinal structure functions in the 
$x$- and 
$y$-directions in midplane for our highest-
$Re$ flow show that they are nearly the same for small and moderate increments, so we average in the two directions. We show in figure 8 the averaged third-order structure function in the normalized form 
$-S_3(r)/(\varepsilon _u r)$ as a function of 
$r/\eta$. The figure shows that the compensated structure function tends to scale with the analytical form of 
$r^2$ in the beginning of the viscous range at 
$r/\eta \sim 1$. It also shows that the normalized structure function exhibits a plateau for an intermediate range of length scales, which implies that 
$S_3(r)$ indeed approximately varies as 
$r$ in the inertial range but the numerical value is slightly larger than 
$4/5$ (Kolmogorov Reference Kolmogorov1941a). One possible reason for this departure is the remnant buoyancy contribution (Yakhot Reference Yakhot1992).
Third-order longitudinal structure function (averaged over the 
$x$- and 
$y$-directions) in the midplane for 
$Pr = 0.001, Ra = 10^7$ varies as 
$r^2$ in the viscous range, whereas it remains nearly a constant in the inertial range. Dashed horizontal line indicates the constant 4/5 appearing in (5.10). Here, 
$\varepsilon _u = \langle \varepsilon _u \rangle _{A,t}(z=H/2)$.
5.4. Dissipative anomaly
The zeroth law of turbulence (or ‘dissipative anomaly’) states that the mean kinetic energy dissipation rate when scaled by a large-scale velocity such as the r.m.s. and the large length scale becomes a constant for sufficiently high Reynolds numbers (Eyink Reference Eyink1994; Frisch Reference Frisch1995; David & Galtier Reference David and Galtier2021). Figure 9 displays this rescaled energy dissipation rate
\begin{equation} \beta=\frac{\varepsilon_u {\mathcal{L}}}{u_{rms}^3}, \end{equation}
where 
${\mathcal {L}}$ is either the height 
$H$ of the convection layer or the integral scale 
$\ell$ calculated from the energy spectrum (5.6) in § 5.2 by
\begin{equation} \ell=2{\rm \pi}\,\dfrac{\int_0^{\infty} k^{{-}1} E(k) \,{\rm d}k}{\int_0^{\infty} E(k) \,{\rm d}k}. \end{equation}
Note that 
$u_{rms}$ is also taken with respect to the midplane. The figure summarizes both versions of 
$\beta$ vs the Taylor microscale Reynolds number
\begin{equation} R_{\lambda}= \left( \frac{25 Ra}{9\varepsilon_u^2 Pr} \right)^{1/4} \, u_{rms}^2 \end{equation}in the bulk region of the flow. The data for different Rayleigh and Prandtl numbers collapse nicely on a curve that saturates at an approximate value of 0.2 (see the inset of the figure). While this is smaller than 0.45 found by Sreenivasan (Reference Sreenivasan1998) for isotropic turbulence, the result implies that the strongly disparate viscous and thermal BL widths (and thus the plume stem widths) do not matter for the driving of the turbulence cascade in the bulk of the flow. However, it also highlights the intrinsic difference in dissipation between convection and isotropic turbulence.
Test of the dissipative anomaly in the bulk of the convection layer using the rescaled mean kinetic energy dissipation rate 
$\beta$, as given by (5.11), vs the Taylor microscale Reynolds number. Asterisks use the integral scale 
${\mathcal {L}}=\ell$ in (5.11), and open circles use 
${\mathcal {L}}=H$. The inset expands the scale for the case 
${\mathcal {L}} = H$. Symbols from left to right in each colour correspond to decreasing 
$Pr$. The Taylor microscale Reynolds number is given by (5.13).
6. Characteristic lengths and times of turbulent superstructures
The characteristic length scale of the dominant energy-containing structures is of the order of the size of the system when 
$\varGamma \approx 1$, e.g. a large-scale circulation covering the entire domain (Schumacher et al. Reference Schumacher, Bandaru, Pandey and Scheel2016; Zürner et al. Reference Zürner, Schindler, Vogt, Eckert and Schumacher2019; Zwirner et al. Reference Zwirner2020). However, for 
$\varGamma \gg 1$, mean circulation rolls with diameters larger than 
$H$ are observed (Emran & Schumacher Reference Emran and Schumacher2015); the resulting large-scale patterns of these rolls are termed turbulent superstructures of convection (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Fonda et al. Reference Fonda, Pandey, Schumacher and Sreenivasan2019; Green et al. Reference Green, Vlaykov, Mellado and Wilczek2020) as we stated already in § 1. Although the superstructures extend all the way from the bottom to the top plate (Pandey et al. Reference Pandey, Scheel and Schumacher2018), they are conspicuous when the vertical velocity 
$u_z$, temperature 
$T$ or the vertical heat flux 
$u_zT$ field is visualized in the midplane (Fonda et al. Reference Fonda, Pandey, Schumacher and Sreenivasan2019; Green et al. Reference Green, Vlaykov, Mellado and Wilczek2020). For instance, figure 1 shows turbulent superstructures for low-
$Pr$ convection in the form of hot upflows and cold downflows.
The characteristic length 
$\lambda$ of the superstructures is the typical distance between two consecutive upwelling or downwelling regions and can be estimated using 1-D spectra of the thermal variance, kinetic energy and convective heat flux in a horizontal plane (Hartlep et al. Reference Hartlep, Tilgner and Busse2003; Pandey et al. Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Green et al. Reference Green, Vlaykov, Mellado and Wilczek2020; Krug et al. Reference Krug, Lohse and Stevens2020). It can also be estimated by computing the two-point auto-correlation function of 
$u_z$ or 
$T$ in a horizontal plane and identifying the location of the first minimum, corresponding to 
$\lambda /2$ (Pandey et al. Reference Pandey, Scheel and Schumacher2018).
We compute the 1-D power spectra of the vertical velocity, temperature and convective heat flux, all averaged with respect to the azimuthal angle. They are given by
$$\begin{gather} S_U(k) = \frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}} |\hat{U}_z(k,\phi_k)|^2 \,{\rm d} \phi_k, \end{gather}$$
$$\begin{gather}S_\varTheta(k) = \frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}} |\hat{\varTheta}(k,\phi_k)|^2 \,{\rm d} \phi_k, \end{gather}$$
$$\begin{gather}S_{U\varTheta}(k) = \frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}} \mathrm{Re}[\hat{U}_z^*(k,\phi_k) \hat{\varTheta}(k,\phi_k)] \,{\rm d} \phi_k, \end{gather}$$
where 
$\hat {U}_z(k,\phi _k)$ and 
$\hat {\varTheta }(k,\phi _k)$ are, respectively, the 2-D Fourier transforms of 
$u_z(x,y,H/2)$ and 
$\theta (x,y,H/2)$. Note that the spectrum 
$S_U(k)$ differs from the energy spectrum 
$E(k)$ defined in § 5.2. We find that the three spectra exhibit a peak at nearly the same wavenumber 
$k^\omega _{max}$ corresponding to the maximum of 
$S_\omega (k)$ (see figure 13), yielding the characteristic spatial scale 
$\lambda _\omega = 2{\rm \pi} /k^\omega _{max}$ of the superstructures, with 
$\omega = \{ U, \varTheta, U\varTheta \}$.
Pandey et al. (Reference Pandey, Scheel and Schumacher2018, Reference Pandey, Schumacher and Sreenivasan2021) found that the characteristic scales 
$\lambda _U$ and 
$\lambda _\varTheta$ do not always agree with each other, the contrast being usually larger at moderate and high Prandtl numbers. Here, we extract the spatial scale of superstructures from the power spectrum of the convective heat flux (see also Hartlep et al. Reference Hartlep, Tilgner and Busse2003; Krug et al. Reference Krug, Lohse and Stevens2020). It was observed in Pandey et al. (Reference Pandey, Scheel and Schumacher2018) that 
$\lambda$ is a function of 
$Pr$ and the maximum of 
$\lambda (Pr)$ was found at 
$Pr \approx 7$ for 
$Ra = 10^5$. In the present simulations (in which the 
$Pr$ range has been extended down to 
$0.001$ and 
$Ra$ by two orders of magnitude), the corresponding values of 
$\lambda$ as a function of 
$Pr$ are plotted in figure 10(a) for all simulations of table 1. We find that 
$\lambda$ decreases with decreasing 
$Pr$ but they seem to level off at 
$\lambda \simeq 3H$ for the lowest Prandtl numbers of 0.005 and 0.001; see also figure 11 where the magnitude of 
$\lambda$ is displayed in the streamline plot. It remains to be studied as to how this scale arises starting from the onset of RBC where 
$\lambda = 2.02 H$ independently of 
$Pr$ (Chandrasekhar Reference Chandrasekhar1981), as indicated in figure 10(a) as a dashed horizontal line. We also refer to Appendix C for further details.
Characteristic spatial (a) and temporal (b) scales of the turbulent superstructures as a function of 
$Pr$. Both temporal and spatial scales decrease with decreasing Prandtl number, but seem to level off for 
$Pr \leq 0.005$.
Plots of streamlines for snapshots at (a) 
$Ra=10^5$, (b) 
$10^6$ and (c) 
$10^7$, all at 
$Pr=0.001$. The view is from the top. Only one half of the cross-section is shown. The corresponding values of 
$\lambda$ are indicated by a horizontal white bar in each panel.
The characteristic temporal scale (
$\tau$) of the superstructures is long compared with the time scale of the turbulent fluctuations. Thus, time averaging the velocity or temperature fields over a time interval 
$\tau$ yields coarse-grained fields, which are nearly devoid of the small-scale fluctuations. The coarse-grained large-scale structures evolve on the time scales which are of the order of 
$\tau$ (Pandey et al. Reference Pandey, Scheel and Schumacher2018; Fonda et al. Reference Fonda, Pandey, Schumacher and Sreenivasan2019), related to the time for a fluid parcel to complete a circulation. It is computed as
\begin{equation} \tau = 3 \frac{{\rm \pi}(\lambda/4 + H/2)}{u_{rms}}, \end{equation}
where the quantity in the numerator is the circumference of superstructure rolls with elliptical cross-section (Pandey et al. Reference Pandey, Scheel and Schumacher2018). The factor of three in the above expression arises from the fact that the circulation time in an extended convection flow is not fixed but exhibits a broad distribution with stretched-exponential tails in the Lagrangian frame of reference along massless tracer trajectories (Schneide et al. Reference Schneide, Pandey, Padberg-Gehle and Schumacher2018; Vieweg et al. Reference Vieweg, Schneide, Padberg-Gehle and Schumacher2021b). The scale computed using expression (6.4) is plotted in figure 10(b) as a function of 
$Pr$. We find that 
$\tau$ increases with 
$Pr$, consistent with the fact that the Reynolds number, and thus the characteristic velocity 
$u_{rms}$, decreases with increasing 
$Pr$, requiring a longer time for fluid parcels to complete a circulation. Note that 
$\lambda$ in (6.4) also increases with increasing 
$Pr$, but the increase is not as significant as the decrease in 
$u_{rms}$.
7. Conclusions and outlook
Our focus here has been the Rayleigh–Bénard convection in a horizontally extended layer for molecular Prandtl numbers as small as 
$Pr=10^{-3}$, which go beyond those accessible in controlled laboratory experiments and approach astrophysical conditions. We extended the parameter space of previous works by Pandey et al. (Reference Pandey, Scheel and Schumacher2018) and Fonda et al. (Reference Fonda, Pandey, Schumacher and Sreenivasan2019) by DNS, both towards lower 
$Pr$ and higher 
$Ra$, and thus determined more conclusively various parameter dependencies such as global heat and momentum transports, temperature fluctuations, as well as the kinetic and thermal dissipation rates. Among others, these results provide a test for existing predictions by the theory of Grossmann & Lohse (Reference Grossmann and Lohse2001) and Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013). Comparisons show that the predictions for the global heat and momentum transports as a function of the Prandtl number at fixed Rayleigh number are in fair agreement.
We also found that the Nusselt number decreases as 
$Nu \sim Pr^{0.19}$, whereas the Reynolds number increases approximately as 
$Re \sim Pr^{-0.6}$ when the Prandtl number decreases from 
$Pr = 0.7$ to 0.001, as detailed in table 2. The dimensionless mean thermal and kinetic energy dissipation rates also decrease with increasing 
$Pr$ and their scaling behaviours are consistent with that of the global heat transport. We studied the depth dependence of these quantities, and found that, due to a high diffusive temperature field in convection at very low 
$Pr$, the bulk fluid is not mixed well and a significant vertical temperature gradient occurs in the bulk region, even for the highest accessible Rayleigh number.
The highly inertial fluid turbulence in the bulk of low-
$Pr$ convection layer was studied by examining the kinetic energy spectra, the 4/5th law, local isotropy and a test of the dissipative anomaly. The results suggest that the fluid turbulence in the bulk for the lowest Prandtl number is close to the classical Kolmogorov turbulence. This implies that the temperature field behaves as a passive and highly diffusive scalar stirred by a highly turbulent flow; the impact of thermal plumes, which are the unstable fragments of the thick thermal BLs, can be considered an efficient large-scale forcing for turbulence, dominantly in the lower wavenumber part of the inertial range (
$k_T$ in figure 7). Some differences do remain. In particular, the asymptotic value of the rescaled mean kinetic energy dissipation rate falls below that in isotropic turbulence (Sreenivasan Reference Sreenivasan1998), which suggests that BLs do matter, despite the nearness to isotropy in the central region.
Our DNS results are fully resolved and can thus have implications for the modelling of small-scale turbulence in coarser-grid simulation studies of mesoscale convection, particularly for the development of subgrid-scale models that go beyond the mixing length theory of Prandtl (Reference Prandtl1925). This class of algebraic turbulence models is still a workhorse in astrophysical simulations; see for example discussions of their limitations and extensions in Miesch (Reference Miesch2005) and Kupka & Muthsam (Reference Kupka and Muthsam2017), or recent extensions by Brandenburg (Reference Brandenburg2016). Finally, we stress that the convection considered here does not incorporate complexities such as rotation, magnetic field, varying molecular transport coefficients or curvature, which are present in geophysical and astrophysical settings. They would have to be included to yield realistic models in these specific instances. The present focus has been the exploration of the effects of low Prandtl numbers, which is an important facet of these flows. A more detailed study of these points in connection with small-scale intermittency in low Prandtl number RBC flows is underway, and will be reported elsewhere.
Acknowledgements
The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (https://www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (https://www.lrz.de).
Funding
A.P. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) within the Priority Programme ‘Turbulent Superstructures’ under Grant No. DFG-SPP 1881. This work was also supported by Grant No. SCHU 1410/30-1 of DFG and NYUAD Institute Grant G1502 ‘NYUAD Center for Space Science.’ D.K. is partly supported by Grant No. KR 4445/2-1 of DFG.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix A. Comparison of the spectral element and finite difference solvers
As mentioned in § 2, we use two different solvers, one based on the spectral element method and one based on the finite difference method (FDM), to perform our simulations due to the demanding requirements at very low Prandtl numbers. Therefore, it is important to ensure that both the solvers yield the same results. To check this, we performed comparison simulations, for 
$Pr = 0.005$ and 
$Ra = 10^5$ and for 
$Pr=0.7$ and 
$Ra=10^7$. In figure 12, we plot the vertical profiles of the convective heat flux, the thermal and the kinetic energy dissipation rates for the former comparison. The profiles from both the spectral element and finite difference solvers agree very well with each other. The globally averaged quantities such as the Nusselt and the Reynolds numbers also agree excellently.
Vertical profiles of (a) the convective heat flux, (b) the thermal dissipation rate and (c) the kinetic energy dissipation rate computed using the spectral element (solid red curves) and the finite difference (dashed green curves) solvers for 
$Pr = 0.005, Ra = 10^5$ agree very well each other.
The global heat flux can also be estimated using the global dissipation rates from the exact relations (Howard Reference Howard1972), and their concurrence is also an indicator of the sufficiency of the spatial and temporal resolutions (Pandey et al. Reference Pandey, Schumacher and Sreenivasan2021). The exact relations yield the Nusselt numbers as
$$\begin{gather} Nu_{\varepsilon_u} = 1 + \sqrt{RaPr} \langle \varepsilon_u \rangle_{V,t} , \end{gather}$$
$$\begin{gather}Nu_{\varepsilon_T} = \sqrt{RaPr} \langle \varepsilon_T \rangle_{V,t}. \end{gather}$$
Further, the heat flux 
$Nu(z)$ in each horizontal plane, see (4.3) in the main text, remains a constant across the convection layer in the statistically steady state and matches with 
$Nu$. We thus compute the averaged heat flux at the top and bottom plates as
\begin{equation} Nu_{\partial_z T} ={-}\left \langle \left( \frac{\partial T}{\partial z}\right)_{z=0,H} \right\rangle_{A,t} \end{equation}
and list it along with 
$Nu, Nu_{\varepsilon _u}, Nu_{\varepsilon _T}$ in table 4. The results in the table show that the agreement between all differently obtained Nusselt numbers is excellent for all the simulations.
The turbulent heat flux, which is computed in four different ways, agrees very well for all reported simulations. The error bars indicate the standard deviation. Here, 
$t_{sim}$ is the total simulation time in the statistically steady state. The time advancement by one free-fall time unit for the biggest simulation took 30 million core hours on 144 000 processor cores on the cluster SuperMUC-NG at Leibniz Rechenzentrum Garching.
A few words on the determination of the kinetic energy dissipation rate from the results of the FDM solver. It has been shown in Viré & Knaepen (Reference Viré and Knaepen2009) that, for the rate of strain tensor 
${\mathsf{S}}_{ij}=(\partial u_i/\partial x_j + \partial u_j/\partial x_i)/2$, the identity 
${\mathsf{S}}_{ij} {\mathsf{S}}_{ij} = -u_{i} \partial {\mathsf{S}}_{ij}/\partial x_j + \partial (u_{i} {\mathsf{S}}_{ij})/\partial x_j$ is not fulfilled in discrete form, which is particularly relevant to finite difference and finite volume methods. The direct computation of 
${\mathsf{S}}_{ij} {\mathsf{S}}_{ij}$ tends to yield underpredicted values, especially if the discretization errors are large. For example, in case of coarse grids, more pertinent to LES however, the difference between direct computation of 
${\mathsf{S}}_{ij} {\mathsf{S}}_{ij}$ and summation by parts can yield a factor of 2–2.5 in the buffer and logarithmic layer regions (Viré et al. Reference Viré, Krasnov, Boeck and Knaepen2011). Albeit this difference scales as 
${\sim}O(h^2)$ (with 
$h$ being the mesh step-size) for the 
$2{\rm nd}$-order approximations, it cannot be completely neglected, even in case of finer resolutions used in DNS. Here, we cannot directly apply the summation-by-parts approach, since it involves the so-called flux variables, prescribed at the cell interface (Viré et al. Reference Viré, Krasnov, Boeck and Knaepen2011). These flux variables are used in the algorithm to secure divergence-free condition and conservative form of the nonlinear terms, but not stored. We have applied a different approach for the FDM results – the 
${\mathsf{S}}_{ij}$ tensor is computed by using 
$6{\rm th}$-order accurate stencils for the velocity gradients to reduce the effect of discretization errors.
Appendix B. Grid sensitivity for our biggest simulation
The Kolmogorov length scale in the flow for 
$Pr = 0.001$ and 
$Ra = 10^7$ becomes very small and, consequently, computational resources required for the numerical investigation of this flow become exorbitant. Therefore, to determine the optimum number of nodes needed to resolve the flow adequately, we performed this simulation on three different grids with 
$15\,360 \times 15\,360 \times 1024$, 
$20\,480 \times 20\,480 \times 1280$ and 
$22\,400 \times 22\,400 \times 1400$ mesh cells, which we denote in the following as mesh-1, mesh-2 and mesh-3, respectively. As summarized in Scheel et al. (Reference Scheel, Emran and Schumacher2013), we compare the horizontally as well as the globally averaged quantities from these simulations to test the effects of grid resolution. The kinetic energy dissipation rate field 
$\varepsilon _u ({\boldsymbol x},t)$, which involves the computation of all the nine terms of the velocity gradient tensor, is very sensitive to the mesh size in a low-
$Pr$ convection. We computed the horizontally averaged kinetic energy dissipation rate 
$\langle \varepsilon _u\rangle _{A,t}(z)$ and investigated the region near the midplane where the computational grid is coarsest. We observed underresolved data for mesh-1, whereas the variation of 
$\langle \varepsilon _u\rangle _{A,t}(z)$ is smooth for simulations with mesh-2 and mesh-3. Furthermore, 
$\langle \varepsilon _u\rangle _{A,t}(z)$ from mesh-2 and mesh-3 agree very well, in particular they both yield the same relative difference between the results of 
$2{\rm nd}$- and 
$6{\rm th}$-order stencils applied for the direct computation of the velocity gradients. This analysis clearly indicates that mesh-2 is able to properly capture the velocity derivatives for these extreme parameters. Thus, we carried out our simulation for 
$Pr = 0.001$ and 
$Ra = 10^7$ with mesh-2.
Appendix C. Estimation of the characteristic lengths and times of superstructures
We show the power spectra 
$S_U(k), S_\varTheta (k)$ and 
$S_{U\varTheta }(k)$ for 
$Pr = 0.001$ at 
$Ra = 10^5$ and 
$Ra = 10^7$ in figure 13 and find that the spectral distributions for the velocity and temperature fields are different. Figure 13 shows that the power first increases with decreasing length scales and attains a maximum before declining sharply with further decrease in the scale size. We find that the decay of the thermal variance spectrum beyond 
$k_{max}$ is rapid compared with that of the squared vertical velocity component. This is because the velocity field in very-low-
$Pr$ convection is vigorously turbulent and possesses larger fine-scale contributions compared with the predominantly large-scale nature of the temperature field (Schumacher et al. Reference Schumacher, Götzfried and Scheel2015). We find, however, that the three spectra exhibit a peak at nearly the same wavenumber 
$k^\omega _{max}$ corresponding to the maximum of 
$S_\omega (k)$, yielding the characteristic spatial scale 
$\lambda _\omega = 2{\rm \pi} /k^\omega _{max}$ of the superstructures, with 
$\omega = \{ U, \varTheta, U\varTheta \}$.
Power spectra 
$S_U$, 
$S_{\varTheta }$ and 
$S_{U\varTheta }$ as defined in (6.1)–(6.3) in the midplane for 
$Pr = 0.001$ at (a) 
$Ra = 10^5$ and (b) 
$Ra = 10^7$. They exhibit peaks at nearly the same wavenumber corresponding to the characteristic spatial scale (or wavelength) of the turbulent superstructures of convection.
Furthermore, the spatial scale does not remain fixed but fluctuates during the evolution of the flow. This can be seen in figure 14, where we plot 
$k_{max}^{U\varTheta }(t)$, extracted from each instantaneous snapshot of our simulations, as a function of time. The figure shows that 
$k^{U\varTheta }_{max}$ is independent of time for the entire duration of simulations for 
$Pr = 0.001$. The same comment also holds for the 
$Pr = 0.005$ and 
$Pr = 0.021$ simulations (not shown in the figure). However, 
$k^{U\varTheta }_{max}$ for the simulations at 
$Pr = 0.7$ and 
$Pr = 7$, though remaining fixed for most of the time, shows occasional excursions. Thus, we determine the characteristic spatial scale of superstructures by time averaging 
$k^{U\varTheta }_{max}(t)$ and finding 
$\lambda = 2{\rm \pi} / \langle k^{U\varTheta }_{max}(t) \rangle _t$.
Temporal evolution of the peak wavenumber 
$k_{max}$ for 
$Ra = 10^5$ (a,c,e) and 
$Ra = 10^6$ (b,d, f). The location of maximum does not change for 
$Pr = 0.001$ (a,b), whereas it occasionally shifts towards the higher wavenumber for 
$Pr = 0.7$ (c,d) and also for 
$Pr = 7$ (e, f).
Open access
