Boundary zonal flows in rapidly rotating turbulent thermal convection

Recently, in Zhang et al. (2020), it was found that in rapidly rotating turbulent Rayleigh-B\'enard convection (RBC) in slender cylindrical containers (with diameter-to-height aspect ratio $\Gamma=1/2$) filled with a small-Prandtl-number fluid ($Pr \approx0.8$), the Large Scale Circulation (LSC) is suppressed and a Boundary Zonal Flow (BZF) develops near the sidewall, characterized by a bimodal PDF of the temperature, cyclonic fluid motion, and anticyclonic drift of the flow pattern (with respect to the rotating frame). This BZF carries a disproportionate amount ($>60\%$) of the total heat transport for $Pr<1$ but decreases rather abruptly for larger $Pr$ to about $35\%$. In this work, we show that the BZF is robust and appears in rapidly rotating turbulent RBC in containers of different $\Gamma$ and in a broad range of $Pr$ and $Ra$. Direct numerical simulations for $0.1 \leq Pr \leq 12.3$, $10^7 \leq Ra \leq 5\times10^{9}$, $10^{5} \leq 1/Ek \leq 10^{7}$ and $\Gamma$ = 1/3, 1/2, 3/4, 1 and 2 show that the BZF width $\delta_0$ scales with the Rayleigh number $Ra$ and Ekman number $Ek$ as $\delta_0/H \sim \Gamma^{0} \Pr^{\{-1/4, 0\}} Ra^{1/4} Ek^{2/3}$ (${Pr<1, Pr>1}$) and the drift frequency as $\omega/\Omega \sim \Gamma^{0} Pr^{-4/3} Ra Ek^{5/3}$, where $H$ is the cell height and $\Omega$ the angular rotation rate. The mode number of the BZF is 1 for $\Gamma \lesssim 1$ and $2 \Gamma$ for $\Gamma$ = {1,2} independent of $Ra$ and $Pr$. The BZF is quite reminiscent of wall mode states in rotating convection.


Introduction
Turbulent convection driven by buoyancy and subject to background rotation is a phenomenon of great relevance in many physical disciplines, especially in geo-and astrophysics and also in engineering applications. In a model system of Rayleigh-Bénard convection (RBC) (Bodenschatz et al. 2000;Ahlers et al. 2009;Lohse & Xia 2010), a fluid is confined in a container where the bottom is heated, the top is cooled, and the vertical walls are adiabatic. The temperature inhomogeneity leads to a fluid density variation which, in the presence of gravity, produces convective fluid motion. When the system rotates with respect to the vertical axis significant modification of the flow occurs owing to the rotational influence including the suppression of the onset of convection (Chandrasekhar 1961;Nakagawa & Frenzen 1955), the enhancement or suppression of turbulent heat transport over different ranges of Ra and Pr (Rossby 1969;Pfotenhauer et al. 1987;Zhong et al. 1993;Julien et al. 1996;Liu & Ecke 1997), the transformation of thermal plumes into thermal vortices with a rich variety of local structure dynamics (Boubnov & Golitsyn 1986, 1990Hart et al. 2002;Vorobieff & Ecke 2002), and the emergence of robust wall modes before the onset of the bulk mode (Buell & Catton 1983;Pfotenhauer et al. 1987;Zhong et al. 1991;Ecke et al. 1992;Kuo & Cross 1993;Herrmann & Busse 1993;Goldstein et al. 1993).
The dimensionless control parameters in rotating RBC (RRBC) are the Rayleigh number Ra ≡ αg∆H 3 /(κν), the Prandtl number Pr ≡ ν/κ, the Ekman number Ek ≡ ν/(2ΩH 2 ), and the diameter-to-height aspect ratio of the container, Γ ≡ D/H. Here α denotes the isobaric thermal expansion coefficient, ν the kinematic viscosity, κ the thermal diffusivity of the fluid, g the acceleration due to gravity, Ω the angular rotation rate, ∆ ≡ T + − T − the difference between the temperatures at the bottom (T + ) and top (T − ) plates, H the distance between the isothermal plates (the cylinder height), and D ≡ 2R the cylinder diameter. The Rossby number Ro ≡ √ αg∆H / (2ΩH) = Ra/Pr Ek is another important non-dimensional parameter that provides a measure of the balance between buoyancy and rotation and is independent of dissipation coefficients.
The important global response parameter in thermal convection is the averaged total heat transport between the bottom and top plates, described by the Nusselt number, Nu ≡ ( u z T z − κ∂ z T z )/(κ∆/H). Here, T denotes the temperature, u is the velocity field with component u z in the vertical direction, and · z denotes the average in time and over a horizontal cross-section at height z from the bottom.
Rotation has various effects on the structure of the convective flow and on the global heat transport in the system. Rotation inhibits convection and causes an increase of the critical Ra c ∼ Ek −4/3 at which the quiescent fluid layer becomes unstable throughout the layer (Chandrasekhar 1961;Nakagawa & Frenzen 1955;Rossby 1969;Lucas et al. 1983;Zhong et al. 1993). In finite containers and sufficiently large rotation rates, however, a different instability occurs at lower Ra w ∼ Ek −1 in the form of anti-cyclonically drifting wall modes (Buell & Catton 1983;Pfotenhauer et al. 1987;Zhong et al. 1991;Ecke et al. 1992;Ning & Ecke 1993;Zhong et al. 1993;Kuo & Cross 1993;Herrmann & Busse 1993;Goldstein et al. 1993Goldstein et al. , 1994Liu & Ecke 1997, 1999Zhang & Liao 2009;Favier & Knobloch 2020). The relative contribution of the wall modes to the total heat transport depends on Γ (Rossby 1969;Pfotenhauer et al. 1987;Zhong et al. 1993;Ning & Ecke 1993;Liu & Ecke 1999) with decreasing contribution -roughly as the perimeter to the area ratio -with increasing Γ .
There are several regions of bulk rotating convection where rotation plays an important role, namely a rotation-affected regime and a rotation-dominated regime. In the former where Ro 1, heat transport varies as a power law in Ra, i.e., Nu = A(Ek )Ra 0.3 and can be enhanced or weakly suppressed by rotation relative to the heat transport without rotation (Rossby 1969;Zhong et al. 1991;Julien et al. 1996;Liu & Ecke 1997;King et al. 2009;Zhong et al. 2009;Liu & Ecke 2009) depending on the range of Ra and Pr. In the latter case in which Ro 1, heat transport changes much more rapidly with Ra in what is known as the geostrophic regime of rotating convection (Sakai 1997;Grooms et al. 2010;Julien et al. 2012;Ecke & Niemela 2014;Stellmach et al. 2014;Cheng et al. 2020).
Despite considerable previous work, the spatial distribution of flow and heat transport in confined geometries has not been well studied for high Ra and low Ro when one is significantly above the onset of bulk convection but still highly affected by rotation. Recently, Zhang et al. (2020) demonstrated in direct numerical simulations (DNS) and experiments that a boundary zonal flow (BZF) develops near the vertical wall of a slender cylindrical container (Γ = 1/2) in rapidly rotating turbulent RBC for Pr = 0.8 (pressurised gas SF 6 ) and over a broad range of Ra (Ra = 10 9 in DNS and for 10 11 Ra 10 14 in experiments) and Ek (10 −6 Ek 10 −5 in the DNS and for 3 × 10 −8 Ek 3 × 10 −6 in experiments). The BZF becomes the dominant mean flow structure in the cell for Ro 1, at which the large scale mean circulation (termed the large scale circulation or LSC) vanishes (Vorobieff & Ecke 2002;Kunnen et al. 2008;Weiss & Ahlers 2011b,a). Further, it contributes a disproportionately large fraction of the total heat transport. Another group (de Wit et al. 2020) also showed the existence of the BZF and its strong influence on heat transport using DNS for Pr = 5 (water) and Γ = 1/5 for Ek = 10 −7 in the range 5 × 10 10 Ra 5 × 10 11 . Thus, the BZF was observed in different fluids, in cells of different aspect ratios, and over a wide range of parameter values. Given the strongly enhanced heat transport in the BZF region (Zhang et al. 2020;de Wit et al. 2020), it is important to explore the BZF in detail. Here we investigate the robustness of the BZF with respect to Pr and to Γ in the geostrophic regime; we do not address here the transition from the low rotation state to the BZF.
Recently, Favier & Knobloch (2020) demonstrated for Ek = 10 −6 through DNS that the linear wall modes of rotating convection (Buell & Catton 1983;Zhong et al. 1991;Ecke et al. 1992;Zhong et al. 1993;Ning & Ecke 1993;Herrmann & Busse 1993;Kuo & Cross 1993;Goldstein et al. 1993;Liu & Ecke 1997, 1999Sánchez-Álvarez et al. 2005;Horn & Schmid 2017;Aurnou et al. 2018) evolve with increasing Ra and appear to be robust with respect to the emergence of bulk convection even with well developed turbulence. They suggested that the BZF may be the nonlinear evolution of wall modes, an idea that we address briefly but that requires significantly more analysis and comparison than can be included here.
In the present work, a series of DNS is carried out to study the robustness and the scaling properties of the BZF with respect to Rayleigh number Ra, Ekman number Ek , Prandtl number Pr , and cell aspect ratio Γ . We explore the extended scalings of the characteristics of the BZF including the width of the BZF, the drift frequency of the BZF, and the heat transport within the BZF in terms of these non-dimensional parameters. We first present our numerical methods, then discuss the results of our calculations, and conclude with our main findings.

Numerical method
We present results of direct numerical simulations (DNS) of RRBC in a cylindrical cell obtained using the goldfish code (Kooij et al. 2018;Shishkina et al. 2015) for Ra up to 5 × 10 9 and Ek down to 10 −7 . In the DNS, the Oberbeck-Boussinesq approximation is assumed as in . Centrifugal force effects are neglected since the Froude number in experiments is typically small, see Zhong et al. (2009) ;.
The governing equations based on the Oberbeck-Boussinesq approximation are Here, u = (u r , u φ , u z ) is the velocity with radial, azimuthal and vertical coordinates, respectively, ρ is the density, p is the reduced pressure, Ω = Ω e z is the angular rotation rate vector, T is the temperature with T 0 = (T + + T − )/2, and e z is the unit vector in the vertical direction. surfaces, constant temperature for the top/bottom plates and adiabatic for the sidewall.
To non-dimensionalise the governing equations, we use ∆ = T + − T − as the temperature scale, the cylinder height H as the length scale, and the free-fall velocity √ αg∆H as the velocity scale (the corresponding time scale is τ f f = H/ (αg∆)).
To evaluate the grid requirements for the simulations, we consider the thermal and velocity boundary layers (BLs) near solid boundaries. The thickness of the BLs near the heated and cooled plates are calculated as (2.4) This is the standard way to define the thermal BL thickness under the assumption of pure conductive heat transport within this layer, cf. Ahlers et al. (2009). The viscous BL thicknesses near the plates (δ u ) and near the sidewall (δ sw ) are defined as the distances from the corresponding walls to the location where the maxima of, respectively < u r 2 > t,φ,r + < u φ 2 > t,φ,r (z) and < u φ 2 > t,φ,z + < u z 2 > t,φ,z (r) are obtained. The velocity components are all averaged in time and over the surface parallel to the corresponding wall. The same criterion was used previously in studies of the sidewall layers in rotating convection, see Kunnen et al. (2011).
The computational grids are set to be sufficiently fine to resolve the mean Kolmogorov microscales (Shishkina et al. 2010) in the bulk and within the BLs (see table 2 in the Appendix). Grid nodes are clustered near the walls to resolve thermal and velocity BLs resulting in grids that are non-equidistant in both the radial and vertical directions. As rotation increases, the viscous BL gets thinner (Kunnen et al. 2008;Stevens et al. 2010;) so more points are required near boundaries: we take at least 7 points within each BL. The details of all simulated parameters and the corresponding grid resolution are listed in Appendix table 2 along with a benchmark comparison between Nu data from these simulations and from experimental data in compressed gases with similar Pr from Wedi et al. (2021); the agreement is excellent. To explore the robustness of the BZF with respect to Ra, Pr and Γ , we conducted simulations in three groups, i.e., in every group we vary only one parameter while keeping the others fixed. The specific parameter ranges are shown in table 1 (also included in several figures with Ra = 10 9 and Pr = 0.8 are data in the range 0.5 1/Ro 5 from Zhang et al. (2020); the calculation details for those values are included in the Appendix).

Boundary Zonal Flow structure
Our goal here is to explore the robustness of the BZF with respect to variations of control parameters. We follow closely the approach and characterization presented in Zhang et al. (2020) but focus on the geostrophic regime where the BZF is well developed. After presenting our main results, we consider the BZF with respect to wall mode structures. We begin with the influence of rotation on the overall temperature and velocity fields in the cell. In figure 1, for particular cases of 1/Ro = 0.5 (weak rotation) and 1/Ro = 10 (fast rotation), 3D instantaneous temperature distributions (figure 1a, d) and 2D vertical cross-sections (figure 1b,c,e,f) of the time-averaged flow fields are shown. The 2D views are taken in a plane P (figure 1 b, e), which in the case of a weak rotation is the LSC plane, and additionally in a plane P ⊥ which is perpendicular to P (figure 1 c, f). For slow rotation, a LSC spanning the entire cell with 2 secondary corner rolls are observed in P whereas a 4-roll structure is seen in P ⊥ , typical of classical RBC at large Ra and for Γ ∼ 1 (e.g., see Shishkina et al. (2014) and Zwirner et al. (2020)). Near the plates, the LSC and the secondary corner flows move the fluid towards the sidewall (figure 1b) so the Coriolis acceleration (−2Ωe z × u) induces anticyclonic fluid motion close to the plates. In the central part of the cell, at z = H/2, the radial component of the mean velocity, u r t , always points towards the cell center (figure 1a, b). Therefore, Coriolis acceleration results in cyclonic fluid motion in the central part of the cell as is also observed in the time-averaged azimuthal velocity field u φ in figure 2a. Cases at higher rotation rates are shown in figures 1d-f and 2 (see also Kunnen et al. (2011)). For both small and large rotation rates, the presence of viscous BLs near the plates implies anticyclonic motion of the fluid there. For strong rotation, the subject of this paper, with high and constant angular velocity Ω, the fluid velocity becomes more uniform along e z owing to the Taylor-Proudman constraint with larger components of lateral velocity compared to the vertical component as in figures 1e, f. Thus, anticyclonic fluid motion is present not only in the vicinity of the plates, but involves more and more fluid volume with increasing Ro −1 . With increasing rotation rate, anticyclonic motion grows from the plates toward the cell center whereas cyclonic motion at z = H/2 remains near the sidewall and becomes increasingly more localized there (figures 2c, d).
As introduced in Zhang et al. (2020), the BZF in rapidly rotating turbulent convection is characterised by an anticyclonic bulk flow, cyclonic vortices clustering near the sidewall, and anticyclonic drift of thermal plumes (see figures 3a,b and figure 4). These structures are associated with the bimodal temperature PDFs obtained in the measurements and DNS near the sidewall (Zhang et al. 2020;Wedi et al. 2021). The radial location r 0 where the mean fluid motion at z/H = 1/2 changes from anticyclonic to cyclonic as indicated by the solid line in figure 3 (see also inset of figure 6a below) is used to describe the width of the BZF δ 0 = R − r 0 . As one might expect, vertical coherence of the BZF is enhanced  by strong rotation. In figure 4, time-angle plots of the temperature at 3 different heights show that the drift frequency ω = 2πR(dφ(r u max φ )/dt)/(2πR/m) = mdφ(r u max φ )/dt is quite constant along z without significant phase differences, i.e., the BZF maintains good vertical coherence. Here, the mode number m equals 1 and dφ(r u max φ )/dt denotes the angular velocity of the temperature drift at r = r u max φ , where the maximum of the time-averaged azimuthal velocity is obtained. In the lower half of the cell, for z = H/4, warm plumes dominate so the warm regions (pink stripes) are wider, whereas in the upper half of the cell, for z = 3H/4, cold plumes dominate resulting in wider cooler regions (blue stripes). Similarly, figures 2c,d and 5a,d show that the zonal flow develops away from the top/bottom plates and extends vertically throughout the bulk. Figure 5 illustrates that owing to the drift, time-averaged fields in the vertical plane average to zero and do not capture important features of the flow motion, in particular, the u z -field. The averaged u 2 z , however, does retain important information about the locations of the Stewartson '1/3' and '1/4' layers (dashed lines) and the BZF (solid line).

Contribution to heat transport
An important and unexpected property of the BZF in rotating RBC is its disproportionately large contribution to the heat transport in the system. Figures 3a and 6 show that the averaged heat flux inside the BZF is much stronger than in the region outside the BZF. To be clear about the averaging we define where r 0 = R − δ 0 . The quantity R f is the ratio of the mean vertical heat flux within the BZF to the vertical heat flux averaged over the whole cell. The quantity R h reflects the portion of the heat transported through the BZF compared to the total transported heat. Especially, in figure 6a, the time-and φ-averaged radial profile at the mid-height for Ra = 10 9 , Pr = 0.8, and Γ = 1/2 shows a significant peak of heat transport inside the BZF, and the peak amplitude increases dramatically as rotation becomes stronger. Thus, although the width of the BZF shrinks with increasing rotation, thereby reducing the effective area of the BZF with respect to the whole domain, as shown in figure 6b, the increasing magnitude of the peak makes the heat transport carried by the BZF quite significant. Note that the annular BZF region of width δ 0 is smaller than the positive contribution to the heat transport as shown in the inset of figure 6a. Figure 6c reveals that the enhancement of the local heat transfer within the BZF increases more rapidly when rotation is very strong (1/Ro 10). As a result of these properties, the heat transport carried by the BZF for these parameter values is always more than 60% of the total heat transport at fast rotation (see figure 6d). Note, however, the effect of the BZF on the heat transport extends over a wider range r < r 0 ; over some range N u is actually negative (see figure 6a) implying an anti-correlation of vertical velocity and buoyancy, i.e., warm fluid going down or cooler fluid moving up. If we modify the annular averaging to take into account the decreased N u region as well as the inner structure of the BZF, i.e., we average over the extended region R − 2δ 0 r R, we get the ratio R * h which is also shown in figure 6d (open symbols) where one sees an even larger fractional contribution.
We also consider the dependence of the heat transport ratio R h as a function of Pr , see inset of figure 6d. Interestingly, for Pr < 1 we find 0.6 < R h < 0.7, whereas for Pr > 1 we have 0.3 < R h < 0.4 with a quite sharp transition for Pr ≈ 1. The origin of this rather sharp change emphasizes the important role Pr plays, perhaps through the competition between thermal and viscous BLs. Finally, comparing our computation of the total N u with increasing rotation with that of Wedi et al. (2021) (see Appendix figure  13), we conclude, given the close agreement, that the contribution of the BZF affects both measures of N u substantially and needs to be taken into account when considering the scaling of geostrophic heat transport in experiments and also in DNS with no-slip sidewall boundary conditions (see also de Wit et al. (2020)).

Dependence on Ra, Pr and Γ
We first discuss the qualitative robustness of the BZF with respect to Ra, Pr and Γ before we consider its quantitative spatial and temporal properties. We demonstrate the character of the BZF with respect to variations of Pr and Γ by considering time-angle plots of temperature T at z = H/2 and r = R. Figure 7a shows that the BZF exists in the flows at different Pr = 0.1, 0.8, 4.38 (also for Pr = 0.25, 0.5, 2, 3, 7, 12.3, not shown), i.e., from small to large Pr . Although there are some quantitative differences among the three cases, they all qualitatively demonstrate the existence of the BZF for more than two decades of Pr . The qualitative dependence of the BZF on the aspect ratio Γ is shown in figure 8 for three different aspect ratios: Γ = 1/2, 1, 2. The BZF is present in all three cases, has the same scaling of BZF width when scaled by H, i.e., δ 0 /H is independent of Γ (see figure 9d inset), and has a drift period (in units of free fall time τ ff = H/ (αg∆) = τ κ Pr −1/2 Ra −1/2 where τ κ = H 2 /κ is the thermal diffusion time) of about 70. The quantitative scaling of the drift frequency is analyzed later, and the data are tabulated in the Appendix (see table 2). The wavelength λ of the traveling BZF mode is independent of Γ for these three values in a straightforward way, as seen in figure 8, namely λ/H = π/2 so that the number of wavelengths around the circumference is m = 2Γ and the wavenumber is k = 2π/λ = 4/H. We note, however, that this relationship is for a limited number of values of Γ and control parameters Ra and Ro. Thus, we make no strong claims to its generality. Indeed, there is already evidence from de Wit et al. (2020) that for Γ = 1/5 one gets m = 1 = 2Γ , and we made additional measurements with Γ = 1/3 and 3/4 that also yield m = 1. We conjecture that owing to periodic azimuthal symmetry m will take on only integer values, similar to the situation for wall mode states (Ecke et al. 1992;Zhong et al. 1993;Goldstein et al. 1993;Ning & Ecke 1993;Liu & Ecke 1999) in cylindrical convection cells. Because of this periodic constraint, one cannot have m < 1 so small aspect ratios with Γ 1 have m = 1. We also note that the mode number dependence on Γ of the BZF is similar to that of the Γ dependence of linear wall state mode number, i.e., m ≈ 3Γ (Goldstein et al. 1993;Kuo & Cross 1993;Herrmann & Busse 1993;Ning & Ecke 1993;Liu & Ecke 1999;Zhang & Liao 2009). Given that our states have values of Ra that are 10 to 100 times greater than the linear wall mode onset Ra w , this difference is not unreasonable and the correspondence is very suggestive. In particular, a range of mode numbers are stable near onset (Zhong et al. 1993;Ning & Ecke 1993;Liu & Ecke 1999) owing to the azimuthally periodic boundary conditions. Significantly above onset there seems to be a selection towards lower mode numbers: e.g., Zhong et al. (1993) figures 3 and 8 with Γ = 2 show stable wall modes with m = 4, 5, 6, 7 near onset but only the m = 4, 5 modes persist for higher Ra which yields m = 2Γ and m = 2.5Γ , respectively, consistent with our results for the BZF (see also Favier & Knobloch (2020)).

Spatial and temporal scales
We next consider the quantitative dependence of the different layer widths on Ra, Ek and Pr, looking for a universal scaling form δ/H ∼ Pr ξ Ra β Ek γ . In figure 9a, the dependence of δ 0 /H on Ek for Ra = 10 9 , Pr = 0.8, and 2 < 1/Ro < 20 is shown to be consistent with a Ek 2/3 scaling whereas the width based on other measures scale closely as Ek 1/3 , i.e., γ takes on values of 2/3 and 1/3 for BZF width and velocity layer widths, respectively. (Because the statistical uncertainty in our reported exponents is of order 5-10%, we report fractional scalings consistent with the data to within these uncertainties; they are not intended to denote exact results.) As mentioned in Zhang et al. (2020), the BZF is characterized by bimodal temperature PDFs near the sidewall. This property was used in both DNS and experimental measurements to identify the BZF over a wide range of Ra. Here, we conduct a more detailed analysis of the DNS data to explore how the width of the BZF changes with Ra. We compute the width at fixed Ro = Ra 1/2 Pr −1/2 Ek so Ek = RoRa −1/2 Pr 1/2 . To determine the scaling with Ra at fixed Ro = 1/10, we have that δ/H ∼ Ra β−γ/2 . By multiplying by Ra γ/2 we obtain the scaling exponent β. In figure 9b, we plot (δ 0 /H)Ra 1/3 and (δ/H)Ra 1/6 corresponding to γ values of 2/3 and 1/3, respectively. From this plot, we obtain values for β of 1/4 and 0, respectively. Similarly for the dependence on Pr , we plot in figure 9c the corrected quantities (δ/H)Pr γ/2 which yields δ 0 /H scalings for ξ of −1/4 for Pr < 1 and 0 for Pr > 1. The other layer widths based on u φ , u z and F z are independent of Pr for Pr < 1 but do not collapse for Pr > 1. The separation of the different widths for Pr > 1 suggests some interesting behaviour not captured by our scaling ansatz. Finally, we can collapse all the data for BZF width onto a single scaling curve by plotting in figure  9d δ 0 /H = δ 0 /H Pr {1/4, 0} Ra −1/4 versus Ek (to compact the different scalings with Pr we denote them as Pr {1/4, 0} for scaling with Pr < 1 and Pr > 1, respectively) so that we can conclude that δ 0 /H ∼ Pr {−1/4, 0} Ra 1/4 Ek 2/3 . The results at one set of parameter values {Ra, Pr , Ro} are independent of Γ , see figure 9d inset, which implies that δ 0 /H ∼ Γ 0 (other dependences on Γ are not ruled out for other parameter values although it is reasonable to assume it to be general in the absence of other data). Thus, we plot in figure 9d all the data with different Γ , Pr , Ra and Ek to obtain scalings δ 0 /H ≈ 0.85Γ 0 Pr −1/4 Ra 1/4 Ek 2/3 for Pr < 1, (3.7) We plot in figure 9e the scaled BZF width (δ 0 /H) / 0.85Γ 0 Pr {−1/4, 0} Ra 1/4 Ek 2/3 . One sees that the data scatter randomly within ±10%, quite good agreement. The BZF drifts anticyclonically, the same as the direction of traveling wall modes of rotating convection (Zhong et al. 1991;Ecke et al. 1992;Kuo & Cross 1993;Herrmann & Busse 1993). We plot in figure 10a  de Wit et al. 2020; Favier & Knobloch 2020) and suggests that there is a correspondence between the states we observe and the nonlinear manifestation of linear wall mode states. The scalings we have determined for ω d with Ek and Pr will be useful in making a more quantitative comparison with the wall mode hypothesis among data sets with different Ek and Pr. Such an analysis is beyond the scope of the present work and will be presented elsewhere. These scalings, of course, depend on the definition of the time unit. Using the free-fall time or the vertical thermal diffusion time, respectively we obtain ω/ αg∆/H ≈ 0.015Γ 0 Pr −5/6 Ra 1/2 Ek 2/3 (3.10) which both show the same Ek scaling as δ 0 , i.e., Ek 2/3 , see figure 9a. For the three choices of time scale, the drift frequency decreases with increasing Pr for all Pr as opposed to the scaling of δ 0 /H which has different scaling for small and large Pr (see figure 9c and figure 12b). As reported in Zhang et al. (2020) and shown here in figure 3, the thermal structures drift anti-cyclonically, opposite to the azimuthal velocity which is cyclonic near the sidewall, as shown in figures 2b-d. We show in figure 11a that the drift frequency decreases as rotation increases with a scaling Ek 2/3 . In figure 11b, we show that the near-plate azimuthal velocity u peak φ is also anticyclonic and shows the same scaling behaviour with Ek (see figure 10b) as the BZF width and drift frequency. Based on this observation, we believe that the drift characteristics of the BZF are determined not only by the presence of the vertical wall but also by the near-plate region.
Finally, we consider the range of Ra and Ek in which the BZF is observed in this study. There are three regions defined by the onset of wall-mode convection Ra w ≈ 32Ek −1 , the onset of bulk convection Ra c = AEk −4/3 , and the transition from geostrophic convection (Grooms et al. 2010;Julien et al. 2012) to buoyancy dominated convection Ra t = Pr Ro 2 t Ek −2 where Ro t ≈ 1 (see Appendix figure 13) is the transition Rossby number out of the geostrophic regime (Julien et al. 1996;Liu & Ecke 2009;King et al. 2009;Weiss & Ahlers 2011b) as indicated in the Ra −Ek phase diagram in figure 12a. Our data fall solely within the geostrophic regime of bulk convection but we include the other zones for context. According to Chandrasekhar (1961) (see also Clune & Knobloch 1993), the critical Rayleigh number for the onset of convection is Ra c ∼ Ek −4/3 with a prefactor A that is weakly dependent on Ek , in the range 6-8.7 (Chandrasekhar 1961;Niiler & Bisshopp 1965); we use a value of 7.5 consistent with our range of Ek . To illustrate one aspect of this range, we consider the BZF width δ 0 /H versus Ek for Ra = 10 9 in figure 12b. A path of constant Ra = 10 9 (figure 12a) yields Ek w ≈ 32Ra −1 = 3.2 × 10 −8 , Ek c = (ARa −1 ) 3/4 = 8 × 10 −7 , and Ek t = Ro t Pr 1/2 Ra −1/2 = 2.8 × 10 −5 . Here the subscripts 'w', 'c' and 't' correspond, respectively, to the onset of wall-mode, bulk convection and transition from rotation to buoyancy dominated regime. These values are indicated by vertical dashed lines in figure 12b. Knowing the dependence of the critical Ra c and Ek and using the relations (3.7)(3.8), we can evaluate the smallest possible δ 0 for any fixed Ek , i.e., δ min 0 ∼ Ra 1/4 c Ek 2/3 ∼ Ek 1/3 (see δ min 0 in figure 12b). Connecting these onset points, we obtain the black line in the diagram, which is parallel to the Stewartson "1/3" layer scaling. The gap between these two black solid lines depends slightly on A but the ratio of the BZF width to the Stewartson layer width is constant (based on δ rms uz ) at the onset of convection (the fixed gap). Thus, although the BZF width decreases faster than the Stewartson layer as rotation increases, there is no crossing of the BZF boundary and the boundary of the Stewartson layer at extreme fast rotation because bulk convection ceases before they can cross. Note that all the data considered here fall within the geostrophic range of rotating convection; what happens in the wall-mode region is not addressed.
The other bound on the BZF scaling depends on when rotation becomes significant. An estimate is made based on a plot of Nu/Nu 0 versus Pr Ro 2 for our DNS and for experimental data from Wedi et al. (2021), see Appendix figure 13, where the data for  Figure 11. For fixed Ra = 10 9 and rotation rates, 1/Ro=5. 6, 6.7, 8.3, 10, 12.5, 16.7, 20: (a) drift frequency ω of BZF, (b) maximum absolute value of u φ near plates (mean value of two maxima). Everywhere Pr = 0.8, Γ = 1/2.
Ra from 10 8 to 10 14 merge onto a single curve. Here Nu 0 is the Nusselt number in nonrotating case. Using an empirical estimate Ro t ≈ 1 for the onset of the rotation dominated regimes, i.e., the geostrophic regime, we get an estimate for the largest possible δ 0 , for any Ek (grey line in figure 12, that is, δ max 0 ∝ Ra 1/4 t Ek 2/3 ∝ Ro 1/2 t Ek 1/6 ≈ Ek 1/6 ). (Pr Ro 2 t ≈ 1 is the onset in figure 13, but in experiment Pr varies from 0.7 to 0.9 and in DNS Pr = 0.8, thus here we take Pr = 1 which gives Ro t ≈ 1 for simplicity.) It is remarkable that the BZF regime is confined by these two critical lines (∼ Ek 1/6 and ∼ Ek 1/3 ) and the range confined in between gets broader for higher Ra. In other words, at low Ra, the BZF is only observed over a small range of rotation rates. At large Ra, the BZF exists over a much broader range of rotation rates (Zhang et al. 2020;Wedi et al. 2021). For any fixed Ra, the BZF exists in a certain Ek -range which is determined by the grey and black lines in figure 12b and the BZF width changes as δ 0 ∼ Ek 2/3 over that range. How the BZF contributes to the heat transport relative to the contribution of the laterally unbounded system in the geostrophic regime remains an open question. Further, the connection between the BZF and linear wall modes requires additional work to understand the relationship between the two convective states.

Conclusion
The BZF is found to be an important flow structure in rapidly rotating turbulent Rayleigh-Bénard convection in the geostrophic regime and is robust over considerable   (Chandrasekhar 1961;Niiler & Bisshopp 1965) (Ek c), and transition to rotation dominated regimes (Ek t) for Pr = 0.8, Ra = 10 9 , Γ = 1/2. ranges of Ra, Ek , Pr and Γ . The main structure, drift of plume pairs, is found to be a m = 2Γ -mode for the choices of Γ = 1/2, 1, 2; additional values of Γ = 1/3, 3/4 yield m = 1 suggesting mode 1 for Γ 1. In addition, the BZF carries a large portion of the total heat; its contribution to the total heat transport is about 60% of the heat transport at fast rotation, Ro < 0.1, and for Pr < 1. For Pr > 1, the BZF heat transport contribution drops to about 35%. Understanding this important contribution to the heat transport is essential in analyzing experiments in rotating convection in the geostrophic regime.
The scaling of the BZF width δ 0 depends on Pr , Ra and Ek as δ 0 /H ∼ Γ 0 Pr {−1/4, 0} Ra 1/4 Ek 2/3 (Pr −1/4 for small-to-moderate Pr and independent of Pr for large Pr ). The universal scaling of the BZF and the sidewall boundary layers is very clean for Pr < 1 but the BZF is less coherent for Pr > 1 and the sidewall boundary layer widths behave differently for those conditions. Further, the sharp decrease in the BZF heat transport contribution similarly marks a transition to a perhaps more complex BZF state for Pr > 1. The drift frequency of the BZF shows scaling ω/Ω ∼ Γ 0 Pr −4/3 RaEk 5/3 , indicating that the drift frequency decreases significantly as Pr increases, is proportional to Ra, and decreases rapidly with increasing rotation (decreasing Ek ). Interestingly, ω seems to be more robust than δ with respect to changes in Pr . Finally, the BZF shares qualitative and some quantitative characteristics with linear wall modes and establishing the connection between these two states will be an important area of future research.

Funding
This work was supported by the Deutsche Forschungsgemeinschaft (X.Z. and O.S., grant number Sh405/8, Sh405/7 (SPP 1881 Turbulent "Superstructures")); and the Los Alamos National Laboratory LDRD program under the auspices of the U.S. Department of Energy (R.E.E.).

Declaration of Interests
The authors report no conflict of interest.

Appendix
We tabulate here a full characterization of the parameters in the DNS (see table 2) and compare our results for N u with experimental data from Wedi et al. (2021). The excellent agreement is a strong indication that our DNS are fully resolved. We also include the resulting numerical values of the BZF drift frequency for different choices of time scale (see table 3), namely ω κ , ω f f , and ω d for different parameters values Γ , Pr , Ra, Ek , 1/Ro.