Convective scale and subadiabatic layers in simulations of rotating compressible convection

(abridged) Context: Rotation is thought to influence the size of convective eddies and the efficiency of convective energy transport in the deep convection zones of stars. Rotationally constrained convection has been invoked to explain the lack of large-scale power in observations of solar flows. Aims: The main aims are to quantify the effects of rotation on the scale of convective eddies and velocity, the depths of convective overshoot, and the subadiabatic Deardorff layers. Methods: Three-dimensional hydrodynamic simulations of rotating convection in Cartesian domains were run. The results were compared with theoretical scaling results that assume a balance between Coriolis, inertial, and buoyancy (Archemedean) forces (CIA balance). Results: The scale of convective eddies decreases as rotation increases, and ultimately reaches a rotationally constrained regime consistent with the CIA balance. Using a new measure of the rotational influence on the system, it is shown that even the deep parts of the solar convection zone are not in the rotationally constrained regime. The simulations capture the slowly and rapidly rotating scaling laws predicted by theory, and the Sun appears to be in between these two regimes. Both, the overshooting depth and the extent of the Deardorff layer, decrease as rotation becomes more rapid. For sufficiently rapid rotation the Deardorff layer is absent. Conclusions: Relating the simulations with the Sun suggests that the convective scale even in the deep parts of the Sun is only mildly affected by rotation and that some other mechanism is needed to explain the lack of strong large-scale flows in the Sun. Taking the current results at face value, the overshoot and Deardorff layers are estimated to span roughly five per cent of the pressure scale height at the base of the convection zone in the Sun.


Introduction
The theoretical understanding of solar and stellar convection was shaken roughly a decade ago when helioseismic analysis suggested that the velocity amplitudes in the deep solar convection zone are orders of magnitude lower than anticipated from theoretical and numerical models (Hanasoge et al. 2012).Significant effort has been put into refining these estimates, but a gaping discrepancy between numerical models and helioseismology remains (e.g., Hanasoge et al. 2016;Proxauf 2021); (see, however, Greer et al. 2015, whose results are at odds with this).This issue is now referred to as the "convective conundrum" (O'Mara et al. 2016).
Several solutions to this conundrum have been proposed, including a high effective Prandtl number (e.g., Karak et al. 2018), rotationally constrained convection (Featherstone & Hindman 2016), and the superadiabatic layer in the Sun being much thinner than thought thus far (Brandenburg 2016); see also Käpylä et al. (2023) and references therein.The two latter ideas are explored further in this current study.The idea that convection is rotationally constrained in the deep convection zone (CZ) has already been borne out through mixing length models of solar convection that imply velocity amplitudes u conv of about 10 m s −1 for the deep convection zone, while the convective length scale conv , which is the mixing length, is of the order of 100 Mm, yielding a Coriolis number Co = 2Ω conv /u conv on the order of ten (e.g., Ossendrijver 2003;Schumacher & Sreenivasan 2020).However, this estimate does not take into account the decreasing length scale due to the rotational influence on convection.Assuming that the Coriolis, inertial, and buoyancy (Archimedean) forces balance, also known as the CIA balance (e.g., Stevenson 1979;Ingersoll & Pollard 1982;King & Buffett 2013;Barker et al. 2014;Aurnou et al. 2020;Vasil et al. 2021;Fuentes et al. 2023), implies that the convective scale is given by conv ∝ Co −1/2 , where Co = 2ΩH/u conv is a global Coriolis number and where Ω is the rotation rate and H is a length scale corresponding to the system size (e.g., Aurnou et al. 2020).This idea has been explored recently by Featherstone & Hindman (2016) and Vasil et al. (2021), who suggested that the largest convectively driven scale in the Sun coincides with that of supergranulation due to rotationally constrained convection in the deep CZ.These studies assumed from the outset that convection is strongly rotationally affected.In this work, a somewhat different perspective is taken in that an attempt is made to independently assess whether this assumption holds for the deep solar CZ.Furthermore, in addition to conv , the scalings of various quantities based on predictions from the CIA balance are studied over a wide range of rotation rates.
Simulations of stratified overshooting convection have revealed that deep parts of CZs are often weakly stably stratified (e.g., Roxburgh & Simmons 1993;Tremblay et al. 2015;Hotta 2017;Bekki et al. 2017;Käpylä et al. 2017;Käpylä 2019).This is interpreted such that convection is driven by cooling at the surface that induces cool downflow plumes that pierce through the entire convection zone and penetrate deep into the stable layers below.This process has been named entropy rain (e.g., Brandenburg 2016) and goes back to ideas presented by Spruit (1997) and the simulations of Stein andNordlund (e.g., Stein &Nordlund 1989, 1998).This picture of convection is a clean break from the canonical view in which convection is driven locally throughout the convection zone by a superadiabatic temperature gradient, an idea which is also encoded into the mixing length concept (e.g., Vitense 1953;Böhm-Vitense 1958).Theoretically, this can be understood such that the convective energy flux that is traditionally proportional to the entropy gradient is supplemented by a non-gradient term proportional to the variance of entropy fluctuations (Deardorff 1961(Deardorff , 1966)).
Analysis of the force balance of upflows and downflows in non-rotating hydrodynamic simulations supports the idea of surface-driven non-local convection (e.g., Käpylä et al. 2017;Käpylä 2019Käpylä , 2021)).Thus far, these studies have mostly concentrated on non-rotating convection (see, however Käpylä et al. 2019;Viviani & Käpylä 2021).In our work, rotation is included to study its impact on the formation and extent of stably stratified Deardorff layers where the convective flux runs counter to the entropy gradient.Another related aspect of interest in astrophysics is convective overshooting (see, e.g., Anders & Pedersen 2023, for a recent review).Numerical studies specifically targeting overshooting have largely concentrated on non-rotating cases (e.g., Singh et al. 1995Singh et al. , 1998;;Saikia et al. 2000;Brummell et al. 2002;Hotta 2017;Käpylä 2019;Anders et al. 2022), and the effects of rotation have received much less attention (e.g., Ziegler & Rüdiger 2003;Käpylä et al. 2004;Brun et al. 2017).It is generally thought that rotation leads to a reduction of overshooting depth (e.g., Ziegler & Rüdiger 2003), but a comprehensive study of this has hitherto been lacking.
The remainder of the paper is organised as follows.The model is described in Sect.2, whereas the results and conclusions of the study are presented in Sects.3 and 4, respectively.The derivations related to the CIA balance are presented in Appendix A.

The model
The model is the same as that used in Käpylä (2019Käpylä ( , 2021)).The Pencil Code (Pencil Code Collaboration 2021) 1 was used to produce the simulations.Convection was modelled in a Cartesian box with dimensions (L x , L y , L z ) = (4, 4, 1.5)d, where d is the depth of the initially convectively unstable layer.The equations for compressible hydrodynamics were solved: 1 https://github.com/pencil-code/ where D/Dt = ∂/∂t + u • ∇ is the advective derivative, ρ is the density, u is the velocity, g = −gê z is the acceleration due to gravity with g > 0, p is the pressure, T is the temperature, s is the specific entropy, ν is the constant kinematic viscosity, and Ω = Ω 0 (− sin θ, 0, cos θ) T is the rotation vector, where θ is the colatitude.The terms F rad and F SGS are the radiative and turbulent subgrid scale (SGS) fluxes, respectively, and C describes cooling near the surface.The term S is the traceless rate-of-strain tensor with The gas is assumed to be optically thick and fully ionised where radiation is modelled via the diffusion approximation.The ideal gas equation of state p = (c P − c V )ρT = RρT applies, where R is the gas constant, and c P and c V are the specific heats at constant pressure and volume, respectively.The radiative flux is given by where K is the radiative heat conductivity where σ SB is the Stefan-Boltzmann constant and κ is the opacity.Assuming that the opacity is a power law of the form κ = κ 0 (ρ/ρ 0 ) a (T/T 0 ) b , where ρ 0 and T 0 are reference values of density and temperature, the heat conductivity is given by The choice a = 1 and b = −7/2 corresponds to Kramers opacity law (Weiss et al. 2004), which was first used in convection simulations by Edwards (1990) and Brandenburg et al. (2000).
Additional turbulent SGS diffusivity is applied for the entropy fluctuations with where s (x) = s(x) − s(z), with the overbar indicating horizontal averaging.The coefficient χ SGS is constant in the whole domain, and F SGS has a negligible contribution to the net energy flux such that F SGS ≈ 0. The cooling at the surface is described by where τ cool is a cooling time, T = e/c V is the temperature, e is the internal energy, and T cool = T top is a reference temperature corresponding to the fixed value at the top boundary.The advective terms in Eqs. ( 1) and (3) are written in terms of a fifth-order upwinding derivative with a hyperdiffusive sixthorder correction with a local flow-dependent diffusion coefficient (see Appendix B of Dobler et al. 2006).

Geometry, initial, and boundary conditions
The  1).The latter corresponds to a marginally stable isentropic stratification.Initially, the uppermost layer above z/d = 1 is isothermal, mimicking a photosphere where radiative cooling is efficient.Convection ensues because the system is not in thermal equilibrium due to the cooling near the surface and due to the inefficient radiative diffusion in the layers above z/d = 0.The velocity field is initially seeded with small-scale Gaussian noise with amplitude 10 −5 dg.
The horizontal boundaries are periodic, and the vertical boundaries are impenetrable and stress free according to A constant energy flux is imposed at the lower boundary by setting where F bot is the fixed input flux and K bot = K(x, y, z bot ).A constant temperature T = T top is imposed on the upper vertical boundary.

Units and control parameters
The units of length, time, density, and entropy are given by where ρ 0 is the initial value of density at z = z top .The models are fully defined by choosing the values of ν, Ω 0 , θ, g, a, b, K 0 , ρ 0 , T 0 , τ cool , and the SGS Prandtl number along with the cooling profile f cool (z).The values of K 0 , ρ 0 , T 0 are subsumed into another constant K 0 = K 0 ρ a+1 0 T b−3 0 , which is fixed by assuming the radiative flux at z bot to equal F bot at t = 0.The cooling profile f cool (z) = 1 above z/d = 1 and f cool (z) = 0 below z/d = 1, connecting smoothly across the interface over a width of 0.025d.The quantity ξ 0 = H top p /d = RT top /gd sets the initial pressure scale height at the surface and determines the initial density stratification.All of the current simulations have ξ 0 = 0.054.
The Prandtl number based on the radiative heat conductivity is where χ(x) = K(x)/c P ρ(x) quantifies the relative importance of viscous to temperature diffusion.Unlike many other simulations, Pr is not an input parameter because of the non-linear dependence of the radiative diffusivity on the ambient thermodynamics.The dimensionless normalised flux is given by where ρ(z bot ) and c s (z bot ) are the density and the sound speed, respectively, at z = z bot at t = 0.At the base of the solar CZ F n ≈ 4 × 10 −11 (e.g., Brandenburg et al. 2005), whereas in the current fully compressible simulations, several orders of magnitude larger values are used (see Table 1).
The effect of rotation is quantified by the Taylor number which is related to the Ekman number via Ek = Ta −1/2 .The Rayleigh number based on the energy flux is given by This can be used to construct a flux-based diffusionfree modified Rayleigh number (e.g., Christensen 2002;Christensen & Aubert 2006) In the current setup, Ra F is given by A reference depth needs to be chosen because Ra F = Ra F (z). Furthermore, H ≡ c P T/g is a length scale related to the pressure scale height.The choice d = H = H p , where H p ≡ −(∂ ln p/∂z) −1 is the pressure scale height at the base of the convection zone, leads to A221, page 3 of 16 Notes.Summary of the runs.Pr SGS = 1 in all runs such that Pe SGS = Re.Runs with rotation were made with θ = 0, corresponding to the north pole of the star.
, where F A n = 4.6 × 10 −6 is the normalised flux in the runs in Set A.

Diagnostics quantities
The global Reynolds and SGS Péclet numbers describe the strength of advection versus viscosity and SGS diffusion where u rms is the volume-averaged rms velocity and where k 1 = 2π/d is an estimate of the largest eddies in the system.The Reynolds and Péclet number based on the actual convective length scale are given by Here, = k −1 mean is chosen, where k mean = k mean (z) is the mean wavenumber (e.g., Christensen & Aubert 2006;Schrinner et al. 2012) and is computed from where E(k, z) is the power spectrum of the velocity field with u 2 (z) = E(k, z)dk.
In general, the total thermal diffusivity is given by A221, page 4 of 16 However, in all of the current simulations, χ χ SGS in the CZ such that the Prandtl and Péclet numbers based on χ eff differ little from Pr SGS and Pe SGS .The Rayleigh number is defined as which varies as a function of height and is quoted near the surface at z/d = 0.85.The Rayleigh number in the hydrostatic, non-convecting, state is measured from a one-dimensional model that was run to thermal equilibrium and where the convectively unstable layer is confined to the near-surface layers (Brandenburg 2016); see also Fig. 1.In the hydrostatic case, χ = χ(z), and χ SGS , which affects only the fluctuations, plays no role.The turbulent Rayleigh number is quoted from the statistically stationary state using the horizontally averaged mean state, where the overbar denotes temporal and horizontal averaging.
The rotational influence on the flow is measured by several versions of the Coriolis number.First, the global Coriolis number is defined as where k 1 = 2π/d is the wavenumber corresponding to the system scale.This definition neglects the changing length scale as a function of rotation and overestimates the rotational influence when rotation is rapid and the convective scale is smaller.A definition that takes the changing length scale into account is given by the vorticity-based Coriolis number where ω rms is the volume-averaged rms value of the vorticity ω = ∇ × u.Another definition of the Coriolis number taking into account the changing length scale is given by where = (k mean ) −1 where the overbar denotes averaging over time and CZ.This is a commonly used choice in simulations of convection in spherical shells (Schrinner et al. 2012;Gastine et al. 2014; see also Aurnou et al. 2020, who considered convection in the limits of slow rotation and rapid rotation).
We further define a flux Coriolis number Co F2 as where u is a reference velocity obtained from where ρ is a reference density, taken here at the bottom of the CZ.The u does not, and does not need to, correspond to any actual velocity, and it rather represents the available energy flux.Therefore, Co F does not depend on any dynamical flow speed or length scale set by complicated interactions of convection, rotation, magnetism, or other relevant physics.On the other hand, Co F depends only on quantities that can either be measured (F tot , Ω 0 ) or deduced from stellar structure models with relatively little ambiguity (H p , ρ ).The significance of Co F is seen when rearranging Eq. ( 20) to yield Identifying the lhs with Co 3 F , Eq. ( 30), gives An often-used phrase in the context of convection simulations targeting the Sun is that while all the other system parameters are beyond the reach of current simulations, the rotational influence on the flow can be reproduced (e.g., Käpylä et al. 2023).Equation ( 33) gives this a more precise meaning in that the solar value of Ra F needs to be matched by any simulation claiming to model the Sun.
The net vertical energy flux consists of contributions due to radiative diffusion, enthalpy flux, kinetic energy flux, and viscous flux, as well as the surface cooling: Here, the primes denote fluctuations and the overbars indicate horizontal averaging.The total convected flux (Cattaneo et al. 1991) is the sum of the enthalpy and kinetic energy fluxes: which corresponds to the convective flux in, for example, mixing length models of convection.Another useful diagnostic is buoyancy, or Brunt-Väisälä frequency, which is given by and describes the stability of an atmosphere with respect to buoyancy fluctuations if N 2 > 0. Finally, the Richardson number related to rotation in the stably stratified layers is defined as Averages denoted by overbars are typically taken over the horizontal directions and time, unless specifically stated otherwise.

Results
Three sets of simulations with varying F n and approximately the same values of Co F are presented in the current study.We refer to these as Sets A, B, and C. The non-rotating runs in these sets correspond respectively to Runs K3, K4, and K5 in A221, page 5 of 16 Käpylä (2019) in terms of F n , although lower values of ν and χ SGS were used in the runs of the present study.We note that when F n is varied between the sets of simulations, the rotation rate Ω 0 and the diffusivities ν and χ SGS are varied at the same time proportionally to F 1/3 n (see, e.g., Käpylä et al. 2020, and Appendix A for more details).Furthermore, the cooling time τ cool is varied proportional to F n .The current simulations have modest Reynolds and Péclet numbers in comparison to astrophysically relevant parameter regimes (e.g., Ossendrijver 2003;Kupka & Muthsam 2017;Käpylä et al. 2023; see Table 1).Earlier studies from non-rotating convection suggest that results obtained at such modestly turbulent regimes remain robust also at the highest resolutions affordable (Käpylä 2021).This is due to the fact that the main energy transport mechanism (convection) and the main driver of convection (surface cooling) are not directly coupled to the diffusivities.However, the current cases with rotation are more complicated because the supercriticality of convection decreases with increasing rotation rate (e.g., Chandrasekhar 1961;Roberts 1968).The effects of decreasing supercriticality are not studied systematically here, but subsets of the runs in Set A were repeated with higher resolutions (576 3 and 1152 3 ) and correspondingly higher Ra F , Re, and Pe (see Sets Am and Ah in Table 1).

Hydrostatic solution
Earlier studies have shown that a purely radiative hydrostatic solution with the Kramers opacity law is a polytrope with index n rad = 3.25 (Barekat & Brandenburg 2014;Brandenburg 2016).Such a solution arises in the case where K = const.and ∇ z T = const.To see if this configuration is recovered with the current setup, Eqs. ( 1) and (3) were solved numerically in a onedimensional z-dependent model with otherwise the same parameters as in the 3D simulations corresponding to the runs in Set A. The resulting temperature profile is shown in Fig. 1a along with corresponding horizontally averaged profiles from convecting Runs A0, A6, and A9.The stratification is consistent with a polytrope corresponding to n rad up to a height of roughly z/d = 0.75.Near the nominal surface of the convection zone, z/d = 1, the temperature gradient steepens sharply because the cooling term relaxes the temperature toward a constant (z-independent) value near the surface.Therefore, neither K nor ∇ z T are constants in this transition region between the radiative and the cooling layers.In the initial state the stratification is isothermal above z/d = 1, but because the cooling profile f cool has a finite width, cooling also occurs below z/d = 1, and the isothermal layer is wider in the final thermally saturated state.This also depends on the value of τ cool .In the convective runs, the stratification is nearly piecewise polytropic, with index n rad near the base of the radiative layer and nearly isentropic with n ad in the bulk of the CZ.
The superadiabatic temperature gradient is defined as where ∇ = d ln T /d ln p is the logarithmic temperature gradient and ∇ ad = 1 − 1/γ is the corresponding adiabatic gradient.
Comparison of the hydrostatic profile and the non-rotating convective model A0 shows that the convectively unstable layer in the former is much thinner than in the latter.This is a direct consequence of the strong temperature and density dependence of the Kramers opacity law.A similar conclusion applies also to the Sun, where the hypothetical non-convecting hydrostatic equilibrium solution has a very thin superadiabatic layer (Brandenburg 2016).The steepness of the temperature gradient near the surface is characterised by the maximum value of (∆∇) hyd , which is 23.4.By comparison, in the convective Run A0, max(∆∇) = 0.12.The Rayleigh number, measured at z/d = 0.85, in the hydrostatic case is Ra = 5.4 × 10 7 , which is about an order of magnitude greater than Ra t in Run A0.

Qualitative flow characteristic as a function of rotation
Figure 2 shows representative flow fields from runs with slow, intermediate, and rapid rotation corresponding to Coriolis numbers 0.13, 1.3, and 16.5, respectively.The effects of rotation are hardly discernible in the slowly rotating case A2, with Co = 0.13.In the run with intermediate rotation, Run A6 with Co = 1.3, the convection cells are somewhat smaller than in the slowly rotating case, and more vortical structures are visible near the surface.For the most rapidly rotating case, Run A9 with Co = 16.5, the size of the convection cells is drastically reduced in comparison to the other two runs, and clear alignment of the convection cells with the rotation vector can be seen.The simulations were run either from the initial conditions described in Sect.2.1 or by continuing from snapshots of already thermally relaxed runs.In a few cases, both approaches were explored with the same system parameters.In all of these cases, the solutions converged to the same final state, suggesting independence from initial conditions.

Convective scale as function of rotation
The power spectra E(k) of the velocity fields for the runs in Set A are shown in Fig. 3 from depths near the surface, at the middle, and near the base of the CZ.As was already evident from visual inspection of the flow fields, the dominant scale of the flow decreases as the rotation rate increases.Quantitatively, the wavenumber k max , where E(k) has its maximum, increases roughly in proportion to Co 1/2 .The mean wavenumber k mean , computed from Eq. ( 23), shows the same scaling for Co 2. This is explained by the broader distribution of power at different wavenumbers at slow rotation in comparison to the rapid rotation cases where fewer, or just a single, convective modes are dominant (see Fig. 3).
A decreasing length scale of the onset of linear instability under the influence of rotation was derived in Chandrasekhar (1961) with k onset ∝ Ta 2/3 .With Ta ∝ Co 2 Re 2 and with Re being approximately constant, k onset ∝ Co 1/3 is obtained.On the other hand, considering the CI part of the CIA balance in the Navier-Stokes equation gives (e.g,.Aurnou et al. 2020, see also Eq. (A.10)) gives or k max ∝ Co 1/2 .This is consistent with the current simulations (see the inset in the left panel of Fig. 3).The same result was obtained in Featherstone & Hindman (2016).Some non-linear convection simulations show scalings that are similar but somewhat shallower than that obtained from the CIA balance (see, e.g., Viviani et al. 2018;Currie et al. 2020).
To estimate the convective length scale in the Sun based on the current results requires that the value of Co F matches that of the deep solar CZ.The quantities on the rhs of Eq. ( 30 2016) at comparable values of Co despite the differences between the model setups.However, the values of Co F in Runs A9, B9, and C9 are at least 16 times higher than in the Sun, suggesting that the simulations of Featherstone & Hindman (2016) were also rotating much faster than the Sun 3 .Therefore, the current results suggest that rotationally constrained convection cannot explain the appearance of supergranular scale as the largest convective scale in the Sun.
Figure 4 shows the velocity power spectra for the most rapidly rotating runs with Co ≈ 17 for Re = 30 . . .142 from Runs A9, A9m, and A9h.There is a marked increase in the power at large scales, which begins to affect k mean at the highest Re (run A9h).This is due to the gradual onset of large-scale vorticity production, most likely due to two-dimensionalisation of the turbulence, which has been observed in various earlier studies of rapidly rotating convection (e.g., Chan et al. 2003;Chan 2007;Chan & Mayr 2013;Käpylä et al. 2011;Guervilly et al. 2014).Despite the rapid rotation with Coriolis numbers exceeding 16, the large-scale vorticity generated in the current simulations is relatively modest, apart from run A9h.A difference from many of the previous studies is that here the relevant thermal Prandtl number (Pr SGS ) is on the order of unity, whereas in many of the earlier studies Pr was lower.Large-scale vorticity produc-3 For example, their run with Ro = 0.011 has Ra F = 6.81 × 10 6 , Ek = 1.91 × 10 −4 , and Pr = 1, and corresponds to Ra F = Ra F Ek 3 /(8Pr) = 5.9 × 10 −6 , or Co F = (Ra F ) −1/3 ≈ 55.We note that here Ek = 2Ta −1/2 .This yields Ω/Ω ≈ [(Ra F ) /Ra F ] 1/3 ≈ 17.6, which is approximate because somewhat different length scales were used.
A221, page 7 of 16 tion was indeed observed in an additional run that is otherwise identical to A9 except that Pr SGS = 0.2 instead of Pr SGS = 1.

Velocity-based Co
The suitability of different measures of rotational influence on the flow has been discussed in various works in the literature (e.g., Käpylä 2023b).A common (and justified) critique regarding the global Coriolis number as defined in Eq. ( 27) is that it does not appreciate the fact that conv = conv (Ω) (e.g., Vasil et al. 2021).The most straightforward way to correct this is to measure the mean wavenumber and use Eq.(29). Figure 5 shows Co as a function of Co for all runs listed in Table 1.For slow rotation, Co 1, Co ∝ Co because u conv and conv are almost unaffected by rotation.For sufficiently rapid rotation this is no longer true because k mean ≈ k max ∝ Co 1/2 , as indicated by Eq. (A.10) and the simulation results (see the inset of Fig. 3).This implies that for rapid rotation Co ∝ Co 1/2 ; see also Eq. (A.12).This is consistent with the numerical results found in the most rapidly rotating cases (see Fig. 5).The higher resolution runs in Set Am have somewhat lower Co than the corresponding runs in Set A because the convective velocities in the higher resolution cases are higher.This shows that the simulations are not yet in an asymptotic regime where the results are independent of the diffusivities.This is further demonstrated by the high resolution runs of Set Ah.Run A5h follows the trend set by Run A5m.Run A9h, with a significantly higher Co than in Runs A9 and A9m, is explained by the increasing k mean due to the large-scale vorticity generation in that case.Aurnou et al. (2020) showed that the dynamical Rossby number is related to the diffusion-free modified flux Rayleigh number Ra F , with different powers for slow and rapid rotation.The corresponding derivations for the Coriolis number Co are presented in Appendix A and show that Co = (Ra F ) −1/3 (slow rotation) and Co = (Ra F ) −1/5 (rapid rotation).Both scalings are also supported by the simulation results (see the inset of Fig. 5).

Vorticity-based Co
Another commonly used definition, Eq. ( 28), is used to take the changing length scale automatically into account.However, Co ω comes with a caveat that has apparently not been discussed hitherto in the astrophysical literature.This is demonstrated by considering a set of rotating systems at asymptotically high Re where u rms is independent of Re.The forcing is assumed fixed by a constant energy flux through the system and the asymptotic value of u rms when Re → ∞ as u ∞ .Furthermore, in this regime the mean kinetic energy dissipation rate where the overbar denotes a suitably defined average, tends to a constant value when normalised by mean length and corresponding rms velocity (e.g., Sreenivasan 1984;Vassilicos 2015).This value is denoted as ∞ .In low-Mach number turbulence, which is a good approximation of stellar interiors, as well as the current simulations with Ma ∼ O(10 −2 ), From the definition of system scale Reynolds number, it follows that and from Eq. ( 45) that Using Eq. ( 28), it is found that This means Co ω → 0 as Re → ∞ at constant Co, while the dynamics at large (integral) scales are unaffected.Therefore Co ω underestimates the rotational influence at the mean scale k mean that dominates the dynamics, as opposed to Eq. ( 27) overestimating it.Equation ( 47) can also be written as For sufficiently large Re, the theoretical prediction is that u rms → u ∞ = const.and k ω ∝ Re 1/2 .This has been confirmed from numerical simulations of isotropically forced A221, page 8 of 16 homogeneous turbulence (e.g., Brandenburg & Petrosyan 2012;Candelaresi & Brandenburg 2013).We show the dependence of k ω on Re in the inset of Fig. 6 for runs with Co ≈ 1.3 and Re ranging between 40 and 174.The results for k ω fall somewhat below the theoretical Re 1/2 expectation.This is likely because the asymptotic regime requires still higher Reynolds numbers.
On the other hand, the mean wavenumber k mean is essentially constant around k mean /k 1 = 7 in this range of Re because the dominating contribution to the velocity spectrum comes from large scales that are almost unaffected by the increase in Re.

Convective velocity as a function of total flux and rotation
The scalings of convective velocity as a function of rotation are derived in Appendix A following the same arguments as in Aurnou et al. (2020).For slow rotation, the convective velocity depends only on the energy flux where u is defined via Eq.( 31).This scaling is altered in the rapidly rotating regime where This result agrees with Eq. (50d) of Aurnou et al. (2020) and Table 2 of Vasil et al. (2021).Therefore the velocity amplitude in the rapidly rotating regime is expected to depend not only on the available flux but also on rotation.Figure 7 shows the corresponding numerical results for the Sets A, B, C, and Am.For slow rotation, Co 1.5, u rms is roughly constant around u rms ≈ 1.15u for Sets A, B, and C and around u rms ≈ 1.23u for Set Am.In the rapid rotation regime, u rms follows a trend that is consistent to that indicated in Eq. ( 51), but the agreement is only obtained for the two highest values of Co.The higher resolution and more supercritical runs from Set Am show a similar trend, although the velocities are higher overall.This demonstrates that the current runs are not yet in an asymptotic regime where the effects of the diffusivities are negligible.

Flow statistics
Compressible non-rotating convection is characterised by broad upflows and narrow downflows (Stein & Nordlund 1989;Cattaneo et al. 1991; see also Fig. 2).This can be described by the filling factor f of downflows as where u z is the mean vertical velocity, whereas u ↑ z and u ↓ z are the corresponding mean upflow and downflow velocities.It was shown in Käpylä (2021) that f is sensitive to the effective Prandtl number of the fluid such that a lower Pr leads to a lower filling factor.In the current paper we performed a similar study focusing on the function of rotation (see Fig. 8).The main result is that f approaches 1/2 in the rapid rotation regime.This is because in rapidly rotating convection, the broad upwellings of non-rotating convection are broken up and the flow consists mostly of smaller scale helical columns where the upflows and downflows are almost invariant.This is due to the Taylor-Proudmann constraint which causes derivatives along the rotation axis to vanish.Hence the tendency for larger structures to appear at greater depths is inhibited, and the average size of convection cells as a function of depth is almost constant (see rightmost panel of Fig. 2).
A221, page 9 of 16 This is also apparent from the probability density functions (PDFs) of the velocity components u i , defined via Figure 9 shows representative examples of PDFs for the extreme cases (Run A0 with Co = 0 and Run A9 with Co ≈ 16.5) and at an intermediate rotation rate (Run A6, Co = 1.3).In nonrotating convection, the PDFs of the horizontal components of the velocity are nearly Gaussian near the surface, whereas for u z the distributions are highly skewed due to the upflow and downflow asymmetry.In deeper parts, the horizontal velocities also deviate from a Gaussian distribution, in agreement with earlier works (e.g., Brandenburg et al. 1996;Hotta et al. 2015;Käpylä 2021).
As the rotation increases, the asymmetry of the vertical velocity decreases such that in the most rapidly rotating cases considered in this work with Co ≈ 17, u z also approaches a Gaussian distribution.Only near the surface (z/d = 0.85) does a weak asymmetry remain.The horizontal components of velocity continue to have a Gaussian distribution as rotation increases, although there is not enough data to say anything concrete concerning the tails of the distributions at high velocity amplitudes.To further quantify the statistics of the flow, skewness S and kurtosis K are computed from: where σ u = (M 2 ) 1/2 , with A221, page 10 of 16 Figure 10 shows S and K for all u i for the same runs as in Fig. 9.The skewness is consistent with zero for the horizontal velocities, which is expected, as there is no anisotropy in the horizontal plane.The negative values of S for u z are a signature of the asymmetry between upflows and downflows.As rotation increases, S approaches zero also for u z .The kurtosis K is a measure of non-Gaussianity, or intermittency.In the non-rotating case K increases from roughly three, indicating Gaussian statistics, to roughly five for horizontal flows as a function of depth within the CZ.For u z , the increase of K is much more dramatic below z/d 0.3.This is because downflows merge at deeper depths such that only a few of them survive deep in the CZ and especially in the overshoot region below roughly z = 0, where K reaches a peak value of roughly 65 for Run A0.A similar, albeit lower, maximum also appears for the horizontal flows.At intermediate rotation (Run A6; Co = 1.3), u z still exhibits strong intermittency below z ≈ 0.1 with max(K) ≈ 54, whereas K for the horizontal flows is significantly reduced in comparison to the non-rotating case.This indicates that the vertical flows in this regime especially are not qualitatively different from those in the non-rotating regime such that the downflows in the overshoot region are rather abruptly decelerated and diverted horizontally.For the most rapidly rotating case (Run A9; Co = 16.5),K ≈ 3 . . . 4 throughout the simulation domain for both vertical and horizontal flows.This is explained by the almost complete wiping out of the upflow-downflow asymmetry also in the deep parts of the CZ and in the overshoot region.The absence of a peak in the kurtosis in the overshoot region in the most rapidly rotating cases is likely due to the deeply penetrating vertical flows in those cases due to the unrealistically small Richardson number.This is discussed in more detail in Sect.3.7.
The average vertical rms velocities from the same runs as in Fig. 9 are shown in Fig. 11.The average rms velocity of the downflows (upflows) is always larger (smaller) than the average total vertical rms velocity.However, the difference between the upflows and downflows and the total rms velocity diminishes monotonically as a function of rotation such that for the most rapidly rotating case the three are almost the same.This is another manifestation of the symmetrisation of upflows and downflows.
Another consequence of the symmetrisation of the vertical flows is that the forces on the upflows and downflows also approach each other (see Fig. 12), where f z = ρDu z /Dt.In accordance with earlier studies (Käpylä et al. 2017; Käpylä 2019), in non-rotating convection, the downflows are accelerated near the surface and decelerated roughly when the stratification turns Schwarzschild stable, whereas the upflows are accelerated everywhere except near the surface.This is interpreted such that the upflows are not driven by buoyancy but by pressure forces due to the deeply penetrating downflow plumes.This qualitative picture remains unchanged for slow and intermediate rotation (see Fig. 12b).For rapid rotation, the forces on the upflows and downflows are nearly identical.However, the situation continues to qualitatively deviate from the mixing length picture also in the rapidly rotating cases in that the downflows are accelerated only near the surface and braked throughout their descent through the superadiabatic CZ (see Fig. 12c).

Overshooting and Deardorff zones
We studied the depths of the overshooting and Deardorff zones as functions of rotation using the same definitions as in previous studies (Käpylä 2019(Käpylä , 2021)).We defined the bottom of the CZ to be the depth z CZ where F conv changes from negative to positive with increasing z.The top of the Deardorff zone (DZ) -or the bottom of the buoyancy zone (BZ)z BZ was taken to be where A221, page 11 of 16 the superadiabatic temperature gradient changes from negative to positive with increasing z.The depth of the DZ is then where ∆t = t 1 − t 0 is the length of the statistically steady part of the time series.We measured a reference value of the kinetic energy flux (F This criterion breaks down in the current models when rotation begins to dominate the dynamics and where F kin → 0. Therefore, the convected flux F conv was also used to estimate the depth of overshooting.The criterion involving F conv assumes the overshoot zone ends at the location (z conv OZ ) where |F conv | falls below 0.02F tot .We computed the corresponding overshooting depth (d conv OZ ) analogously to Eq. ( 57).The layer below the OZ is the radiative zone (RZ).Figure 13 shows the energy fluxes from representative runs at different Coriolis numbers from Set A. For slow and moderate rotation up to Co ≈ 1, the situation is qualitatively similar: the positive (upward) enthalpy flux exceeds F bot in the bulk of the CZ, and it is compensated by a negative (downward) kinetic energy flux F kin .As rotation increases, the maxima of F enth and |F kin | decrease monotonically.Similarly, the extents of the overshoot and Deardorff zones diminish with rotation.For the case with the most rapid rotation, Run A9 with Co = 16.5, the kinetic energy flux is almost zero, and F conv ≈ F enth .This is yet another manifestation of the decreasing asymmetry between the upflows A221, page 12 of 16 The positions of the boundaries of the different layers and their depths are summarised for all runs in Table 2, and Fig. 14 shows a summary of the overshooting and Deardorff zone depths as a function of rotation from Sets A, B, and C. The main difference between the sets of simulations is the applied flux F n .The overshooting depth measured from the kinetic energy flux decreases with increasing rotation as in earlier studies (e.g., Ziegler & Rüdiger 2003;Käpylä et al. 2004).However, the lowermost panel of Fig. 13 shows that the upper part of the radiative zone is mixed far beyond the regions where F kin is non-negligible in the rapidly rotating cases.This is confirmed when the convected flux is used to estimate the overshooting depth.Furthermore, d conv OZ increases with rotation for Co 1.This is explained by the fact that the Mach number, and therefore also the rotation rate Ω 0 , in the current simulations are much larger than in real stars.This means that the convective, rotation, and Brunt-Väisälä frequencies are closer to each other in the simulations in comparison to the overshoot region of the Sun.For example, in the most rapidly rotating runs, the Richardson number Ri Ω based on the rotation rate is smaller than unity (see the 12th panel of Table 1).This, in addition to the smooth transition from the convective to the radiative region, can lead to gravity waves breaking in the radiative zone, thus contributing to the burrowing of the flows into the RZ (e.g., Lecoanet & Quataert 2013).As a comparison, Ri Ω in the upper part of the solar radiative zone is expected to be on the order of 10 4 .Another possibility is that shear caused by the rotationally constrained convective columns lowers the corresponding shear Richardson number close to the limit where turbulence can occur also in stable stratification.
Lowering the luminosity in Sets B and C shows that both measures of d OZ decrease with F n in qualitative accordance with earlier results (e.g., Käpylä 2019).Even though Ri Ω modestly increases in these runs (see the 12th column in Table 1), the most rapidly rotating cases, even in the runs with the lowest luminosities, continue to show deep mixing, which is most likely due to the still unrealistically low Ri Ω .It is numerically very expensive to increase the Richardson number in fully compressible simulations much further, at least without accelerated thermal evolution methods (e.g., Anders et al. 2018Anders et al. , 2020)).Comparing the overshooting depths between Runs [A,B,C]5 with solar Co F and the non-rotating Runs [A,B,C]0 shows a reduction between about a third to a half (see the seventh and eight columns in Table 2).In Käpylä (2019) the overshooting depth extrapolated to the solar value of F n was found to be roughly 0.1H p , and the current results including rotation reduce this to 0.05 . . .0.07H p .However, the dependence of the overshooting depth on F n is steeper here ([ dkin OZ , dconv OZ ] ∝ F 0.15 n ) than in the nonrotating cases where Käpylä (2019) found d OZ ∝ F 0.08 n .On the other hand, the thickness of the Deardorff zone d DZ decreases monotonously as a function of Co.In the most rapidly rotating cases, the Deardorff zone vanishes altogether and even reverses such that at the base of the CZ, the stratification is unstably stratified but the convective flux is inward (see the lowermost panel of Fig. 13).This does not change significantly in more supercritical runs, such as A9m and A9h.In the entropy rain picture (e.g., Brandenburg 2016), cool material from the surface A221, page 13 of 16 is brought down deep into the otherwise stably stratified layers.This is mediated by relatively few fast downflows with a filling factor of f (z) < 1/2 that also produce a strong net downward kinetic energy flux, as seen in the top panel of Fig. 13 (see also Fig. 8 and Table 1 and Sect . 3.3 in Brandenburg 2016).If, on the other hand, the upflows and downflows are symmetrised such that f (z) = 1/2 and their velocities are nearly the same, F kin vanishes, and non-local transport due to downflows is no longer significant.Therefore, the kinetic energy flux is a proxy of the non-local transport due to downflows, and its absence signifies the absence of a Deardorff zone.The depth of the Deardorff zone is independent of the energy flux F n .This further illustrates that the DZ is caused by surface effects that are kept independent of F n in the current simulations.A reduction of d DZ of about a third between the non-rotating runs [A,B,C]0 and the runs with the solar value of Co F (Runs [A,B,C]5) was found (see the sixth column of Table 2).

Conclusions
Simulations of compressible convection were used to study the convective scale and scalings of quantitites such as the Coriolis number and convective velocity as functions of rotation.The results were compared to those expected from scalings obtained for incompressible convection with slow and fast rotation (Aurnou et al. 2020).The actual length scale is almost unaffected by rotation for Co 1 and decreases proportionally to Co 1/2 for rapid rotation.Correspondingly, the dynamical Coriolis number Co is proportional to Co for slow rotation and to ∝ Co 1/2 for rapid rotation.Furthermore, Co is proportional to (Ra F ) −1/3 for slow and ∝ (Ra F ) −1/5 for rapid rotation, where Ra F is the diffusion-free flux-based modified Rayleigh number.Finally, the convective velocity is compatible with proportionality to (F tot /ρ) 1/3 for slow and ∝ (F tot /ρ) 1/3 Co −1/6 for rapid rotation.All of these scalings are consistent with those derived by Aurnou et al. (2020) andVasil et al. (2021).
In an earlier work (Käpylä 2023b), several measures were used to characterise the rotational influence on convection.A commonly used definition where the changing length scale of convection is taken into account is Co ω = 2Ω/ω rms .We have shown that this quantity cannot be used to characterise the effects of rotation on the mean scale because ω rms is expected to increase with the Reynolds number as Re 1/2 .Therefore the only reliable way to account for the changing convective length scale as a function of rotation is to compute the mean wavenumber.This was not correctly identified in Käpylä (2023b), and it is now clear that Co ω will diverge as Re increases.On the other hand, Käpylä (2023b) introduced a stellar Coriolis number Co that depends on luminosity and rotation rate, which are observable, and a reference density, which is available from stellar structure models, but not on any dynamical length or velocity scale.In our work, we renamed this quantity as Co F , and we showed that with a suitable choice of length scale, Co F = (Ra F ) −1/3 .Matching Co F (or equivalently Ra F ) with the target star gives more concrete meaning to the often-expressed idea that it is possible to match the Coriolis number of, for example, the Sun with 3D simulations, while most other dimensionless parameters are out of reach (cf.Käpylä et al. 2023).
The current simulations suggest that convection, even in the deep parts of the CZ in the Sun, is not strongly rotationally constrained and that the CIA balance is therefore inapplicable there.The latter has been assumed to be the case by Featherstone & Hindman (2016) and Vasil et al. (2021) to argue that the largest convectively driven scale in the Sun is the supergranular scale.The current results seem to refute this conjecture and we instead argue that the actual scales may be larger.
Finally, the effects of rotation on convective overshooting and subadiabatic Deardorff zones were studied.The effects of rotation are relatively mild such that for the case with the solar value of Co F , the overshooting depth and the extent of the Deardorff zone are reduced by between 30% and 50% in comparison to the non-rotating case.Therefore the current results suggest an overshooting depth of about 5% of the pressure scale height at the base of the solar CZ.Taking the current results at face value, a similar depth is estimated for the Deardorf zone.However, the latter is still subject to the caveat that the current simulations do not capture the near-surface layer very accurately (e.g., Käpylä 2023a) and that the driving of entropy rain can be significantly stronger in reality.Another aspect that needs to be revisited in the future is the effect of magnetic fields.

Fig. 1 .
Fig. 1.Comparison of hydrostatic and convecting models.Panel a: Temperature as a function of height from a 1D hydrostatic model (black solid line) and convective runs A0 (red dashed line), A6 (blue dashed line), and A9 (orange dashed line).The black (red) dotted line shows a polytropic gradient corresponding to index n rad = 3.25 (n ad = 1.5) for reference.Panel b: Absolute value of the superadiabatic temperature gradient ∆∇ from the same runs as indicated by the legend.Red (blue) indicates convectively unstable (stable) stratification.The dotted vertical lines at z = 0 and z/d = 1 denote the base and top of the initially isentropic layer.

Fig. 3 .
Fig. 3. Normalised velocity power spectra ẼK (k) = E K (k)/ E K (k)dk near the surface (left panel), middle (middle), and base (right) of the CZ from runs in Set A with Co varying between 0 and 16.5.The inset in the left panel shows the mean scale kmean = k mean /k H and wavenumber where E(k) has its maximum ( kmax = k max /k H ) as functions of Co for z/d = 0.85.The error bars indicate the standard deviation.The grey dashed line shows a power law proportional to Co 1/2 .

Fig. 5 .
Fig. 5. Dynamical Coriolis number Co as a function of Co for all of the runs in Table1.Power laws proportional to Co (slow rotation; Co < 1) and Co 1/2 (rapid rotation; Co > 2) are shown for reference.The inset shows Co as a function of Ra F with power laws proportional to (Ra F ) −1/5 (rapid rotation; Ra F < 3 × 10 −3 ), and (Ra F ) −1/3 (slow rotation; Ra F > 0.03).

Fig. 6 .
Fig. 6.Normalised velocity power spectra near the surface of simulations with Co ≈ 1.3 and Re = 40 . . .174 (Runs A6, A6m, and A6h).The dotted line shows Kolmogorov k −5/3 scaling for reference.The inset shows kmean = k mean /k H (black symbols) and kω = k ω /k H (red) as functions of Re.The dotted lines are proportional to powers 0, and 1/2 of Re.

Fig. 7 .
Fig. 7. Volume-averaged rms velocity normalised by u from Sets A, B, C, and Am.The dotted lines for Co < 1.5 denote constants 1.15 and 1.23, respectively, whereas for Co > 6, they are proportional to Co −1/6 as indicated by the theoretical CIA scaling (see Eq. (51)).

Fig. 10 .
Fig.10.Skewness (S, dashed lines) and kurtosis (K, solid) from the same runs as in Fig.9.The black, blue, and red colours indicate data corresponding to u x , u y , and u z , respectively.We note the difference in scale between each of the panels.The insets show a zoom-in of the region z/d ≥ 0.

Fig. 11 .
Fig. 11.Horizontally averaged vertical rms velocity for the same runs as in Fig. 9.The overall vertical velocity (ũ rms z ) is shown in black, and the corresponding quantities for upflows (ũ ↑rms z ) and downflows (ũ ↓rms z ) are shown in red and blue, respectively.The tildes refer to normalisation by gd.

Fig. 12 .
Fig. 12. Horizontally averaged total force (black) and the force acting on upflows (red) and downflows (blue).The dotted red and blue lines show the superadiabatic temperature gradient.Data are shown for (a) a non-rotating Run A0, (b) an intermediate rotation rate (Co = 1.3, Run A6), and (c) for rapid rotation (Co = 16.5, Run A9).
ref kin ) at z CZ .The base of the overshoot zone (OZ) is taken to be the location (z kin OZ ) where |F kin | falls below 0.01F ref kin and

Fig. 13 .
Fig. 13.Time-averaged mean energy fluxes as defined in Eqs.(34) and (38) (apart from the negligibly small viscous flux F visc ) as functions of z from Runs A0, A6, and A9.The red circles indicate (from left to right) the bottoms of OZ, DZ, BZ, and the top of the BZ.The grey (orange) shaded areas indicate mixed (radiative) regions.

Fig. 14 .
Fig. 14.Depth of the overshoot zone from kinetic energy (d kin OZ ; black lines) and convective fluxes (d conv OZ ; red), and depth of the Deardorff zone (d DZ ; blue) as functions of rotation measured by Co.All quantities have been normalised by the pressure scale height at the base of the CZ.The different lines correspond to the three different values of F n or to Sets A (solid lines), B (dashed), and C (dotted).

Table 1 .
Summary of the runs.

Table 2 .
Summary of the buoyancy, Deardorff, and overshoot zones.Run z BZ /d z DZ /d z kin OZ /d z conv OZ /d d DZ /d d kin OZ /d d conv