# Bridging the Rossby number gap in rapidly rotating thermal convection

Adrian van Kan, Keith Julien, Benjamin Miquel, Edgar Knobloch

► To cite this version:

Adrian van Kan, Keith Julien, Benjamin Miquel, Edgar Knobloch. Bridging the Rossby num-
ber gap in rapidly rotating thermal convection. Journal of Fluid Mechanics, 2025, 1010, pp.A42.
10.1017/jfm.2025.290 . hal-05368049

HAL Id: hal-05368049

https://hal.science/hal-05368049v1

Submitted on 17 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Banner appropriate to article type will appear here in typeset article

# Bridging the Rossby number gap in rapidly rotating thermal convection

Adrian van Kan<fnref n="1" />*, Keith Julien<fnref n="2" />†, Benjamin Miquel<fnref n="3" />, Edgar Knobloch<fnref n="1" />

(Received xx; revised xx; accepted xx)

Geophysical and astrophysical fluid flows are typically driven by buoyancy and strongly constrained at large scales by planetary rotation. Rapidly rotating Rayleigh-Bénard convection (RRRBC) provides a paradigm for experiments and direct numerical simulations (DNS) of such flows, but the accessible parameter space remains restricted to moderately fast rotation rates (Ekman numbers $Ek \gtrsim 10^{-8}$), while realistic $Ek$ for geo- and astrophysical applications are orders of magnitude smaller. On the other hand, previously derived reduced equations of motion describing the leading-order behaviour in the limit of very rapid rotation ($Ek \to 0$) cannot capture finite rotation effects, and the physically most relevant part of parameter space with small but finite $Ek$ has remained elusive. Here, we employ the rescaled rapidly rotating incompressible Navier-Stokes equations (RRRiNSE) – a reformulation of the Navier-Stokes-Boussinesq equations informed by the scalings valid for $Ek \to 0$, recently introduced by Julien et al. (2025b) – to provide full DNS of RRRBC at unprecedented rotation strengths down to $Ek = 10^{-15}$ and below, revealing the disappearance of cyclone-anticyclone asymmetry at previously unattainable Ekman numbers ($Ek \approx 10^{-9}$). We also identify an overshoot in the heat transport as $Ek$ is varied at fixed $\overline{Ra} \equiv RaEk^{4/3}$, where $Ra$ is the Rayleigh number, associated with dissipation due to ageostrophic motions in the boundary layers. The simulations validate theoretical predictions based on thermal boundary layer theory for RRRBC and show that the solutions of RRRiNSE agree with the reduced equations at very small $Ek$. These results represent a first foray into the vast, largely unexplored parameter space of very rapidly rotating convection rendered accessible by RRRiNSE.

**Key words:** to be added in proof

MSC Codes 76U05, 76F35

† Deceased on 14th April 2024.

* Email address for correspondence: adrianvankan@gmx.de

<footnote n="1">¹Department of Physics, University of California at Berkeley, Berkeley, CA 94720, USA.</footnote>

<footnote n="2">²Department of Applied Mathematics, University of Colorado, Boulder, CO 80309, USA.</footnote>

<footnote n="3">³CNRS, École Centrale de Lyon, INSA Lyon, Université Claude Bernard Lyon 1, Laboratoire de Mécanique des Fluides et d'Acoustique, UMR5509, F-69134 Écully, France.</footnote>
# 1. Introduction

The universe abounds with examples of turbulent flows that are driven by buoyancy and constrained by rotation. These highly complex flows shape our environment across scales: in protoplanetary disks, they mediate the formation of planetesimals (Lesur & Ogilvie 2010); in planetary atmospheres, including on Earth, they carry heat, momentum and moisture, crucially impacting the climate (Emanuel 1994; Dauxois et al. 2021; Siegelman et al. 2022); in the Earth's oceans (Marshall & Schott 1999; Cheon & Gordon 2019), they transport heat, carbon dioxide, and nutrients, thereby crucially influencing biotopes (Severin et al. 2014); in planetary and stellar cores, they are responsible for the generation of large-scale magnetic fields via the dynamo instability (King et al. 2010; Jones 2011; Schaeffer et al. 2017); in the subsurface oceans of the gaseous giants' moons, they shape the icy crusts of Europa, Enceladus, Titan, etc. (Pappalardo et al. 1998; Mitri & Showman 2008; Nimmo & Pappalardo 2016; Soderlund et al. 2024); in stars like the Sun (Miesch 2000; Fan 2021) they are responsible for large-scale magnetic fields as well as rich dynamics such as the 22-year cycle in solar activity.

Such geo- and astrophysical flows, in spite of their differences, typically share two important properties. Owing to their turbulent nature, their dynamics results from a broad range of interacting spatial and temporal scales. Accordingly, advection by the flow dominates over viscous forces when the Reynolds number is large, $Re_H = UH/ν \gg 1$, based here on the fluid layer depth $H$, a typical velocity of the flow $U$, and kinematic viscosity $ν$. Simultaneously, advection by the flow occurs on time scales much longer than the planetary rotation period, as measured by small Rossby numbers $Ro_H = U/(2ΩH) \ll 1$, where $Ω$ is the planetary rotation rate.

Also of particular importance is the ratio of the viscous and Coriolis forces providing an *a priori* external parameter referred to as the Ekman number $Ek = Ro_H/Re_H = v/(2ΩH^2)$, as well as the Prandtl number $σ = v/k$ characterising the fluid under consideration, where $κ$ is the thermal diffusivity. The strength of the thermal forcing is controlled by the Rayleigh number $Ra$, which is proportional to the temperature drop $Δ_T$ across the fluid layer. This number can be used to define an alternate *a priori* parameter, the convective Rossby number $Ro_{conv} = \sqrt{Ra/σ Ek}$, measuring the strength of the rotation relative to thermal forcing.

Celestial bodies typically feature extreme values for the non-dimensional parameters defined above, as summarised in table 1: rapid planetary rotation is reflected in very low Ekman and Rossby numbers, whereas large Reynolds numbers indicate the highly turbulent character of the flow. Such rapidly rotating, highly turbulent flows are challenging to analyze and specifically to simulate numerically owing to the sheer number of spatial and temporal degrees of freedom which they involve. One common approach for obtaining theoretical predictions at geophysically relevant parameters (specifically, $Ek \to 0$) relies upon identifying transitions between different flow regimes in a region of more moderate parameters ($Ek \gtrsim 10^{-8}$) that is amenable to direct numerical simulations (DNS) or laboratory experiments, together with the scaling laws for transport coefficients valid within each regime. These results are then extrapolated to extreme parameters on the assumption that the regimes and scaling laws observed for accessible parameter values extend to geo-/astrophysically relevant parameter values. Assuming that no nontrivial transition occurs outside the observed moderate parameter interval is a particularly strong assumption, akin to a leap of faith due to the current paucity of data at the most extreme parameter values.

The dynamics of the fluid flows found in the celestial bodies listed in table 1 is further rendered complex by additional ingredients such as spherical geometry, multiple contributors to the density, compressibility, and the presence of magnetic fields. Absent such complexities, the quintessential paradigm for investigating rotationally influenced buoyant flows is provided
<table><caption>Table 1: Non-dimensional parameter estimates for planetary (Schubert & Soderlund 2011), satellite (Soderlund 2019) and stellar interiors (Miesch 2000; Garaud & Garaud 2008). Estimates of the Rossby number are derived from the relation $Ro_H = Re_H Ek$.</caption><thead><tr><th>Celestial body</th><th>Ek</th><th>σ</th><th>Ro<sub>H</sub></th><th>Re<sub>H</sub></th></tr></thead><tbody><tr><td>Earth's outer core</td><td>10<sup>-15</sup></td><td>0.1</td><td>10<sup>-7</sup></td><td>10<sup>8</sup></td></tr><tr><td>Mercury (core)</td><td>10<sup>-12</sup></td><td>0.1</td><td>10<sup>-4</sup></td><td>10<sup>8</sup></td></tr><tr><td>Jupiter (core)</td><td>10<sup>-19</sup></td><td>0.1</td><td>10<sup>-10</sup></td><td>10<sup>9</sup></td></tr><tr><td>Europa (ocean)</td><td>10<sup>-12</sup></td><td>11.0</td><td>10<sup>-2.5</sup>–10<sup>-1.5</sup></td><td>10<sup>9.5</sup>–10<sup>10.5</sup></td></tr><tr><td>Ganymede (ocean)</td><td>10<sup>-10</sup>–10<sup>-13</sup></td><td>10.0</td><td>10<sup>-3.5</sup>–10<sup>1.5</sup></td><td>10<sup>9.5</sup>–10<sup>11.5</sup></td></tr><tr><td>Saturn (core)</td><td>10<sup>-18</sup></td><td>0.1</td><td>10<sup>-9</sup></td><td>10<sup>9</sup></td></tr><tr><td>Enceladus (ocean)</td><td>10<sup>-10</sup>–10<sup>-11</sup></td><td>13.0</td><td>10<sup>-3.5</sup>–10<sup>-1</sup></td><td>10<sup>7.5</sup>–10<sup>9</sup></td></tr><tr><td>Titan (ocean)</td><td>10<sup>-11</sup>–10<sup>-12</sup></td><td>10.0</td><td>10<sup>-3</sup>–1</td><td>10<sup>9</sup>–10<sup>11</sup></td></tr><tr><td>Neptune (core)</td><td>10<sup>-16</sup></td><td>10.0</td><td>10<sup>-6</sup></td><td>10<sup>10</sup></td></tr><tr><td>Uranus (core)</td><td>10<sup>-16</sup></td><td>10.0</td><td>10<sup>-6</sup></td><td>10<sup>10</sup></td></tr><tr><td>Sun (convection zone)</td><td>10<sup>-15</sup></td><td>10<sup>-6</sup></td><td>10<sup>-3</sup></td><td>10<sup>12</sup></td></tr></tbody></table>

by rapidly rotating Rayleigh-Bénard convection (RRRBC). A large number of studies have been published on this model system, which is very well suited for detailed experimental, numerical and theoretical studies, see for example Chandrasekhar (1953); Nakagawa & Frenzen (1955); Veronis (1959); Rossby (1969); Boubnov & Golitsyn (1986); Zhong et al. (1991); Julien et al. (1996); Knobloch (1998); Hart et al. (2002); Vorobieff & Ecke (2002); Boubnov & Golitsyn (2012); King et al. (2012); Stevens et al. (2013) and the recent review by Ecke & Shishkina (2023). In RRRBC, gravity is often taken to be antiparallel to the rotation axis on the assumption that the Froude number $Fr_Ω ≡ Ω^2H/g$, where g is the gravitational acceleration, is small. The resulting configuration is translation-invariant and believed to be relevant to the polar regions of planets, moons and stars as well as to their interiors, although the case of misaligned gravity and rotation axis has also been studied (Hathaway & Somerville 1983; Julien & Knobloch 1998; Julien et al. 2006; Miquel et al. 2018; Novi et al. 2019; Currie et al. 2020; Cai 2020; Julien et al. 2025a).

A major obstacle in the study of rapidly rotating convection, even in the highly idealised setting of RRRBC, is the difficulty, to date, of attaining the extremely low Ekman numbers, typically $10^{-12}$ or smaller (cf. table 1), that characterize geo- and astrophysical flows. By contrast, current state-of-the-art DNS of the full Navier-Stokes-Boussinesq equations that describe RRRBC have been restricted to $Ek \gtrsim 10^{-8}$ (Kunnen et al. 2016; Guervilly et al. 2019; Song et al. 2024a,b), despite enormous computational effort. Inclusion of magnetic fields in the simulations to capture the capability for dynamo action, as well as spherical shell geometry, further limits the reported range of nearly all studies to $Ek \gtrsim 10^{-7}$, $Re_H \lesssim 10^3$ (Schaeffer et al. 2017; Cooper et al. 2020; Mason et al. 2022; Kolhey et al. 2022; Majumder et al. 2024; Gastine & Aurnou 2023); simulations down to $Ek = 10^{-10}$ are available for linearised dynamics only (He et al. 2022). State-of-the-art laboratory experiments are likewise restricted to $Ek \gtrsim 10^{-8}$ (Shew & Lathrop 2005; Aurnou et al. 2015; Rajaei et al. 2017; Cheng et al. 2020; Lu et al. 2021; Hawkins et al. 2023; Pothérat & Horn 2025) although the Eindhoven TROCONVEX experiment is in principle capable of achieving yet lower $Ek \sim 5 \times 10^{-9}$, cf. Cheng et al. (2018). These accessible parameter values in numerical or laboratory experiments are orders of magnitude larger than the estimates for
realistic geo- and astrophysical applications. On the other hand, existing reduced equations (Julien et al. 1998; Sprague et al. 2006), known as the non-hydrostatic quasi-geostrophic (NHQG) equations, describe the dynamics in the limit $Ek \to 0$, but their applicability is, by construction, confined to asymptotically small values of $Ek$. This leaves an unexplored gap corresponding to small but finite Ekman numbers, precisely the regime that is most relevant for geo- and astrophysical applications.

It is extremely rare to have a numerical or laboratory approach that allows one to span the parameter range from prototypical laboratory and DNS settings all the way to realistic geophysical and astrophysical settings. Remarkably, taking advantage of the characteristics of rapidly rotating thermal convection leads to a new numerical approach that relies on appropriately rescaled equations, and enables exploration of this very wide parameter range.

In the following section we describe the rescaled rapidly rotating Navier-Stokes equations (RRRiNSE) and sketch the numerical method used to solve them, followed in § 3 by a description of the results obtained, focusing both on bulk transport properties and the structure of the boundary layers at the top and bottom. The paper concludes in § 4 with a discussion of the significance of the RRRiNSE reformulation and prospects for future work. An Appendix provides additional details on the flow statistics and structure.

## 2. Methodology

In this section, we succinctly present the RRRiNSE, an equivalent reformulation of the full Navier-Stokes-Boussinesq equations, informed by the scaling behaviour valid in the limit of vanishing $Ek$ (Julien et al. 1998, 2012b). The RRRiNSE formulation is described in detail in a recent paper (Julien et al. 2025b), where it was first introduced and validated against existing results in the literature, and where the full details of the numerical implementation are given.

In short, the RRRiNSE retain all terms in the Navier-Stokes-Boussinesq equations while rescaling each variable according to its leading-order asymptotic scaling behaviour valid for $Ek \to 0$ in terms of the natural small parameter of rapidly rotating convection $\epsilon \equiv Ek^{1/3}$ (Chandrasekhar 1953; Sprague et al. 2006). In particular, the rapid rotation generates a high degree of anisotropy, which is explicitly incorporated in the RRRiNSE formulation: vertical length scales are non-dimensionalised by the layer height $H$, while horizontal length scales are significantly smaller, given by $\epsilon H$. Consequently, time is measured in units of the viscous diffusion time based upon $\epsilon H$. As shown by Julien et al. (2025b), this procedure ensures an adequate preconditioning of the equations of motion at small but finite $Ek$, thus facilitating DNS in this highly relevant but unexplored regime.

### 2.1. Rescaled Navier-Stokes equation

The convective fluid layer confined by horizontal surfaces separated by a distance $H$ is subjected to rotation around the vertical axis $\Omega = \Omega \hat{e}_z$ and uniform gravity $-g \hat{e}_z$. The temperature difference between the bottom and the top plate is $\Delta_T > 0$. In the Oberbeck-Boussinesq limit, the fluid density varies linearly with temperature $\rho = \rho_0[1-\alpha(T-T_0)]$ near a reference state $(\rho_0, T_0)$, with both the thermal expansion coefficient $\alpha$, the fluid kinematic viscosity $\nu$, and the temperature diffusivity $\kappa$ considered constant. The dimensional governing equations for the solenoidal velocity field $U$, the temperature $T$ and the pressure $P$ are:

$$ \partial_t U + U \cdot \nabla U + 2\Omega \times U = -\frac{1}{\rho_0} \nabla P + \alpha g T \hat{e}_z + \nu \Delta U, \quad \nabla \cdot U = 0, \quad (2.1) $$

$$ \partial_t T + U \cdot \nabla T = \kappa \Delta T. \quad (2.2) $$
As described in detail by Julien et al. (2025b), the power behind the RRRiNSE formulation lies in their use of rotation-influenced characteristic scales for non-dimensionalisation. For convenience, we introduce the small parameter $\epsilon \equiv Ek^{1/3} = (\nu/2\Omega H^2)^{1/3}$ and measure lengths in units of the layer depth $H$ in the vertical, but in units of $\ell_* = \epsilon H$ in the horizontal, so that $(X, Y, Z) = (\epsilon H \tilde{x}, \epsilon H \tilde{y}, H \tilde{z})$. In these coordinates, the rescaled anisotropic gradient $\nabla = \epsilon^{-1}H^{-1}\tilde{\nabla}_{\epsilon} \equiv \epsilon^{-1}H^{-1}(\partial_{\tilde{x}}, \partial_{\tilde{y}}, \epsilon \partial_{\tilde{z}})$ and the Laplace operator becomes $\Delta = \epsilon^{-2}H^{-2}\tilde{\nabla}_{\epsilon}^2$, where we have introduced the anisotropic diffusion operator $\tilde{\nabla}_{\epsilon}^2 = \partial_{\tilde{x}\tilde{x}} + \partial_{\tilde{y}\tilde{y}} + \epsilon^2 \partial_{\tilde{z}\tilde{z}}$. All tilde will be dropped in the following. Time is non-dimensionalised using a diffusive time-scale based on the horizontal scale rather than the layer depth alone: $t_* = \ell_*^2/\nu = \epsilon^2 H^2/\nu$. Accordingly, the velocity scale is $U_* = \ell_*/t_* = \nu/(\epsilon H)$, giving the non-dimensional velocity $\mathbf{u} = \mathbf{U}/U_* \equiv (u, v, w)$; the pressure scale $P_* = \rho_0 U_*^2/\epsilon = \rho_0 v^2/(\epsilon \ell_*^2)$ is used to define the non-dimensional pressure $p = P/P_*$. This isotropic velocity scale, combined with our anisotropic coordinate system, results in horizontally-dominated material derivatives:

$$D_t^{\perp, \epsilon} = \partial_t + u \partial_x + v \partial_y + \epsilon w \partial_z = \partial_t + \mathbf{u}_{\perp} \cdot \nabla_{\perp} + \epsilon w \partial_z. \quad (2.3)$$

Here, the perpendicular symbol $\perp$ in the subscript denotes the projection upon the horizontal plane: specifically, $\nabla_{\perp} = (\partial_x, \partial_y, 0)^{\top}$ and $\mathbf{u}_{\perp} = \mathbf{u} - w \hat{e}_z$. The incompressibility condition becomes:

$$\nabla_{\perp} \cdot \mathbf{u}_{\perp} + \epsilon \partial_z w = 0. \quad (2.4)$$

The (dimensional) temperature $T$ is decomposed as the sum of three terms: an unstable linear stratification, a vertical profile $\Theta(z)$, and horizontal fluctuations $\theta(x, y, z)$ with zero horizontal mean:

$$T - T_0 = \Delta_T [1 - z + \Theta(z) + \epsilon\theta(x, y, z)], \quad (2.5)$$

where $0 \le z \le 1$. We thus formally embed in our decomposition one hallmark of rapidly rotating convection: horizontal temperature fluctuations are dominated by vertical variations. With these elements, the dimensionless governing equations for the horizontal fluctuations $(\partial_x, \partial_y) \ne (0, 0)$ of momentum and temperature are found to be:

$$D_t^{\perp, \epsilon} \mathbf{u} + \epsilon^{-1} \hat{e}_z \times \mathbf{u} = -\epsilon^{-1} \nabla_{\epsilon} p + \frac{\tilde{R}a}{\sigma} \theta \hat{e}_z + \nabla_{\epsilon}^2 \mathbf{u}, \quad (2.6a)$$

$$D_t^{\perp, \epsilon} \theta + w (\partial_z \Theta - 1) = \frac{1}{\sigma} \nabla_{\epsilon}^2 \theta, \quad (2.6b)$$

where $\tilde{R}a = Ek^{4/3}Ra$, and the Rayleigh number $Ra$ is given by $Ra = \alpha g H^3 \Delta_T / (\nu K)$ in terms of the thermal expansion coefficient $\alpha$ of the fluid and gravitational acceleration $g$. By denoting the horizontal average of any quantity using an overline

$$\bar{q}^{xy} \equiv \frac{1}{L_x L_y} \int_0^{L_x} dx \int_0^{L_y} dy \, q(x, y, z, t), \quad (2.7)$$

we write separately the equations governing the vertical profiles of horizontal velocity $(\bar{u}^{xy}, \bar{v}^{xy})$ and temperature $\Theta$:

$$\partial_t \bar{u}^{xy} + \epsilon \partial_z \bar{w} \bar{u}^{xy} - \epsilon^{-1} \bar{v}^{xy} = \epsilon^2 \partial_{zz} \bar{u}^{xy}, \quad (2.8a)$$

$$\partial_t \bar{v}^{xy} + \epsilon \partial_z \bar{w} \bar{v}^{xy} + \epsilon^{-1} \bar{u}^{xy} = \epsilon^2 \partial_{zz} \bar{v}^{xy}, \quad (2.8b)$$

$$\epsilon^{-2} \partial_t \Theta + \partial_z w \theta^{xy} = \frac{1}{\sigma} \partial_{zz} \Theta. \quad (2.8c)$$

Equations (2.4), (2.6) and (2.8) form the RRRiNSE system, which is analytically equivalent to the dimensional equations Eqs. (2.1), (2.2) as we have only rescaled, but not discarded any
terms in the equations. Numerically, however, these equations differ drastically in the sense that (2.4), (2.6) and (2.8) can be solved accurately for very small values of Ekman number, whereas the poor conditioning of the unscaled equations (2.1), (2.2) renders their solution impossible in practice whenever $Ek \lesssim 10^{-9}$ for reasons explained by Julien et al. (2025b).

## 2.2. Asymptotically reduced non-hydrostatic quasi-geostrophic equations

In the limit of rapid rotation $Ek, Ro \to 0$, one can derive a reduced set of governing equations for the flow (presented below) in a systematic and controlled way using asymptotic methods (Julien et al. 1998; Sprague et al. 2006; Julien & Knobloch 2007; Julien et al. 2012b). Here, we emphasize, in a quick sketch, how the rescaled equations (2.6) naturally lead to the reduced equations when the formal limit $\epsilon \to 0$ is considered. Upon inspection, if the scaling employed for deriving the rescaled equations holds and all dimensionless variables remain $O(1)$, the leading-order terms at order $O(\epsilon^{-1})$ in equation (2.6a) indicate geostrophic balance:

$$ \hat{\mathbf{e}}_z \times \mathbf{u}_{\perp} = -\nabla_{\perp} p, \qquad (2.9) $$

so that the pressure serves as the geostrophic streamfunction:

$$ \mathbf{u} = -\nabla \times (p\hat{\mathbf{e}}_z) + w\hat{\mathbf{e}}_z + O(\epsilon), \qquad (2.10a) $$

$$ = -\partial_y p \hat{\mathbf{e}}_x + \partial_x p \hat{\mathbf{e}}_y + w \hat{\mathbf{e}}_z + O(\epsilon). \qquad (2.10b) $$

To obtain prognostic equations for $p, w, \theta, \Theta$ we substitute expression (2.10) into the rescaled momentum equation and project the result upon $\hat{\mathbf{e}}_z$, obtaining at leading order in $\epsilon$:

$$ \partial_t w + J[p, w] + \partial_z p = \frac{\widetilde{Ra}}{\sigma} \theta + \nabla_{\perp}^2 w, \qquad (2.11a) $$

where the Jacobian $J[p, w] = \partial_x p \partial_y w - \partial_y p \partial_x w$ indicates horizontal advection of $w$ and $\nabla_{\perp}^2 = \partial_{xx} + \partial_{yy}$ emphasizes that horizontal diffusion dominates. We complement the vertical velocity equation with a governing equation for the streamfunction $p$ obtained by projecting the curl of the momentum equation on the vertical. Using incompressibility, one proves $\hat{\mathbf{e}}_z \cdot \nabla \times (\hat{\mathbf{e}}_z \times \mathbf{u}) = -\epsilon \partial_z w$, so that the leading-order contribution of $\hat{\mathbf{e}}_z \cdot \nabla \times (2.6a)$ yields:

$$ \partial_t \nabla_{\perp}^2 p + J[p, \nabla_{\perp}^2 p] - \partial_z w = \nabla_{\perp}^4 p. \qquad (2.11b) $$

Finally, one readily obtains the leading-order contribution to the governing equation for the temperature fluctuation $\theta$ and the mean temperature $\Theta$:

$$ \partial_t \theta + J[p, \theta] + w (\partial_z \Theta - 1) = \frac{1}{\sigma} \nabla_{\perp}^2 \theta, \qquad (2.11c) $$

$$ \epsilon^{-2} \partial_t \Theta + \overline{\partial_z w \theta^{xy}} = \frac{1}{\sigma} \partial_{zz} \Theta. \qquad (2.11d) $$

Equations (2.11) constitute the reduced Non-Hydrostatic Quasi-Geostrophic (NHQG) equations previously derived using a multiple-scale asymptotic expansion (Julien et al. 1998; Sprague et al. 2006). These equations filter out fast inertial waves and thin Ekman layers and govern the leading order geostrophically balanced flow.

The Ekman number dependence that remains in the governing equation (2.11d) for the mean temperature reflects the fact that the timescale $\tau = \epsilon^2 t$ associated with the evolution of $\Theta$ is much slower than the convective timescale. Previous studies (Sprague et al. 2006; Julien et al. 2012b) have established that, in the limit of horizontally extended domains where a large number of convective cells or plumes contribute to the horizontal spatial averaging, this term can be omitted for the computation of global diagnostics of a statistically stationary state $\partial_t \Theta \approx 0$. This assumption can be justified based on the following rationale: as the horizontal
extent $L_\perp = L_x = L_y$ of the domain increases, so does the number of independent plumes $N_p$ that feed back into the temperature profile. The central limit theorem then suggests that, in the ergodic limit (time and ensemble averages behave similarly), the amplitude of temporal fluctuations of the vertical temperature profile $\Theta(z)$ should decrease like $1/\sqrt{N_p} \approx 1/L_\perp$, justifying the neglect of the term $\partial_t \Theta$. In the context of the RRRiNSE formulation, this assumption has been successfully validated and has proved to shorten the transient that leads to the stationary state (Julien et al. 2025b). Thus in all simulations of the RRRiNSE (2.4), (2.6) and (2.8) and the NHQG equations (2.11) described in this paper, $\partial_t \Theta$ is set to zero, as is expected to be the case for a sufficiently large horizontal extent $L_\perp$ of the domain.

## 2.3. Details on the numerical method

We compute solutions to both the rescaled equations (2.4), (2.6) and (2.8) and the reduced equations (2.11) using the pseudo-spectral code Coral (Miquel 2021), a flexible platform for solving partial differential equations in a Cartesian domain. We assume periodicity in $(x, y)$ and thus decompose the problem variables using a Fourier basis in the horizontal, whereas a Chebyshev expansion is used in the vertical in order to impose boundary conditions on the bounding surfaces $z = 0, 1$ using a Galerkin basis recombination. Differential operators are discretised using the quasi-inverse method (Clenshaw 1957; Greengard 1991; Julien & Watson 2009) to optimize memory usage and accuracy. The resolution was varied from 192 grid points in each direction (corresponding to 128 alias-free modes) at the lowest values of $\overline{Ra}$, up to 768 grid points (512 Chebyshev modes) in the vertical direction and 576 grid points (384 Fourier modes) in the horizontal directions at larger $\overline{Ra}$. The time-marching relies on the family of stable semi-implicit Runge-Kutta schemes introduced by Ascher et al. (1997), specifically their four-stage, third-order scheme, known as ARS443.

For some of the runs presented here, we have used low-amplitude noise as the initial condition, and waited for the convective turbulent flow to settle into a statistically stationary state after a transient period. However, for strongly supercritical flow, this transient may prove much more turbulent than the stationary state itself. Hence, in some cases, we resort to an alternative: instead of low-amplitude noise, the initial condition for a given run was taken from the solution computed at a neighboring position in parameter space. This procedure facilitates the computation of the transient.

Once the transient has subsided, we compute time averages of different diagnostic quantities reported here. For instance, the barotropic vertical vorticity displayed in figure 2 is defined as:

$$ \overline{\omega}_z(x, y) = \int_0^1 \omega_z \, dz, \qquad (2.12) $$

where the vertical vorticity is given by $\omega_z = \partial_x v - \partial_y u$. The Nusselt number is computed in terms of a volume average as $Nu(t) - 1 = \sigma \langle w\theta \rangle_{\text{vol}}$, where the volume average is given by

$$ \langle \cdot \rangle_{\text{vol}} = \overline{(\cdot)^{xy}} = \frac{1}{L_x L_y} \int_0^{L_x} dx \int_0^{L_y} dy \int_0^1 (\cdot) \, dz, \qquad (2.13) $$

and the small-scale Reynolds number based on the vertical velocity and the horizontal length scale is given, in the non-dimensional RRRiNSE formulation, by $Re(t) = \sqrt{\langle w^2 \rangle_{\text{vol}}}$. In addition, we denote the combined space and time average by:

$$ \langle \cdot \rangle = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} dt \frac{1}{L_x L_y} \int_0^{L_x} dx \int_0^{L_y} dy \int_0^1 (\cdot) \, dz. \qquad (2.14) $$
**3. Results**

Here, we present the first DNS of RRRBC at realistically small Ekman numbers for planetary, satellite and stellar interiors, using the RRRiNSE formulation, and analyze the physics of the flow in this extreme parameter range. Our extensive set of simulations was performed using the pseudospectral solver Coral (Miquel 2021) for a periodic domain of horizontal extent $10\ell_c$, where $\ell_c \approx 4.82Ek^{1/3}$ is the critical wavelength in the limit $Ek \to 0$ (Chandrasekhar 1953). Specifically, we consider a fluid with $\sigma = 1$ subject to stress-free boundaries (for simplicity) at top and bottom where constant temperatures are imposed, and periodic boundary conditions in the horizontal directions. We explore unprecedentedly small but finite Ekman numbers, as low as $Ek = 10^{-15}$, values that are more than six orders of magnitude smaller than those previously attainable in DNS or in the laboratory.

**3.1. Overview of parameter space**

Figure 1a shows the parameter space explored in our simulations. The top panel shows the estimates of these parameters for planetary interiors and oceans on icy satellites given in table 1 in the plane spanned by the Ekman number $Ek$ and the bulk Reynolds number $Re_H$, together with blue circles showing the values reached in our simulations, indicating that we attain realistically small $Ek$ equal to the estimated values for the outer core of the Earth ($Ek \approx 10^{-15}$) and the convective zone of the Sun (not shown, cf. table 1). We further indicate the approximate range of parameters reached in previous numerical and laboratory studies, which are restricted to higher $Ek$. Dashed lines highlight important regime boundaries where the Rossby number $Ro$ equals unity in different parts of the flow. The green dashed line is an estimate of where the bulk of the flow is expected to become rotationally supported as the rotation rate increases, namely $Ro_H = 1$, i.e., $Re_H = 1/Ek$. In contrast, the rotation rate required to generate rotational support in the thermal boundary layers at the top and bottom is much higher. This rotation rate is determined by the condition that the local Rossby number in the thermal boundary layer $Ro_{BL} = 1$ (red dashed line). As shown by Julien et al. (2012a) this condition is satisfied when the convective Rossby number associated with the bulk flow $Ro_{conv} = Ek^{1/5}$. With velocity following the rotational free-fall scale, $U = \alpha g \Delta_T / (2\Omega)$, we see that $Re_H = Ro_{conv}^2/Ek$, cf. Julien et al. (2025b), and hence that $Re_H = Ek^{-3/5}$ (red dashed line in figure 1a).

In rapidly rotating flows, higher values of $Ra$ are required to trigger convection as a consequence of the Taylor-Proudman constraint (Proudman 1916). Specifically, the critical Rayleigh number for convection in the limit of small $Ek$ is proportional to $Ek^{-4/3}$ as shown by Chandrasekhar (1953). Therefore, it is natural to introduce $\tilde{Ra} \equiv RaEk^{4/3}$, the reduced Rayleigh number, which measures the supercriticality of rapidly rotating convection. Panel b of figure 1 shows an alternative representation of the parameter space in terms of the reduced Rayleigh number $\tilde{Ra}$ and the Taylor number $Ta = Ek^{-2}$. We reach Taylor numbers up to $10^{30}$ and reduced Rayleigh numbers up to 300, corresponding to rapidly rotating and highly turbulent flows. The associated bare Rayleigh number $Ra$, indicated by colors in the background, ranges over a similarly wide interval from $O(10^3)$ to very large values in excess of $10^{20}$.

This wide array of numerical simulations allows us to probe, for the first time in DNS, the physics of rotating convection from moderate Ekman numbers down to the fully rotationally constrained regime, previously only accessible in asymptotically reduced equations (Sprague et al. 2006; Julien et al. 2012a) within the $Ek \to 0$ limit.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Figure 1: Overview of the parameter space. Panel a: estimates of the non-dimensional parameters for different celestial objects (see table 1) shown in the plane spanned by the Ekman and Reynolds numbers. The parameters reached in previous laboratory experiments and simulations are indicated by shaded regions. The dashed lines indicate regime boundaries discussed in the text: the green dashed line corresponds to the bulk Rossby number $Ro_H = 1$ while the red dashed line corresponds to the local Rossby number $Ro_{BL} = 1$ in the thermal boundary layers (BLs). These lines are given by $Re_H = 1/Ek$ and $Re_H = Ek^{-3/5}$, respectively. The regime bounded by these two dashed lines is characterised by unbalanced boundary layers and a balanced bulk flow. The symbols refer to the simulations summarised in panel b: $\tilde{Ra} = \text{const.}$ (white diamonds), $Ra = \text{const.}$ (black triangles) and $Ek = \text{const.}$ (blue plus signs). Panel c: overview of simulations (with symbols identical to panel a) in the plane spanned by the reduced Rayleigh number $\tilde{Ra} = Ra Ek^{4/3}$ and the Taylor number $Ta = Ek^{-2}$. Dashed lines represent regime boundaries corresponding to those in panel a.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Top row: vertically averaged (i.e., barotropic) vertical vorticity field $\overline{\omega}_z(x, y)$ at $Ek = 10^{-15}$ (panel a), $10^{-8}$ (panel b), $10^{-7}$ (panel c), $10^{-6}$ (panel d) and $\overline{Ra} = 120$ (all panels). Panel e: skewness of the barotropic vorticity field versus $Ek$ indicating a strong cyclone-anticyclone asymmetry for $Ek \gtrsim 10^{-9}$, and an approximate (statistical) cyclone-anticyclone symmetry for $Ek \lesssim 10^{-9}$.</figcaption></figure>

## 3.2. Flow morphology – bulk and boundary layers

An important feature of rapidly rotating flows is the formation of large-scale flow structures (vortices or jets), associated with spectral condensation of kinetic energy at the domain scale as a result of an inverse energy cascade (Favier *et al.* 2014; Guervilly *et al.* 2014; Rubio *et al.* 2014). This behaviour, which is robust for rotating convection in different domain geometries (Guervilly & Hughes 2017; Julien *et al.* 2018; Lin & Jackson 2021), is also observed more broadly in anisotropic turbulent flows, such as stably stratified rotating turbulence, turbulence subject to strong magnetic fields or geometric confinement to a thin-layer geometry, among others (Alexakis & Biferale 2018; Alexakis 2023; van Kan 2024). In the case of a periodic, horizontally square Cartesian domain, the condensate in the limit of asymptotically large rotation $Ek \to 0$ takes the form of a large-scale vortex dipole, first observed in the reduced equations (Julien *et al.* 2012b; Rubio *et al.* 2014), resembling findings in two-dimensional turbulence (Boffetta & Ecke 2012). In contrast, direct numerical simulations of rotating convection at previously attainable moderately small Ekman and Rossby numbers reveal an asymmetry between cyclones and anticyclones, namely a preference for strong cyclonic structures and diffuse, weaker anticyclones (Favier *et al.* 2014). While a dominance of cyclones over anticyclones is also observed in stably stratified flows at moderately small Rossby numbers (Roullet & Klein 2010; Gallet *et al.* 2014), the transition from the moderately small Rossby number regime to very small $Ek$ (relevant to the celestial bodies listed in table 1), where the cyclone-anticyclone symmetry is expected to be restored, could not previously be reached in RRRBC.

To address this question, we illustrate and analyze in figure 2 how the vertically averaged (i.e., barotropic) vertical vorticity $\overline{\omega}_z(x, y)$ (see Sec. 2) changes as $Ek$ varies for a fixed supercriticality $\overline{Ra} = 120$. At moderate rotation rates ($Ek = 10^{-6}$ and $Ro \approx 0.4$, panel (d)), a striking asymmetry exists between the strong, coherent cyclone and a weak, incoherent
anticyclonic region. Here, the *a posteriori* Rossby number $Ro$ is computed as the root mean square of the ratio $||u \cdot \nabla u_{\perp}||/(2\Omega||u_{\perp}||)$, where $u_{\perp} = (u, v, 0)$. As rotation increases (top row of figure 2, right to left), this asymmetry decreases but remains visible for $Ek = 10^{-8}$ ($Ro \approx 0.07$, panel b), a value close to the current state of the art for RRRBC studies. The anticipated cyclone-anticyclone symmetry present in the asymptotic rapidly rotating limit $Ek \to 0$ based on the NHQG equations is in fact observed for $Ek = 10^{-15}$ ($Ro \approx 8 \times 10^{-5}$, panel a). To quantify this qualitative observation, we display in figure 2(e) the skewness of the barotropic vertical vorticity:

$$ \text{Skew}(\bar{\omega}_z) = \left\langle (\bar{\omega}_z - \langle \bar{\omega}_z \rangle)^3 \right\rangle / \left\langle (\bar{\omega}_z - \langle \bar{\omega}_z \rangle)^2 \right\rangle^{3/2}. \qquad (3.1) $$

A finite positive value, indicating the prevalence of strong cyclones, is obtained for $Ek \ge 10^{-9}$, whereas an approximately vanishing skewness is observed at smaller Ekman numbers, indicating a restored (statistical) cyclone-anticyclone symmetry, with a relatively sharp transition between the two regimes.

In addition to the barotropic component of the flow, which informs us about the bulk flow structure, the flow inside the thermal boundary layer near the stress-free, constant-temperature boundaries at the top and bottom is of great importance. Figure 3(a) shows visualizations of the flow inside the bottom thermal boundary layer (at height $z = \delta_{\omega_z}$, defined thereafter, cf. figure 5) at $\tilde{Ra} = 120$ for four distinct Ekman numbers $Ek = 10^{-15}, 10^{-8}, 10^{-7}, 10^{-6}$ in terms of three fields: the vertical vorticity $\omega_z$, the temperature perturbation $\theta$ away from the horizontally averaged temperature profile and the vertical velocity $w$ at the top of the (bottom) boundary layer, where the flow transitions between qualitatively distinct regimes. Panel b shows histograms of each field at the same Ekman numbers, while panel c displays the associated statistical skewness versus $Ek$. At very low Ekman numbers, $Ek = 10^{-15}$, where the bulk flow features two large-scale vortices and displays an approximate statistical cyclone-anticyclone symmetry, the same symmetry is obeyed in the boundary layer as evidenced by the approximately vanishing skewness not only of $\omega_z$, but also $\theta$ and $w$ (see also Supplementary Movies 1-3). A near perfect correlation can be seen between the vertical vorticity and $\theta$. We mention that the RRRiNSE formulation remains numerically stable even when the Ekman number is reduced further, as far as $Ek = 10^{-24}$; the corresponding flow fields obtained from DNS are shown in Appendix A.

As the Ekman number increases, the approximate symmetry of the histograms persists up to $Ek \approx 10^{-9}$ (not shown). At $Ek \approx 10^{-8}$, where cyclone-anticyclone asymmetry and cyclone dominance emerges in the barotropic mode, the boundary layer flow morphology is also drastically altered: in addition to the trace of the strong large-scale cyclone, similarly strong, albeit short-lived, anticyclonic vortices emerge in the boundary layer, whose anticyclonic core is seen in some cases to be surrounded by a cyclonic shield (see also Supplementary Movie 4). These strong vortices with anticyclonic cores have a clear signature in the temperature fluctuation $\theta$ as strong, cold perturbations (in the bottom boundary layer shown here, see Supplementary Movie 5), and also in the vertical velocity $w$ which shows descending flow inside the vortex core, surrounded by rising fluid motion (see Supplementary Movie 6). As a result, the statistics of $\theta$ and $w$ are skewed towards negative values at this $Ek$. It is important to stress that the structures displayed in figure 3 and discussed above are confined within the boundary layer, and do not extend into the bulk flow. To the best of our knowledge, the strong, partially shielded, short-lived anticyclonic vortices in the boundary layer with stress-free boundaries are described here for the first time, although shielded structures have been previously observed in various studies of rotating convection (Sprague et al. 2006; Julien et al. 2012b; Stellmach et al. 2014; Song et al. 2024a). We also note that similar structures were predicted based on Ekman layer theory for no-slip boundary conditions by Julien et al.
<figure><img src="image_5.png" /></figure>

Figure 3: Panel a: overview of boundary layer flow morphology in terms of the vertical vorticity $\omega_z$, temperature perturbation $\theta$, and vertical velocity $w$ at different Ek and $\bar{Ra} = 120$. Panel b: histograms of $\omega_z, \theta, w$ at different Ek. Near-Gaussian statistics are observed at $Ek = 10^{-15}$ with larger skewness at larger Ek. Panel c: skewness of $\omega_z, \theta, w$ versus Ek. All results shown are obtained at $z = \delta\omega_z$ (top of the momentum boundary layer near the bottom boundary) defined in figure 5(a).
(1996) in response to a Gaussian temperature perturbation. As *Ek* is increased to 10<sup>-7</sup> and 10<sup>-6</sup>, cyclones clearly dominate the flow in the boundary layer (as in the bulk), while the strong localised anticyclones become more diffuse, leading to a significant skewness in the statistics of ω<sub>z</sub>. The responses in θ and *w* remain visible in the boundary layer, particularly outside the large-scale cyclone, although their skewness decreases with further increase in *Ek*.

3.3. Flow statistics

Going beyond the visualizations shown in figure 3, we display in figure 4 key quantitative
information on the flow in terms of important non-dimensional numbers. We define the non-
dimensional heat flux across the layer, measured by the Nusselt number Nu − 1 ≡ σ⟨wθ⟩
where ⟨·⟩ is the combined volume and time average, and denote by Nu(t) the instantaneous
Nusselt number obtained by averaging only over the volume, not time. In addition, we
measure the small-scale Reynolds number Re, based on the horizontal length scale and the
root-mean-square vertical velocity, by taking the time average of the instantaneous Reynolds
number Re(t) = √⟨w²⟩<sub>vol</sub> (see Sec. 2). Both Nu and Re are commonly used to quantify
turbulent convection since their statistics are known to converge quickly in the nonlinear
regime following the initial exponential growth phase of the convective instability, see e.g.
Maffei et al. (2021). This remains true even in the presence of an inverse cascade leading
to a slow growth in the horizontal kinetic energy. Figure 4(a) displays the time series
of the instantaneous Nusselt number at Ek = 10<sup>-10</sup> for three different reduced Rayleigh
numbers Ra = 40, 60, 80, showing increasing turbulent heat flux (mean and fluctuations)
with increasing Ra. Rescaling the Nusselt numbers by the well-known turbulent scaling law
(Julien et al. 2012b) for rapidly rotating convection Nu − 1 ~ Ra<sup>3/2</sup> (see panel b), the data
collapse satisfactorily for Ra = 60, 80, while the least turbulent simulation at Ra = 40 shows
a small mismatch. The corresponding theoretical scaling prediction for the Reynolds number
is Re ~ Ra, which we use to rescale our Reynolds number data.

Averaging over time series like those shown in the top row of figure 4, we obtain the Nusselt and Reynolds number statistics in panels *c* and *d* of figure 4, shown versus the Taylor number *Ta* = *Ek*⁻². Specifically, four sets of simulations at constant *Ra* = 40, 60, 80, 120 are shown for a wide range of 10² ≤ *Ta* ≤ 10³⁰. Triangles in different orientations for each *Ra* show the observed average *Nu* and *Re*, while error bars indicate the measured standard deviation about those mean values. Horizontal dashed lines indicate, for each *Ra*, the values corresponding to the reduced NHQG equations, with a shaded region indicating the corresponding standard deviation measured in the NHQG equations. At small *Ta* (large *Ek*), for the given *Ra*, the flow is only weakly supercritical and therefore features only small *Nu* and *Re*. A break is seen in the curves around *Ta* = 10⁶ (*Ek* = 10⁻³), where the bulk convective Rossby number equals one, marked by faint vertical dashed lines, indicating that the bulk flow begins to be rotationally constrained. The inverse cascade sets in around *Ta* = 10⁸ for *Ra* = 40, and closer to *Ta* = 10¹⁰ for larger *Ra* = 60, 80, 120 (hollow circular markers). This is in agreement with the literature, see e.g. Favier *et al.* (2014). As *Ta* increases (*Ek* decreases) further from this point, an overshoot in *Nu* is observed leading to a maximum in the convective heat transport for the given control parameter *Ra*, whose amplitude and location depend on *Ra*. This is followed by convergence within one standard deviation to the value obtained in the NHQG equations at yet larger *Ta*, beyond a threshold indicated by a second set of vertical dashed lines, one for each *Ra*, corresponding to convective Rossby number in the boundary layer of order unity, a condition equivalent to *Ek* = σ³/² *Ra*⁻¹⁵/⁴ (Julien *et al.* 2012a), i.e., the rotation rate where the boundary layer also becomes rotationally constrained. These regime
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

Figure 4: Panel a: sample time series of $Nu - 1$ at $Ek = 10^{-10}$ for $\tilde{Ra} = 40, 60, 80$ in the statistically stationary state, revealing turbulent fluctuations about a well-defined mean value. Both mean and amplitude of fluctuations increase with $\tilde{Ra}$. Panel b: same data, collapsed by rescaling the y axis by $\tilde{Ra}^{3/2}$. Panel c: temporally averaged $Nu - 1$ in the stationary state, compensated by $\tilde{Ra}^{3/2}$, with error bars indicating standard deviation, vs $Ta = Ek^{-2}$. $Nu$ increases with $Ta$ at fixed $\tilde{Ra}$ up to a maximum, beyond which it decreases to converge to a $Ta$-independent constant, which agrees with the value obtained in the NHQG equations at the same $\tilde{Ra}$ (horizontal dashed line), within one standard deviation (shaded area). Vertical dashed lines (same color scheme) correspond to $Ro_{conv} = 1$ and $Ro_{BL} = 1$. Panel d: corresponding data for the small-scale Reynolds number $Re \equiv \sqrt{\langle w^2 \rangle}$, compensated by $\tilde{Ra}$, show a similar structure with an overshoot and eventual convergence to the NHQG prediction at large $Ta$. Panel e: a posteriori Rossby numbers versus $Ek^{1/3}$ for the bulk and the boundary layer, both of which exhibit collapse in panel f when shown versus $Ek^{1/3} \tilde{Ra}^{5/4}$.
boundaries are equivalent to those shown in the top panel of figure 1. A similar structure is observed in the Reynolds number data, which collapse slightly less satisfactorily under the rescaling (this is known to be a consequence of the presence of the large-scale vortex (Maffei et al. 2021; Oliver et al. 2023)). At small $Ta$, the flow is only weakly supercritical and therefore the Reynolds number is small. It increases with $Ta$ up to a maximum, beyond which it decreases to converge (within one standard deviation) to the value obtained from the NHQG equations. Since Taylor numbers in most laboratory experiments or unrescaled DNS have been limited to $Ta \lesssim 10^{14} - 10^{16}$, this $Nu$ and $Re$ vs $Ta$ dependence has not previously been described. In particular, the presence of a maximum in the heat transport and turbulence intensity at a finite Taylor (or equivalently, Ekman) number as well as the statistical convergence to the NHQG equations are reported here for the first time.

In addition to the *a priori* (bulk/boundary layer) convective Rossby numbers shown by the dashed vertical lines in panels c and d of figure 4, panels e and f show the *a posteriori* Rossby number $Ro$ computed as the root-mean-square of the ratio $||u \cdot \nabla u_{\perp}|| / (2\Omega ||u_{\perp}||)$ in the bulk and in the boundary layer, where it takes distinct values (see Appendix A). When shown versus $\epsilon = Ek^{1/3}$, the standard small parameter of rotating convection, the data, while scattered, clearly show that the bulk Rossby number increases monotonically with $Ek$ as expected, crossing unity near the observed onset of the inverse cascade. It is interesting that this occurs at Ekman numbers which are two orders of magnitude smaller than what may naively be expected based on the convective Rossby number $Ro_{conv} = 1$. When the same data is plotted against the rescaled parameter $Ek^{1/3} \overline{Ra}^{5/4}$, it collapses to a remarkable degree, both in the bulk and in the boundary layer. We stress that this is somewhat surprising, since this effective small parameter is obtained explicitly from the derivation of the reduced equations for the boundary layer (Julien et al. 2012a,b). However, we note that the bulk flow develops cyclone-anticyclone symmetry at similar Ekman numbers as the boundary layer, which is another empirical link between the bulk and boundary layer flows observed here. The collapse of the data in figure 4e,f implies in particular that when $Ek^{1/3} \overline{Ra}^{5/4} = 1$, below which the average Nusselt and Reynolds numbers converge to the NHQG data, there is a well-defined associated critical Rossby number in the boundary layer, which is found to be approximately 0.2. This is where the boundary layer loses rotational constraint and becomes unstable to strong ageostrophic motions with pronounced vertical gradients. These additional vertical gradients are associated with additional dissipation contributions. In steady state, the dissipation rates of kinetic energy and temperature variance are directly linked to the Nusselt number, via the so-called power integrals. In the RRRiNSE formulation, these read

$$ \langle |\nabla_{\perp} \theta|^2 \rangle + \epsilon^2 \langle (\partial_z \theta)^2 \rangle + \langle (\partial_z \Theta)^2 \rangle = Nu - 1, \quad (3.2a) $$

$$ \langle |\partial_x u|^2 + |\partial_y u|^2 \rangle + \epsilon^2 \langle |\partial_z u|^2 \rangle = \frac{\overline{Ra}}{\sigma^2} (Nu - 1), \quad (3.2b) $$

where $\langle \cdot \rangle$ is the combined time and volume average over the whole domain. Directly measuring all contributions in our simulations, we find that both balances are approximately satisfied (indicating that nonstationarity is weak, see Appendix A), with important contributions to the dissipation of temperature variance originating in the boundary layers near $Ek^{1/3} \overline{Ra}^{5/4} = 1$. The increase in the Nusselt number away from the small $Ek$ NHQG limit can therefore be understood as a consequence of increased dissipation due to enhanced (vertical and horizontal) gradients, at least in part due to the strong anticyclonic vortical structures appearing in the boundary layer at $Ek = 10^{-8}$ (which is near $Ek^{1/3} \overline{Ra}^{5/4} = 1$ for $\overline{Ra} = 120$), cf. figure 3.
### 3.4. Boundary layer structure

From the above analysis, it is apparent that the boundary layers play a nontrivial role in the transition from small to large $Ek$, controlling the departure from the NHQG branch. To further characterize the boundary layer structure, we measure and analyze vertical root-mean-square profiles of different quantities as a function of $z$. Figure 5(a) shows a log-log plot of the thermal boundary layer depth $\delta_\theta$ computed as the location of the maximum in the root-mean-square vertical profile of $\theta$ (shown in the inset for different $Ek$ at $\tilde{Ra} = 120$). At small Ekman numbers, $\delta_\theta$ converges to a constant value for a given reduced Rayleigh number, corresponding to the NHQG limit. Panel c confirms that the thermal boundary layer depth scales with $\tilde{Ra}^{-15/8}$, which is the scaling obtained from reduced NHQG equations for the boundary layer (Julien et al. 2012a,b). As $Ek$ is increased close to the point where the boundary layer loses rotational support ($Ek = \sigma^{3/2} \tilde{Ra}^{-15/4}$, not shown), the boundary layer depth undershoots before increasing with a power law not far from $Ek^{1/2}$ (and steeper than $Ek^{1/3}$, not shown). Figure 5(b) shows an alternative measure of the boundary layer thickness, based on the root-mean-square vertical profile of $\partial_z\omega_z$. The inset shows sample profiles at $\tilde{Ra} = 120$ for different Ekman numbers. A local maximum is seen in these profiles, which we use to define a momentum boundary layer thickness $\delta_{\omega_z}$. Since we are considering stress-free boundaries, where no linear Ekman layer is present, the existence of this momentum boundary layer is nontrivial. We interpret $\delta_{\omega_z}$ as a measure of how far the effect of the stress-free boundary penetrates into the bulk of the flow.

We have furthermore verified that the momentum boundary layer seen here in the RRRiNSE formulation is also present in the reduced NHQG equations (not shown), although to our knowledge this has not been described previously. Like $\delta_\theta$, at sufficiently small $Ek$, $\delta_{\omega_z}$ becomes constant, while at larger values of $Ek$, there is again an undershoot followed by a power-law scaling regime where $\delta_{\omega_z}$ increases approximately as $Ek^{1/2}$. The approximate $Ek^{1/2}$ scaling, which is observed both in terms of $\delta_\theta$ and in $\delta_{\omega_z}$ and is also observed for no-slip boundary conditions (Song et al. 2024a), is surprising in the present case, given that we consider stress-free boundaries, which do not lead to the formation of linear Ekman layers. While a boundary layer thickness proportional to $Ek^{1/2}$ is theoretically predicted for the transitional thermal boundary layer (Julien et al. 2012b), the presence of such a scaling at significantly larger Ekman numbers, where the boundary layer flow is no longer rotationally dominated, goes beyond this prediction and the precise origin of this scaling law remains to be elucidated.

### 3.5. Alternative cuts through parameter space

The results presented thus far were obtained along a particular cut through the physical control parameter space, namely varying $Ek$ at fixed $\tilde{Ra}$. This is the natural approach in the framework of the RRRiNSE. By contrast, keeping $\tilde{Ra}$ fixed in an experimental setup or in an unrescaled DNS requires tuning two parameters. In such contexts, it is instead more natural to vary $\tilde{Ra}$ at fixed $Ek$ (equivalent to varying $\tilde{Ra}$ at fixed $Ek$) or vary $Ek$ at fixed $\tilde{Ra}$ (equivalent to varying $\tilde{Ra}$ at fixed $\tilde{Ra}$). To make contact with existing experimental and numerical studies, the corresponding cuts through parameter space are shown in figure 6. Panel a shows $Nu - 1$ versus $\tilde{Ra}$ for different sets of RRRiNSE simulations, each with $\tilde{Ra}$ fixed to a different value between $\tilde{Ra} = 10^6$ and $\tilde{Ra} = 10^{10}$, as well as the corresponding values from the NHQG equations ($\tilde{Ra} \to \infty$) at the same $\tilde{Ra}$. The NHQG data approximately follow the turbulent $\tilde{Ra}^{3/2}$ scaling law, cf. figure 4, while at finite $\tilde{Ra}$, there is a well-defined transitional value $\tilde{Ra}_t$ of $\tilde{Ra}$ where $Nu - 1$ intersects the $\tilde{Ra} \to \infty$ limit, and subsequently
<figure><img src="image_12.png" /><figcaption>Figure 5: Boundary layer structure for a range of Ek and $\tilde{Ra}$. Panel a: log-log plot of the non-dimensional widths of the thermal boundary layer $\delta_\theta/H$, given by the location in z of the maxima in the vertical root-mean-square (r.m.s.) profiles of $\theta$ shown in the inset, versus Ek for different $\tilde{Ra}$. Dashed line indicates a $Ek^{1/2}$ power law. Panel b: log-log plot of the momentum boundary layer thickness $\delta_{\omega_z}/H$ based on the location of the maxima of the r.m.s. profile of $\partial_z\omega_z$ shown in the inset. Panel c: log-log plot of the non-dimensional thermal boundary layer thickness $\delta_\theta/H$ at $Ek = 10^{-15}$ compared with the power law $\tilde{Ra}^{-15/8}$ predicted by the NHQG equations (Julien et al. 2012b).</figcaption></figure>

follows a shallower slope. This structure of the Nusselt number flattening out beyond a threshold $\tilde{Ra} = \tilde{Ra}_t$ is robust among the datasets of runs with different $Ra$. Qualitatively similar behavior was recently reported for RRRBC with no-slip boundaries (Song et al. 2024a), although the eventual low $Ek$ regime remains elusive in that case. Panel b shows the corresponding Reynolds number data for the same simulations. As discussed earlier, the corresponding turbulent scaling law $Re \propto \tilde{Ra}$ has proved to be more elusive owing to the impact of the inverse cascade, see Maffei et al. (2021), Oliver et al. (2023) and figure 4. Panel b shows a similar structure to panel a, in that the $Re$ curves intersect the curve corresponding to the NHQG equations at a (numerically slightly different) threshold value for each $\tilde{Ra}$. Since these runs are performed at constant $Ra$, each $\tilde{Ra}$ is equivalent to an Ekman number $Ek$. Therefore $Ra_t$ is equivalent to a transitional Ekman number $Ek_t$.

Panel c shows $Ek_t$, corresponding to the vertical dashed lines in panels a and b, indicating that the data are quantitatively compatible with the scaling law $Ek_t \propto \tilde{Ra}^{-5/8}$. This scaling law corresponds to the theoretical prediction (Julien et al. 2012a) for the Ekman number at which the boundary layer loses rotational support. Here, we recover this scaling from a simple measurement of the Nusselt and Reynolds numbers. An equivalent representation of
<figure><img src="image_13.png" /><figcaption>Figure 6: Nusselt number (panel a) and Reynolds number (panel b) versus the reduced Rayleigh number $\tilde{Ra}$ at fixed bare Rayleigh number $Ra$. Dashed lines indicate intersection of a given $Ra = const.$ set with the asymptotic limit. Transitional Ekman number (panel c) and reduced Rayleigh number (panel d) versus the bare Rayleigh number $Ra$. $Nu - 1$ (panel e) and $Re$ (panel f) versus $\tilde{Ra}$ at fixed Ekman number.</figcaption></figure>

the same data is shown in panel *d* in terms of the reduced Rayleigh number $\tilde{Ra}$ versus the bare Rayleigh number $Ra$, which also reveals satisfactory agreement with the corresponding theoretical scaling prediction $\tilde{Ra}_t \propto Ra^{1/6}$.

Panels *e* and *f* of figure 6 show a third possible cut through the parameter space of rotating convection (qualitatively resembling results presented in Song *et al.* (2024a)), holding *Ek* fixed and changing *Ra* (or equivalently changing $\tilde{Ra}$). The turbulent scaling laws with $\tilde{Ra}$ compare similarly well with these datasets and exhibit similar changes in scaling in the Nusselt and Reynolds numbers with $\tilde{Ra}$ at finite *Ek*, up to $Ek = 10^{-6}$, while for $Ek = 10^{-9}$
and $Ek = 10^{-10}$ the values of $\tilde{Ra}$ attainable in our simulations were not sufficiently high to observe a departure from the curve corresponding to the NHQG equations.

The flattening of the Nusselt number curves with increasing $\tilde{Ra}$ at finite Ekman numbers is a robust signal that could potentially be reproduced in laboratory experiments. Our results in figure 6 qualitatively resemble similar findings recently obtained for RRRBC with no-slip boundaries (Song *et al.* 2024a).

# 4. Discussion

The results described above show that the RRRiNSE formulation, based on low Ekman number asymptotics, allows efficient direct numerical simulations (DNS) far beyond the current state of art in laboratory experiments or unrescaled DNS, down to $Ek = 10^{-15}$ and below, which is comparable to the estimated values of $Ek$ in the outer core of the Earth and the convection zone of the Sun. Specifically, this indicates that the RRRiNSE enable DNS at Ekman numbers which are over six orders of magnitude smaller than the smallest value obtained previously, in an impressive effort and at great computational cost, by means of unrescaled DNS (Song *et al.* 2024a,b). Using the RRRiNSE formulation, we revealed rich new physics in this previously inaccessible parameter regime. First, we uncovered a transition at $Ek \approx \sigma^{3/2} \tilde{Ra}^{-15/4}$ ($Ek \approx 10^{-9}$ for $\tilde{Ra} = 120$) from a regime at larger $Ek$ where the depth-averaged flow features a strong large-scale cyclone and a weak diffuse anticyclone towards a regime characterised by cyclone-anticyclone symmetry and a large-scale vortex dipole, in agreement with the predictions of the NHQG equations. Second, we identified a nontrivial transition in the boundary layer dynamics corresponding to the loss of rotational support in the boundary layer (confirming a previously untestable theoretical prediction of Julien *et al.* (2012b)), a transition accompanied by the emergence of strong, albeit short-lived, anticyclonic structures near the $Ek$ threshold that weaken as $Ek$ increases, leading to dominance of strong cyclonic structures. It is interesting to note that the cyclone-anticyclone symmetry-breaking in the bulk approximately coincides in $Ek$ with the boundary layer transition, although a theoretical explanation of this observation remains unavailable. While flows in the Earth's atmosphere are characterised by only moderately small Rossby numbers and are therefore dominated by strong cyclones, our results suggest that in the interior of the Earth or other celestial bodies, cyclones and anticyclones may be statistically of equal strength.

We showed quantitatively that the time-averaged Nusselt and Reynolds numbers in steady state reflect these regime transitions, taking values consistent with the NHQG limit for $Ek \lesssim \sigma^{3/2} \tilde{Ra}^{-15/4}$, while overshooting as $Ek$ is increased at fixed $\tilde{Ra}$, reaching maximum heat transport near the value of $Ek$ where the bulk Rossby number reaches unity, and decreasing to close to zero as $Ek$ increases further at fixed $\tilde{Ra}$ due to reduced supercriticality. The increase in the Nusselt number was shown to be tied to increased dissipation associated with a change in the boundary layer flow morphology. The boundary depth was quantified in terms of the thermal fluctuations and vertical vorticity, revealing that, at $Ek \lesssim \sigma^{3/2} \tilde{Ra}^{-15/4}$, the boundary layer depth becomes $Ek$-independent, taking the value associated with the NHQG equations, but undershooting for $Ek \gtrsim \sigma^{3/2} \tilde{Ra}^{-15/4}$ before increasing approximately with an Ekman layer-like scaling $Ek^{1/2}$ as $Ek$ increases. Owing to the presence of stress-free boundary conditions, this behavior cannot be explained by a linear Ekman layer and is therefore a nonlinear effect. Finally, we have considered alternative cuts through the parameter space, with one set of runs varying $Ek$ at fixed $\tilde{Ra}$ and another set of runs varying $\tilde{Ra}$ at fixed $Ek$. This procedure revealed that at finite Ekman numbers, the Nusselt and Reynolds
numbers remain close to the values found in the NHQG equations as $\tilde{Ra}$ increases, but start to deviate from them at a value of $\tilde{Ra}$ that depends on $Ra$ in a way consistent with boundary layer theory (Julien et al. 2012a). We mention that the NHQG equations apply to systems with no-slip boundaries since in the limit $Ek \to 0$ the no-slip boundaries become effectively stress-free (Julien et al. 2012b). At finite $Ek$, however, Ekman boundary layers and the associated Ekman pumping necessarily modify the boundary layer structure, an effect that can be included in the NHQG equations following Julien et al. (2016).

The results presented here are based on the RRRiNSE reformulation of the equations and explain why previous state-of-the art DNS and laboratory experiments could not reach the parameter regime for observing the transition to fully rotationally constrained dynamics. The RRRiNSE formulation therefore opens the door to the further exploration of the parameter regime of very small but finite Ekman and Rossby numbers, highly relevant to planetary, satellite and stellar interiors. This is all the more timely since current and future observational missions prompted by the last Planetary Decadal Survey (National Academies of Sciences, Engineering, and Medicine 2023), including the recently launched ESA's JUICE (JUPiter Icy Moons Explorer) mission (Grasset et al. 2013) and NASA's Europa Clipper mission (Howell & Pappalardo 2020) will provide new data requiring interpretation based on faithful model simulations at parameters as close as possible to realistic values. Our results are a first step in this direction.

The RRRBC as studied here, with antiparallel gravity and rotation axis is an appropriate, albeit idealised, model of the North Pole regions of the celestial bodies listed in table 1. Owing to the latitudinal non-uniformity of rotating convection in spherical shells (Gastine & Aurnou 2023), the RRRiNSE formulation cannot be straightforwardly generalised to that setting, but an investigation of the RRRiNSE with misaligned gravity and rotation axes could provide important insights in this regard. While RRRiNSE enables simulations at realistic values of $Ek$, many challenges remain: the increasing resolution constraint for increasing levels of turbulence cannot be circumvented, and for simplicity our simulations considered stress-free boundaries as well as Prandtl number $\sigma = 1$ while ignoring other important ingredients such as magnetic fields. To further test the applicability of our findings to planets, icy moons and stars, future work must also investigate how the regime transitions described in this study are influenced by additional effects, including fluids with Prandtl numbers different from one, compressibility, internal heating (Barker et al. 2014; Bouillaut et al. 2021; Hadjerci et al. 2024), and strong magnetic fields, as well as alternative choices of boundary conditions. The doors to these explorations are now open thanks to the RRRiNSE formulation. Comparison of the RRRiNSE numerical results with future laboratory experiments would also be highly desirable.

**Acknowledgments.** K. J. passed away before this manuscript was finalised. We have attempted to present the results of our collaboration in accordance with his high standards. Any errors or misinterpretations remain our own. Part of the writing of this manuscript was done while one of the authors (AvK) was a Staff Member at the Woods Hole Oceanographic Institution Geophysical Fluid Dynamics summer program. This research used the Savio computational cluster resource provided by the Berkeley Research Computing program at the University of California Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information Officer). This research also utilised the Alpine high performance computing resource at the University of Colorado Boulder. Alpine is jointly funded by the University of Colorado Boulder, the University of Colorado Anschutz, and Colorado State University. Data storage for this project was supported by the University of Colorado Boulder PetaLibrary. This project was also granted access to computational resources of TGCC under the allocation 2024-A0162A10803 made by GENCI, and to resources of PMCS2I (Pôle de Modélisation et de Calcul en Sciences de l'Ingénieur et de l'Information) of Ecole Centrale de Lyon. We thank Daria Holdenried-Chernoff for pointing out useful references on dynamo simulations at low Ekman number and Rudie Kunnen for helpful comments. This research was funded, in whole or in part, by Agence Nationale de la Recherche (Grant ANR-23-CE30-0016-01). A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to
the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access
conditions.

**Funding.** This work was supported by the National Science Foundation (Grants DMS-2009319 and DMS-2308338 (KJ), Grants DMS-2009563 and DMS-2308337 (EK)), by the German Research Foundation (DFG Projektnummer: 522026592) and by Agence Nationale de la Recherche (Grant ANR-23-CE30-0016-01).

**Supplementary movies.** Supplementary movies are available at https://doi.org/10.1017/jfm.2025.290.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.**

A. van Kan, https://orcid.org/0000-0002-1217-3609;

K. Julien, https://orcid.org/0000-0002-4409-7022;

B. Miquel, https://orcid.org/0000-0001-6283-0382;

E. Knobloch, https://orcid.org/0000-0002-1567-9314.

**Appendix A. Additional details on flow statistics and structure**

Supplementary Movies 1 to 3 show the evolution of the vertical vorticity $\omega_z$, the temperature perturbation $\theta$ and the vertical velocity $w$, respectively, at the top of the momentum boundary layer $z = \delta\omega_z$ for $Ek = 10^{-15}$ and $\tilde{Ra} = 120$. Supplementary Movies 4 to 6 provide the same information for $Ek = 10^{-8}$. At $Ek = 10^{-15}$, all fields display close to zero skewness and cyclones are of approximately equal strength and structure as anticyclones. In contrast, at $Ek = 10^{-8}$, the boundary layer is characterised by the presence of strong, albeit short-lived, anticyclonic vortical structures typically with a shielded structure that impose a clear signature on $\theta$ and $w$, associated with cold spots in the centre where the fluid descends while rising within a surrounding ring. Figure 7 shows an analysis of the dissipation budget predicted by the power-integral relations given in (3.2) in the main text, as a function of the Taylor number $Ta \equiv Ek^{-2}$. Panel a shows that in the statistically stationary state at $\tilde{Ra} = 120$ the sum of the dissipation contributions from the temperature variations due to gradients in the horizontal and vertical directions sums up approximately to $Nu - 1$ as predicted by (3.2a). The departure from the approximately constant Nusselt number at small $Ek$ is seen to be associated with an increase in dissipation due to both horizontal and vertical gradients. Panel b shows the contributions from vertical gradients, which primarily stem from the boundary layers near the departure from the NHQG limit. Panel c shows the kinetic energy dissipation budget as a function of the Taylor number at $\tilde{Ra} = 120$. The sum of the dissipation of kinetic energy due to vertical and horizontal gradients shows a satisfactory overall agreement with $(Nu - 1)\tilde{Ra}/\sigma^2$, as predicted by (3.2b). However, while for $Ta \lesssim 10^8$, where no inverse energy cascade is observed, the agreement is close to perfect, it is notable that for $Ta \gtrsim 10^{10}$, the measured kinetic energy is slightly below the value predicted by (3.2b). This is because the large-scale vortex is still slowly growing in amplitude via the condensation process. Panel d shows the contributions from vertical gradients, which show a similarly drastic transition from negligibly small to finite values, contributing to the observed overshoot in the Nusselt number.

The RRRiNSE formulation remains numerically stable at even lower Ekman numbers than those discussed in the main text. Figure 8 illustrates this in terms of the flow field obtained for $\tilde{Ra} = 120$ at $Ek = 10^{-24}$. The flow morphology is indistinguishable from that at $Ek = 10^{-15}$. Specifically, in the boundary layer, there is no cyclone-anticyclone asymmetry and the $\omega_z$ and $\theta$ fields are nearly perfectly correlated. Moreover, the barotropic vertical vorticity displays a counter-rotating vortex dipole at the domain scale, similar to what is seen in two-dimensional turbulence. While the observed flow structure and the associated flow statistics at $Ek = 10^{-24}$ are close to $Ek = 10^{-15}$ for the case considered here, more analysis
<figure><img src="image_14.png" /><figcaption>Figure 7: Panel a shows the dissipation budget of the temperature variance arising from vertical flow variations at $\tilde{R}a = 120$ versus $Ta = Ek^{-2}$, with $\epsilon_{\theta}^{\perp} = \langle |\nabla_{\perp}\theta|^2 \rangle$, $\epsilon_{\theta}^{\parallel} = \epsilon^2 \langle (\partial_z\theta)^2 \rangle$ and $\epsilon_{\Theta}^{\parallel} = \langle (\partial_z\Theta)^2 \rangle$, showing that the power-integral equation (3.2a) of the main text is well satisfied and hence that a statistically steady state has been reached. The overshoot in $Nu - 1$ is primarily due to increased $\epsilon_{\theta}^{\perp}$ and $\epsilon_{\theta}^{\parallel}$. Panel b shows the dissipation due to vertical variations in $\theta$ at $\tilde{R}a = 80, 120$ versus $Ta = Ek^{-2}$ together with the boundary layer contributions (marked by $\times$, integrated over a depth of $2\delta_{\theta}$). Panel c shows the kinetic energy budget showing that the contributions from the horizontal $(\epsilon_{kin}^{\perp} \equiv \langle |\nabla_{\perp}u|^2 \rangle)$ and vertical gradients $(\epsilon_{kin}^{\parallel} \equiv \epsilon^2 \langle |\partial_zu|^2 \rangle)$ approximately add up to $(Nu - 1)\tilde{R}a/\sigma^2$, as predicted by (3.2b) of the main text. A small mismatch is seen at $Ek \lesssim 10^{-5}$ since the emerging large-scale vortex has not fully saturated in amplitude. Panel d shows the contribution from vertical variation to the kinetic energy dissipation versus $Ta$ at $\tilde{R}a = 80, 120$.</figcaption></figure>

is needed to determine the lowest Ekman numbers attainable with RRRiNSE given finite machine precision, and to characterize the breakdown of the method.

At a given Ekman number, the importance of the Coriolis force compared to nonlinear acceleration is not uniform in the vertical direction. To illustrate this, figure 9 shows the vertical r.m.s. profiles of the Rossby number, computed based on the ratio $|u \cdot \nabla u_{\perp}| / (2\Omega|u_{\perp}|)$. This quantity has a well-defined value in the bulk of the flow, but assumes a distinct value
<figure><img src="image_15.png" /><figcaption>Figure 8: Snapshots of the flow state obtained via direct numerical simulation of the RRRiNSE at $Ek = 10^{-24}$ and $Ra = 120$ in terms of $\omega_z$, $w$, $\theta$ at the top of the boundary layer $z = \delta\omega_z$, and the depth-averaged vorticity field $\overline{\omega_z}$ revealing a symmetric large-scale vortex dipole.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 9: Vertical profiles of the r.m.s. local Rossby number $|u \cdot \nabla u|/(2\Omega|u_\perp|)$ versus non-dimensional height at $Ra = 120$ for Ekman numbers $Ek = 10^{-15}, 10^{-12}, 10^{-9}$. Inset shows same data in log-log representation, showing that bulk and boundary layer have distinct values of the Rossby number. At $Ek = 10^{-9}$, there is a local maximum in the Rossby number observed near the top of the thermal boundary layer. The values shown in the main text are obtained by averaging over the thermal boundary layer volume.</figcaption></figure>
in the boundary layer, as can be seen in the doubly logarithmic representation in the inset. These two distinct values characterising the bulk flow and the boundary layer are shown in panels *e* and *f* of figure 4 in the main text.

## REFERENCES

ALEXAKIS, A. 2023 Quasi-two-dimensional turbulence. *Rev. Mod. Plasma Phys.* **7**, 31.

ALEXAKIS, A. & BIFERALE, L. 2018 Cascades and transitions in turbulent flows. *Phys. Rep.* **767**, 1–101.

ASCHER, U. M., RUUTH, S. J. & SPITERI, R. J. 1997 Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. *Appl. Numer. Math.* **25**, 151–167, Special Issue on Time Integration.

AURNOU, J. M., CALKINS, M. A., CHENG, J. S., JULIEN, K., KING, E. M., NIEVES, D., SODERLUND, K. M. & STELLMACH, S. 2015 Rotating convective turbulence in Earth and planetary cores. *Phys. Earth Planet. Inter.* **246**, 52–71.

BARKER, A. J., DEMPSEY, A. M. & LITHWICK, Y. 2014 Theory and simulations of rotating convection. *Astrophys. J.* **791**, 13.

BOFFETTA, G. & ECKE, R. E. 2012 Two-dimensional turbulence. *Annu. Rev. Fluid. Mech.* **44**, 427–451.

BOUBNOV, B. M. & GOLITSYN, G. S. 1986 Experimental study of convective structures in rotating fluids. *J. Fluid Mech.* **167**, 503–531.

BOUBNOV, B. M. & GOLITSYN, G. S. 2012 *Convection in Rotating Fluids*. Springer Science & Business Media.

BOUILLAUT, V., MIQUEL, B., JULIEN, K., AUMAÎTRE, S. & GALLET, B. 2021 Experimental observation of the geostrophic turbulence regime of rapidly rotating convection. *Proc. Natl. Acad. Sci. U.S.A.* **118**, e2105015118.

CAI, T. 2020 Penetrative convection for rotating Boussinesq flow in tilted *f*-planes. *Astrophys. J.* **898**, 22.

CHANDRASEKHAR, S. 1953 The instability of a layer of fluid heated below and subject to Coriolis forces. *Proc. R. Soc. A* **217**, 306–327.

CHENG, J. S., AURNOU, J. M., JULIEN, K. & KUNNEN, R. P. J. 2018 A heuristic framework for next-generation models of geostrophic convective turbulence. *Geophys. Astrophys. Fluid Dyn.* **112**, 277–300.

CHENG, J. S., MADONIA, M., AGUIRRE GUZMÁN, A. J. & KUNNEN, R. P. J. 2020 Laboratory exploration of heat transfer regimes in rapidly rotating turbulent convection. *Phys. Rev. Fluids* **5**, 113501.

CHEON, W. G. & GORDON, A. L. 2019 Open-ocean polynyas and deep convection in the Southern Ocean. *Sci. Rep.* **9**, 6935.

CLENSHAW, C. 1957 The numerical solution of linear differential equations in Chebyshev series. *Math. Proc. Camb. Philos. Soc.* **53**, 134–149.

COOPER, R. G., BUSHBY, P. J. & GUERVILLY, C. 2020 Subcritical dynamos in rapidly rotating planar convection. *Phys. Rev. Fluids* **5**, 113702.

CURRIE, L. K., BARKER, A. J., LITHWICK, Y. & BROWNING, M. K. 2020 Convection with misaligned gravity and rotation: simulations and rotating mixing length theory. *Mon. Not. R. Astron. Soc.* **493**, 5233–5256.

DAUXOIS, T., PEACOCK, T., BAUER, P., CAULFIELD, C. P., CENEDESE, C., GORLÉ, C., HALLER, G., IVEY, G. N., LINDEN, P. F., MEIBURG, E. & OTHERS 2021 Confronting grand challenges in environmental fluid mechanics. *Phys. Rev. Fluids* **6**, 020501.

ECKE, R. E. & SHISHKINA, O. 2023 Turbulent rotating Rayleigh-Bénard convection. *Annu. Rev. Fluid Mech.* **55**, 603–638.

EMANUEL, K. A. 1994 *Atmospheric Convection*. Oxford University Press, USA.

FAN, Y. 2021 Magnetic fields in the solar convection zone. *Living Rev. Sol. Phys.* **18**, 5.

FAVIER, B., SILVERS, L. J. & PROCTOR, M. R. E. 2014 Inverse cascade and symmetry breaking in rapidly rotating Boussinesq convection. *Phys. Fluids* **26**, 096605.

GALLET, B., CAMPAGNE, A., CORTET, P.-P. & MOISY, F. 2014 Scale-dependent cyclone-anticyclone asymmetry in a forced rotating turbulence experiment. *Phys. Fluids* **26**, 035108.

GARAUD, P. & GARAUD, J.-D. 2008 Dynamics of the solar tachocline–II. The stratified case. *Mon. Not. R. Astron. Soc.* **391**, 1239–1258.

GASTINE, T. & AURNOU, J. M. 2023 Latitudinal regionalization of rotating spherical shell convection. *J. Fluid Mech.* **954**, R1.

GRASSET, O., DOUGHERTY, M., COUSTENIS, A., BUNCE, E., ERD, C., TITOV, D., BLANC, M., COATES, A., DROSSART, P., FLETCHER, L., HUSSMANN, H., JAUMANN, R., KRUPP, N., LEBRETON, J.-P., PRIETO-
BALLESTEROS, O., TORTORA, P., TOSI, F. & VAN HOOLST, T. 2013 JUpiter ICy moons Explorer (JUICE): An ESA mission to orbit Ganymede and to characterise the Jupiter system. *Planet. Space Sci.* **78**, 1–21.

GREENGARD, L. 1991 Spectral integration and two-point boundary value problems. *SIAM J. Numer. Anal.* **28**, 1071–1080.

GUERVILLY, C., CARDIN, P. & SCHAEFFER, N. 2019 Turbulent convective length scale in planetary cores. *Nature* **570**, 368–371.

GUERVILLY, C. & HUGHES, D. W. 2017 Jets and large-scale vortices in rotating Rayleigh-Bénard convection. *Phys. Rev. Fluids* **2**, 113503.

GUERVILLY, C., HUGHES, D. W. & JONES, C. A. 2014 Large-scale vortices in rapidly rotating Rayleigh-Bénard convection. *J. Fluid Mech.* **758**, 407–435.

HADJERCI, G., BOUILLAUT, V., MIQUEL, B. & GALLET, B. 2024 Rapidly rotating radiatively driven convection: experimental and numerical validation of the ‘geostrophic turbulence’ scaling predictions. *J. Fluid Mech.* **998**, A9.

HART, J. E., KITTELMAN, S. & OHLSEN, D. R. 2002 Mean flow precession and temperature probability density functions in turbulent rotating convection. *Phys. Fluids* **14**, 955–962.

HATHAWAY, D. H. & SOMERVILLE, R. C. J. 1983 Three-dimensional simulations of convection in layers with tilted rotation vectors. *J. Fluid Mech.* **126**, 75–89.

HAWKINS, E. K., CHENG, J. S., ABBATE, J. A., PILEGARD, T., STELLMACH, S., JULIEN, K. & AURNOU, J. M. 2023 Laboratory models of planetary core-style convective turbulence. *Fluids* **8**, 106.

HE, J., FAVIER, B., RIEUTORD, M. & LE DIZÈS, S. 2022 Internal shear layers in librating spherical shells: the case of periodic characteristic paths. *J. Fluid Mech.* **939**, A3.

HOWELL, S. M. & PAPPALARDO, R. T. 2020 NASA’s Europa Clipper—a mission to a potentially habitable ocean world. *Nat. Comm.* **11**, 1311.

JONES, C. A. 2011 Planetary magnetic fields and fluid dynamos. *Annu. Rev. Fluid Mech.* **43**, 583–614.

JULIEN, K., AURNOU, J. M., CALKINS, M. A., KNOBLOCH, E., MARTI, P., STELLMACH, S. & VASIL, G. M. 2016 A nonlinear model for rotationally constrained convection with Ekman pumping. *J. Fluid Mech.* **798**, 50–87.

JULIEN, K., ELLISON, A., MIQUEL, B., CALKINS, M. & KNOBLOCH, E. 2025a Quasi-geostrophic convection on the tilted *f*-plane. (to be submitted to *J. Fluid Mech.*).

JULIEN, K., VAN KAN, A., MIQUEL, B., KNOBLOCH, E. & VASIL, G. 2025b Rescaled equations for well-conditioned direct numerical simulations of rapidly rotating convection. *J. Comput. Phys.* **541**, 114274.

JULIEN, K. & KNOBLOCH, E. 1998 Strongly nonlinear convection cells in a rapidly rotating fluid layer: the tilted *f*-plane. *J. Fluid Mech.* **360**, 141–178.

JULIEN, K. & KNOBLOCH, E. 2007 Reduced models for fluid flows with strong constraints. *J. Math. Phys.* **48**, 065405.

JULIEN, K., KNOBLOCH, E., MILLIFF, R. & WERNE, J. 2006 Generalized quasi-geostrophy for spatially anisotropic rotationally constrained flows. *J. Fluid Mech.* **555**, 233–274.

JULIEN, K., KNOBLOCH, E. & PLUMLEY, M. 2018 Impact of domain anisotropy on the inverse cascade in geostrophic turbulent convection. *J. Fluid Mech.* **837**, R4.

JULIEN, K., KNOBLOCH, E., RUBIO, A. M. & VASIL, G. M. 2012a Heat transport in low-Rossby-number Rayleigh-Bénard convection. *Phys. Rev. Lett.* **109**, 254503.

JULIEN, K., KNOBLOCH, E. & WERNE, J. 1998 A new class of equations for rotationally constrained flows. *Theoret. Comput. Fluid Dyn.* **11**, 251–261.

JULIEN, K., LEGG, S., MCWILLIAMS, J. & WERNE, J. 1996 Rapidly rotating turbulent Rayleigh-Bénard convection. *J. Fluid Mech.* **322**, 243–273.

JULIEN, K., RUBIO, A. M., GROOMS, I. & KNOBLOCH, E. 2012b Statistical and physical balances in low Rossby number Rayleigh-Bénard convection. *Geophys. Astrophys. Fluid Dyn.* **106**, 392–428.

JULIEN, K. & WATSON, M. 2009 Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods. *J. Comp. Phys.* **228**, 1480–1503.

VAN KAN, A. 2024 Phase transitions in anisotropic turbulence. *Chaos* **34**, 122103.

KING, E. M., SODERLUND, K. M., CHRISTENSEN, U. R., WICHT, J. & AURNOU, J. M. 2010 Convective heat transfer in planetary dynamo models. *Geochem. Geophys. Geosyst.* **11**, Q06016.

KING, E. M., STELLMACH, S. & AURNOU, J. M. 2012 Heat transfer by rapidly rotating Rayleigh-Bénard convection. *J. Fluid Mech.* **691**, 568–582.

KNOBLOCH, E. 1998 Rotating convection: recent developments. *Int. J. Eng. Sci.* **36**, 1421–1450.
KOLHEY, P., STELLMACH, S. & HEYNER, D. 2022 Influence of boundary conditions on rapidly rotating convection and its dynamo action in a plane fluid layer. *Phys. Rev. Fluids* **7**, 043502.

KUNNEN, R. P. J., OSTILLA-MÓNICO, R., VAN DER POEL, E. P., VERZICCO, R. & LOHSE, D. 2016 Transition to geostrophic convection: the role of the boundary conditions. *J. Fluid Mech.* **799**, 413–432.

LESUR, G. & OGILVIE, G. I. 2010 On the angular momentum transport due to vertical convection in accretion discs. *Mon. Not. R. Astron. Soc.* **404**, L64–L68.

LIN, Y. & JACKSON, A. 2021 Large-scale vortices and zonal flows in spherical rotating convection. *J. Fluid Mech.* **912**, A46.

LU, H.-Y., DING, G.-Y., SHI, J.-Q., XIA, K.-Q. & ZHONG, J.-Q. 2021 Heat-transport scaling and transition in geostrophic rotating convection with varying aspect ratio. *Phys. Rev. Fluids* **6**, L071501.

MAFFEI, S., KROUSS, M. J., JULIEN, K. & CALKINS, M. A. 2021 On the inverse cascade and flow speed scaling behaviour in rapidly rotating Rayleigh-Bénard convection. *J. Fluid Mech.* **913**, A18.

MAJUMDER, D., SREENIVASAN, B. & MAURYA, G. 2024 Self-similarity of the dipole-multipole transition in rapidly rotating dynamos. *J. Fluid Mech.* **980**, A30.

MARSHALL, J. & SCHOTT, F. 1999 Open-ocean convection: Observations, theory, and models. *Rev. Geophys.* **37**, 1–64.

MASON, S. J., GUERVILLY, C. & SARSON, G. R. 2022 Magnetoconvection in a rotating spherical shell in the presence of a uniform axial magnetic field. *Geophys. Astrophys. Fluid Dyn.* **116**, 458–498.

MIESCH, M. S. 2000 The coupling of solar convection and rotation. *Sol. Phys.* **192**, 59–89.

MIQUEL, B. 2021 Coral: A parallel spectral solver for fluid dynamics and partial differential equations. *J. Open Source Softw.* **6**, 2978.

MIQUEL, B., XIE, J.-H., FEATHERSTONE, N., JULIEN, K. & KNOBLOCH, E. 2018 Equatorially trapped convection in a rapidly rotating shallow shell. *Phys. Rev. Fluids* **3**, 053801.

MITRI, G. & SHOWMAN, A. P. 2008 Thermal convection in ice-I shells of Titan and Enceladus. *Icarus* **193**, 387–396.

NAKAGAWA, Y. & FRENZEN, P. 1955 A theoretical and experimental study of cellular convection in rotating fluids. *Tellus* **7**, 2–21.

NATIONAL ACADEMIES OF SCIENCES, ENGINEERING, AND MEDICINE 2023 *Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023-2032*. Washington, DC: The National Academies Press.

NIMMO, F. & PAPPALARDO, R. T. 2016 Ocean worlds in the outer solar system. *J. Geophys. Res. E Planets* **121**, 1378–1399.

NOVI, L., VON HARDENBERG, J., HUGHES, D. W., PROVENZALE, A. & SPIEGEL, E. A. 2019 Rapidly rotating Rayleigh-Bénard convection with a tilted axis. *Phys. Rev. E* **99**, 053116.

OLIVER, T. G., JACOBI, A. S., JULIEN, K. & CALKINS, M. A. 2023 Small scale quasigeostrophic convective turbulence at large Rayleigh number. *Phys. Rev. Fluids* **8**, 093502.

PAPPALARDO, R. T., HEAD, J. W., GREELEY, R., SULLIVAN, R. J., PILCHER, C., SCHUBERT, G., MOORE, W. B., CARR, M. H., MOORE, J. M., BELTON, M. J. S. & OTHERS 1998 Geological evidence for solid-state convection in Europa's ice shell. *Nature* **391**, 365–368.

POTHÉRAT, A. & HORN, S. 2025 Seven decades of exploring planetary interiors with rotating convection experiments. *Comptes Rendus. Physique* **25**.

PROUDMAN, J. 1916 On the motion of solids in a liquid possessing vorticity. *Proc. R. Soc. A* **92**, 408–424.

RAJAEI, H., KUNNEN, R. P. J. & CLERCX, H. J. H. 2017 Exploring the geostrophic regime of rapidly rotating convection with experiments. *Phys. Fluids* **29**, 045105.

ROSSBY, H. T. 1969 A study of Bénard convection with and without rotation. *J. Fluid Mech.* **36**, 309–335.

ROULLET, G. & KLEIN, P. 2010 Cyclone-anticyclone asymmetry in geophysical turbulence. *Phys. Rev. Lett.* **104**, 218501.

RUBIO, A. M., JULIEN, K., KNOBLOCH, E. & WEISS, J. B. 2014 Upscale energy transfer in three-dimensional rapidly rotating turbulent convection. *Phys. Rev. Lett.* **112**, 144501.

SCHAEFFER, N., JAULT, D., NATAF, H.-C. & FOURNIER, A. 2017 Turbulent geodynamo simulations: a leap towards Earth's core. *Geophys. J. Int.* **211**, 1–29.

SCHUBERT, G. & SODERLUND, K. M. 2011 Planetary magnetic fields: Observations and models. *Phys. Earth Planet. Inter.* **187**, 92–108.

SEVERIN, T., CONAN, P., DE MADRON, X. D., HOUPERT, L., OLIVER, M. J., ORIOL, L., CAPARROS, J., GHIGLIONE, J. F. & PUJO-PAY, M. 2014 Impact of open-ocean convection on nutrients, phytoplankton biomass and activity. *Deep Sea Res. Part I: Ocean Res.* **94**, 62–71.
SHEW, W. L. & LATHROP, D. P. 2005 Liquid sodium model of geophysical core convection. *Phys. Earth Planet. Inter.* **153**, 136–149.

SIEGELMAN, L., KLEIN, P., INGERSOLL, A. P., EWALD, S. P., YOUNG, W. R., BRACCO, A., MURA, A., ADRIANI, A., GRASSI, D., PLAINAKI, C. & OTHERS 2022 Moist convection drives an upscale energy transfer at Jovian high latitudes. *Nat. Phys.* **18**, 357–361.

SODERLUND, K. M. 2019 Ocean dynamics of outer solar system satellites. *Geophys. Res. Lett.* **46**, 8700–8710.

SODERLUND, K. M., ROVIRA-NAVARRO, M., LE BARS, M., SCHMIDT, B. E. & GERKEMA, T. 2024 The physical oceanography of ice-covered moons. *Annu. Rev. Mar. Sci.* **16**, 25–53.

SONG, J., SHISHKINA, O. & ZHU, X. 2024a Direct numerical simulations of rapidly rotating Rayleigh-Bénard convection with Rayleigh number up to $5 \times 10^{13}$. *J. Fluid Mech.* **989**, A3.

SONG, J., SHISHKINA, O. & ZHU, X. 2024b Scaling regimes in rapidly rotating thermal convection at extreme Rayleigh numbers. *J. Fluid Mech.* **984**, A45.

SPRAGUE, M., JULIEN, K., KNOBLOCH, E. & WERNE, J. 2006 Numerical simulation of an asymptotically reduced system for rotationally constrained convection. *J. Fluid Mech.* **551**, 141–174.

STELLMACH, S., LISCHPER, M., JULIEN, K., VASIL, G., CHENG, J. S., RIBEIRO, A., KING, E. M. & AURNOU, J. M. 2014 Approaching the asymptotic regime of rapidly rotating convection: Boundary layers versus interior dynamics. *Phys. Rev. Lett.* **113**, 254501.

STEVENS, R. J. A. M., CLERCX, H. J. H. & LOHSE, D. 2013 Heat transport and flow structure in rotating Rayleigh-Bénard convection. *Eur. J. Mech. B Fluids* **40**, 41–49.

VERONIS, G. 1959 Cellular convection with finite amplitude in a rotating fluid. *J. Fluid Mech.* **5**, 401–435.

VOROBIEFF, P. & ECKE, R. E. 2002 Turbulent rotating convection: an experimental study. *J. Fluid Mech.* **458**, 191–218.

ZHONG, F., ECKE, R. & STEINBERG, V. 1991 Asymmetric modes and the transition to vortex structures in rotating Rayleigh-Bénard convection. *Phys. Rev. Lett.* **67**, 2473–2476.