On the inverse cascade and flow speed scaling behavior in rapidly rotating Rayleigh-B\'enard convection

Rotating Rayleigh-B\'enard convection is investigated numerically with the use of an asymptotic model that captures the rapidly rotating, small Ekman number limit, $Ek \rightarrow 0$. The Prandtl number ($Pr$) and the asymptotically scaled Rayleigh number ($\widetilde{Ra} = Ra Ek^{4/3}$, where $Ra$ is the typical Rayleigh number) are varied systematically. For sufficiently vigorous convection, an inverse kinetic energy cascade leads to the formation of a depth-invariant large-scale vortex (LSV). With respect to the kinetic energy, we find a transition from convection dominated states to LSV dominated states at an asymptotically reduced (small-scale) Reynolds number of $\widetilde{Re} \approx 6$ for all investigated values of $Pr$. The ratio of the depth-averaged kinetic energy to the kinetic energy of the convection reaches a maximum at $\widetilde{Re} \approx 24$, then decreases as $\widetilde{Ra}$ is increased. This decrease in the relative kinetic energy of the LSV is associated with a decrease in the convective correlations with increasing Rayleigh number. The scaling behavior of the convective flow speeds is studied; although a linear scaling of the form $\widetilde{Re} \sim \widetilde{Ra}/Pr$ is observed over a limited range in Rayleigh number and Prandtl number, a clear departure from this scaling is observed at the highest accessible values of $\widetilde{Ra}$. Calculation of the forces present in the governing equations shows that the ratio of the viscous force to the buoyancy force is an increasing function of $\widetilde{Ra}$, that approaches unity over the investigated range of parameters.


Introduction
Convection is common in the fluid regions of planet and stars. In particular, convection is the primary energy source for the generation of large-scale planetary and stellar magnetic fields (Jones 2011;Gastine et al. 2014;Aurnou et al. 2015), and it is thought to be a source of energy for the observed zonal flows in the atmospheres of the giant planets (e.g. Heimpel et al. 2016). The flows in many of these natural systems are considered turbulent and strongly influenced by rotation; previous studies have shown that the combination of these physical ingredients can lead to an inverse kinetic energy cascade (Smith & Waleffe 1999). The inverse cascade transfers kinetic energy from small-scale convection up to domain-scale flows, and, in a planar geometry of square cross section, results in the formation of depth-invariant, large-scale vortices (LSVs) (Julien et al. 2012b;Rubio et al. 2014;Guervilly et al. 2014;Favier et al. 2014;Stellmach et al. 2014). Such vortices tend to be characterized by flow speeds that are significantly larger than the underlying convection, and can have an influence on heat transport and magnetic field generation (Guervilly et al. 2014(Guervilly et al. , 2015. However, the convective vigor required to excite LSVs, how fluid properties (via the thermal Prandtl number) influence their formation, and the scaling of their saturated amplitude with buoyancy forcing still remain poorly understood.
In planetary and astrophysical fluid systems, rapid rotation is thought to play an essential role in shaping the dynamics of convection. The importance of rotation on the dynamics of such systems is quantified by the Ekman and Rossby numbers, respectively defined as where ν is the kinematic viscosity, Ω is the rate of rotation, H is the spatial scale of the system (i.e. the depth of the fluid region) and U is a characteristic flow speed. Ek and Ro quantify, respectively, the ratio of viscous forces to the Coriolis force, and the ratio of inertia to the Coriolis force. Systems in which (Ek, Ro) 1 are said to be rapidly rotating and rotationally-constrained. For the Earth's outer core, for instance, estimates suggest that Ek 10 −15 and Ro 10 −5 (Finlay & Amit 2011;de Wijs et al. 1998;Rutter et al. 2002). It is currently impossible to use such extreme values of the governing parameters with direct numerical simulation (DNS). Quasi-geostrophic (QG) models have helped to overcome this deficiency, and have been critical for elucidating several convective phenomena that are thought to be of significant interest for planets, including identification of the primary flow regimes of rotating convection, heat transport behavior, and the convection-driven inverse cascade (Julien et al. 2012b;Rubio et al. 2014). QG models accurately capture the leading order dynamics in systems characterized by small values of Ek and Ro.
Previous work has shown that the structure of the LSVs is dependent on the relative importance of rotation: for sufficiently small values of Ek and Ro the vortices are dipolar in structure with zero net (spatially-averaged) vorticity (Julien et al. 2012b;Rubio et al. 2014;Stellmach et al. 2014); whereas for larger values of Ro the vortices tend to be cyclonic with a net vorticity that is parallel to the rotation axis (Chan & Mayr 2013;Guervilly et al. 2014;Favier et al. 2014). QG models find only dipolar LSVs since they capture asymptotically-small values of Ek and Ro only, in which there is no preferred sign for the vorticity. DNS studies have found dipolar LSVs when Ek ≈ 10 −7 (Stellmach et al. 2014), and cyclonic LSVs when Ek 10 −6 (Chan & Mayr 2013;Guervilly et al. 2014;Favier et al. 2014). This distinction in LSV structure may have consequences for the influence of LSVs on heat transport and flow speeds. Moreover, for Ro = O(1), the presence of LSVs appears to be dependent upon the initial condition used in the simulations (Favier et al. 2019). In contrast, for the rapidly rotating regime studied here, we find that LSV formation is not dependent upon initial conditions. Natural systems are characterized by a broad range of fluid properties and, as a result, the Prandtl number P r = ν/κ (where κ is the thermal diffusivity) can take on a wide range of values, ranging from P r = O(10 −6 ) in stellar interiors (Ossendrijver 2003) to P r = O(10 −2 ) for the liquid metals characteristic of planetary interiors (Pozzo et al. 2013). More generally, the density heterogeneities that lead to buoyancy-driven convection can also result from compositional differences, as is expected to be the case within terrestrial planetary interiors, for instance; under such circumstances the thermal Prandtl number in the governing equations is replaced by the compositional Schmidt number Sc = ν/D where D is a chemical diffusivity. For the Earth's outer core, studies suggest Sc = O(100) for representative chemical species (e.g. Pozzo et al. 2013). This wide range of diffusivities that characterize geophysical and astrophysical fluids motivates the need for additional investigations that explore the influence of the Prandtl number on the dynamics, since all previous numerical calculations investigating the inverse cascade have focussed on fluids with P r = 1. QG simulations (Julien et al. 2012b) have characterized the flow regimes that occur when P r 1. In general, it is found that low P r fluids reach turbulent regimes for a lower value of the thermal forcing (as measured by the Rayleigh number Ra) than higher P r fluids. More turbulent flows can be characterized by an increase in the Reynolds number which is a fundamental output parameter for convection studies. Many of the results found in QG studies have been confirmed by DNS calculations (Stellmach et al. 2014).
In both approaches, the formation of LSVs has been tied to the geostrophic turbulence regime, which, prior to the present work, has not been observed for P r > 3. Laboratory experiments are an important tool for exploring the dynamics of rapidly rotating convection (e.g. Vorobieff & Ecke 2002;King et al. 2009;Kunnen et al. 2010;Stellmach et al. 2014;Ecke & Niemela 2014) and, much like natural systems, can access a wide range of P r values. Liquid metals (P r ≈ 0.025) (Aurnou & Olson 2001;Aubert et al. 2001;King & Aurnou 2013;Zimmerman et al. 2014;Adams et al. 2015;Aurnou et al. 2018;Vogt et al. 2018), water (P r ≈ 7) (Sumita & Olson 2002;Vorobieff & Ecke 2002;King et al. 2009;Cheng et al. 2015), and gases (P r ≈ 1) (Niemela et al. 2010;Ecke & Niemela 2014), have all been utilized. As for numerical calculations, fluids with smaller values of P r reach higher values of Re than higher P r fluids for equivalent Ra. From this perspective, such fluids are of significant interest for exploring the geostrophic turbulence regime of rapidly rotating convection (e.g. Julien et al. 2012a). However, the use of such low P r fluids is not always practical and can reduce or eliminate flow visualization opportunities. Furthermore, whereas small P r fluids can lead to more turbulent flows, lower Ekman numbers must be used to provide a sufficiently large parameter regime over which the fluid remains rotationally-constrained (e.g. King & Aurnou 2013). It would therefore be of use to identify the general dynamical requirements for observing inversecascade-generated LSVs for a variety of fluid properties.
One of the basic goals in convection studies is to determine the functional dependence of Re on the input parameters, namely, determination of the functional form Re = f (Ra, P r). Power-law scalings of the form Re = c 1 Ra c2 P r c3 are often sought, where each c i is a constant. A well-known example is the so-called 'free-fall' scaling of the form Re ∼ (Ra/P r) 1/2 observed in non-rotating convection (e.g. Ahlers et al. 2009;Orvedahl et al. 2018). The free-fall scaling arises from a balance between nonlinear advection and the buoyancy force in the momentum equation, and represents a 'diffusion-free' scaling in the sense that the flow speeds are independent of both the thermal and viscous diffusion coefficients. Motivated by this free-fall form of the scaling, and the assumption that natural systems are expected to be highly turbulent in the sense that Re 1, rotating convection studies have also sought to find diffusion-free scalings for the flow speeds. For instance, the recent work of Guervilly et al. (2019) observed Re ∼ RaEk/P r in spherical convection simulations, which is also a diffusion-free scaling for the flow speeds.
In the present work we show that this scaling is equivalent to Re ∼ Ra/P r, where Re = Ek 1/3 Re and Ra = Ek 4/3 Ra are, respectively, a Reynolds number and a Rayleigh based on the small convective scale . The Re ∼ RaEk/P r scaling appears to be present for our larger P r cases over a finite range in Ra; as Ra is increased a significant departure in this scaling is observed.
In the present work we investigate the properties of the inverse cascade for varying Ra and P r in the rapidly rotating asymptotic (QG) limit for thermal convection in a plane layer geometry (rotating Rayleigh-Bénard convection). This choice allows for comparison with previous results from QG (Sprague et al. 2006;Julien et al. 2012b;Rubio et al. 2014) and DNS (Stellmach et al. 2014;Guervilly et al. 2014) plane-layer calculations. The paper is organised as follows: in section 2 we describe the governing equations and diagnostic quantities; in section 3 the results of the simulations are presented and analysed; and a discussion is given in section 4.

Governing equations
We consider rotating Rayleigh-Bénard convection in a plane layer Cartesian geometry of depth H, with constant gravity vector g = −g z pointing vertically downward, perpendicular to the planar boundaries. The fluid is Boussinesq with thermal expansion coefficient α. The top boundary is held at constant temperature T 1 and the bottom boundary is held at constant temperature T 2 such that ∆T = T 2 − T 1 > 0. The system rotates about the vertical with rotation vector Ω = Ω z. In the limit of strong rotational constraint (i.e. small Rossby and Ekman numbers), the governing equations, can be reduced to the following set of equations Sprague et al. 2006) where ζ is the vertical component of the vorticity, ψ is the dynamic pressure and also the geostrophic streamfunction, and w is the vertical component of the velocity. The vertical component of vorticity and the streamfunction are related via ζ = ∇ 2 ⊥ ψ, where ∇ 2 ⊥ = ∂ x 2 + ∂ y 2 . The non-dimensional temperature θ is decomposed into mean and fluctuating components Θ and ϑ, respectively, such that θ = Θ + Ek 1/3 ϑ. Here the mean is defined as an average over the small spatial (x, y, z) and the fast temporal (t) scales. The Jacobian operator J[ψ, f ] = ∂ x ψ∂ y f − ∂ y ψ∂ x f = u ⊥ · ∇ ⊥ f describes advection of the generic scalar field f by the horizontal velocity field u ⊥ = (u, v, 0) = (−∂ y ψ, ∂ x ψ, 0). The reduced Rayleigh number Ra is defined by where the standard Rayleigh number is (2.6) Equations (2.1)-(2.4) have been nondimensionalized by the small-scale viscous diffusion time 2 /ν, where the small horizontal convective length scale is = HEk 1/3 . The derivation of the reduced system relies on the assumption that the Coriolis force and pressure gradient force are balanced at leading order, i.e. geostrophic with z×u = −∇ ⊥ ψ. This force balance implies the Taylor-Proudman constaint is satisfied on small vertical scales such that ∂ z (u, ψ) = 0 (e.g. Stewartson & Cheng 1979). Therefore, along the vertical direction all fluid variables vary on an O(H) dimensional scale associated with the coordinate Z = Ek 1/3 z. As a result, fast inertial waves with dimensional frequency O(Ω) are filtered from the above equations, allowing for substantial computational savings. However, slow, geostrophically-balanced inertial waves are retained (Julien et al. 2012b).
The mean temperature Θ evolves on the slow time-scale τ = Ek 2/3 t associated with the vertical diffusion time H 2 /ν. However, Julien et al. (1998);Plumley et al. (2018) found that spatial averaging over a sufficient number of convective elements on the small scales is sufficiently accurate to (i) omit fast-time averaging and (ii) assume a statistically stationary state where the slow evolution term ∂ τ Θ that would appear in (2.4) is omitted.
Finally, we note that three-dimensional incompressibility is invoked through the solenoidal condition for the ageostrophic, sub-dominant horizontal velocity, where u ag ⊥ = O(Ek 1/3 u ⊥ ). The equations are solved using impenetrable, stress-free mechanical boundary conditions, and constant temperature boundary conditions. However, it should be noted that the specific form of the thermal boundary conditions are unimportant in the limit of rapid rotation , and the present model can be generalized to no-slip mechanical boundary conditions . Each variable is represented with a spectral expansion consisting of Chebyshev polynomials in the vertical (Z) dimension, and Fourier series in the horizontal (x, y) dimensions. The resulting set of equations are truncated and solved numerically with a pseudo-spectral algorithm that uses a 3rd-order implicit/explicit Runge-Kutta time-stepping scheme (Spalart et al. 1991). The code has been benchmarked successfully and used in many previous investigations Maffei et al. 2019;Yan et al. 2019).
Spatial and temporal resolutions are given in table 1. The horizontal dimensions of the domain are periodic and scaled by the critical horizontal wavelength λ c = 2π/k c ≈ 4.8154, measured in small-scale units . Most of the simulations use horizontal dimensions of 10λ c × 10λ c , though some additional simulations with different domain size were also carried out to quantify the influence of the geometry. We find that a domain size of 10λ c × 10λ c is sufficient for accurate computation of statistical quantities, though the role of LSVs appears to become increasingly important with increasing domain size; we discuss this effect in our results.

Depth-averaged dynamics and energetics
For the purpose of investigating the inverse energy cascade, we decompose the vertical vorticity into a depth-averaged (barotropic) component, ζ , and a fluctuating (baroclinic) component, ζ , such that where, by definition, ζ = 0. The depth-averaged (barotropic) vorticity equation is then found by vertically-averaging equation (2.1), and is given by Thus, the barotropic dynamics are governed by a two-dimensional vorticity equation in which the sole forcing comes from convective dynamics, represented by the first term on the righthand side of the above equation. The barotropic, time-dependent kinetic energy density is defined as follows: (2.10) where · V indicates an average over the small, horizontal spatial scales, consistent with the notation employed in Plumley et al. (2018). In Fourier space, the barotropic kinetic energy equation is derived by multiplying the Fourier representation of equation (2.9) by the complex conjugate of − ψ k exp (ik · x), the spectral representation of ψ , and integrating over physical space to obtain where the box-normalized horizontal wavenumber vector is k = (k x , k y , 0), and k = |k|. This equations describes the evolution of the kinetic energy contained in the barotropic mode of wavenumber k that is due to (1) the interaction with the other barotropic modes, (2.12) (2) the interaction with the baroclinic, convective modes, ( 2.13) and (3) the viscous dissipation of the barotropic mode, (2.14) In the above definitions, the superscript * denotes a complex conjugate, F k [·] indicates the horizontal Fourier transform of the argument in square brackets, the symbol • indicates a Hadamard (element-wise) product, Re {·} is the real part of the argument in curly brackets and the sum is taken over all horizontal wavenumbers. The barotropic-tobarotropic and baroclinic-to-barotropic transfer functions T k and F k can be explicitly expressed in terms of a triadic interaction due to the Jacobian (i.e. non-linear) terms . The formation of LSVs is due to a positive contribution from T k and F k in equation (2.11) at the domain-scale wavenumber k = 1. As the LSV forms, the kinetic energy grows in time until dissipation balances the positive transfer at k = 1. Eventually a statistically stationary state is reached where D k ≈ −(T k + F k ), where we notice that, for these quantities, · is equivalent to an average over the fast temporal scale only; in contrast with previous work, all of the simulations presented here have reached this stationary state. Hereafter, in order to simplify notation we omit the averaging operator and refer only to the time-averaged values of K bt , T k , F k , and D k , unless otherwise stated.

Diagnostic quantities
Here we define several diagnostic quantities that will be used to characterize the dynamical state of the convective system. The heat transfer across the fluid layer is quantified by the non-dimensional Nusselt number N u = 1 + P r wθ . (2.15) In the present study the small-scale, or convective, Reynolds number is defined as where W rms = (W 2 ) 1/2 and w rms = (w 2 ) 1/2 are the rms values of the dimensional and non-dimensional vertical velocity component, respectively. The above definition is particularly useful for characterizing the amplitude of the convective motions, rather than the large-amplitude horizontal motions that occur in the presence of a strong inverse cascade. We also find it useful to refer to instantaneous values of the Nusselt and Reynolds number, and denote these by N u(t) and Re(t), respectively. Together with the barotropic kinetic energy (2.10) we will also consider the timeaveraged baroclinic, vertical and total kinetic energy densities, respectively defined as: With the above definitions, the Reynolds number can be expressed as Re = √ 2K z . As for N u and Re, we find it useful to refer to the instantaneous values of the total kinetic energy density as K(t).

Flow morphology: two-scale flows
The details of the simulations performed for this study are given in table 1. The choice of parameters allows us to refine the results of previous QG calculations (Julien et al. 2012b) in the range 1 P r 3 and for P r = 7, of particular relevance for laboratory experiments. The temporally averaged values of Re and N u displayed in table 1 are calculated over a temporal window in which the system reached a statistically stationary state. If an LSV is present in the domain, Re and N u might reach stationary values only when the barotropic kinetic energy has saturated, in accordance with previous studies (Julien et al. 2012b;Favier et al. 2014;Guervilly et al. 2014;Rubio et al. 2014) where N u has been shown to evolve over the time needed for the total kinetic energy to saturate. The interested reader is directed to supplemetary figure 1 for an example of this behavior. Figures 1a and 1c shows Re and N u as functions of Ra and P r. The continuous lines in figure 1a indicate the least-square fit to a power law of the kind Re = α r ( Ra− Ra c ) βr P r γr with α r = 0.1883, β r = 1.1512 and γ r = −1.2172. In figure 1b we illustrate the collapse of the Re data points to the law Re ∼ ( Ra − Ra c )P r −1 , empirically found and consistent with the the coefficients β r and γ r . Figure 1b suggests that the reduced Grashof number, RaP r −1 plays a key role in controlling the dynamics. Figure 1d shows the collapse of the N u data points to a power law of the kind N u ∼ ( Ra − Ra c ) 3/2 P r −1/2 , distinctive of the ultimate regime of thermal convection. Further details concerning these fits are given in section 3.5, but are shown here to summarize the cases that were computed.
Following Julien et al. (2012b), inspection of volumetric renderings of the fluctuating temperature (figure 2) suggests that we can qualitatively classify the flows into: the cellular regime (C); the convective Taylor column regime (CTC); the plume regime (P); and the geostrophic turbulence regime (G) . Regime C is only obtained close to the onset of thermal convection, i.e. for Rayleigh numbers not much larger than the critical value of Ra c 8.7; the CTC regime is characterized by columns that stretch across the fluid layer, surrounded by "sleeves" of oppositely signed vorticity (also visible in the fluctuating temperature) that prevent columns from interacting with each other; in the P regime the insulation mechanism weakens and column-column interaction shortens these structures, (a) Ra = 40, P r = 2 (CTC/P)  Table 1: Details of the numerical simulations. P r is the Prandtl number; Ra is the reduced Rayleigh number; N x , N y and N Z are, respectively, the number of Fourier modes in the x and y directions and the number of Chebyshev modes in the Z direction; ∆t is the timestep size used during the simulation; Re = w rms is the time-averaged, reduced Reynolds number based on the vertical component of the velocity; N u is the time-averaged Nusselt number; σ Re and σ N u are the standard deviations of Re and N u, respectively. The superscript † indicates that the horizontal box size for the simulation is taken to be 20λ c × 20λ c , where λ c = 2π/k c is the critical wavelength for the onset of thermal convection; for all other cases the box size is 10λ c × 10λ c . The superscript * indicated cases for which the influence of the lack or presence of an LSV in the initial condition on the saturated state has been checked (see section 3.3 for details).
transforming them into plumes; finally, geostrophic turbulence prevails at sufficiently large Rayleigh numbers where no obvious coherence in the fluctuating temperature field is observed. Although distinct transitions in the flow statistics can sometimes be used to separate these flow regimes (Nieves et al. 2014), an obvious distinction cannot always be made, e.g. cases ( Ra = 40, P r = 2) and ( Ra = 60, P r = 3) shown in figure 2, where plumes generated at each horizontal boundary seem to coexist with columns spanning the whole vertical extension of the computational domain. For a given value of P r, we observe the formation of LSVs as Ra is increased. These depth-invariant, dipolar vortices are readily identified from visual inspection of the geostrophic streamfunction ψ. Some representative cases are shown in figure 3. Crucially, we observe LSV formation for all P r values reported in table 1, including, for the first time to our knowledge, P r > 3. We find that for an LSV to be present in the domain, convection does not need to be in the geostrophic turbulent regime, as was previously suggested (Julien et al. 2012b;Stellmach et al. 2014).

LSV characterization
To quantify the presence of an LSV in the domain we analysed the time-averaged barotropic kinetic energy spectra K bt (k). We define flows in which the LSV is energetically dominant by the two conditions: K bt > K bc ; and K bt (k = 1) K bt (k > 1). As examples, figures 4a and 4b show the barotropic kinetic energy spectra for P r = 1 and P r = 2 over a range in Ra. The transition to LSV-dominant states occurs within the ranges 20 < Ra < 30 and 40 < Ra < 45 for the P r = 1 and the P r = 2 cases, respectively. As Ra is further increased beyond the transition, the LSV becomes increasingly dominant, as shown by the spectra. Note that for the ( Ra = 20, P r = 1) case, the barotropic spectra has a maximum at k = 1. However, for this case, the barotropic kinetic energy is not dominant, rather, we find that K bc 3K bt for this case. Therefore, there is no energetically dominant LSV in the domain for this particular case. A similar process is observed for all P r considered in the present study, although the threshold Rayleigh number for an LSV-dominant state increases with P r. However, we observe that the LSV becomes energetically dominant provided Re 5.812, independent of P r. This value of Re (indicated by the horizontal dashed line in figures 1a and 1b) corresponds to the case ( Ra = 45, P r = 2) and it is the lowest Re value for which LSV formation has been observed. The only exception is the ( Ra = 55, P r = 2.5) case for which no energetically dominant LSV is observed, although Re = 6.067 ± 0.134. This discrepancy can be explained by noting that these two values of Re are (considering their temporal fluctuations) consistent with each other and by admitting that the transition to the LSV-dominated regime is not abrupt. A more detailed exploration of the parameter space around the transition could reveal other exceptions to the threshold we identified and possibly a subtle P r dependence. In addition to the transition shown in the barotropic kinetic energy spectra, we also find (with increasing Ra) a distinct transition in the character of the three terms present in the A black vertical line is drawn in correspondence of T k , F k = 0. The insets highlight the behaviour for Ra 50. All quantities have been time-averaged over a statistically stationary state for which T k + F k ≈ −D k . spectral kinetic energy equation (2.11). In figure 4c we illustrate how the time-averaged, barotropic energy transfer T k +F k evolves with Ra for the specific case of P r = 2. We note that [T k + F k ] k=1 > 0 for all of the cases investigated, indicating that energy is always being transferred to the k = 1 mode, regardless of the value of Ra. However, we find that T k +F k changes from possessing a peak at k > 1, to then peaking at k = 1 for a sufficiently Analyzing all of our simulations shows that this transition occurs when Re 6.491, for any value of P r. This threshold value of corresponds to the simulation ( Ra = 135, P r = 7) and it is noted that no simulation with Re < 6.491 satisfies The only exception is the case ( Ra = 40, P r = 1.5), for which Re = 6.7529 ± 0.227 and the energy transfer at k = 1 is subdominant. Again, given the finite fluctuations in the Re values, we argue that a transition region may exists for which a simple threshold rule may not always work; a more detailed exploration of the parameter space may reveal subtle P r dependencies in the transition into the regime for which the barotropic energy transfer at k = 1 is dominant. Closer inspection of F k and T k separately (figure 4d) indicates that the baroclinic, convective dynamics primarily transfers energy to the barotropic dynamics around k 5 (notice the positive peak in F k for k 5). Energy transferred to the barotropic dynamics is then transferred upscale by the barotropic non-linear interactions, as indicated by negative values of T k for wavenumbers k > 4, and positive values at the largest scales. The inverse cascade that leads to LSV formation is therefore directly driven by the barotropic non-linear interactions in (2.9), whereas the energy is provided by the interaction of the baroclinic dynamics with the barotropic flows. The fact that k 0 T k = 0 confirms that the non-linear barotropic interactions do not inject or extract barotropic energy and that the saturation of K bt is controlled by the balance between energy injected from the baroclinic dynamics and energy dissipated through viscosity: This mechanism is in agreement with the formation of large-scale condensates in 2D calculations (Chertkov et al. 2007;Laurie et al. 2014) where the transfer of energy is due to the non-linear interactions between different scales of the 2D flow (equivalent to the barotropic-to-barotropic energy transfer described by T k ). This view is also confirmed by recent 3D studies (Buzzicotti et al. 2018). As mentioned in Section 2.2, in the saturated state D k ≈ −(T k + F k ). Examples are shown in figure 5 for Ra = 60 and P r = 7, P r = 2.5 and P r = 1. These cases are representative of, respectively, a low-turbulence case at the edge of the CTC and P regimes where Re 5.812 and, consequently, an inverse cascade that is not strong enough to drive LSV-formation; a case in the P regime with Re 5.812, slightly above the critical value for LSV formation and an inverse cascade; and a case in the G regime, with a robust inverse cascade possessing a strong peak at k = 1 driving an energetically dominant LSV. This figure illustrate how the pattern observed in figure 4 for fixed P r and increasing Ra can be discerned for fixed Ra and decreasing P r.
Following Guervilly et al. (2014), we can characterize the kinetic energy of the barotropic flow using the ratio of total kinetic energy to vertical kinetic energy, (3. 2) The factor of 3 in the denominator ensures that Γ → 1 if the kinetic energy is equipartitioned between the horizontal and vertical components of the velocity. Conversely, when the barotropic kinetic energy dominates, we expect this ratio to become significantly larger than unity. Figures 6a and 6b show Γ for all of the simulations as a function of Ra and Re, respectively. In agreement with the DNS calculations of Guervilly et al. (2014), Γ ≈ 1 for small values of Ra, then increases rapidly once LSVs begin to form. We find that Γ reaches a maximum of Γ max ≈ 5.5, that appears to be independent of the particular value of P r, though only the P r = 1 and P r = 1.5 simulations show a maximum value. For P r = 1, Γ reaches a maximum value at Ra = 80, whereas for P r = 1.5, Γ max occurs at Ra = 120, suggesting that the value of Ra at which Γ max is observed increases rapidly with Prandtl number.
Since the value of Ra at which LSVs begin to form is P r-dependent, Γ is also plotted as a function of Re in figure 6b. The data suggests that the evolution of Γ is uniquely determined by Re (or ( Ra − Ra c )P r −1 according to figure 1b) since all curves show selfsimilar behaviour, independent of the particular value of P r. The dashed vertical line denotes Re = 5.812, the approximate value at which the (k = 1) LSV becomes dominant. For cases in which Re is below this threshold value, Γ is close to 1 for all values of P r, and the convective pattern, or flow regime, for all of these cases can be qualitatively classified as cellular or convective Taylor columns. Above this threshold value of the Reynolds number, we find both plumes and eventually geostrophic turbulence as Re grows. For both the P r = 1 and P r = 1.5 cases, Γ max is reached for Re 24.
We emphasize that since all of the calculations in the present work were carried out with the QG model, the observed decrease of Γ for large values Ra is not due to a loss of rotational constraint. Although the DNS study of Guervilly et al. (2014) also report a decrease in Γ for sufficiently large forcing, their observed decrease might be caused by an increase in the Rossby number with increasing forcing. In the present study, Ro remains asymptotically small, regardless of the thermal forcing. Also, as pointed out previously, the LSV observed in Guervilly et al. (2014) is cyclonic, whereas the LSV observed in the present simulations is dipolar. Figures 6c and 6d show the barotropic kinetic energy versus Ra and Re, respectively. In are shown as reference, along with the vertical dashed line denoting the threshold Reynolds number Re = 5.812. The 's-shaped' behavior of the data, along with Γ , suggests that the barotropic mode is growing at an ever-decreasing rate as Ra is increased. Taking the barotropic kinetic energy K bt scaling with Re as illustrated in figure 6d and with K bc ∼ Re 2 (see supplementary figure 2), we can derive the expected evolution of Γ in the growing (5.812 Re 24) and decaying ( Re > 24) regimes. In the former, under the assumption that K bt K bc , K z , we obtain Γ ∼ Re 3/2 ; in the latter, taking Re → ∞, we obtain Γ ∼ Re −1/2 . These slopes are illustrated in figure 6b, for reference.
To better understand the change in scaling behavior of the barotropic kinetic energy with increasing Rayleigh number, we examine the nonlinear convective forcing term in the barotropic vorticity equation (2.9). In particular, the nonlinear baroclinic term can be written as baroclinic vorticity, defined as and analogously for the cross correlation for the y-component of the baroclinic velocity field and the baroclinic vorticity, C(v , ζ ). We note that this definition leads to C = 0.5 for perfect correlation between one component of the baroclinic velocity vector and the baroclinic vorticity, since the statistics are isotropic in the horizontal plane when sufficiently time-averaged. The coefficients were computed over the entire investigated range of Ra for the case P r = 1. Figure 7 shows the average value (3.5) We observe that C(u , ζ ) decays as Ra is increased from Ra > 20, suggesting one possible reason for the reduced rate of growth of the barotropic kinetic energy with increasing Rayleigh number.

The influence of initial conditions
For the cases indicated by the superscript * in table 1, additional simulations were carried out to test the influence of initial conditions on the occurrence of LSVs. In particular, for cases capable of forming an LSV ( Re > 5.812), we checked that the kinetic energy of the saturated state is independent of the presence of an LSV in the initial condition. Our results indicate that both baroclinic (or convective) amplitude (measured by Re) and the barotropic kinetic energy in the saturated state do not depend on the initial condition, but only depend on P r and Ra. As an example, in figure 8 we show this by illustrating the time evolution of the kinetic energy per unit volume, K(t), and the vertical Reynolds number, Re(t), for the case ( Ra = 50, P r = 1.5). Two simulations  Figure 8: Influence of initial conditions on large-scale vortex (LSV) formation. Kinetic energy (a) and reduced vertical Reynolds number (b) for the parameters P r = 1.5 and Ra = 50 from two different initial conditions: the case marked as "random IC" has a random initial condition with no initial LSV present; "LSV IC" marks an initial condition with a well-developed LSV in the system.
were run for this case: one started from a random initial condition that does not contain a pre-existing LSV (labeled as "random IC" in figure 8), and one started from an initial condition with an LSV already present in the domain, given from the saturated state of the ( Ra = 40, P r = 1) case (labeled as "LSV IC"). In the former case, the initial growth of kinetic energy is due to the formation of the LSV due to the imbalance T k + F k > |D k | for k = 1. In the latter case, the LSV initially present in the system was formed at a higher Re and due to a stronger inverse cascade than the one developed for ( Ra = 50, P r = 1.5). Therefore, initially the imbalance T k + F k < |D k | leads to a kinetic energy decay as the LSV cannot be energetically sustained. For both cases, a new state is eventually reached for which, statistically speaking |D k | = T k + F k . Similarly, we also found that for cases in which Re < 5.812, an LSV would eventually decay if it was present in the initial conditions. This result is in contrast with DNS calculations where large-intensity, domain-scale cyclonic vortices appear to be long lived when injected in a convective system in which large-scale structure would not spontaneously form (Favier et al. 2019).

The influence of box dimensions
The horizontal dimensions of the simulation domain are represented in terms of integer multiples of the critical wavelength λ c . We indicate the horizontal size of the computational domain by n c λ c × n c λ c , with n c being an integer. Most of the simulations were carried out with horizontal dimensions of 10λ c ×10λ c (i.e. n c = 10), which represents a trade-off between using a box size that is large enough to allow for computing converged statistics, and using a horizontal spatial resolution that is computationally feasible for an extensive exploration of the parameter space. Previous work has used values up to n c = 20, but, to our knowledge, no systematic investigations of the box size on key quantities such as the Nusselt number and Reynolds number have been reported for rotating convection. For non-rotating convection, however, Stevens et al. (2018) showed that surprisingly large box dimensions are needed to obtain convergence in all statistical quantities; in contrast, the same authors found that globally averaged quantities such as the Nusselt number converged with relatively small box dimensions. For the present work we have carried out simulations for fixed Rayleigh number and Prandtl number of Ra = 40 and P r = 1. A robust LSV is present with this parameter combination. Timeseries of these simulations are available in the supplementary material (see supplementary figure 1). Figure 9 shows the convective Reynolds number and Nusselt number, and the barotropic and baroclinic kinetic energy for a range of box sizes. The solid lines in figure  9(a) show the Nusselt and Reynolds number for a simulation in which the barotropic mode was set to zero at each timestep. We observe a nearly 23% increase in the heat transport when the barotropic mode is not present. This result might be interpreted in terms of the horizontal mixing that is induced by the barotropic mode; the vertical transport of heat is reduced when horizontal motions sweep heat laterally. In addition, we find that the Reynolds number is reduced by ≈ 4.4% with respect to the n c = 20 case. This observation suggests that the inverse cascade plays a relatively small role in influencing the amplitude of the convective flow speeds.
An estimate for the intensity of the LSV based on the domain size can be made from the following simple argument. When a well developed LSV is present, the dominant component of the kinetic energy spectra (k 5) scales approximately as K bt (k * ) ∼ k * −3 (Kraichnan 1967;Smith & Waleffe 1999;Rubio et al. 2014), where k * = k k box , with k box = 2πL −1 box and L box = n c λ c , is the dimensional box-scale wavenumber. Calculating the total kinetic energy we obtain: where K bt K bt (k * = k box ) when a robust LSV is observed in the system. In particular, by doubling the linear size of the (squared) domain the kinetic energy of the LSV is  Table 2: Least-squares fits to the Reynolds number, Re = α r ( Ra − Ra c ) βr P r γr (for data encompassing multiple P r) or Re = α r ( Ra − Ra c ) βr (when a single P r is considered).
allowed to quadruple in magnitude. The DNS study of Favier et al. (2014) also observed an increase in the barotropic kinetic energy with increasing box size. Our QG data shown in figure 9b is supportive of this quadratic dependence on box size.

Scaling laws for the baroclinic dynamics
Here we discuss least-squares fits to the baroclinic quantities Re and N u. Power law scalings were computed from all data collected in this study (see table 1 and figures 1a) for various subsets of Ra and P r. For Re with varying P r, we used power law fits of the form where α r , β r and γ r are all numerically computed constants. For constant P r, we used The numerically computed constants are denoted by α r , β r and γ r and given in table 2. Fitting to all available data reported in this study (including the P r = 10 dataset from Calkins et al. (2016)) we obtain (α r , β r , γ r ) = (0.1883, 1.1512, −1.2172). We notice that these values are not too different from a linear scaling of the form Re ∼ RaP r −1 , again suggesting that the reduced Grashof number plays a key role in controlling the dynamics. For many of the cases we find that β r is closer to unity when a single value of P r is used. Figure 10a shows the compensated Reynolds number ReP r/ Ra, where we see that there is a range of Ra values over which this scaling provides a reasonably good fit. However, significant departure from this linear Grashof number scaling is observed for the lower values of P r, i.e. those simulations that are characterized by the largest values of Re. Interestingly, this departure seems to be correlated with the behavior of the kinetic energy ratio Γ ; the largest departures from the linear Grashof number scaling are observed for cases that possess the peak Γ max , i.e. those cases in which Re 24.
We note that because the QG model employed here is asymptotically reduced, the Ekman number does not appear explicitly in the governing equations. However, we can relate our small-scale Reynolds number to the large-scale Reynolds number typically employed in DNS studies by noting that the convective length scale and fluid depth are and with values of α r , β r and γ r reported in (3.11) The above scaling is consistent with the recent spherical convection study of Guervilly et al. (2019). However, we note that although the above large-scale Reynolds number scaling is diffusion-free, the corresponding small-scale scaling is not. Analogous least-squares fits to the Nusselt number (N u) are given by (3.12) or for fixed values of Prandtl number. Results for various subsets of the explored parameters space are given in table 3. Figure 11b shows the compensated N u according to (3.12) using all available data from the present study. From table 3 we see that the same fit using only the 1 P r 2 cases suggests a fit that is roughly consistent with the ultimate scaling (3.14) in agreement with Julien et al. (2012b). Cases characterised by a lower Re ( e.g. P r 3) lead to a lower value for the exponent β n while using only the highest Re cases (P r = 1, Ra 120) leads to a higher value of β n . Compensated N u values, based on the ultimate  Table 3: Least-squares fits to N u = α n ( Ra − Ra c ) βn P r γn (for data encompassing multiple P r) or N u = α n ( Ra − Ra c ) βn (when a single value of P r is considered). scaling (3.14), are illustrated in figure 11a. This plot suggests that cases that reach the ultimate N u regime are the same as those that are characterised by an increase in reduced Re values in figure 10, and a corresponding drop in Γ in figure 6b.

Balances
Vertical profiles of the horizontal rms of each term present in the baroclinic vertical vorticity, vertical momentum and fluctuating heat equations are shown in figures 12a, 12b and 12c, respectively, for the most extreme calculation of Ra = 200 and P r = 1 ( Re ≈ 84). All of the quantities shown have been time-averaged. As shown previously (Julien et al. 2012b), within this high-Ra regime, the dominant terms in the governing which shows that horizontal advection of all these quantities is a key characteristic of this regime. Close inspection of the first of these balances reveals that, as Re grows, the advection of vorticity is increasingly dominated by the advection due to the barotropic flow, i.e. J[ ψ , ζ] J[ψ , ζ] for Re 0. On their own, the 'balances' given above reveal little about the resulting dynamics. Higher-order, or subdominant, effects are necessary in the dynamics, especially with regard to heat transport. Figure 12b suggests that small differences between the rms values of ∂ t w and J[ψ, w] are necessary to balance the vertical pressure gradient, ∂ Z ψ. This perturbative effect repeats again at even higher order, as figure 12b shows that the buoyancy force and vertical viscous force are approximately balanced, i.e.
Moreover, we find a subdominant balance in the fluctuating heat equation between the advection of the mean temperature and horizontal thermal diffusion, To better understand the role of the subdominant balance between viscosity and buoyancy, we show in figure 13 the ratio of of the vertical components of the rms viscous force, F v,z = ∇ 2 ⊥ w, to the rms buoyancy force, F b,z = RaP r −1 ϑ, both as a function of Ra and Re. Surprisingly, this ratio is an increasing function of Ra (and therefore also of Re). This result is in stark contrast to non-rotating convection in which viscous forces become ever smaller (relative to other forces) with increasing Rayleigh number. Indeed, the socalled 'free-fall' scaling for convective flow speeds, characterized by a balance between buoyancy and inertia, relies on the influence of viscosity being weak (e.g. see Yan et al. 2019). All of the different Prandtl number cases appear to show qualitatively similar behavior, and, as figure 13(b) shows, a reasonable collapse of the data can be obtained when the force ratio is plotted versus the Reynolds number.

Discussion and conclusions
A systematic investigation of rapidly rotating convection was carried out to determine the necessary conditions under which large-scale vortices (LSVs) form, and how the amplitude of such vortices and associated convective flow speeds scale with the input parameters. To achieve the extreme parameter regimes that are thought to be representative of natural systems such as planetary and stellar interiors, we have made use of an asymptotic description of the governing equations that rely on the assumption of a leading-order geostrophic balance. Varying the thermal Prandtl number has allowed us to determine the influence of fluid properties on the convective dynamics, and has also allowed for a more detailed control of the convective Reynolds numbers over our investigated range of Rayleigh numbers.
The LSVs form as a consequence of an inverse cascade that transports kinetic energy from small scale, convective motions up to the system-wide scale, characterized by a box-normalized wavenumber of k = 1. These LSVs grow in time until the energy input from the convection is balanced by large-scale viscous dissipation. All of the simulations presented show evidence of this equilibration process, regardless of the particular values of P r and Ra. We find that LSV-dominant convection can be characterized by a critical convective Reynolds number Re ≈ 6 across the range of investigated Prandtl numbers, in agreement with low-Ek DNS simulations performed at P r = 1 (Favier et al. 2014;Guervilly et al. 2014). Although an increase in P r leads to a concomitant increase in viscous dissipation for a fixed value of Ra, we find, for the first time to our knowledge, evidence of LSV-dominant convection in the "plume" regime. In particular, we observed the formation of barotropic vortices with a Prandtl number as high as P r = 7, a value that is representative of water at typical laboratory conditions. This finding suggests that LSVs may be detectable in laboratory experiments that use water as the working fluid. From the data reported in this study we can estimate a threshold value of Ra t 120 for LSV to form at P r = 7 which can be translated into large-scale Ra t for a given Ek via (2.5). State-of-the-art laboratory experiments can reach Ek = 10 −8 (Cheng et al. , 2019 giving Ra t 5.6 · 10 12 , a value for which heat transfer data suggest convection to be in a transitional regime between rotationally-dominated and non-rotating dynamics. We note that the presence of no-slip boundaries (not used in the present study) has been shown to suppress the formation of LSVs (Plumley et al. 2016), so additional studies are needed to determine the threshold for LSV-dominant convection with no-slip boundary conditions.
Several properties of LSVs have been studied. In agreement with the DNS study of Guervilly et al. (2014), we find that the relative size of the kinetic energy of the barotropic flow to that of the convection reaches a maximum value at a particular value of the Rayleigh number. However, with DNS studies, there is a corresponding increase in the Rossby number with increasing Rayleigh number. In contrast, the asymptotic model used here only captures the asymptotically small Rossby number limit, showing that this change in the growth of the LSV must be present in the rapidly rotating regime. When data from the entire range of P r is plotted as a function of Reynolds number, the peak in the barotropic kinetic energy occurs near Re ≈ 24. Therefore, the growth of the barotropic kinetic energy slows as the Rayleigh number is increased, suggesting that there is an optimum forcing level. We find that this change in behavior is related to a decrease in the velocity and vorticity correlations that are necessary to drive the inverse cascade. Our findings suggest that additional regimes, beyond the accessible limits of the present investigation, may be present in the convective dynamics as the Rayleigh number is increased further.
The horizontal dimensions of the simulation domain are shown to have a direct influence on the energy present in the LSV. It is found that the LSV energy grows quadratically with the horizontal dimension of the simulation domain (assuming domains of square cross-section), in agreement with DNS calculations (Favier et al. 2014). This finding is likely linked to the total available convective kinetic energy, which also grows quadratically with the horizontal dimensions of the simulation domain. Although a detailed investigation of the dynamical effect of this scaling was beyond the scope of the present investigation, this geometry-dependent effect may nevertheless have implications on the resulting dynamics.
The simulations suggest that there is no obvious scaling regime in the convective flow speeds with increasing Rayleigh number. A linear scaling of the form Re ∼ Ra/P r appears to collapse the data over a limited range in Ra, but the highest Re cases diverge from this scaling at the highest accessible values of Ra. We note that this linear scaling can be translated to an equivalent large-scale Reynolds number scaling of the form Re ∼ EkRa/P r, which has been noted in previous studies of rotating convection . Although this scaling is independent of the diffusion coefficients ν and κ when viewed on the largest scales of the system, the small-scales remain influenced by viscosity. Indeed, the simulations have revealed that the ratio of the rms viscous force to the rms buoyancy force in the vertical component of the momentum equation is an increasing function of Ra (or equivalenty Ra). This observation may simply be a result of the energetics of the Boussinesq system that requires the net heat transport to be balanced by viscous dissipation. In this regard, it might be argued that viscosity is fundamental to rotating convective dynamics.