Wall modes and the transition to bulk convection in rotating Rayleigh–B´enard convection

,


I. INTRODUCTION
Thermal convection driven by buoyancy and subject to rotation [1,2] is a phenomenon of great relevance in many physical disciplines, especially in geo-and astrophysics and also in engineering applications.In the Rayleigh-Bénard (RBC) geometry a fluid with thermal expansion coefficient α, kinematic viscosity ν, and thermal diffusivity κ is bounded above and below by isothermal boundaries separated by a height H and heated from below so as to produce a temperature difference ∆.In a gravitational field with acceleration g, the strength of thermal forcing is determined by the Rayleigh number Ra ≡ αg∆H 3 /(κν), and the fluid type by the Prandtl number Pr = ν/κ.In the presence of rotation about an axis parallel to gravity, the influence of the Coriolis force is a balance of angular rotation rate Ω d and viscous diffusion time H 2 /ν.There are many ways to represent this balance namely non-dimensional Ω ≡ Ω d H 2 /ν and the Taylor number Ta ≡ (2Ω) 2 .Another measure popular in geophysics that we use here is the Ekman number Ek ≡ ν/ 2Ω D H 2 that nicely captures the limit of rapid rotation as a small parameter, i.e., Ek → 0. A similar parameter that measures the balance of buoyancy and rotation is the convective Rossby number Ro ≡ Ra/ (PrTa) = (Ra/Pr) 1/2 Ek = gα∆/H/ (2Ω d ) where the latter is the ratio of buoyancy frequency and rotation frequency.Note that Ro is independent of dissipation parameters.Unless otherwise specified, we use length scale H, time scale H/(gα∆), and temperature scale ∆ for non-dimensional variables.
In physical realizations of convection, there are solid horizontal boundaries, usually of cylindrical geometry with aspect ratio Γ ≡ D/H (D is the cylinder diameter), that have a finite thermal conductivity and diffusivity.In direct numerical simulations (DNS), these sidewall boundaries can be perfectly insulating or, less frequently, perfectly conducting.Experimental boundaries are somewhere in between these limits and tend towards the insulating case so as to minimize heat transport through the sidewalls.We consider here the case of perfectly insulating sidewalls; results for perfectly conducting boundaries will be presented elsewhere.
The global response of the system is measured by the molecular-diffusion-normalized heat transport, 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.
The multitude of states of rotating Rayleigh-Bénard convection (RRBC) in a fluid layer bounded above and below by isothermal boundaries have been elucidated over many decades going back to early linear stability calculations [3,4] and preliminary experiments [5] (see also [2]).They found that the onset of bulk convection was strongly suppressed by the rotation-induced Coriolis force which acts to inhibit the vertical motions associated with thermal convection via mechanisms associated with the Taylor-Proudman constraint on slow, inviscid flows [6].In particular, the critical Rayleigh number Ra c was found to vary approximately as Ek −4/3 for rapid rotation.Whereas the form of convection without rotation is in the form of convection rolls with characteristic length λ ≈ 2H (for a roll pair), rotation induces vortical motions near onset [7] with a corresponding characteristic length λ ∼ Ek 1/3 [3].Experimental heat transport measurements [8][9][10][11] confirmed the suppression of the convective onset by rotation, noted that some manner of convection ensued prior to strong bulk flow of linear theory, and revealed that rotation could, over some ranges in Ra and Pr, enhance heat transport over non-rotating turbulent RBC.The nature of the "subcritical instability" identified in the heat transport measurements was partially explained by linear stability calculations of stationary wall-mode states [11,12] but simultaneous heat transport and shadowgraph visualization experiments [13,14] together with theory [15,16] showed that these wall-mode states precess in the rotating frame and have precession frequency ω and integer (resulting from the cylindrically periodic azimuthal symmetry) mode number m.Further experiments [17][18][19] and linear stability calculations [16,20,21] showed beautiful correspondence between theory and experiment -here we denote Ra w ∼ Ek −1 as the critical value for the transition from the diffusive conduction state to one of wall-mode convection.A phase diagram showing the boundaries of different convective states in RRBC is presented in Fig. 1 where we plot Ra/Ra c (Ek = 0) versus Ek.
The enhancement of heat transport (Nu) in RRBC compared to its value in non-rotating RBC (Nu 0 ) has been observed in many experimental studies [8,14,18,[25][26][27][28][29] and DNS [30][31][32][33][34] and has been attributed to Ekman pumping from the thermal boundary layers.This enhancement occurs over a range of Ra, Pr and Ek, with states at low Pr or small Ek not showing an enhancement but rather a decrease in Nu for all Ra [9,11,24,[35][36][37][38]. We denote the value Ra t at fixed Ek as the transition from buoyancy-dominated flow to rotation-influenced behavior.
Another well-defined region of rapidly rotating convection is the geostrophic regime where the flow maintains geostrophic balance and is rotation-dominated.The importance of this regime to geophysical and astrophysical situations was articulated through the derivation [39] and simulation [40][41][42] of a reduced set of equations valid in the limit Ek → 0 and Ra → ∞ such that the scaled variable Ra ≡ RaEk 4/3 remains finite.The evolution of states from the onset of bulk convection at Ra c ≈ 8.7 depends on Pr and includes cellular patterns, convective Taylor vortices, plume states, and geostrophic turbulence.The range 1 ≤ Ra/Ra c ≲ 20 was explored using these reduced equations and a prediction was made that in the geostrophic turbulence regime Nu ∼ Ra RaEk 4/3 /A0 (dashed black line) where A0 = 8.7 is the asymptotic Ek → 0 value versus Ek.The lines for Rac (dashed blue line) [4,22,23] and Raw (solid blue line) [20] are from linear stability analysis with idealized boundary conditions on laterally infinite domains.Rag marks the approximate boundary of the geostrophic, rotation-dominated region and Rat denotes the transition to buoyancy-dominated flow [1,2,24].The precise details of this diagram depend on parameters including Pr and Γ and are not included here.The data points show the parameter values reported in this paper.
the interior quasi-geostrophic (QG) flow rather than being boundary layer controlled, a result that was demonstrated in the simulations only for Pr = 1 over a short range in RaEk 4/3 .The experimental and DNS test of these predicted states [26,28,29,36] and their Nu scalings is challenging [43] owing to the need for very small Ek < 10 −6 , large Ra, and the experimental resolution to measure rather small deviations from Nu ∼ 1, i.e., over a range 1 ≤ Nu ≲ 20.
The flow structure of rotating convection for sufficiently small Ek ≲ 0.02 starts with wall modes that arise as a supercritical bifurcation from the no-flow base state [13-16, 20, 21].These initial flow states persist to higher Ra and serve as the nonlinear base state for the onset of fully three-dimensional bulk convection.In the idealized case usually considered in DNS, the sidewall is perfectly adiabatic.In experiments, however, this is not the case.The details of the onset and evolution of wall modes depend on the ratio of thermal conductivity of the sidewalls to that of the non-convecting fluid.Often the sidewalls are acrylic or polycarbonate with thermal conductivity of order 0.2 W/(m•K) compared to, for example, water with 0.6 W/(m•K).These conditions yield small (order 10%) modification of the asymptotic critical parameters for insulating sidewall boundary conditions [21] with higher Rayleigh number for the onset of the wall modes Ra w (74.4 versus 63.6), lower procession frequency at the onset ω κc = ω c H 2 /κ = (57.7 versus 66.1), and lower wave number k w (5.55 versus 6.07).Here we consider the idealized case of perfectly insulating sidewalls.

A. Direct Numerical Simulations
We consider rotating Rayleigh-Bénard convection (RRBC) in a cylindrical container with radius R and height H (Γ ≡ 2R/H).The governing equations for RRBC in the Oberbeck-Boussinessq approximation include the continuity equation ∇ • u = 0 and the momentum equation and heat equations: Here, p denotes the reduced kinematic pressure, e z is the unit vector pointing upward, and T 0 ≡ (T + + T − )/2.We do not include any centrifugal acceleration because, in laboratory realizations of RRBC, it is negligible compared to the gravitational acceleration, i.e., the Froude number is very small, F r ≡ Ω 2 d R/g ≪ 1.We apply the standard boundary conditions for RBC which are no-slip for the velocity (u = 0) at all walls, isothermal at the bottom plate (T = T + at z = 0) and top plate (T = T − at z = H, where T − < T + ) and adiabatic (∂T /∂n = 0) at the sidewall for the temperature.Here n denotes a unit vector orthogonal to the boundary surface.
Taking as the reference quantities τ = H/ √ αg∆H for time, ∆ for temperature, H for length, √ αg∆H for velocity, and αg∆H for the reduced kinematic pressure, from equations ( 1)-( 2) we obtain the following dimensionless governing equations The resulting set of equations ( 3)-( 4) together with ∇ • u = 0 is solved numerically, using the direct numerical solver goldfish [44,45] in its latest version [46,47], which is optimized for efficient computation with massive parallelization.The computational code uses a fourth-order finite-volume discretization on staggered grids and a third-order Runge-Kutta time integration scheme.The direct numerical solver goldfish utilizes a two-dimensional pencil decomposition employing the 2DECOMP library [48] and a hdf5 file management.Unless otherwise specified, all quantities are made dimensionless by τ , H, and ∆ (velocity by H/τ ).The exception is the radial distance which we explicitly scale with the radius R, i.e., r/R = (2/Γ)r/H.

B. Linear Stability Theory
One can establish using linear stability calculations the onset of RRBC of a laterally-infinite fluid layer [4,22,23] and of precessing wall modes [16,20,21,25] that arise owing to the existence of sidewalls.Here we summarize the asymptotic results (see Appendix for empirical fits for larger Ek) for Pr ≳ 0.68 and Ek → 0. For the bulk instability the onset Rayleigh number Ra c , critical wave number k c and characteristic length λ c are: where k c and λ c are normalized by H −1 and H, respectively, and λ c k c = 2π.
For the precessing wall modes with perfectly insulating boundaries, one has for a planar wall (no curvature) and to second order in Ek (see also Appendix) [20,21,25]: where the spatial scale is H and the precession frequency ω w is made dimensionless with the time scale H 2 /κ for ω κw , H/(gα∆) for ω f fw , and Ω −1 for ω κc and ω dw , respectively.Recall that for realistic sidewall boundary conditions [21] these asymptotic values are modified by about 10%; something that needs to be considered when comparing experimental and DNS results.For finite cylindrical geometries, especially for small Γ, the values also shift slightly with smaller Ra w for m = 1, 2 [16,20] compared to their planar wall values.
The paper is organized to generally track the evolution of the RRBC state for Ek = 10 −6 and Pr = 0.8 over the range 3 × 10 7 ≤ Ra ≤ 5 × 10 9 .In Section II, we describe the properties of wall modes including (II A) effects of aspect ratio, (II B) their spatial and temporal structure in terms of linear state eigenfunctions, and (II C) the evolution of the nonlinear state through a subcritical bifurcation.We also consider (II D) the transition to the bulk state and (II E) the associated spatial and temporal properties of the coexisting wall mode and bulk mode states where we use the term boundary zonal flow (BZF) to label the wall-localized state in the presence of bulk flow.In subsection II F, we focus on the heat transport (Nu).We provide conclusions in Section III and end with acknowledgements in Section IV.The Appendix has details of the data presented here, some extra features of wall modes, and the empirical fits of critical parameters for larger Ek.

II. WALL MODES
The goal of this paper is to address the evolution of the states of rotating convection at small Ek from the onset of wall modes through the onset of bulk convection and up to the upper threshold of the geostrophic regime (Fig. 1).We begin with a description of pure wall modes for Pr = 0.8, Ek = 10 −6 and 10 −4 , and for a range of 1/5 ≤ Γ ≤ 5 (details of these parameters are presented in Tables II in the Appendix).Within this range we identify the critical onset Ra w , the mode number m (and corresponding azimuthal wavenumber k), and different measures of temporal and spatial degrees of freedom.Owing to the previous lack of data for wall modes at small Ek, we systematically consider aspects of the onset of wall modes with a particular emphasis on the effects of small Γ < 1.We then address the nonlinear evolution of the wall modes including a subcritical transition to time dependence and the mechanisms for that transition.We identify three distinct regions.First, one has steady (in the precessing frame) wall mode states for 3 × 10 7 ≲ Ra ≲ 3 × 10 8 .For 4 × 10 8 ≲ Ra ≤ 9 × 10 8 , there is a subcritical bifurcation to time dependence followed by increasingly nonlinear states.The mechanism for this bifurcation is the lateral ejection of plumes originating near the azimuthal zero-crossings of temperature and u z and propagating into the interior, see also [49].Finally, one gets a transition at Ra ≈ 10 9 to bulk convection where the wall mode states are in a dynamic balance with the bulk states in the form of a BZF.We characterize this later state in the range 10 9 ≤ Ra ≤ 5 × 10 9 .Some aspects of these states emphasizing the connection between wall modes and BZF were reported earlier by us [50].
The critical Ra w for Γ = 1/5, 1, 5 are determined from the DNS with Pr = 0.8 at Ek = 10 −4 and Γ = 1/2, 2 at Ek = 10 −6 , see Table I.One can understand the dependence of Ra w on Γ by considering the discrete wavenumber k m of the corresponding mode number m (k m = 2m/Γ) which is selected owing to azimuthal periodicity and comparing it to the critical wavenumber k w for a planar (flat) sidewall boundary.The marginal stability curve, Fig. 3, is defined by the k dependence of the onset Ra designated as the marginal stability Rayleigh Number Ra M (k).For small differences k − k w the dependence of Ra M is quadratic, i.e., ϵ M ≡ Ra M (m)/Ra w − 1 = ξ 2 0 (k − k w ) 2 [17,19,21].The quadratic approximation is valid only for small k − k w and ϵ whereas for larger values there are terms higher order in k − k w (see, e.g., [55]).Nevertheless, this representation gives a reasonable semi-quantitative description of the variation of Ra M for different parameter values.The difference k − k w is particularly important for small Γ as we will see below.The  coefficient ξ 2 0 was evaluated for experimentally realized [14,17] sidewall thermal conductivity as a function of Ek [21] (see Appendix).In Fig. 3(a), the curve is evaluated for Ek = 10 −6 with k w = 5.75 and ξ 0 = 0.18, and the labeled data points are the discrete values k m for Γ = 1/5, 1/3, 1/2, 3/2 [49], and 2 which have the lowest values of Ra M /Ra w .One can see the large and unintuitive dependence on Γ ≲ 1/2 where m = 1: for Γ = 1/5, Ra M /Ra w = 1.6 whereas Γ = 1/3 yields Ra M /Ra w = 1.003 and Γ = 1/2 has Ra M /Ra w = 1.1.The increase in Ra M /Ra w is non-monotonic in Γ and depends sensitively on k w .Another observation concerns the predicted lowest value for Γ = 3/2 which is m = 4 whereas m = 6 was observed in DNS [49] upon increasing Ra but m = 4 when decreasing Ra.This would be consistent with the m = 6 mode resulting from starting the DNS for Ra ≳ Ra w from random initial conditions at Ra = 0 so that the state would evolve to an m ̸ = m w where m w is the value corresponding to the smallest k m − k w .The effect of a sudden increase in Ra was used in experiments on wall modes [14,17,19] to access a wide range of available m.In Fig. 3(b), we consider Ek = 10 −4 where k w = 4.9 and ξ 0 = 0.2 for Γ = 1/5, 1, 5.Here the shift in Ra M /Ra w is even larger for Γ = 1/5 than for Ek = 10 −6 but as Γ increases the density of discrete k m increases as ∆k ≡ k m+1 − k m = 2/Γ so that they increasingly cluster near Ra M /Ra w ≈ 1.For example, for Γ = 5, we observe modes m = 14, 15 compared to the predicted values with the lowest Ra M , m = 12, 13 (a small shift in k w would account for that difference).More importantly, we again appear to find that the selected m can be different than the one closest to critical m w owing to starting from random initial conditions; m = 14 near onset whereas m = 15 above onset and similarly for Γ = 1, m = 2 near onset and m = 3 above onset as shown in Fig. 4. When we slowly raise Ra from near onset, the selected m near onset is preserved.
To directly evaluate the observed and predicted onset values Ra M , we plot in Fig. 5(a), Ra M /Ra w = ϵ M + 1 for different observed values of m at corresponding Γ.The dramatic effect for Γ = 1/5 where m = 1 noted above is well captured by the marginal stability model with an observed factor of 2.6 in the onset Ra M compared to a predicted ratio of about 2.1.For Γ = 1/2, 1, there is actually a decrease relative to Ra w which is captured in calculations for finite Γ [16], see Fig. 5(b).As Γ is further increased, the observed m values yield k ≈ k w so there is little or no shift in the onset Ra M .Given that ϵ M is not precisely quadratic except for small ϵ, ξ 0 is computed for imperfectly insulating sidewall boundary conditions, and the computations for Ra w and k w do not include rigid top/bottom boundary conditions, the agreement is quite satisfactory.

B. Wall mode spatial and temporal structure
We now turn to the spatial and temporal structure of the wall mode states.In previous work [50,53,54], we considered the radial boundary length scales of wall modes and BZF states defined by δ uz = max ⟨u 2 z ⟩ 1/2 ϕ,t and the zero crossing δ 0 of ⟨u ϕ ⟩ ϕ,t , i.e., δ 0 is the distance to the sidewall where ⟨u ϕ ⟩ ϕ,t = 0, where ⟨•⟩ ϕ,t means averaging in the azimuthal direction and in time.As we discuss below, δ 0 may not be a good measure of radial extent owing to the midplane symmetry of u ϕ .Here we consider characteristics of the stationary and subcritical wall modes as represented in the fluid fields: temperature T , vertical velocity u z , azimuthal velocity u ϕ , radial velocity u r , and vertical vorticity ω z .
One feature of wall modes that is perhaps not appreciated is the decoupling of the radial length scales of the temperature field T from the vertical velocity u z and vertical vorticity ω z fields with decreasing Ek [20].In Figs.6(a,b), we show, respectively, the radial eigenfunctions of u z (r) and T (r) (interpolated from [20] for Ek = 10 −6 ) where a Fourier-Bessel expansion is an excellent representation of u(z) (anticipating the results in a cylindrical cell) (1/e decay length) being dependent on Ek only through its weak dependence k w (Ek) whereas the u z length scale (first zero crossing) δ uz ∼ Ek 1/3 has an explicit Ek scaling.This separation of scales has important implications for RRBC in small aspect ratio containers with rapid rotation and for understanding the heat transport properties of wall modes and the BZF.
In Fig. 6(c), we consider the different length scales associated with the eigenfunctions of u z and T .We take the geometric length scale L R of the container as R so that L R = Γ/2 and compare with the length scales δ T , and δ uz 0 evaluated at Ek = 10 −6 and 10 −4 .One obtains δ T ≈ 0.13 (consistent with [20] which predicts scaling of order k −1 w = 0.17), and δ uz = 2.5Ek 1/3 which yields values of 0.12 for Ek = 10 −4 and 0.025 for Ek = 10 −6 .Thus, for Ek = 10 −4 , δ T ≈ δ uz whereas for Ek = 10 −6 , δ T ≈ 5δ uz .For Γ ≈ 1/3, δ T ≈ L R .Thus, the temperature field is highly constrained by the geometry for small Γ ≲ 1.On the other hand, for Ek = 10 −6 , δ uz 0 ≪ δ T and only becomes of order R for very small Γ ∼ 1/20.Thus, u z should not be affected by geometry and should have a form in cylindrical geometry similar to that of a planar boundary.We expand on the disparity in length scales by considering in Figures 7(a-c), the wall-mode linear temperature eigenfunction of a planar wall mode state with perfectly insulating sidewall boundary conditions plotted versus r/R for different Γ: (a) 1/2, (b) 1, and (c) 2, 5.The positive and negative representations corresponding to opposite sides of the cylinder and the sum of the two as an approximation of the effects of finite geometry (see also [16]) are illustrated.For Γ = 1/2, there is considerable overlap of the two functions.For Γ = 1, one has just about achieved separation into individual contributions but one needs Γ ≳ 2 for a clean separation of the thermal field.For comparison, we have plotted vertical dashed lines showing δ uz .
Whereas the analysis above was for a planar wall [20], our finite cylindrical convection cell suggests that the appropriate representation in that geometry are radial eigenfunctions using sums of Bessel functions of the first kind, J m (kr) of order m [16].It is therefore natural to consider u z (r) fields as Fourier-Bessel transforms with m = 1: where u z (r) is the radial component of the vertical velocity field and j 1n are the n-th order zeros of J 1 (x).Empirically, we find that N = 50 yields excellent fits to u z (r) which is primarily consequential in the advection of heat and the (a determination of the global heat transport Nu.In Fig. 8(a), we show the radial dependence of mean u z for Pr = 0.8, Γ = 1/2, and Ek = 10 −6 evaluated at its azimuthal maximum (in the precessing frame) and at the mid-plane (z = 1/2) height for 3 × 10 7 ≤ Ra ≤ 3 × 10 8 .The inset is the normalized Fourier-Bessel decomposition coefficients c n which has a simple form and is quite independent (modulo an overall scale factor that increases with Ra for Ra < 3×10 8 , i.e., for steady wall modes).An expanded view in Fig. 8(b) shows the details of the profiles near the sidewall and the definition of the radial extent of the main wall peak -here we take the first zero crossing of u z (r) as the characteristic width δ uz .The inset shows the radial profile obtained by averaging the normalized c n values over the Ra range.The data points are from the amplitude-normalized planar eigenfunction of u z shown in Fig. 6(a) and match almost perfectly, indicating again that cylindrical curvature is not important here for u z .For larger Ra ≥ 4 × 10 8 , the profiles change slightly as shown in Fig. 8(b) with a larger spread in δ uz .This change is associated with a secondary bifurcation of the wall modes as discussed below.It is quite surprising that the radial profile of u z maintains a close correspondence with the linear eigenfunction over such a large range of Ra.In Fig. 8(d), we compare the linear eigenfunctions (dashed lines) for u z , u ϕ , and u r with data for Ra = 3 × 10 7 (solid lines) where one sees very close correspondence between the eigenfunctions and the data for u z and u ϕ with a slightly smaller first zero crossing for u ϕ .u r has a very different shape with a maximum value at a radius r max (we define the characteristic length of u r as δ ur = 1 − r max ) near the first zero crossing of u z , and there is a distinct difference between the eigenfunction and the data indicating that u r is more affected by finite wall curvature than the other convective fields.Figure 8(e) shows the peak amplitudes of u z , u ϕ , and u r plotted as a function of Ra − Ra w .Close to onset u ∼ (Ra − Ra w ) 1/2 as expected for a supercritical bifurcation.The monotonic increase in u z and u r appears to saturate near the onset of bulk convection at Ra ≈ 10 9 whereas u rmax drops rapidly after the onset of the subcritical instability described below.Finally, there is a small drop in u ϕmax at Ra ≈ 9 × 10 7 which signals the end of the weakly nonlinear wall mode regime.This signature shows up in measures of T .[20] for u ϕ (red) uz (blue), and ur (black).The first zero crossing of uz is at about the same r as the maximum of ur.All profiles are scaled such that the maximum values are 1.(e) Variation of radial maximum values uz max (blue solid circles), uz rms,max (light blue inverted triangles), u ϕmax (red solid squares), and ur max (black solid triangles) vs Ra − Raw.Linear scaling u ∼ (Ra − Raw) 1/2 (blue, red, and black dashed lines, respectively).Vertical dashed lines from smaller to larger Ra: end of quasi-linear regime (long-dashed) Ra ≈ 9 × 10 7 -note decrease in u ϕ , subcritical instability to lateral jet ejection (dashed), transition to aperiodic (chaotic) jet ejection (long-dashed), and bulk onset (dashed).The other dominant field contributing to the heat transport is the temperature field which, as noted above, has a much broader radial distribution and is strongly affected by small Γ as was suggested by considerations of the temperature eigenfunction [20] shown in Fig. 7(a).To get an overall flavor of the wall mode states for Ra ≤ 3 × 10 8 , we show in Fig. 9 some representative azimuthal mid-plane contours and corresponding fields; in all of them T extends significantly over most of the radial domain.Fig. 9(a) shows the nearly sinusoidal mean contour for Ra = 3 × 10 7 (ϵ = Ra/Ra w − 1 = 0.07) with counter clockwise rotation and clockwise precession.For Ra = 2 × 10 8 (ϵ = 6.1) in Fig. 9(b), the contour is more nonlinear but the precession remains steady whereas for Ra = 5 × 10 8 (ϵ = 17) there is increased distortion of the mean contour and the state is time dependent.By going to the precessing frame one can obtain well-converged statistical averages of the radial profiles of different fields.We show in Fig. 10, the radial distribution of the max value of T (averaged in precessing frame) in a mid-plane (z = 1/2) cross section.For Ra ≤ 10 8 , the profiles are remarkably linear, more so than the eigenfunction, indicating that the finite geometry is having a significant effect on T (r).For larger Ra > 10 8 , T (r) becomes increasingly nonlinear with radial undulations owing to the transition to time-dependent states and bulk/BZF convection described below.In Fig. 10(c), we show the Ra variation of the max value of T which indicates the weakly nonlinear growth of the wall mode amplitude but for higher Ra ≈ 9 × 10 7 this reverses and T decreases slowly owing to the decreasing vertical gradient of T in the cell center as discussed below.This feature is correlated with the small drop in u ϕmax in Fig. 8(e).Also shown is T rms for Ra ≳ 10 9 owing to the unsteady and intermittent precession once the bulk mode is present which makes averaging in the precessing frame much slower to converge.
With the insights provided by the structure of u z (r) and T (r), one can understand the radial structure of the heat transport Nu(r).In Fig. 11(a), the normalized (peak values set to 1) fields u z (r), T (r), u z (r) T (r), and Nu(r) are shown for an average over time-independent wall modes 3 ≤ Ra/10 7 < 30 .T (r) varies slowly with r and the fields are time independent in the precessing frame.If the azimuthal phase difference between u z and T is small, i.e., the place where u z is maximally upward (downward) is the same as where T is maximally positive (negative) with respect to the mean, then one has Nu(r) ∼ u z (r) to within about 10%.Taking Nu ∼ ⟨T ⟩⟨u z ⟩ improves the agreement to about ±5%.Averaging over the time-dependent wall mode states in the range 40 ≤ Ra/10 7 ≤ 90 yields similar results shown in Fig. 11(b).More surprising is that averaging over the bulk/BZF states 100 ≤ Ra/10 7 ≤ 500, Fig. 11(c), also produces a fairly good correspondence where here we subtract out the bulk contribution in the cell interior.There is a widening zone of boundary zone influence in the BZF but the form of Nu near the wall is driven predominately by the vertical velocity u z of the wall mode which remains highly localized to within a radial width of order Ek 1/3 characteristic of the linear wall-mode eigenfunctions.A more quantitiative test of this comparison requires detailed amplitude averaging of the u z (ϕ) and T (ϕ) for which we would need better converged statistical quantities in the bulk phase with Ra > Ra c ; we will consider this in future work.
In addition to the radial structure of the fields, we consider their azimuthal structure.In Figs. 12 (a-d), we show instantaneous vertical temperature fields T (r = 0.98R, ϕ, z) for values of Ra in the interval 3 × 10 7 ≤ Ra ≤ 5 × 10 9 .A very nonlinear profile develops with increasing Ra with only the smallest Ra = 3 × 10 7 being well approximated by a single-mode sinusoidal function.To characterize this nonlinearity for Ra = 5 × 10 8 , we show in Figs. 12 (e-h) azimuthal profiles of T , u z , u ϕ , and u r , respectively, averaged in the precessing frame (ϕ p = ϕ 0 + ωt) at r = 0.98R and at three different z = 0.5, 0.8, and 0.95 (horizontal dashed lines in Fig. 12 (c)).For T and u z , the mid-plane profiles are rather close to the linear eigenfunction, i.e., to the lowest Fourier mode.The eigenfunction solutions for the ϕ dependence of T and u are discrete Fourier series N n=1 a n sin(nmϕ)e iωt with the first term being sin(mϕ + ϕ p ).For u ϕ and u r , the mid-plane is problematic because it corresponds to a zero-crossing.Considering z = 0.8, one also sees reasonable correspondence with the lowest Fourier mode for u ϕ but u r is very different with almost no weight in the lowest Fourier mode.Close to the top (or bottom) boundary with z = 0.95, the nonlinearity of the wall mode is revealed where there is a sharp front-like feature appearing in T , u z , and u ϕ and a pulse-like structure for u r .The emerging contribution of bulk-like excitations are visible in the azimuthal oscillations of order 2λ c at the leading edge of u r shown in the inset of Fig. 12(h).The shock-like features for T , u z , and u ϕ vanish at the upper boundary consistent with rigid and isothermal boundary conditions as shown in the insets of (f-g): T approaches zero consistent with global heat transport constraint, u z because the zero normal velocity BC, and u ϕ because of the zero tangential BC.The latter should be controlled by an Ekman boundary layer condition of order δ E = Ek 1/2 = 0.001 but more quantitatively there is a coefficient relating δ E to the viscous boundary layer of order 3 [56].The blue dashed line in (g) is proportional to e (1−z)/0.003consistent with an Ekman boundary layer.More detail on the z-dependence is provided in the Appendix (Fig. 28) where we compute for Ra = 5 × 10 8 : ⟨X⟩ rms and ⟨X − X (1) ⟩ rms where X is the convective field and X (1) is its lowest Fourier mode.
We characterize the transition to an oscillatory state through the time dependence of the heat transport which is time-independent for precessing wall modes.For 10 8 ≤ Ra ≤ 3 × 10 8 , Nu(t) is a decaying oscillation of the form Nu ≈ Nu ∞ + Nu 0 cos (ω 2 (t − t 0 ))e −t/τ0 starting from some random initial condition as shown in Fig. 13(a-b).For 4×10 8 ≤ Ra ≤ 6×10 8 , the waveform is approximately stationary (or slightly increasing saturation, see Appendix V C) with amplitude ∆Nu (mean peak-to-peak) and an oscillation frequency ω 2 as in Fig. 13 (c).The nature of the onset of time dependence is a subcritical bifurcation as shown in Fig. 14 where we plot the max-min envelope of Nu(t), the difference Nu max − Nu min = ∆Nu, the oscillation frequency ω d2 /ω d , and the decay frequency ω d0 /ω d versus Ra where ω d is the precession frequency in units of Ω −1 .
We show in Fig. 14 (a), the Nu(t) envelope defined by its maximum and minimum values as a function of Ra.For Ra ≤ 3 × 10 8 , Nu is time-independent and Nu max = Nu min .There is a rapid increase in Nu max for Ra ≥ 4 × 10 8 whereas the minima stay roughly constant at Nu min ≈ 6.3, very close to its value at Ra = 3 × 10 8 .The solid (dashed) blue line is an 3 parameter adjustable fit to the stable (unstable) portions of a subcritical bifurcation.This schematic illustration of the nonlinear time-dependent amplitude (solid line) and the unstable subcritical region (dashed line) is presented in a conventional bifurcation diagram in Fig. 14 (b).In Fig. 14 (c), we show ω d2 /ω d where the black circles and corresponding black dashed line indicate the oscillatory part of the decaying solutions associated with the time-independent state whereas the blue circles and dashed line show ω d2 /ω d ∼ Ra ≈ 1.7.The decay frequency ω d0 /ω d ∼ (Ra w2 − Ra) (magenta diamonds and dashed line), where its zero intercept indicates the Rayleigh number Ra w2 at which the wall mode state becomes unstable to infinitesimal perturbations.For Ra ≥ 7 × 10 8 , the waveform becomes chaotic with increasing fluctuations and irregular oscillation frequency.The inference from Fig. 14 (a) is that the minimum value Nu min is when the system is close to the pure wall mode state whereas Nu max occurs when the lateral jet makes the maximal contribution to the heat transport.We discuss this in more detail below when we consider the full range of Nu.The time dependence of the wall mode was noted earlier [49] for Pr = 1, Γ = 3/2, and Ek = 10 −6 where the emission of horizontally-propagating thermal plumes originating within the wall-mode region and moving into the interior were observed for Ra ≳ 5 × 10 8 , consistent with our results (see also [57,58]).The mechanism for these fluctuations was hypothesized to result from a shear instability of a Stewartson layer of width ∼ Ek 1/4 [49,59] that develops a net mean flow for Ra ≫ Ra w .The criterion for instability in a differentially sheared rotating layer (not a wall-bounded flow) is Re ≈ 10 to 20 [60] which is qualitatively consistent with the computational results [49].For comparison we get Re values of 6 and 12 for Ra = 2 × 10 8 and 5 × 10 8 , respectively, based on the wall-bounded layer and 15 and 30 based on the outer free shear layer.Although the qualitative Re are about right for shear instability of a zonal barotropic mode [61], a more detailed analysis was beyond their scope [49].Our results show that small scale fluctuation structures appear in the wall-bounded zone of order Ek 1/3 , see Figs. 15 (d,g,h), rather than in the outer shear layer of order Ek 1/4 (see also below for length scale of striations).Further analysis and characterization of this instability is also beyond the scope of our work.Here we elucidate the onset of this time-dependent degree of freedom for Pr = 0.8, Γ = 1/2, and Ek = 10 −6 .In Fig. 15, we show instantaneous fields of horizontal and vertical profiles, respectively: (a,e) T , (b,f) u z , (c, g) u ϕ , and (d, h) u r for Ra = 5 × 10 8 , above the subcritical onset of time dependence.Figs.15(a-d) show horizontal cross sections at (a) z = 1/2 and (b-d) z = 0.8 whereas (e-f) illustrate the wall mode region r/R = 0.98 with vertical 0 ≤ z ≤ 1 and horizontal 0 ≤ ϕ ≤ 2π.The horizontal fields show the lateral plume emission into the interior near the crossing point of positive/negative u ϕ and u z whereas the vertical fields show a spatially oscillatory signature of order 2λ c that generates in the build-up to the emission.The oscillations are strongest for u ϕ and u r compared to less obvious features in T and u z .In the bulk region represented in (a-d), u z is quite small compared to its value near the wall whereas the other fields have significant contributions in the interior.Thus, although the contribution to Nu ∝ u z seems small; the perturbations in the interior region probably act as the foundation for the nucleation of bulk convection.At the same time, the vertical striations of order 2λ c seem to be related to bulk instability implying that the underlying wall mode acts as the foundation for the nucleation of bulk instability.Thus, the jet instability and the 2λ c bulk structures seem inextricably linked.The striation length scale (proportional to Ek 1/3 ) suggests that the origins of the instability arise within the wall shear zone rather than in the outer free shear zone of order Ek 1/4 [49].Other characteristics of the sub-critical state include unsteady precession and the qualitative change in the shape of the contours of T : compare Figs.9(b, c).

D. Transition to bulk modes
The transition to finite amplitude bulk convection occurs for Ra ≈ 10 9 as compared to the theoretical prediction of Ra c ≈ 7.8 × 10 8 for an infinite layer without sidewall effects.For comparison with the theoretical predictions of the NHQGS equations (see also below), the range 8 × 10 8 < Ra ≤ 5 × 10 9 corresponds to 1 ≲ Ra/Ra c ≲ 6 (8 ≲ RaEk 4/3 ≲ 60) where the plume state is predicted to be prevalent for Pr = 1 [40,41].In Fig. 16, representative horizontal cross sections for instantaneous fields of T , ω z , u z , u ϕ , and u r are shown (z = 1/2 for T and ω z and z = 0.8 for u).At Ra = 10 9 , there are strong lateral jets of contrasting thermal signature for T in (a) which are generated near the intersection of cold and hot zones, i.e., where T ≈ ⟨T ⟩ and u z ≈ 0 (see Fig. 12).There are corresponding regions of cyclonic vorticity generation as shown for ω z in (b).Around the sidewall boundary, there are azimuthal variations of T and ω z with some spatially-intermittent structure.When bulk convection is well established at Ra = 5 × 10 9 as in Figs.16(c,d) for T and ω z , respectively, the BZF remains visible in the T field but there are significant small scale temperature and vorticity fluctuations with spatial scales of the order of the linear length scale λ c = 4.8Ek 1/3 (0.051 for Ek = 10 −6 , see Eq. 7).In Figs.16(e-g), u z , u ϕ , and u r , evaluated for Ra = 4 × 10 9 and at the same time, show that u z displays little remnant of the horizontal jet flow but u ϕ and u r show that even at Ra/Ra c ≈ 4, there are large scale motions suggesting residual emanations from the BZF region -note in particular the in/out radial velocity structure in Fig. 16(g).
The vertical profiles of the same fields evaluated for r/R = 0.98, Figs.17(a-g) further characterize interactions of the boundary flow and the bulk interior flow.The vertically-coherent structures of T and ω z in (a,b), respectively, at Ra = 10 9 have a characteristic length of order 2λ c shown in both plots and the m = 1 BZF remains quite coherent.At Ra = 5 × 10 9 , the BZF is instantaneously less coherent in T (c) but with vertical structures with remaining vertical correlation.In the non-hydrostatic, quasi-geostrophic description [40,41], the state at Ra = 10 9 corresponds to the cellular regime of bulk convection whereas 2 × 10 9 ≤ Ra ≤ 5 × 10 9 corresponds to the plume regime (ending at about Ra/Ra c ≈ 6, i.e., for Ra = 6 × 10 9 for our data at Ek = 10 −6 ).The vertical profiles of the velocity fields (e-g) for Ra = 4 × 10 9 show similar spatial structure with some vertical coherence and a range of horizontal spatial scales.The BZF structure is again most visible in u z (e) compared to u ϕ (f) and u r (g).The spatial structure of u ϕ is similar to u z whereas u r has much finer spatial structures.The u r field in particular shows the strong radially inward and outward exchange with the bulk interior flow.To make these statements involving spatial length scales more quantitative, we compute the different measures of spatial length scale covering the range from the wall modes up to full bulk convection.

E. Length scales of wall modes and BZF
There are multiple different length scales to consider when spanning the regimes of steady wall modes, 3 × 10 7 ≤ Ra ≤ 3 × 10 8 , unsteady sub-critical wall modes 4 × 10 8 ≤ Ra ≤ 8 × 10 8 , and bulk modes 9 × 10 8 ≤ Ra ≤ 5 × 10 9 .For our Γ = 1/2, the largest length scale is the m = 1 wall mode state with λ 1 = π/2 which persists over the full range considered here.The radial length scale of the wall-mode/BZF states can be found in different ways [53,54].Here we take a radial length defined either by the first zero-crossing r 0 of each field in the radial direction starting from the sidewall boundary (r/R = 1), for u r , the radial position of the peak maximum (see Fig. 8(d)), and for ω z we take the second zero crossing.Note that δ u /H = (1−r 0 /R)Γ/2 = (1−r 0 /R)/4.We compute this in the precessing frame at the maximum of the field in the azimuthal direction; the approach we used above for comparing with linear eigenfunctions [20].Previously, we measured the radial length scale of the wall-mode/BZF states by the first zero-crossing of ⟨u ϕ ⟩ ϕ,t at the midplane [53,54].This does not work well for the wall mode state because u ϕ (z) has an approximate 0 at the midplane, a radial length scale δ u ϕ is about 1/2 of its value at, for example, z = 0.8.This explains the apparent scaling ∼ Ek 2/3 of δ 0 found earlier for a range of Ek [53,54] instead of the expected Ek 1/3 scaling.This approach works well even for unsteady precession in the sub-critical regime but there are insufficient data to get statistical convergence for the highly-intermittent BZF state in this study where strong bulk flow is present.In that case, we use the first minimum of the root-mean-square (rms) field ⟨u z ⟩ ϕ,t which should be at the approximate location of the first zero crossing used for the wall-mode regime.The results are shown in Fig. 18 for u z , u ϕ , u r , and ω z where the widths δ are normalized by Ek 1/3 for comparison with the wall-mode linear eigenfunctions.In the weakly nonlinear regime with Ra ≤ 7 × 10 7 , δ u ϕ ≈ 2.25Ek 1/3 consistent with the linear eigenfunction whereas for Ra ≥ 10 8 , we have δ u ϕ ≈ δ uz ≈ 2.5Ek 1/3 .In contrast, δ ωz ≈ 3.3Ek 1/3 ≈ 1.3δ uz for Ra ≤ 10 9 .Both δ u ϕ and δ uz start to increase slightly in the region of subcritical instability.For Γ = 2 (blue, open circle), δ uz ≈ δ u ϕ at Ra = 5 × 10 8 .For Ra ≥ 10 9 , δ uz rms increases rapidly with Ra. δ ur shows more variability, perhaps owing to the smaller amplitude of u r which might require better statistical averaging.
Other length scales are the horizontal length scales in the interior of the cell as measured in horizontal crosssections and the horizontal and vertical length scales in the azimuthal and vertical direction around the circumference in the sidewall boundary region.The first is characteristic of bulk convection whereas the second is reflective of BZF-bulk interactions.For horizontal cross sections, we compute the 2D auto-correlation function of the field F as C F F (r) ≡ ⟨F (x)F (x + r)⟩ A,t /⟨F 2 (x)⟩ A,t in a centered rectangular domain with dimensions 0.6D × 0.6D, where ⟨•⟩ A,t means averaging over a horizontal cross-section and in time.We define the length scale δ such that C F F (δ) = 0.25 with uncertainty defined by dδ = ±(dC(r)/dr) −1 ∆C = ±0.05(dC(r)/dr)−1 (as compared to [57] who use δ = ∞ 0 C(r)dr so our lengths are a bit larger).In Fig. 19(a), we show δ T Ek −1/3 versus Ra with corresponding representative images (left) in (b) for (1) Ra = 2 × 10 8 , (2) Ra = 1 × 10 9 , and (3) Ra = 5 × 10 9 .The size of the interior domain of area 0.36D 2 is shown.In the steady wall-mode region, δ T ≈ 13 which corresponds roughly to the exponential decay shown in Fig. 7(a).In the unsteady subcritical regime, δ T decreases owing to an increasing fraction of lateral spatial structures with smaller length scale.In the bulk regime for Ra ≥ 1 × 10 9 , δ T Ek −1/3 ≈ λ c Ek −1/3 = 5.1 characteristic of bulk rotating convection.We also compute the 1D correlation functions for the anisotropic directions in a vertical azimuthal slice T (r/R = 0.98, ϕ, z): C T T (Rϕ) = ⟨⟨T (r/R = 0.98, θ, z)T (r/R = 0.98, (θ + ϕ, z))⟩ z /⟨T 2 ⟩ z ⟩ t and C T T (z) = ⟨⟨T (r/R = 0.98, ϕ, y)T (r/R = 0.98, ϕ, y + z)⟩ ϕ /⟨T 2 ⟩ ϕ ⟩ t .The results are shown in Fig. 20(c) with corresponding representative images in (b) (right).The azimuthal length scale is dominated by the m = 1 sinusoidal character of the wall mode and is about 30 up to the onset of bulk convection -a m = 1 cosine function has a correlation length of about 20 in these units.In the bulk regime, it drops to about 17.The vertical length scale is roughly independent of Ra which is consistent with the main vertical variation of T being a linear gradient with a corresponding slope independent correlation length of about 25 (the correlation length of a linear gradient is independent of its slope).Instead, the value is somewhat lower, about 18-20.To obtain a more complete description of the length scales defined above, we apply the auto-correlation analysis to the velocity and vertical vorticity fields for horizontal cross sections and for the horizontal and vertical directions in a azimuthal slice near the sidewall.In particular, we use for the horizontal sections, z = 0.8 for u z , u ϕ , u r , and ω z to avoid the approximate vertical zero crossing of u ϕ at z = 1/2.We show in Figs.20(a-c) the correlation lengths δEk −1/3 versus Ra/10 7 for (a) horizontal cross sections, (b) horizontal variations in a vertical azimuthal profile, and (c) vertical variations in a vertical azimuthal profile.The length scales for horizontal cross sections all show the same overall trend with a decrease from about 10 starting at Ra ≈ 5 × 10 8 to between 2 and 5 corresponding to λ c /2 and λ c , respectively, for higher Ra.Whereas, δ T decreases gradually over the subcritical regime, the velocities and ω z abruptly decrease and maintain roughly constant values through the subcritical and bulk regimes 5 × 10 8 ≤ Ra ≤ 5 × 10 9 .This suggests that the subcritical wall-mode instability is related to the bulk instability with structures of the order of the linear stability length scale.The slower decline of δ T is related to the longer (exponential) width of the temperature field relative to the radially-confined velocity and ω z fields.
The horizontal and vertical lengths in the sidewall boundary region, Figs.20 (b,c), are more complicated to interpret.The length scales of T , u and ω z are consistently about 28 for Ra ≤ Ra c ≈ 1 × 10 9 with the exception of δ ur ≈ 8 in the subcritical regime.For 1 × 10 9 ≤ Ra ≤ 5 × 10 9 , δ T ≈ δ uz ≈ δ u ϕ ≈ δ ωz ∼ F (Ra/Ra c − 1) (assuming a dependence on Ra/Ra c − 1).An empirical approximation to the data is F (x) ∼ x 1/3 .The exception is u r ≈ 1.5 which reflects the small scale nature of the boundary-interior interactions, see Fig. 17 (g).The vertical correlations are somewhat different although δ uz and δ u ϕ are reasonably described by F (x) ∼ x γ with γ = 1.2, 0.7, respectively.Comparatively, δ ur starts to decrease in the subcritical regime, similar to the horizontal component.The quantitative trends shown here are consistent with the qualitative impressions provided by representative images in Figs. 17.

F. Heat transport Nu
There are heat transport contributions from the wall modes and their remnants in the BZF and from the bulk flow.As we have seen above, the two states are not independent.The wall-mode jet instability for Ra ≥ 4 × 10 8 injects temperature and velocity perturbations into the interior prior to and continuing into the region of bulk instability, see Figs. 16.In the other direction, the finite convective amplitude near the sidewall boundary appears to be the nucleation site for bulk-like perturbations with spatial separation of order λ c and with strong vertical coherence, see Figs. 17.Thus, cleanly separating the individual contributions, especially for convections cells with Γ < 1, is a  It appears from the analysis presented in Figs.21(a, b) that the lateral jets from the wall region contribute substantially to the heat transport at the onset of bulk convection at Ra ≈ 1 × 10 9 .With that in mind, we look at the average Nu and consider how best to separate the wall and bulk contributions.Previously, we have done so in several ways.First, owing to the rather unique radial structure of the wall modes and their remnants, we defined the BZF contribution as 2π 1 r0 Nu(r)rdr where r 0 is the first zero crossing of Nu(r) [53,54].But this misses the negative contribution in Nu(r) as discussed above, see Figs. 11(a, b).Recently, we generalized this approach by dividing the total Nu into inner and outer portions (2/r 2 s ) rs 0 Nu(r)rdr and 2/(1 − r 2 s ) 1 rs Nu(r)rdr, respectively, where r s was varied [50].Here we choose r s = 0.8 (in units of R) which is a good approximation to where the radial signature of the wall-mode/BZF seems to have become small, see Fig. 11(b).So from this perspective there are 3 contributions to the heat transport: pure wall mode (with jet instability), BZF, and bulk.In Fig. 22(a), these different pieces are plotted versus Ra/Ra c .For r s = 0.8, Nu contributions of BZF and bulk are about the same and both vary quite linearly with Ra/Ra c .For comparison, Nu from the non-hydrostatic quasi-geostrophic model (NHQGS) [40,41] for Ek → 0 and from Nu from a laterally-periodic DNS for Pr = 1 and Ek = 10 −7 [42] which includes Ekman pumping contributions (see also [62,63]) are plotted.The NHQGS is smaller owing to ignoring the contributions of Ekman pumping whereas the DNS is a bit higher, perhaps owing to larger Ek, i.e., Ek > 0. Another possible separation strategy is to define r s (Ra) based on some feature of, for example, u z (r) or N u(r) rather than having it be a constant value.Varying r s slightly for different Ra introduces a net upward (downward) curvature to Nu versus Ra for BZF (bulk) pieces.We choose instead to evaluate the separation for several constant r s , recognizing that there is no unambiguous, testable separation strategy using the data we have here.In Fig. 11(c), the features of Nu, u z , and T suggest 0.7 < r s /R < 0.8.In Fig. 22(b), we show results using r s /R = 0.7 to contrast with r s /R = 0.8 in Fig. 22(a).The trends are both consistent with a linear increase of Nu with Ra/Ra c but with a greater (lesser) slope for BZF (bulk) contributions.Here the bulk contribution is quite comparable to the NHQGS data with little Ekman pumping correction whereas the seeming correspondence of the BZF contribution with the DNS results is purely coincidental.
Although the procedure above is a tempting one that offers a more or less clean separation, there are concerns about this approach that can be elucidated by considering the vertical temperature profile ⟨T (r = r max , ϕ, z)⟩ ϕ,t , where r max is the value of, and its variation with, Ra.In Fig. 23(a), profiles of ⟨T ⟩ are shown as functions of z for 3 × 10 7 ≤ Ra ≤ 7 × 10 8 and in the inset for 10 9 ≤ Ra ≤ 5 × 10 9 .There is a continual increasing slope at the top and bottom boundaries but the interior gradient is not monotonic with a smaller slope at Ra = 5×10 8 compared to that at Ra = 7 × 10 8 .This reversal is also observed in the inset where higher Ra have steeper interior slope.These variations are features associated with the different states of wall modes for Ra ≤ Ra c ≈ 9 × 10 8 and the development of the bulk instability.An important point here is that the slope is continuously varying (see Appendix Fig. 23) so that the definition of a thermal BL thickness is quite problematic.Perhaps there is localized Ekman pumping associated with the BZF -a topic for further investigation.
In Fig. 23(b), we present the slope at the lower boundary −dT /dz as −dT /dz(z = 0, r = r max )−1 for which there are 3 distinct regions: weakly nonlinear growth ∼ (Ra − Ra w ), the nonlinear state including the periodic time-dependent state up to Ra = 6 × 10 8 , and the combined BZF/bulk state with increasing slope ∼ Ra.If one averaged over r as well as over ϕ and t, one would expect an exact correspondence between global heat transport and local temperature slope Nu − 1 = −dT /dz(z = 0, r = r max ) − 1.In Fig. 23(c), there is very close agreement between these two quantities except in the weakly nonlinear regime.Returning to the interior slope, we plot in Fig. 23(b), dT /dz(z = 1/2) + 1 which is 0 at Ra w and would approach 1 with increasing Ra for non-rotating convection but has a finite slope for bulk rotating convection [50].Here we have between 0.8 and 0.9 compared to a value of 0.6 for pure bulk geostrophic convection [40].The difference is not surprising given the coexistence of bulk fluctuations on top of the wall-localized state as previously documented above.
Perhaps the most interesting feature of these profiles of ⟨T ⟩ is the close agreement of the localized (in r) measure of Nu local = −dT /dz(z = 0, r max ) and the global heat transport Nu that averages over the whole domain.If there was a clean separation between wall and bulk modes, Nu local would not increase rapidly and in close correspondence to Nu.One can then infer that the bulk-like fluctuations in the wall zone continue to increase in strength and dominate the heat transport with the average wall mode foundation playing a minor role.One cannot escape the realization that there are strong interactions among the bulk and coexisting BZF that plays a major role in the heat transport in some transition region from wall modes to bulk rotating convection.

III. CONCLUSION
We employed DNS using the Goldfish code of RRBC in a cylindrical geometry with aspect ratio Γ = 1/2, insulating sidewall boundary conditions, Ek = 10 −6 , and Pr = 0.8 to characterize the progression of states over a range 3 × 10 7 ≤ Ra ≤ 5 × 10 9 .We also computed a smaller set of properties for the same parameters but with Γ =1, 2 and for Ek = 10 −4 with Γ = 0.3, 1/2, 1, 2, 5, primarily to determine the mode number m.The set of values for Ek = 10 −6 are tabulated in the Appendix Table II).We elucidated how the mode number m observed for different Γ was consistent between wall modes and the BZF state, how the observed wall mode critical Ra w varies with Ek, Γ, and m, and how small Γ ≲ 0.7 always has m = 1 owing to the azimuthal periodicity.We further demonstrate the decoupling of radial length scales for T ∼ (k w ) −1 and u ∼ Ek −1/3 in the steady wall mode regime 3 × 10 7 ≤ Ra ≤ 3 × 10 8 and show that the eigenfunctions of the planar wall linear solutions [20] provide an excellent representation of the data, particularly for the radial dependence of u(r) where, for Ek = 10 −6 , the radial localization near the sidewall of about 0.1 makes the sidewall curvature effects very small, i.e., 0.1/(2π) ≈ 0.02.On the other hand, for T the linear eigenfunction extends over a significant fraction of R and deviations owing to finite curvature and growing nonlinearity are more keenly felt.Similarly, the azimuthal and vertical eigenfunctions are also more effected by nonlinearity and curvature but give insight into the development of the steady wall mode state.
In the steady wall-mode regime, Nu is constant despite the traveling wave nature of the wall mode.For 4 × 10 8 ≤ Ra ≤ 8 × 10 8 , the steady wall mode undergoes a subcritical bifurcation to a state of time-dependent Nu through a nonlinear mechanism of lateral jet ejection from the wall mode into the bulk interior (see [49], [57]).Throughout the linear (steady) and nonlinear (time dependent) regimes for Ek = 10 −6 , Nu(r) ∼ ⟨u z (r)⟩⟨T (r)⟩ owing to the rapid radial variation of u z compared to the very slow spatial variation of T and the fixed small phase difference between them.The jet instability strengthens as Ra → Ra c ≈ 1 × 10 9 and fine scale spatial (and temporal oscillations, see Appendix) structures form around the azimuth in the wall localized region.We show that the radial length scale δ uz and δ u ϕ defined by first zero-crossing is about 2.5 Ek 1/3 whereas δ ωz ≈ 3.3Ek 1/3 using its second zero crossing.We also use auto-correlation analysis to extract correlation lengths for different cross sections of fields, namely horizontal cross sections and vertical surfaces of {ϕ, z} at constant r/R = 0.98 (the maximum value of u z (r)).The correlation lengths are fairly constant in the wall mode regime and decrease rapidly once the bulk mode sets in for Ra ≳ 10 9 .In horizontal regions in the cell center and Ra ≳ 10 9 , δ T , δ u ϕ ≈ λ c whereas δ uz , δ ur , δ ωz ≈ λ c /2.In the wall localized regime, the azimuthal correlation lengths are are roughly constant at about a half wavelength of the m = 1 wall mode but decrease to 2-4 λ c in the bulk mode region.The vertical correlation lengths are more varied among the different fields but systematically decrease in the bulk mode region.
When considered in the context of heat transport scaling, the analysis of the structure of wall modes and their remnant robust features in the presence of bulk convection -the BZF -play a key role in attempting to separate bulk influence from that of the wall-localized portion.The negative feature in Nu is understood as arising from the oscillatory radial structure of u z (r) and the much slower radial dependence of T .Once bulk convection begins in the interior, the BZF spreads out radially and different features of Nu, u z and T including the lateral jet instability makes separation of the 2 components complicated.The resulting separation using r s /R = 0.8 and 0.7 and subtracting the wall mode contribution [50,53], yields Nu ∼ Ra/Ra c but with different slopes: one that gives agreement with the NHQGS theory and the other with Ekman pumping corrections.There are also subtle changes in the Ra dependence for the two r s values.The average vertical temperature profiles reinforce the understanding that the BZF carries strong bulk fluctuations that contribute heavily to the total heat transport.There seems to be no unambiguous way to disentangle the two components of Nu in small aspect ratios such as studied here.A closer inspection of the interplay of thermal boundary layers in the transitioning wall mode and for bulk convection with coexisting BZF states is certainly warranted.

IV. ACKNOWLEDGEMENT
The authors acknowledge the support from the German Research Foundation (DFG), grants Sh405/20 and Sh405/22, and the Leibniz Supercomputing Centre (LRZ) for providing computing time.REE acknowledges support from the Los Alamos National Laboratory LDRD Program.

V. APPENDIX
In this section, we provide tabulated data from our DNS, empirically expand the set of critical parameters to higher order in Ek, and show some additional features of wall modes in the BZF regime.In addition, we tabulate movies that we present in supplemental material.

B. Wall Mode Critical Parameters and Aspect Ratio Effects
The asymptotic approximations for the critical values of Ra w , ω w , and k w are valid over varying ranges of Ek but the approximations are not adequate to accurately estimate these critical values for larger Ek.Using the numerical data presented in [20], one can empirically extend the asymptotic results to larger Ek using fits to the data and expansions in higher powers of Ek 1/3 .In Fig. 24, we fit the data for the critical values assuming insulating sidewall boundary conditions at Pr = 0.7, 7 (see Fig. 3 [20]).Only Ra w has significant variation with Pr.Finite wall curvature  (Γ = 1 [16]) is shown to slightly decrease Ra w whereas finite wall conductivity slight increases Ra w [19] as shown in Fig. 24(a).We also fit the curvature parameter ξ 0 of the marginal stability boundary computed for Pr = 7 and for partially insulating sidewalls [21].The results are shown in Fig. 24.
with the resulting fit parameters:

C. Transients
Whereas the form of the time dependence of Nu is very precisely sinusoidal (e.g., Fig. 25(b), the subcritical nonlinear state has some interesting additional features.At Ra = 4 × 10 8 , there is a very clear large/small Nu asymmetry in Nu(t) as shown in Fig. 25(a), i.e., the system prefers to be in a larger Nu state relative to a smaller one.The solid line is a fit with a − b[cos 4 ((ω 2 /2)t − ϕ 0 )] to the normalized Nu(t) where Nu n is its approximate mean value.In Fig. 13(c) (Ra = 5 × 10 8 ), there is some remaining asymmetry with both cos (ωt) and cos ((ω 2 /2)t) 4 contributions.Whereas the amplitude ∆Nu at Ra = 4 × 10 8 is very constant over several oscillation periods, the amplitude increases slowly for Ra = 6 × 10 8 as a saturation of the form 1 − e −t/τ0 ); this is also observed for Ra = 5 × 10 8 , see Fig. 13(c), but the single oscillation period makes the saturation difficult to measure accurately.The time series for Ra = 7 × 10 8 in Fig. 25(c) is quite chaotic with an aperiodic oscillation and an intermittent amplitude.

D. Γ = 1/3 prograde dynamics
Most of the BZF states that we explored had retrograde (anti-cyclonic) precession.The exception is for Γ = 1/3 with Ra = 10 8 , and Ek = 1.1 × 10 −5 , where we find prograde (cyclonic) precession and alternations between prograde and retrograde directions, see Figs. 26 (a,b).At early t ≲ 400 (Fig. 26(b) shows an expanded interval) and late times t ≳ 700, the precession is prograde but for intermediate times, an interval of retrograde precession is observed.It is not determined from our data whether the prograde state is stable at long times.Uniform cyclonic precession dynamics

F. Movies
The dynamics and structure of the wall modes and the bulk/BZF states are hard to fully characterize through static single frame images.In the Supplementary Material, we provide representative movies that illustrate a variety

FIG. 1 .
FIG.1.Schematic phase diagram of states of rotating convection (perfectly insulating sidewall boundary conditions): RaEk 4/3 /A0 (dashed black line) where A0 = 8.7 is the asymptotic Ek → 0 value versus Ek.The lines for Rac (dashed blue line)[4,22,23] and Raw (solid blue line)[20] are from linear stability analysis with idealized boundary conditions on laterally infinite domains.Rag marks the approximate boundary of the geostrophic, rotation-dominated region and Rat denotes the transition to buoyancy-dominated flow[1,2,24].The precise details of this diagram depend on parameters including Pr and Γ and are not included here.The data points show the parameter values reported in this paper.

6 FIG. 6 .
FIG. 6. Linear eigenfunctions (discretized from data in Figures 4 and 5 of [20] for Ek = 10 −6 and Pr = 7): (a) Eigenfunction of uz versus Ek −1/3 x, where x reflects the distance from a flat sidewall.Characteristic lengths are indicated by the location of the peak δu zp , the first zero crossing δu z , and the minimum δu zm .(b) Eigenfunction of T versus x.The dashed line is an exponential fit to the data 1.05e −x/0.13 that fits extremely well over the whole range except very near the sidewall on the order of the uz radial width (see Inset where short red (long) dashed vertical line is rmax (r0) of uz).(c) The geometric length scale LR = R = Γ/2 is plotted versus values of Γ considered (solid black circles) and compared to the temperature length scale δT and the vertical velocity length scale δu z evaluated at Ek = 10 −4 (dashed blue) and Ek = 10 −6 (solid blue).

10 FIG. 8 .
FIG. 8. Mid-plane uz (averaged in precessing frame at maximum in ϕ) versus r/R for (a) Ra/10 7 = 3, 4, 5, 7, 10, 20, 30 where dark (light) blue corresponds to the smallest (largest) Ra.Inset of (a) shows Fourier-Bessel coefficients cn, normalized by cn max (Ra) and averaged over the set of Ra in (a).(b) Expanded version of (a) with 0.7 ≤ r/R ≤ 1.0.Inset of (b) is uz(r) reconstructed from the average of cn/cn max in inset of (a).(c) uz versus r/R of time-dependent wall-mode state for 4 × 10 8 ≤ Ra ≤ 9 × 10 8 (colors per legend) with saturating peak amplitude and wider spread in δu z .(d) Ra = 3 × 10 7 average radial structure (solid line) compared to linear eigenfunctions (dashed)[20] for u ϕ (red) uz (blue), and ur (black).The first zero crossing of uz is at about the same r as the maximum of ur.All profiles are scaled such that the maximum values are 1.(e) Variation of radial maximum values uz max (blue solid circles), uz rms,max (light blue inverted triangles), u ϕmax (red solid squares), and ur max (black solid triangles) vs Ra − Raw.Linear scaling u ∼ (Ra − Raw) 1/2 (blue, red, and black dashed lines, respectively).Vertical dashed lines from smaller to larger Ra: end of quasi-linear regime (long-dashed) Ra ≈ 9 × 10 7 -note decrease in u ϕ , subcritical instability to lateral jet ejection (dashed), transition to aperiodic (chaotic) jet ejection (long-dashed), and bulk onset (dashed).

2 ω d 2 1 0∼ (Raw 2 −
FIG. 14.(a) Extreme maximum and minimum values of Nu(t) (minima: black solid circles, maxima: blue solid circles and solid red square for chaotic state), where the error bars depict the range of less extreme max/min values, (b) ∆Nu (mean peak-to-peak), and (c) normalized oscillation frequency ω d 2 /ω d (solid circles: black -stable state and blue -unstable state) and normalized decay frequency (scaled by 1/2 for comparison) (1/2)(2π/τ d 0 )/ω d (magenta diamonds) vs Ra/10 8 where ω d is the wall mode precession frequency.Error bars denote variability of ∆Nu, ω d 2 /ω d , and (2π/τ d 0 )/ω d owing to shorter or longer lengths of time series and/or unsteady oscillations.The blue shaded region denotes an approximate zone of subcritical instability.The solid (dashed) lines are schematic suggestions for stable (unstable) behavior.The square (red) indicates a state where the oscillations have become chaotic with increased fluctuations and variable frequency.The dashed magenta line represents the stable state transient inverse time constant which is expected to vary as τ −1 0 ∼ (Raw 2 − Ra).

FIG. 15 .
FIG.15.For Ra = 5 × 10 8 , Ek = 10 −6 , and Pr = 0.8: horizontal cross sections of (a) T (r, ϕ) and (b-d) uz (dashed lines are at radii r/R = 0.98 and 0.90, corresponding to the peak and the first zero crossing, respectively), u ϕ , and ur, respectively, for z = 0.8.Red (blue) corresponds to hot (cold) and up (down), AC -anticyclonic (C -cyclonic), and radially out (radially in), respectively.Vertical profiles (e-h) for corresponding fields (ϕ, z) evaluated at r/R = 0.98 (see (b) outer dashed line).Arrows and ± symbols indicate qualitative motion directions.Spatial oscillations of order 2λc are visible in horizontal cross-sections of uz and ur and in vertical azimuthal profiles of uz, u ϕ , and ur.Note that the major ejections of jets (and return flow) represented in ur happen near the top and bottom boundaries.
FIG. 21.(a) Numax (blue solid circles, red solid square and diamonds) and Numin (black, solid circles and triangles) vs. Ra/10 8 .Regions are labeled with vertical dashed lines separating them.Red and black dashed lines are linear fits (of the form a + bRa) to Nu vs. Ra for Ra ≥ 10 9 .(b) T fields for (top) Ra = 5 × 10 8 and (bottom) Ra = 10 9 showing example at the minimum (left image) and maximum (right image) of Nu.Color bar of T is shown.

Fig. 26 (
Fig. 26 (b) the boundary layer profiles dT /dz(z) near the bottom boundary in the bulk/BZF region 1 ≤ Ra/10 9 ≤ 5.The curve (blue) for Ra = 1 × 10 9 has a form still dominated by the wall-mode whereas the higher Ra curves fall off much faster.For comparison the the boundary layer profiles of dT /dz for non-rotating convection for the same total Nu are shown for Ra/10 9 = 4 and 5.The RRBC profiles are quite different than those of non-rotating RBC.
FIG. 28.(a) ⟨⟨T (z)⟩ ϕ,t ⟩rms (blue, solid circles) and ⟨⟨T (z)⟩ ϕ,t − T (1) (z)⟩rms (red, solid squares) vs. Ra − Raw; T (1)(z) is the first Fourier component of ⟨T (z)⟩ ϕ,t .Inset shows the ratio of ⟨T − T (1)(z) ⟩rms/⟨T ⟩rms.Vertical dashed lines show boundaries of different regions.(b) (dT /dz)/(dT /dz(z = 0)) vs. z for values of Ra/10 7 in the range of bulk rotating convection: 100, 200, 300, 400, 500.There is no clear demarcation of a thermal boundary layer with an approximately constant slope only for z values within an Ekman boundary layer δ Ek = Ek 1/2 = 0.001 indicated by dashed black vertical line.The dashed orange and black curves are the corresponding profiles for non-rotating convection with the same Nu as for the corresponding solid orange and black curves.

TABLE I .
Data from DNS with Pr = 0.8 showing Ek, Γ, Raw, Row, and m.