On the Penetration of Large-scale Flows into Stellar Radiative Zones

The propagation of meridional circulation below the base of the convection zone (CZ) of low-mass stars may play a crucial role in the transport of angular momentum and also significantly contribute to the transport of chemical species and magnetic fields within their stable radiative zone (RZ). We systematically study these large-scale mean flows by performing three-dimensional global numerical simulations in a spherical shell that consists of a convective region overlying a stably stratified region. We find that the meridional flows can penetrate distances as large as ∼0.21r o (where r o is the outer radius) below the base of the CZ, provided that the Eddington–Sweet timescale t ES is much shorter than the viscous timescale t ν , as measured by the parameter σ=(tES/tν)1/2 . In the solar-like regime, where σ ≲ 1 in the upper RZ, we find that the angular momentum transport in the deep RZ is determined primarily by the action of the Coriolis force on meridional flows. In contrast, in models run in the σ > 1 regime, the meridional flows become weaker and the viscous effects dominate. We find that the penetration lengthscale δ MC of these mean flows when σ ≲ 1 is proportional to σ −0.22. Our findings may provide a better understanding of the role of the meridional flows in the dynamics of the solar interior and inform future numerical studies that are focused on capturing solar-like dynamics self-consistently.


INTRODUCTION
Helioseismic observations have provided us with the Sun's rotational profile and have revealed that its radiative interior rotates almost uniformly while the outer convective region rotates differentially, with the equator rotating faster than the poles (e.g., Brown et al. 1989;Thompson et al. 2003;Howe 2009).Differential rotation such as that observed in the Sun (i.e., with a rapidly-rotating equator) is thought to be maintained through anisotropic angular momentum transport arising in response to the tendency of convective cells to organize into columnar-like structures (e.g., Zhang 1992;Busse 2002;Aurnou et al. 2007;Camisassa & Featherstone 2022).
The action of the Coriolis force on the differential rotation is such that near the rapidly-rotating equator, fluid motions are driven away from the rotation axis whereas near the slowly rotating poles, the Coriolis force acts instead to push the flow toward the rotation axis.The result of this process, known as "gyroscopic pumping", is to drive axisymmetric motions (the nominal meridional circulation) in the north/south and radially inward/outward directions.Gyroscopic pumping had been initially studied in the context of the Earth's circulation and later discussed in the astrophysical context for its role in the solar interior dynamics (McIntyre 2007;Garaud & Acevedo Arreguin 2009;Garaud & Bodenheimer 2010;Wood et al. 2011;Miesch & Hindman 2011).
In the near-surface layers, a variety of helioseismic and surface-feature-tracking techniques have revealed the Sun's meridional flow to be predominantly poleward within each hemisphere, with a characteristic speed of 10-20 m/s (e.g., Komm et al. 1993;Švanda et al. 2007;Ulrich 2010;Hathaway et al. 2022).Unlike the well-known solar rotational profile, however, there is no clear consensus on the structure of meridional circulation in the deep convection zone.In particular, substantial disagreement regarding the depth of the return flow, as well as the multi-or mono-cellular nature of the meridional flow, persists between measurements made using different instruments and helioseismic techniques (e.g., Zhao et al. 2013;Jackiewicz et al. 2015;Gizon et al. 2020).This disagreement increases with depth and, as a result, the question of the penetration of these mean flows below the base of the convection zone into the underlying radiative zone remains unanswered.Meridional flows are nevertheless thought to play an important role in the dynamics that arise near the interface between the convection and radiative zones, mediating the mixing of chemical species (e.g., Pinsonneault 1997), the transport of angular momentum and the transport of magnetic fields (Gough & McIntyre 1998;McIntyre 2007;Wood et al. 2011;Acevedo-Arreguin et al. 2013;Wood & Brummell 2018).
Direct numerical simulations offer an alternative means to explore meridional flows below the base of the convection zone.As with all such studies, however, care must be taken to ensure that the simulated system is operating in a parameter regime relevant to the physical system under consideration.The sense of differential rotation realized in a spherical convection zone, for example, is known to depend on the ratio of the buoyancy and Coriolis forces.When buoyancy is dominant, a so-called antisolar differential rotation with rapidly-rotating poles and a slowly-rotating equator develops (e.g., Gastine et al. 2014).For this reason, models designed to study the Sun are constructed such that the Coriolis force in dominant over buoyancy.
In a similar vein, the structure and penetration depth of meridional flows are thought to be determined by the competition of two timescales, namely the timescale associated with viscous diffusion and the timescale related to the advection of the gyroscopically-pumped meridional flows.When viscous effects dominate, the meridional flows are exponentially damped in depth and confined to a small region below the inner convective boundary.When viscous effects are negligible, gyroscopic pumping causes meridional flows to penetrate long distances below the base of the convection zone (Garaud & Brummell 2008;Garaud & Acevedo Arreguin 2009;Garaud & Bodenheimer 2010;Wood & Brummell 2012;Acevedo-Arreguin et al. 2013).
In solar-type stars, viscous effects are negligible near the convection zone-radiative zone interface, and in the case of the Sun, the timescale associated with the gyroscopically-pumped meridional flows is much faster than that related to viscous diffusion in the upper part of the radiative region (see Section 2, and e.g.Wood & Brummell 2012).A number of studies have examined angular momentum transport and the driving of meridional flow in this low-diffusivity regime.Notably, the work of Garaud & Acevedo Arreguin (2009) and Wood & Brummell (2012) (see, also Garaud & Brummell 2008;Garaud & Bodenheimer 2010) demonstrated that when the proper ordering of timescales was respected, meridional flows penetrated large distances below the base of the convection zone.These calculations were carried out, however, either in axisymmetric spherical shells or in 3D Cartesian boxes.That restricted geometry prevented the self-consistent development of differential rotation and meridional circulation, which require a fully 3D, spherical domain to self-consistently capture the mean flows and the convective motions that drive them.On the other hand, those models that do account for these nonlinear and geometric effects were run in a regime where viscous effects were dominant.For instance, Brun & Zahn (2006) and Strugarek et al. (2011) found in fully nonlinear simulations employing 3D, spherical geometry that the angular momentum transport was dominated by viscosity within the bulk of the RZ and the meridional circulation did not propagate substantially beyond the base of the convection zone.Recently, Matilsky et al. (2022) also ran 3D magnetohydrodynamic (MHD) global simulations to study the self-consistent formation of the tachocline due to dynamo action within the stable radiative zone.Their simulations were run in the non-solar regime and consequently they also found that viscous diffusion played a dominant role in their dynamics.
In this paper, we seek to bridge the gap between these two approaches.Our goal is to understand the dynamics occurring within the convection zone-radiative zone (CZ-RZ) interface, in a fully 3D, spherical geometry, by exploring a range of parameter regimes in which either viscous or gyroscopic pumping effects dominate.In Section 2, we describe our two-zone spherical shell formulation and provide the set of anelastic Navier-Stokes equations along with the initial conditions, input parameters, and the boundary conditions used in our simulations.In Section 3, we present our results related to the differential rotation and meridional circulation profiles within and below the CZ and their dependence on the input parameters.We also discuss angular momentum transport within the different parameter regimes and compare the time evolution of the corresponding fluxes among the different cases.Finally in Section 4, we provide a summary of our findings along with a discussion of their relation to other studies and the implications of these results for the solar interior dynamics.

Dimensional Equations
We are interested in exploring the dynamics related to the global, large-scale flows arising in rotating overshooting convection by accounting for a two-layered system with a stably stratified region underlying a convective region.We assume a fixed aspect ratio r i /r o = 0.45 where r i is the inner radius and r o is the outer radius of the spherical shell while the depth of the convective region is given by L = r o − r c = 0.2408r o , where r c = 0.7592r o is the radius at the inner boundary of the CZ.The depth L is chosen such that it corresponds to the solar convection zone from ∼ 0.7187R ⊙ to ∼ 0.9041R ⊙ − 0.9695R ⊙ , since we aim to focus on the deep interior dynamics and not the radiative transfer processes occurring near and at the outer solar convective boundary.Consequently, we employ the anelastic approximation which assumes small perturbations of the thermodynamic variables compared with their mean while it also filters out the sound waves (Gough 1969;Gilman & Glatzmaier 1981).We solve the 3D Navier-Stokes equations under the anelastic and Lantz-Braginsky-Roberts approximations (Lantz 1992;Braginsky & Roberts 1995) with the latter being exact where the reference state is adiabatic.
Then, the dimensional Navier-Stokes equations become where u = (u r , u θ , u φ ) is the velocity field, P is the pressure, ρ is the reference density, T is the reference temperature, d S/dr is the background entropy gradient, S corresponds to the entropy perturbations about the reference state, Ω o is the stellar mean frame rotation rate (for the Sun, this corresponds to Ω 0 = 2.87 • 10 −6 rad/s), g(r) is the gravity (where g ∝ 1/r 2 ), and c p is the specific heat at constant pressure.We assume that the viscosity ν and the thermal diffusivity κ are constant for simplicity.The viscous stress tensor D is given by and the viscous heating is denoted by where e ij is the strain rate tensor.We adopt an internal heating term Q that is similar to the one formed in Model S (see e.g., Christensen-Dalsgaard et al. 1996) by making it proportional to the reference pressure P (for more details see Korre & Featherstone 2021).The internal heating satisfies Q(r) = −∇ • F rad , where F rad is the radiative flux in the system given by The equation of state is ρ ρ = assuming the ideal gas law P = ℜρ T , where ρ and T characterize respectively the density and temperature perturbations about the reference state, ℜ is the gas constant, and γ = c p /c v , where c v is the specific heat at constant volume.
To create our two-zone spherical shell, we choose a profile of d S/dr that satisfies an adiabatic polytropic solution with γ = 5/3 in the convective region r c ≤ r ≤ r o and that forms a stably stratified region for r i ≤ r < r c .
We also designate the number of density scale-heights in the convection zone through the parameter N ρ = ln(ρ(r c ))/ρ(r o )) (see Appendix of Korre & Featherstone (2021) for details).Then, we integrate d S/dr from the bottom of the convection zone r c inward which results in a solution for S with an integration constant that matches the solution to the value of S at r c .
Given that the entropy for a monatomic ideal gas is S = ln( P 1/γ /ρ), we can differentiate S with respect to r, and apply the hydrostatic balance equation ∂ P /∂r = −ρg.Then, we obtain This allows us to numerically solve for the reference density profile ρ(r) of any chosen background entropy gradient profile.

Nondimensional Equations
We nondimensionalize equations ( 1)-(3) using the depth of the convection zone [l] = L as the lengthscale, the velocity scale [u] = ν/L and the viscous timescale [t] = L 2 /ν.Throughout the paper, we define the volume average of a quantity q(r, θ, φ) in the CZ as the time and spherical average of a quantity q as q(r) = 1 4π(t 2 − t 1 ) and the time and azimuthal average of a quantity q as q(r, θ) = 1 2π(t 2 − t 1 ) For the density and temperature scales, we use [ρ] = ρcz , and [T ] = Tcz , respectively.For the entropy perturbations S, we choose a scaling associated with the thermal energy flux F such that [S] = LF cz /(ρ cz Tcz κ), where F (r) = r r i Q(r ′ )r ′2 dr ′ /r 2 .We determine our selected nondimensional d S/dr profile in the RZ using the function where we vary A S , r b and d.In Figure 1, we show the profiles of d S(r)/dr versus r/r o for five typical runs.We note that for larger values of A S , we obtain a more stably stratified RZ while the different values of r b and d control the smoothness of the transition width between the two layers.
The nondimensional reference density ρ, reference temperature T and gravity are equivalent to the dimensional reference density, reference temperature and gravity divided by ρcz , Tcz , and g cz , respectively.The nondimensional anelastic Navier-Stokes equations become  13)).The black dashed vertical like marks the transition from the adiabatic convection zone (where d S(r)/dr = 0) to the subadiabatic layer (where d S(r)/dr > 0).
where Q nd = LQ/F cz is the nondimensional internal heating function (for details, see Korre & Featherstone 2021).We note that u = u m + u f , where u m is the mean, axisymmetric part of the velocity corresponding to the large-scale motions related to rotation and u f is the fluctuating component of the velocity field associated with the convective motions.In all that follows, all the variables and parameters are now nondimensional and this nondimensionalization introduces the flux Rayleigh number Ra, the Prandtl number Pr, the Ekman number Ek and the dissipation number Di, which are defined respectively as We note that although Ra, Pr and Ek are free input parameters, Di is not as it depends on the chosen N ρ and accounts for the fact that we have both a thermal scale and a kinetic scale.Thus, Di appears in the thermal equation to account for the viscous heating term, which is intrinsically kinetic.
We have run a series of 3D numerical simulations solving equations ( 14)-( 16) using the open-source convection code Rayleigh (Featherstone & Hindman 2016;Matsui et al. 2016;Featherstone 2018).We vary the background entropy gradient d S/dr as described above, the Rayleigh number, the Ekman number, the number of density scale-heights in the convection zone N ρ while we keep the Prandtl number fixed at Pr = 1 in order to avoid the excitation of unphysical modes that are known to appear in rapidly rotating anelastic convection systems for values of Pr < 1 (e.g., Calkins et al. 2015).
We provide the parameters used in our simulations in Table 1.We also report on the convective Rossby number Ro c = RaEk 2 /(4Pr) and the output Rossby number defined as Ro = U/(2LΩ 0 ) = ReEk/2, which is the ratio of inertial forces to the Coriolis force.Here, the first expression involves dimensional quantities where U = u rms ν/L is a typical dimensional velocity in the CZ and where for the value of u rms , we use the total (fluctuating part and mean part) velocity averaged over the CZ and weighted by density.The second expression includes nondimensional quantities where Re= UL/ν is the Reynolds number extracted from the simulations while the Péclet number Pe is Pe = PrRe = Re in our Pr = 1 runs (see Table 1).
We use impenetrable and stress-free boundary conditions for the velocity while for the entropy perturbations, we assume ∂S/∂r| r i = 0 at the inner boundary and S| ro = 0 at the outer boundary.Each simulation is evolved from a zero initial velocity and small-amplitude perturbations in the entropy field until a statistically stationary and thermally equilibrated state is achieved.

Choice of Input Parameters
As discussed in Section 1, in order for the meridional circulation to penetrate below the base of the CZ more deeply into the stable region, the correct ordering of timescales needs to be achieved in numerical simulations even though the latter cannot account for real stellar values.The gyroscopicallypumped meridional circulation operates on a timescale of the order of the Eddington-Sweet timescale Spiegel & Zahn 1992;Wood & Brummell 2012) which is shorter than the viscous timescale t ν = L 2 /ν in the upper solar radiative zone (Fig. 2).The ordering of timescales, namely t ES < t ν , can be expressed via the parameter σ = (t ES /t ν ) 1/2 = ν/κN/Ω 0 (Garaud & Acevedo Arreguin 2009;Garaud & Bodenheimer 2010;Wood & Brummell 2012).In our nondimensional formulation, σ is written as where N(r) = (Ra/Pr)g(d S(r)/dr) is the Brunt-Väisälä frequency.We are interested in systematically investigating the dependence of the meridional flows on σ within a 3D spherical shell that self-consistently accounts for large-scale mean flows.To achieve this, we vary σ(r) below the base of the convection zone.In our Pr = 1 simulations, we can do so by varying Ek and/or N(r) through Ra and/or d S(r)/dr.We note, however, that in our runs the variation in σ(r) is mostly a result of the different d S/dr profiles chosen.
In Table 1, we show all the cases considered in this study by reporting σ ov which is σ evaluated at the radius down to which the convective motions overshoot in each case.That is at r c − δ G r o , where δ G = δ Gf /r o is the overshoot lengthscale (for more details on how we calculate this, see Section 3 from Korre & Featherstone (2021)).In Korre & Featherstone (2021), we found that a convective Rossby number Ro c ≈ 0.5 leads to solar-like rotation with the equator rotating faster than the poles in the convective region.Since we are interested in low-mass stars, we vary our input parameters such that we are in the Ro c ≈ 0.5 regime and study the influence of σ on the mean flow dynamics.
In Figure 2, we show the function σ(r) versus the radius r/r o below the base of the convection zone for the same five representative runs illustrated in Figure 1.We also show the solar σ(r) profile calculated from Model S (Christensen-Dalsgaard et al. 1996), where we have moved the solar convection zone by 0.046r o to match the base of the CZ of our spherical shell so that all the σ(r) profiles can be directly comparable in the same plot.For the calculation of the solar σ(r), we have used the density, temperature, pressure, Brunt-Väisälä frequency, specific heat at constant pressure and gravity of Model S while we have calculated the thermal diffusivity and viscosity (and as a result the Prandtl number) using the formulae given in Garaud & Garaud (2008) (see also Gough 2007).We see that case 9 has a smaller than one σ(r) profile within the whole RZ, while in case 10, σ(r) is overall small but smaller than unity only near the upper part of the radiative zone, similarly to the actual solar case.Case 12 possesses a σ < 1 profile near the upper stable layer unlike cases 3 and 5 which both have σ >> 1 across the whole RZ.By illustrating different diagnostics of these five typical cases that span a wide range of σ ov values, in the following sections, we will explore the dynamics associated with the large-scale mean flows, and their dependence on σ.

Differential Rotation Profiles
We begin our investigation by focusing on the rotational profiles of our simulations and the mean flows associated with the differential rotation within and below the CZ.In Figure 3, we illustrate profiles of the differential rotation u φ /(r sin θ), for the five representative cases shown in Figure 2. We verify that for Ro c ≈ 0.5, the simulations have a solar-like rotation whereby the equatorial region rotates more rapidly than the polar regions throughout the convection zone.
We do note that, while these models incorporate a region of overshoot, contours of isorotation are largely parallel to the rotation axis, as opposed to tilted in the radial direction as observed in the Sun (e.g., Howe et al. 2005).This "tilting" of the isorotation contours has been posited to arise in response to strong latitudinal temperature variations established in the tachocline region (Rempel 2005;Miesch et al. 2006).As in Korre & Featherstone (2021), we observe only very weak deviations from cylindrical isorotation contours, in spite of the fact that a region of overshoot has been included in these simulations.This is possibly because the thermal-wind balance required to induce deviations from cylindrical isorotation is instead produced in response to uniform thermal flux at the upper boundary, which is not imposed for the models considered here (e.g., Matilsky et al. 2020).None of the models illustrated in Figure 3 possess a true tachocline of shear at the base of the convection zone.For all values of σ ov , differential rotation profiles established in the convection zone are reflected to some extent in the radiative zone.This behavior is largely in accord with the widely-accepted consensus that additional physics, such as magnetic fields (either primordial, or dynamo, or both) are needed to reproduce the Sun's rotation and the self-consistent emergence of a tachocline (Gough & McIntyre 1998;Acevedo-Arreguin et al. 2013;Wood & Brummell 2018;Matilsky et al. 2022).
We note, however, that for the low-σ ov cases, differential rotation in the equatorial region is considerably more confined to the upper radiative zone than in the high-σ ov regime runs.This can be seen more clearly in Figure 4 where we illustrate the radial dependence of the differential rotation at different colatitudes for the case with σ ov = 0.1 and the case with σ ov = 9.8.We verify that the profiles are quite different between the two cases, with the low-σ ov run presenting a smaller variation of differential rotation below the base of the CZ and a profile at θ ≈ 90 o that drops more rapidly beyond r c , unlike the high-σ ov case.As we discuss in Section 3.3, in the systems with σ ov 1.1, it is advection of angular momentum by meridional circulation that halts the spread of differential rotation that is otherwise observed in the equatorial regions of models with σ ov >> 1, where viscous effects dominate.To see why this might be the case, we now examine the accompanying meridional circulations established in these different models.

Meridional Flows Below the Base of the Convection Zone
We present profiles of the meridional circulation in Figure 5. There, we plot the timeand azimuthally-averaged meridional circulation streamlines overlying the mass flux function  It is below the base of the convection zone where significant differences in the meridional flow morphology arise.For the low-σ ov models 9 and 10, the mean flows penetrate deeply into the radiative zone, even reaching the lower boundary of the simulation domain, although their amplitude is much diminished there.Models 3 and 5 possess values of σ that are larger than unity everywhere in the radiative zone, and their resulting meridional flows are confined to a small depth just below the base of the CZ.The intermediate case 12, which possesses σ(r) < 1 only in the upper radiative zone (and a σ ov of 1.1) has meridional flows that travel well below the base of the CZ, beyond the overshoot region but their amplitude rapidly decays in the deeper stable zone.
These results qualitatively support the earlier theoretical arguments and restricted-geometry models discussed in e.g.Garaud & Brummell (2008), Garaud & Acevedo Arreguin (2009), Wood & Brummell (2012), and Acevedo-Arreguin et al. (2013).When σ >> 1, in our runs N(r) is large, namely the thermal stratification in the RZ is stronger and the meridional flows are viscouslydamped before they can extend to large depths below the base of the CZ unlike in the cases with σ < 1 where the meridional flows are stronger and are able to propagate longer distances into the RZ and even span the whole stable layer (Garaud & Brummell 2008;Garaud & Bodenheimer 2010).
A more quantitative picture of the meridional circulation can be seen in Figure 6.There, the timeand azimuthally-averaged velocity u θ is shown as a function of radius at different colatitudes for the solar-like case 10 with σ ov = 0.1 and the non-solar case 5 with σ ov = 9.8.We observe that the meridional circulation is roughly similar in amplitude and structure at all θ in the CZ for both cases, with the multicellular profiles seen by the change of sign of u θ in the convective region.We note that the meridional circulation profiles are unlike those obtained from observations near the solar surface due to diffusive effects in the boundary layer associated with the outer boundary condition employed in our numerical simulations (for a more detailed discussion on this topic, see, Fuentes et al. 2024).
These trends in the meridional flow behavior are also evident in the radial profiles of the time-and spherically-averaged kinetic energy related to the meridional circulation, ẼMC , which we define as We present profiles of ẼMC for each one of our representative cases in Figure 7a.We observe that ẼMC (r) is approximately the same in the convection zone for all five cases but quite different within the RZ.Indeed, ẼMC (r) decays faster beyond the inner convective boundary as σ becomes larger there for the cases with σ ov > 1.More specifically, in the intermediate case 12 where σ ov = 1.1, ẼMC (r) decreases in the RZ but the damping is slower compared with cases 3 and 5 where the amplitude of ẼMC (r) becomes weaker more rapidly below the base of the CZ.In contrast, both cases 9 and 10, which have σ ov < 1, possess a larger kinetic energy in the meridional flows in the RZ as well as a profile that gradually decays below r c but is not changing by much overall.This is similar to what we qualitatively observed in Figure 5.We can use these profiles to develop a more quantitative description of the penetration lengthscale of the meridional flows in the RZ.For each simulation, we fit an exponential function of the form to the ẼMC (r) data within the region spanning from the base of the CZ (r = r c ) down to the point where the profile of ẼMC is well fit by Equation ( 20) as depicted in Figure 7b for case 10.
Once we have measured the decay rate λ for each case, we define the characteristic penetration lengthscale of the meridional flow as δ M C = δ M Cf /r o = 2/λ.We note that the prefactor 2 comes from the fact that we fit the exponential function to the kinetic energy of the meridional flows, but we are interested in the lengthscale associated with the meridional flows themselves, namely ẼMC ∝ e (−λ(r−rc)) 1/2 = e (−λ(r−rc)/2) .The relationship between the meridional penetration lengthscale δ M C and the parameter σ ov is illustrated in Figure 8.We find that this relationship is well described by a broken power law.For those models with σ ov 1.1, we find that δ M C ∝ σ −0.22 ov .This scaling law suggests a weak scaling of the penetration lengthscale of the meridional circulation with respect to σ ov for values of σ ov 1.1.It also provides a predictive tool that could be used to estimate the penetration depth of the meridional circulation in stellar radiative zones that have σ < 1 (see Section 4).For values of σ ov >> 1.1, we find that δ M C ∝ σ −1 ov .This latter result is similar to what has been shown in e.g.Garaud & Acevedo Arreguin (2009) and Wood & Brummell (2012), associated with the viscous damping of the meridional circulation below the base of the CZ when σ >> 1.

Angular Momentum Transport
As discussed in Section 1, meridional flows can be driven in response to the redistribution of angular momentum, an effect known as gyroscopic pumping.Here we explore this effect by examining the   20)) is fitted to the kinetic energy data below the base of the CZ down to the point where ẼMC is well fit.The computed decay rate λ for each simulation provides a penetration lengthscale δ M C = 2/λ for all runs shown in Table 1.The inner dashed black line corresponds to the base of the CZ in both panels.balances struck in our models between angular momentum transport due to meridional advection, convective transport, and viscous stresses.We can obtain an expression for the conservation of angular momentum by taking the zonal component φ of the momentum equation, multiplying it by the momentum arm µ = r sin θ, averaging over both time and longitude and finally assuming a steady state.When contributions due to the Lorentz force are neglected, we arrive at where u mc = (u r , u θ ) and where L = µ 2 Ω, which can be nondimensionally expressed as The nondimensional transport of angular momentum due to Reynolds stresses and due to viscous stresses is given respectively by We expect that in the thermally-relaxed, statistically-steady state, the time-averaged term associated with the Coriolis force ρ u mc • ∇L will be balanced primarily by the Reynolds stresses −∇ • F RS in the convection zone when viscous stresses can largely be neglected, as is typically true for high values of Ra (see e.g., Featherstone & Hindman 2016).Deeper in the stable region, where the Reynolds stresses are weaker, we expect the viscous stresses −∇ • F V S to contribute significantly to the angular momentum balance.For all models presented in Table 1, we have verified that in the thermally-relaxed and statistically stationary state, ρ u mc •∇L ≈ −(∇•F RS +∇•F V S ) within the CZ and the RZ.In all cases, we notice that −∇ • F RS is more dominant in the bulk of the CZ while −∇ • F V S is larger within the equatorial region, especially near and somewhat below the bottom of the CZ.These results are illustrated in Figure 9, where we show the balances achieved in two extreme cases, with σ ov = 0.1 and σ ov = 9.8.We observe that in the CZ, the balance is achieved between the fluxes associated with both the Reynolds and viscous stresses and the Coriolis force.The Reynolds stresses are somewhat stronger than the viscous stresses as they seem to span the whole CZ unlike the viscous stresses which are mostly present near the equatorial region.This is a result of the level of turbulence dictated by the chosen Ra of our runs shown here, namely, even higher Ra simulations (which are harder to achieve computationally in our two-zone spherical shell) need to be considered for the resulting Reynolds stresses to completely balance the Coriolis term and the viscous stresses to be negligible across the whole CZ.Nonetheless, deeper in the RZ (below the overshoot region), where convective motions are much weaker or non-existent, viscosity becomes more important and counteracts ρ u mc • ∇L.That is in fact expected in these simulations where there is no magnetism and reinforces what has already been suggested in the literature, namely that a magnetic field, either primordial or dynamo in nature, needs to be accounted for in solar-like simulations for the dynamics to be more consistent with the Sun (see, e.g.Gough & McIntyre 1998;Acevedo-Arreguin et al. 2013;Wood & Brummell 2018;Matilsky et al. 2022).
We are interested in understanding the role that each one of the terms in Equation ( 21) plays in the angular momentum transport over time, hence we now investigate the time evolution of their absolute value.We expect to notice differences in the amplitudes of the fluxes that will depend on σ ov for the cases demonstrated in Figure 10 which are all run at the same Ra.That is due to the fact that in the solar-like regime where σ ov 1.1, the meridional flows have a strong presence deeper in the RZ below the overshoot region while for cases with σ ov >> 1, the mean flows are mostly confined to a small region below the base of the CZ and their amplitude decays much faster beyond the convective boundary.
To examine the angular momentum transport more quantitatively and understand the dependence of the fluxes on σ ov , in Figure 10, we present the profiles of the absolute values of the fluxes, namely , volume-averaged from their corresponding overshoot region down to the bottom of the RZ (same calculation as in Equation ( 10) but radially integrated from r i to r c − δ G r o ).We find that | − ρ u mc • ∇L| is much larger than | − ∇ • F V S | for cases 9 and 10 with σ ov << 1 while its amplitude gradually decreases as σ ov becomes larger from panel a) to panel c).In Figure 10d, for case 3, | − ρ u mc • ∇L| ≈ | − ∇ • F V S |, while for case 5 where σ ov >> 1 (Figure 10e), the flux associated with the Coriolis term is now smaller than the flux associated with viscosity, namely | − ρ u mc • ∇L| < | − ∇ • F V S | which indicates that the angular momentum transport is ultimately dictated by viscosity unlike what is expected to occur in the solar interior.Finally, we notice that as expected (and discussed above), the Reynolds stresses are very weak deeper in the RZ below the overshoot region.
These findings are consistent with the differential rotation profiles observed in Figure 3.For the cases with lower σ ov values, we noticed that the differential rotation is more confined within the equatorial region compared with the high-σ ov cases.In Figure 10, we found that in the low-σ ov cases, the Coriolis term is dominant in the angular momentum transport in the RZ while as σ ov becomes larger, the flux related to the mean flows decreases and as a result the viscous stresses eventually control angular momentum transport.Consequently, in the cases with low σ ov , the propagation of  21) for two σ ov cases at the same Ra.The gyroscopic pumping balance is achieved in both cases.a) In case 10 with σ ov = 0.1, the meridional flows penetrate deeply into the RZ.The gyroscopic pumping occurs across the whole shell but is more prominent within the CZ.b) In case 5 with σ ov = 9.8, the mean flows are confined to a small region below the base of the CZ.The gyroscopic pumping balance holds everywhere across the shell for this run as well, however it is negligible below the base of the CZ where the meridional flows rapidly decay.
the differential rotation in the RZ is due to the advection of the mean flows while in the cases with high σ ov , it is due to viscous diffusion.
To better understand this behavior, in Figure 11, we show the the time-and azimuthally-averaged radial fluxes related to the Reynolds stresses, the mean flows, the Coriolis force, and the viscous stresses as well as the total radial flux, all plotted at the base of the CZ (at r = r c ) and given respectively by We notice that in the lower σ ov cases, the total flux F r,T at the base of the CZ is positive near the equatorial region, similarly to the flux related to the Coriolis force F r,C , while the viscous flux F r,V < 0 there.On the other hand, in the cases with σ ov >> 1, F r,T is negative around the equator following the behavior of F r,C and F r,V .This finding suggests that the suppression of the propagation of the differential rotation within the equator for the low-σ ov runs is a result of the radially outward flux of angular momentum, which derives primarily from the interaction of meridional flow with the Coriolis force.This lies in contrast with the high-σ ov cases where F r,T << 0 leads to the inward angular momentum transport and the deepening of the differential rotation in the RZ.We note that the amplitudes and profiles of these fluxes do not only depend on σ(r) but also on Ra and the density stratification N ρ .However, we found that in all of our runs, the differential rotation for the lower σ ov runs is more contained within the equator and that F r,C > 0 there, unlike in the higher σ ov cases where both F r,C and F r,T are negative.Our aim in this work has been to examine the dynamical balances associated with axisymmetric mean flows in the radiative interior of a low-mass star such as the Sun.As discussed in Wood & Brummell (2012), in the absence of magnetism, the balance that expresses the relative importance of advective and viscous transport in the radiative interior is well-characterized by a single parameter, σ (Eq.18).Viscous effects are negligible when σ < 1, as is appropriate for a stellar interior, and are dominant when σ > 1.
To date, models in the stellar-relevant regime have only been run in axisymmetric spherical or Cartesian geometries, which cannot account for latitudinal variations of the mean flows.Spherical models incorporating a region of overshoot have in turn been run only in the σ > 1 regime.In order to bridge this gap, we have analyzed the results from a series of 3D global simulations spanning both the σ ≤ 1 and σ > 1 regimes.These models are constructed in spherical geometry, incorporate rotation, and consist of an outer convective region overlying a stably stratified layer.
Much as in the models of Garaud & Acevedo Arreguin (2009), Wood & Brummell (2012), and Acevedo-Arreguin et al. ( 2013), we find that in the σ < 1 regime appropriate for the Sun, meridional flows can penetrate deeply into the RZ.In the σ > 1 regime, they decay exponentially in depth with a lengthscale proportional to 1/σ (see Fig. 8).The deeply-and shallowly-penetrating flows arising in these two regimes lead to markedly different balances in the angular momentum transport.For increasing values of σ ov , the angular momentum flux arising from the Coriolis force becomes weaker, and for cases with σ ov >> 1, the angular momentum transport is dictated primarily by the viscous flux in the long term.In contrast, for simulations with σ ov << 1, the flux associated with the Coriolis term is much larger than the viscous flux indicating that in the solar-like regime, the mean flows are responsible for the transport of angular momentum.
The dominance of the Coriolis term in the σ ov << 1 regime is particularly notable near the equator and leads to differences in the differential rotation that develops below the convection zone.The differential rotation established in all models consists of a rapidly-rotating equator and slowly-rotating polar regions.This general profile was found to extend below the base of the CZ for all models due to either viscous diffusion (σ ov >> 1) or advection by the mean flows (σ ov ≤ 1.1).We note, however, that for models with σ ov < 1, the rapid equatorial differential rotation does not extend as deeply as it does in those with σ ov > 1 (see Fig. 3).These results suggest that meridional flows may play a role in tachocline confinement at low latitudes, though a definitive statement cannot be made here owing to the fact that this study did not consider the effects of magnetism.
Currently, even the most state-of-the-art simulations cannot capture real solar (and stellar) values due to computational constraints arising from the required spatio/temporal resolution.As such, it is important to ensure that since solar values are unattainable, they should not be used for certain input parameters (e.g.N, and Ω) and ignored for others (e.g.Pr) as this may lead to non-solar/stellar parameter regimes.Consequently, the resulting dynamics cannot effectively capture the dynamical processes taking place in the solar interior and at the same time, any extrapolations to solar/stellar parameters through derived scaling laws are dubious when they result from simulations that operate in non-solar/stellar regimes.

Implications for Dynamical Processes in the Solar Interior
As we showed in Section 2, the solar σ profile is smaller than one in the upper radiative zone and approximately equal to 0.5 around the base of the solar tachocline.Since we are reporting our results based on the location of the overshoot lengthscale, we will use our estimate of the solar overshoot depth from Korre & Featherstone (2021) and attempt to estimate the penetration lengthscale of the meridional flows in the solar radiative interior.We have to note, however, that even though our scaling law is derived from our runs that were operated within the solar-like regime where σ < 1, we still need to be cautious when extrapolating our numerical results to the Sun and stars in general.That is due to the fact that most of the other parameters have non-solar values, for instance our largest Ra = 10 6 while Ra ⊙ = 10 20 and we used Pr = 1 while Pr ⊙ = 10 −6 (e.g., Ossendrijver 2003).
With that in mind, we may proceed to provide an estimate of δ M C,⊙ .In Korre & Featherstone (2021), we assessed that the solar overshoot lengthscale is approximately equal to 0.1H p , where H p = 5.7e9 cm is the pressure scale-height measured at the base of the solar CZ (using Model S).Then the location of the overshoot depth is at r ov ≈ 0.7r ⊙ hence σ ov,⊙ ≈ 0.42 (where r ⊙ = 0.9983R ⊙ is the radius at the base of the granulation layer and R ⊙ = 6.957e10 cm).From the scaling law derived in Section 3, we can then obtain δ M C,⊙ = 0.1r ⊙ /(σ 0.22 ov,⊙ ) ≈ 0.12r ⊙ which corresponds to δ M C,⊙ ≈ 1.46H p .This result indicates that the meridional circulation might be penetrating much more deeply in the solar interior than accounted for or anticipated by previous numerical studies and observations.In fact, we find that the meridional flows might travel as far as down to r = 0.59r ⊙ = 4.098e10 cm, which is even below the overshoot region and the tachocline.
If that is indeed the case, the meridional circulation may have a greater role to play in the transport of chemical species within the RZ, as well as in the dynamical balances taking place within and beyond the tachocline region.For instance, Gough & McIntyre (1998) suggested that a primordial magnetic field in the RZ could confine the solar tachocline and enforce uniform rotation there by interacting with the large-scale mean flows gyroscopically-pumped from the CZ deeper into the RZ in a way that they could halt the magnetic field from entering the CZ.Also, our results may hold implications for flux transport models where, in tandem with the tachocline, the meridional flow plays an important role in modulating the magnetic cycle (e.g., Dikpati & Charbonneau 1999;Charbonneau 2010).We plan to explore these topics in a future study that incorporates magnetism into the models presented here.

Figure 1 .
Figure 1.Profiles of the nondimensional background entropy gradient d S(r)/dr versus the radius r/r o for five representative runs with different A S , and/or r b and/or d (see Eq. (13)).The black dashed vertical like marks the transition from the adiabatic convection zone (where d S(r)/dr = 0) to the subadiabatic layer (where d S(r)/dr > 0).

Figure 2 .
Figure 2. Profiles of σ(r) versus r/r o for the five runs with the d S/dr profiles presented in Figure 1.The horizontal dashed black line corresponds to σ(r) = 1, while the black line corresponds to the computed solar σ(r) profile with the base of the solar convection zone moved by 0.046r o to coincide with the base of the CZ of the spherical shell (see text for details).

Figure 3 .
Figure 3. Differential rotation profiles u φ /(r sin θ) for the five representative runs spanning a wide range of σ ov from a) σ ov = 0.032 up to e) σ ov = 9.8.For decreasing values of σ ov , the differential rotation is more confined to the upper radiative zone within the equatorial region.

Figure 4 .
Figure 4. Profiles of the time-and azimuthally-averaged differential rotation profiles versus the radius r/r o at different colatitudes θ for a) case 10 with σ ov = 0.1 and b) case 5 with σ ov = 9.8.

Figure 5 .
Figure 5. Time-and azimuthally-averaged profiles of the meridional circulation streamlines with underlying contours of the mass flux where blue color corresponds to clockwise motion (CW), and red color to counterclockwise motion (CCW) for the five typical runs spanning values of σ ov from a) 0.032 to e) 9.8.The meridional flows present multiple cells within the convective region and for the cases where σ ov < 1 (a) and b)), the flows penetrate more deeply into the stable layer below.In panel c) where σ ov = 1.1, there is still substantial propagation of the meridional circulation below the base of the CZ while in panels d) and e) where σ ov > 1, the mean flows only travel a small distance beyond the inner convective boundary.

Figure 6 .
Figure 6.Meridional circulation profiles at different θ.Time-and azimuthally-averaged profiles of u θ plotted from 0.65r o to r o at different colatitudes.The profiles illustrate the meridional flows in and below the base of the CZ for a) σ ov = 0.1 (case 10) and b) σ ov = 9.8 (case 5).

Figure 7 .
Figure 7. Time-and spherically-averaged profiles of the kinetic energy in the meridional flows versus r/r o .a) Profiles for the five representative runs.ẼMC is almost the same for all runs within the CZ, but changes for the different σ ov cases in the RZ.For increasing values of σ ov , ẼMC decreases much faster below the base of the CZ compared with the smaller σ ov cases.b) ẼMC profile for case 10 with σ ov = 0.1.An exponential function (see Eq. (20)) is fitted to the kinetic energy data below the base of the CZ down to the point where ẼMC is well fit.The computed decay rate λ for each simulation provides a penetration lengthscale δ M C = 2/λ for all runs shown in Table1.The inner dashed black line corresponds to the base of the CZ in both panels.

Figure 9 .
Figure 9. Time-and azimuthally-averaged gyroscopic pumping terms of Equation (21) for two σ ov cases at the same Ra.The gyroscopic pumping balance is achieved in both cases.a) In case 10 with σ ov = 0.1, the meridional flows penetrate deeply into the RZ.The gyroscopic pumping occurs across the whole shell but is more prominent within the CZ.b) In case 5 with σ ov = 9.8, the mean flows are confined to a small region below the base of the CZ.The gyroscopic pumping balance holds everywhere across the shell for this run as well, however it is negligible below the base of the CZ where the meridional flows rapidly decay.

Figure 10 .
Figure 10.Time evolution of the absolute value of the terms in Equation (21) volume-averaged from their overshoot region down to the inner RZ boundary for increasing values of σ ov spanning values from a) σ ov = 0.032 to e) σ ov = 9.8.For smaller values of σ ov , the angular momentum transport is dictated by the flux associated with the Coriolis force, however as σ ov increases, the Coriolis term becomes weaker and viscosity plays a more important role and eventually dictates the angular momentum transport in the long term.

Figure 11 .
Figure 11.Profiles of the time-and azimuthally-averaged fluxes given in Equations (25)-(29) computed at the base of the CZ for increasing values of σ ov spanning values from a) σ ov = 0.032 to e) σ ov = 9.8.For smaller (larger) values of σ ov , both F r,T and F r,C are positive (negative) near the equator indicating that the angular momentum transport is outward (inward) there.
Summary and Discussion of Our Results

Table 1 .
Input and output parameters of the simulations Case N ρ Ra Ek Ro c N r × N θ × N φ Note-Columns 2 − 5 indicate the input parameters, column 6 provides the resolution and columns 7 − 11 report on the output parameters of the simulations.