Issue 
A&A
Volume 655, November 2021



Article Number  A79  
Number of page(s)  14  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202040052  
Published online  23 November 2021 
Generation of mean flows in rotating anisotropic turbulence: The case of solar nearsurface shear layer
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: barekat@mps.mpg.de
^{2}
Department of Computer Science, Aalto University, 15400, 00076 Aalto, Finland
^{3}
Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
^{4}
Institut für Astrophysik, GeorgAugustUniversität Göttingen, 37077 Göttingen, Germany
^{5}
Department of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08544, USA
^{6}
Princeton Plasma Physics Laboratory, Princeton University, Princeton, New Jersey 08543, USA
Received:
2
December
2020
Accepted:
2
July
2021
Context. Results from helioseismology indicate that the radial gradient of the rotation rate in the nearsurface shear layer (NSSL) of the Sun is independent of latitude and radius. Theoretical models using the meanfield approach have been successful in explaining this property of the NSSL, while global direct or largeeddy magnetoconvection models have so far been unable to reproduce this.
Aims. We investigate the reason for this discrepancy by measuring the mean flows, Reynolds stress, and turbulent transport coefficients under conditions mimicking those in the solar NSSL.
Methods. Simulations with as few ingredients as possible to generate mean flows were studied. These ingredients are inhomogeneity due to boundaries, anisotropic turbulence, and rotation. The parameters of the simulations were chosen such that they matched the weakly rotationally constrained NSSL. The simulations probe locally Cartesian patches of the star at a given depth and latitude. The depth of the patch was varied by changing the rotation rate such that the resulting Coriolis numbers covered the same range as in the NSSL. We measured the turbulent transport coefficient relevant for the nondiffusive (Λeffect) and diffusive (turbulent viscosity) parts of the Reynolds stress and compared them with predictions of current meanfield theories.
Results. A negative radial gradient of the mean flow is generated only at the equator where meridional flows are absent. At other latitudes, the meridional flow is comparable to the mean flow corresponding to differential rotation. We also find that the meridional components of the Reynolds stress cannot be ignored. Additionally, we find that the turbulent viscosity is quenched by rotation by about 50% from the surface to the bottom of the NSSL.
Conclusions. Our local simulations do not validate the explanation for the generation of the NSSL from meanfield theory where meridional flows and stresses are neglected. However, the rotational dependence of the turbulent viscosity in our simulations agrees well with theoretical predictions. Moreover, our results agree qualitatively with global convection simulations in that an NSSL can only be obtained near the equator.
Key words: hydrodynamics / turbulence / Sun: rotation
© A. Barekat et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
The convection zone (CZ) of the Sun, despite being highly turbulent, shows a wellorganized largescale axisymmetric rotation profile depending on both depth and latitude. The entire CZ rotates faster at the equator than at the poles, and the rotation rate decreases mildly with depth, except near the radial boundaries, where there are regions of strong shear (Thompson et al. 1996; Schou et al. 1998). Additionally, a largescale circulation in the meridional plane, known as the meridional flow (MC), is also present. The amplitude of the MC is about 15–20 ms^{−1}, which is two orders of magnitude smaller than that of the rotational velocity (Duvall 1979; Hathaway 1996).
The nearsurface shear layer (NSSL) occupies about 17% of the CZ, or roughly 35 Mm in depth, from the photosphere. Recently, two further properties of it have been reported. First, the value of the logarithmic radial gradient of the rotation rate is reported to be
in the upper 13 Mm of the NSSL, independent of latitude up to 60° (Barekat et al. 2014). Second, the gradient evolves over time by an amount between 5–10% of its timeaveraged value. This closely follows the magnetic activity cycle (Barekat et al. 2016). On the other hand, the MC maintains its poleward motions throughout the cycle (Hathaway & Upton 2014).
Shear flows play an important role in generating and maintaining the solar magnetic field and its activity cycle (e.g., Krause & Rädler 1980). In particular, radial shear is important in the αΩ dynamo model for explaining the equatorward migration of the magnetic activity (Parker 1955; Yoshimura 1975). In this model, negative radial shear in combination with positive α is required to produce the correct equatorward migration of the activity. This negative shear exists prominently in the NSSL. The effect of the NSSL has been tested numerically in meanfield dynamo models by Käpylä & Korpi (2006), where it was found to aid equatorward migration. More observational and theoretical arguments that the NSSL strongly shapes the solar dynamo process were presented by Brandenburg (2005). The role of NSSL can be easily investigated in meanfield models, where it can be added or removed by hand. In contrast, global 3D convection simulations typically fail to selfconsistently generate a realistic NSSL (e.g., Guerrero et al. 2013; Hotta et al. 2015), and thus its role in the resulting dynamo solutions is unclear. Understanding the role that the NSSL plays for the dynamo therefore requires that we understand its formation mechanism and why global simulations do not capture it.
The equations governing the generation of largescale flows in the solar CZ are the following: First, the azimuthally averaged angular momentum equation describes the time evolution of the differential rotation. This equation is obtained using the Reynolds decomposition, where each physical quantity, A, is decomposed into its mean and fluctuations around the mean, a, and where averages are taken over the azimuthal direction. Then, we obtain the equation
where , , , ν, μ_{0}, and are density, meridional flow, angular velocity, molecular viscosity, the vacuum permeability, and the magnetic field, respectively. Furthermore, ϖ = r sin θ, where θ is the latitude, Q_{ϕi} and M_{ϕi} are the Reynolds and Maxwell stresses, and is the mean rate of the strain tensor. The Reynolds and Maxwell stresses are the correlations of fluctuating components and , respectively, where i denotes r and θ. Density fluctuations are omitted, which corresponds to an anelastic approximation.
Second, the azimuthally averaged equation for the azimuthal component of vorticity describes the time evolution of the MC:
where is the vorticity, s is the specific entropy, T is the temperature, and ∂/∂z is the derivative along the rotation axis. The first and second terms describe the centrifugal and baroclinic effects, respectively. From these two equations it becomes clear that meridional flow can drive differential rotation and vice versa, and additionally, any misalignment of the density and temperature gradients can drive meridional circulation through the baroclinic term, while turbulent stresses are important in driving both flows.
Theoretical studies have shown that the main forces generating stellar differential rotation are described by the first two terms in Eqs. (2) and (3) (Rüdiger 1989; Kitchatinov et al. 2013). Additionally, the Coriolis number Ω_{⋆}, describing the degree of rotational influence on the flow, defined as
where is the rotation rate of the star and τ is the turnover time of the turbulence, has been found to be a key parameter. Ω_{⋆} describes the role of rotation in different parts of the CZ. In particular, it leads to a completely different rotation profile within the NSSL in comparison to the rest of the CZ.
In the solar structure model of Stix (2002), Ω_{⋆} changes from the surface to the bottom of the CZ as . Helioseismic measurements (Greer et al. 2015) have revealed that in the region of the NSSL, the Rossby number (the inverse of Ω_{*}) is 0.11 at a depth of 30 Mm and 2.16 in the nearsurface layers. This shows that the rotational influence is weak in the surface regions and grows to a significant level in the bottom part of the NSSL. The accuracy of the measurements is generally very good at such small depths.
Nonrotating densitystratified convection is dominated by radial motions, in which case the vertical anisotropy parameter , where and are the fluctuating horizontal and radial velocities. u_{ϕ} and u_{θ} are the longitudinal and latitudinal fluctuating velocities. Rotation tends to suppress convection (e.g., Chandrasekhar 1961), reducing the radial motions, and typically A_{V} decreases when Ω_{⋆} increases such that the maximum of A_{V} is achieved for Ω_{⋆} = 0 (e.g., Chan 2001; Käpylä et al. 2004). A_{V} is also expected to depend on latitude due to the suppression of and the rotational influence on the horizontal components.
Rotation also introduces an anisotropy between the latitudinal and azimuthal directions. This is described by the horizontal anisotropy parameter . As has been empirically shown using numerical simulations in both local convection (Käpylä et al. 2004) and forced turbulence simulations (Käpylä & Brandenburg 2008), A_{H} is typically positive and increases with Ω_{⋆}. Its value also increases toward the equator due the increasingly deviating effect of the Coriolis force on the horizontal flows. Furthermore, A_{H} → 0 as Ω_{⋆} → 0. Thus, Ω_{⋆} in the solar CZ also reflects the anisotropy of turbulence, which arises due to the presence of the Coriolis force and density stratification. Consequently, rotation and gravity vectors define the necessary two misaligned preferred directions for nonzero offdiagonal Reynolds stress (Rüdiger 1989).
A theoretical model that reproduces the entire rotation profile of the Sun including the NSSL was presented in Kitchatinov & Rüdiger (2005, hereafter KR05). They derived an extension to the quasilinear turbulence model of Kichatinov & Rüdiger (1993) that did not take the slowrotation regime into account that is required to treat the NSSL. In the new model (KR05), the vertical anisotropy parameter A_{V} is crucial to produce strong inward transport of the angular momentum in the slowly rotating surface layers, resulting in the formation of the NSSL. The horizontal anisotropy parameter is assumed to be small, hence A_{V}≫A_{H} for Ω_{⋆} ≲ 1 applies in the surface layers, although the turbulent transport coefficients depend on latitude. They used a hydrodynamic meanfield (MF) model and studied the value of the anisotropy parameter required to reproduce the NSSL. They concluded that values that agreed reasonably well with the local magnetoconvection models were sufficient to reproduce the NSSL.
The remarkable agreement of the recently observed latitudinal independence of the angular velocity gradient with their model motivated a further development of the theory, in which the effect of the magnetic field in the NSSL was included (Kitchatinov 2016). This led to the prediction that the angular velocity gradient varies during the solar cycle, which qualitatively agrees with the observations. As the variations are caused by the magnetic field, Kitchatinov (2016) suggested that measurements of the rotational properties of the NSSL can be used as an indirect probe to measure the subsurface magnetic field. In their model, however, the Reynolds stresses were computed using a secondorder correlation approximation (SOCA), the validity of which in astrophysical regimes with high Reynolds numbers is questionable.
To avoid the necessity of using such simplifications, it is desirable to build numerical simulations of stellar convection, directly solving for the relevant, either hydro or magnetohydrodynamic, equations in spherical geometry. Such models have been developed and used since the 1970s (e.g., Gilman 1977, 1983; Glatzmaier 1985), but reproducing the NSSL is a serious challenge for these models. These global convection simulations can generate a shear layer close to the equator, mostly confined outside the tangent cylinder, where rotationaligned, elongated largescale convection cells form (see, e.g., Robinson & Chan 2001; Käpylä et al. 2011a; Guerrero et al. 2016; Warnecke et al. 2016; Matilsky et al. 2019). Only for a higherdensity stratification was a shear layer found that extended to latitudes of 60° (Hotta et al. 2015). In this case, however, the gradient of the radial shear was positive in the range 0° < θ < 45°, which contradicts the helioseismic inferences of the NSSL. Hotta et al. (2015) concluded that the meridional Reynolds stress, originating from the radial gradient of the poleward meridional flow, is the most important driver of the NSSL. In their model, the luminosity was decreased to obtain an accelerated equator, therefore the influence of rotation on convection (Coriolis number) was overestimated, and they also speculated about an unfavorable effect of the boundary conditions on their results. It is therefore unclear whether these results are indeed applicable to the NSSL. Overall, it is cumbersome to use global direct numerical simulations (GDNS) to study the origin of the NSSL because the computational cost is high, many effects are present, and it is difficult to reliably separate them from each other. For such an approach, a simpler modeling strategy is required. This is what we attempt here.
In addition to MF and GDNS models, the NSSL has also been studied from the point of view of different types of equilibria. The most recent of these models, Gunderson & Bhattacharjee (2019), considered the formation of the NSSL in a magnetohydrostatic equilibrium model, where it is driven by a poleward meridional flow near the surface. In addition to the assumption of stationarity, which does not apply because the magnetic field of the Sun is oscillatory, the model considers only nonturbulent states; nevertheless, a largescale poloidal flow, when inserted in addition to the equilibrium configuration, is seen to reduce the rotational velocity near the surface and therefore leads to an NSSLlike condition there. Miesch & Hindman (2011) accounted for the Reynolds and Maxwell stresses in the governing equations, hence allowing for turbulent effects. The authors considered a case in which an equilibrium condition exists for the angular momentum transport, Eq. (2), in which case the meridional circulation and the relevant stresses must balance. Any imbalance in the term encompassing the stresses was then postulated not only to drive differential rotation, but more importantly, to induce a meridional flow. Similarly, the azimuthal vorticity equation, Eq. (3), in a steady state was postulated not only to drive meridional flow, but more importantly, to contribute to maintaining the differential rotation profile. In the earliest scenario explaining the NSSL, Foukal & Jokipii (1975) proposed that the reason for the existence of the NSSL would be the local angular momentum conservation from rising and falling convective fluid parcels, which would lead to inward angular momentum transport. In the scenario of Miesch & Hindman (2011), however, this angular momentum transport is not a sufficient condition to sustain the NSSL. Another necessary ingredient is the meridional force balance in between the turbulent stresses and centrifugally driven circulation within the NSSL, however. In the bulk of the convection zone, a meridional force balance would be rather provided by the baroclinic effect, and the bottom of the NSSL would be determined by the transition point from baroclinic to Reynolds stress balancing. Some agreement with this scenario was found in the study by Hotta et al. (2015), whose models showed that in the region of the NSSL, the force caused by the turbulent stresses was balanced by the Coriolis force.
Our approach here is entirely different from the approaches reviewed above. We formulate a model with as few ingredients as possible to generate largescale flows to study the role of rotationinduced Reynolds stress specifically in a rotational regime relevant for the NSSL. This involves replacing convection with anisotropically forced turbulence and omitting density stratification, magnetic fields, and spherical geometry. The simplicity of the model allows an unambiguous identification of the drivers of the mean flows that can be used to assess the generation mechanisms of the solar NSSL. Additionally, this simplicity facilitates measuring the turbulent transport coefficients that can be compared with those in the MF model. This can provide insights regarding the discrepancy between MF and GDNS in obtaining the NSSL.
2. NSSL in terms of meanfield hydrodynamics
In this section we briefly explain the theory of the Λeffect and its relevance for the formation of the NSSL (Kitchatinov et al. 2013; Kitchatinov 2016). We refer to Rüdiger (1989) for a thorough treatise. In this theory, rotating and anisotropic turbulence contribute to diffusive and nondiffusive transport of angular momentum. The nondiffusive part is known as the Λeffect (Lebedinski 1941). Therefore the Reynolds stress consists of two parts,
where N_{ijkl} and Λ_{ijk} are fourth and thirdrank tensors describing the turbulent viscosity and Λeffect, respectively. In spherical geometry, Q_{rϕ}, Q_{θϕ}, and Q_{rθ} are the vertical, horizontal, and meridional stresses, respectively. We note here that the meridional stress appears only in the vorticity equation, and in the model by KR05, it does not play any role in the generation of the NSSL.
Ignoring magnetic fields, the vertical and horizontal stresses are given by
where ν_{∥} and ν_{⊥} are the diagonal and offdiagonal components of the turbulent viscosity tensor N_{ijkl}, respectively. Component ν_{⊥} is caused by the effect of the rotation on the turbulent motions (Rüdiger et al. 2019). V and H are the vertical and horizontal Λeffect coefficients, which to the lowest order are proportional to A_{V} and A_{H} (Rüdiger 1980). These coefficients are typically expanded in latitude in powers of sin^{2}θ as
In the NSSL, Ω_{⋆} ≤ 1 and A_{H} ≈ 0, such that due to the Λeffect vanishes. The offdiagonal viscosity ν_{⊥} is nonzero, but small, such that its effect is negligible (Rüdiger et al. 2019). It has been shown analytically (Kitchatinov & Rüdiger 2005) and numerically (Käpylä 2019b) that in the slowrotation regime, only the first term in the expansion of the vertical coefficient V^{(0)} survives and tends to a constant. Furthermore, when a stressfree boundary condition is applied at the radial boundaries, Q_{rϕ} = Q_{rθ} = 0. Using this in Eq. (7) and equating the diffusive and nondiffusive stresses, we obtain
which agrees reasonably well with observational results, where the radial rotational gradient is independent of latitude (Barekat et al. 2014).
3. Model
We investigated the role of Reynolds stresses in the generation of mean flows in a quite simplified model that encompasses all physics required by the MF model to generate the NSSL. The physics are A_{V}≫A_{H}, Ω_{⋆} ≤ 1, and a stressfree boundary condition in the radial direction. The first simplification comes from the geometry, in which we used cubic boxes to map the NSSL at different depths and latitudes instead of using spherical geometry, see Fig. 1. This approach certainly brings specific limitations to this study. They are the absence of the MC and the adoption of stressfree boundary conditions in models that are located at different depths of the NSSL. The consequences of the latter limitation is discussed in detail in Sect. 5.5.
Fig. 1. Schematic representation of the geometry of the current models and their relation to the NSSL. The depth of the layer is magnified for clarity. The simulation boxes are located at nine depths (not all shown) and seven latitudes. Ω_{⋆} increases gradually from the surface to the bottom of the NSSL. 
The second simplification comes from replacing selfconsistently generated convection with an artificially forced turbulent medium. We used a forcing that produces anisotropic turbulence, which has two important advantages. First, it allowed us to control the turbulent parameters, for example, the size of the eddies. Second, the desired anisotropy of the NSSL can be provided artificially without gravity. Consequently, the effect of stratification can be simply eliminated from the medium because it is difficult to separate this effect from the effect caused by Reynolds stresses in generating the mean flows.
The third simplification comes from the treatment of thermodynamic properties of the medium. We assumed the medium to be isothermal with a constant sound speed obeying , where p, ρ, and c_{s} are the pressure, density, and sound speed, respectively. We also note here that full compressibility of the equations is retained, but the effect is weak (Ma = u_{rms}/c_{s} = 0.04). This choice leads to a simpler MF analysis in which the density fluctuation can be ignored. This model was used successfully to study the properties of the Reynolds stresses in slow and fastrotating regimes by Käpylä & Brandenburg (2008) and Käpylä (2019b).
The hydrodynamic equations that we solved directly are
where D/Dt = ∂/∂t + U⋅∇ is the advective derivative and Ω = Ω_{0}(cos θ, 0, sin θ)^{T} is the rotation vector. The viscous force is given by
where is the traceless rate of the strain tensor, δ_{ij} is the Kronecker delta, and the commas denote differentiation. The forcing function is given by
where x, k, and ϕ are the position, wave vector, and a random phase, respectively. The desired vertical (z) anisotropy can be enforced using a tensorial normalization factor of the forcing, where f_{0} and f_{1} are the amplitudes of the isotropic and anisotropic parts, respectively. δt and Θ are the time step and the angle between the vertical direction z and k, respectively, and k = k determines the dominant size of the eddies. In the forcing, f_{k} is given by
which results in the forcing transversal waves; is an arbitrary unit vector. The details of the forcing can be found in Brandenburg (2001).
4. Simulation setup
We used the PENCIL CODE^{1} (Brandenburg et al. 2021) to run the simulations. We considered a cubic box with size (2π)^{3} discretized over 144^{3} grid points. z corresponds to the vertical, x to the latitudinal, and y to the azimuthal direction, respectively, the latter two being referred to as the horizontal directions. Horizontal boundaries are periodic, and stressfree conditions are imposed at vertical boundaries with
where z_{bot} and z_{top} represent the bottom and top of the domain. The box size is represented by the wave number k_{1} = 2π/L, and we chose a forcing wave number k_{f}/k_{1} = 10. The units of length, time, and density are , , and ρ_{0}, respectively, where ρ_{0} is the initially uniform value of the density. The forcing parameters f_{0} = 10^{−6} and f_{1} = 0.04 were chosen such that the effects of compressibility are weak with a Mach number Ma = u_{rms}/c_{s} ≈ 0.04 in all simulations. Moreover, with f_{1} ≫ f_{0}, we fulfilled the NSSL condition in which A_{V}≫A_{H}; see Table 1. The vigor of turbulence is quantified by the Reynolds number,
Summary of the runs in which we varied the Taylor number and latitude.
where is the root mean square of the fluctuating velocity field. When a fixed value of the kinematic viscosity is used, , the Reynolds number is about 13 for all simulations. We placed the box at seven equidistant latitudes from the equator to the pole by setting the angle θ between the rotation vector and the vertical direction as shown in Fig. 1. The vertical placement was determined by the value of Ω_{0}, which was varied such that the range of Ω_{⋆} from Eq. (4) is relevant for the NSSL. The turnover time was defined as τ = ℓ/u_{rms}, where ℓ is the size of the eddies. In our simulations the energycarrying scale of turbulence is the forcing scale ℓ = 2π/k_{f}. The Coriolis number in the simulations is therefore given by
The corresponding input parameter is the Taylor number,
The values of Ta, Ω_{⋆}, and the anisotropy parameters are given in Table 1. An additional run with Ω_{0} = 0 was made. This run was used to obtain the reference anisotropy; see Sect. 5.2. It was also used to remove a contribution to the Reynolds stress that appears in the nonrotating case; see Sect. 5.4
Mean quantities were defined as horizontal (xy) averages. The local Cartesian quantities are related to their counterparts in spherical polar coordinates via (r, θ, ϕ)→(z, x, y), , Q_{θϕ} → Q_{xy}, Q_{θr} → Q_{xz}, and Q_{rϕ} → Q_{yz}. We normalized the quantities such that and ; the tilde denotes this operation. Additionally, the error on the measured physical quantities, which were obtained directly from the simulations, was estimated by dividing the time series into three parts and comparing their timeaveraged values with the value obtained from the whole time series. The maximum deviation from the latter was considered to be the error of the measurement.
5. Results
We start our analysis by studying the properties of the velocity field. Then we apply horizontal averages to separate the mean and fluctuating velocities. Using the velocity fluctuations obtained in this way, we measure the diagonal components of the Reynolds stresses to ensure that the flow fulfills the NSSL condition A_{V}≫A_{H} in Sect. 5.2. In Sect. 5.3, by measuring the mean velocity field at all depths and latitudes, we search for the NSSLlike flow. Next, we focus on the properties of the offdiagonal components of the Reynolds stresses in Sect. 5.4, which are the only transporters of the angular momentum in our model. We also compare their properties with those in previous works of both local and relevant global convection simulations. These comparisons show the differences between our artificially generated turbulent flow with more realistic flows. This can also indicate the relevance of the excluded physics from our model. In Sect. 5.5 we investigate whether these stresses are the main driver of the observed mean flows measured in Sect. 5.3.
5.1. Velocity field
A statistically stationary turbulent state appears after about a few τ independent of Ω_{⋆} everywhere except at the equator, where the statistically stationary state is reached between a few to about 300 τ from the lowest to highest Ω_{⋆}. As an example, we show snapshots of the zonal flow normalized by the sound speed at about 1000 τ for set C46 at the equator and at a latitude of 30° in panels A and B of Fig. 2, respectively. The other components of the velocity field are very similar to the zonal component shown in panel B. The dominant scale of the turbulence is the forcing scale k_{f}/k_{1} = 10. The expected largescale zonal flow similar to the actual NSSL is generated only at the equator, as shown in panel A. All other sets show a similar behavior.
Fig. 2. Streamlines of the velocity field. The color table shows the amplitude of the azimuthal component of the velocity field normalized by the sound speed. Panels A and B show U_{y}/c_{s} at the equator and at θ = 30° for set C46, respectively. 
5.2. Anisotropy of the flow
We used the diagonal components of the Reynolds stress tensor to measure the anisotropy parameters, which are given by
We show the timeaveraged diagonal stresses of a nonrotating run in panel A of Fig. 3. All of our runs show similar profiles, regardless of θ and Ω_{⋆}. The stresses are almost constant in the entire domain, except at the boundaries, where and the horizontal components rise to twice higher values. In the interior, the values of are about twice as high as the other two components, reflecting the fact that A_{V} ≈ −0.5. This is shown in panel B of Fig. 3, where we show volume averages of A_{V} excluding the boundaries (−2 ≤ zk_{1} ≤ 2) to avoid boundary effects. As expected, A_{V} has its maximum value in the absence of rotation at Ω_{⋆} = 0.
Fig. 3. Panel A: timeaveraged and normalized diagonal components of the Reynolds stress as functions of z. The dotted (solid) line shows ( and ) of set C0. The vertical dashed lines mark the part of the domain from which A_{V} and A_{H} were measured. Anisotropy parameters A_{V}panel B and A_{H}panel C are shown as functions of Ω_{⋆} at the latitudes indicated in the legend. The diamonds in panel B show the values of A_{V} at the equator from the runs from which the mean flow was removed. 
In the slowrotation regime, Ω_{⋆} ≲ 0.2, the effect of rotation on the turbulence anisotropy is weak at all latitudes, in agreement with the theory. In deeper layers of the NSSL (Ω_{⋆} > 0.2), the effect of rotation in generating horizontal motions and suppressing vertical motions becomes more prominent, which is reflected in the two anisotropy parameters. A_{H} stays positive at all Ω_{⋆} and latitudes, which shows that independent of the two parameters. We note that the value of increases toward the equator and the bottom of the NSSL. In contrast to , increases toward higher latitudes and at bottom of the NSSL, except at 15°, where it decreases by increasing Ω_{⋆}. A_{H} increases significantly by about one order of magnitude from the top to the bottom of the NSSL at low latitudes. In contrast, A_{V} decreases by only about 15% at low latitudes and stays constant at the equator. This shows that the latitudinal dependence of A_{V} is much weaker than that of A_{H} in Ω_{⋆} < 1.
Although the latitudinal dependence of A_{V} is weak, it shows a nonmonotonous profile. For a specific Ω_{⋆}, the largest magnitude of A_{V} is obtained at the equator. A_{V} decreases significantly at a latitude of 15° and then starts to monotonically increase toward higher latitudes, but it never reaches the highest values, which are obtained at the equator. The numbers listed in Table 1 reflect the differences measured between the pole and the equator. They indicate an overall decrease in magnitude of A_{V} as a function of latitude due to this nonmonotoicity. Käpylä & Brandenburg (2008) applied fully periodic boundary conditions and retrieved a monotonic profile. We performed an additional set of runs at the equator without the mean flows, in which case the monotonic behavior was recovered. The results of these runs are shown in panel B of Fig. 3 and are indicated with a diamond. This shows that the nonmonotonicity is caused by the strong equatorial mean flows generated in our model that are not permitted by the periodic boundary conditions. We note that the GDNS models commonly exhibit thermal Rossby waves (also known as Busse columns or banana cells) close to the equator (e.g., Käpylä et al. 2011b; Guerrero et al. 2013; Hotta et al. 2015; Matilsky et al. 2019), and the mean flows seen in our stressfree boundary condition model might reflect their appearance. Studying this phenomenon is beyond the scope of this study, however and will be addressed at a later time.
Although the anisotropy of the flow changes with increasing rotation rate, we see from comparing the absolute values of A_{V} and A_{H} that the model fulfills the NSSL condition A_{V}≫A_{H}. Therefore the model has all ingredients that were suggested by the MF model as necessary to generate the NSSL.
5.3. Mean flows
The development of mean flows in rotating cases means that it takes significantly longer to reach a statistically steady state than in nonrotating runs. Furthermore, long time averages are needed for a statistical convergence of the turbulent quantities. We ran all the simulations for at least 1100 turnover times. As an example, we show a subset of the time evolution of the three components of the normalized mean velocity field for about 1200 turnover times for set C24 at the equator and at θ = 15° in Fig. 4. At the equator, a large zonal flow with a negative vertical gradient developed gradually over 100τ, as shown in panel A. All other sets show a similar zonal flow profile at the equator, but the amplitude and steepness of the gradient increase with increasing Ω_{⋆}. Farther away from the equator, the amplitude of the mean zonal flow is significantly reduced and the negative gradient disappears, as shown in panel B of Fig. 4. The dependence of the mean zonal flow on rotation is shown in the panel A of Fig. 5, where we show the timeaveraged at selected Ω_{⋆} at a latitude of 15°. When Ω_{⋆} increases, the gradient of changes sign and becomes steeper up to Ω_{⋆} = 0.46, then it becomes shallower, and it slowly vanishes in the middle at Ω_{⋆} = 1. The latitudinal dependence of is shown for sets C06 and C46 in panels C and E of Fig. 5, respectively. We find that decreases as a function of latitude, vanishes at the poles, and the amplitude is lower than 5% of u_{rms} everywhere except at the boundaries.
Fig. 4. Normalized mean components of the velocity field vs. time in terms of turnover time in representative runs in set C24. The rows from top to bottom show , , and . The left and right columns show the mean velocities at the equator and at a latitude of 15°, respectively. To make the comparison of the velocity components feasible, we clipped the values of the color table of panel A at 50% of the maximum value. 
The timeaveraged meridional component of the mean flow is consistent with zero at the equator for all runs similar to the run shown in panel C of Fig. 4. In contrast to the zonal flow, its value increases away from the equator; see panel D of Fig. 4. The timeaveraged value of this component at 15° is shown in panel B of Fig. 5 for selected values of Ω_{⋆}. The negative gradient persists up to Ω_{⋆} = 0.24. Above this Ω_{⋆}, the shear slowly vanishes at the center of the box and becomes slightly positive at increasing Ω_{⋆}. However, the strong shear persists only near the boundaries. We show the latitudinal dependence of for the two sets C06 and C46 in panels D and F of Fig. 5. The amplitude of decreases as a function of latitude. The amplitudes of and are comparable everywhere except at the equator, and the negative gradient of for Ω_{⋆} < 0.1 persists at all latitudes. These results clearly show that in contrast to the MF model, the NSSLlike flow is generated only at the equator. It is also interesting that although our model is significantly simpler than the GDNS models and there is no connection between latitudes, it reproduces results similar to the GDNS.
Fig. 5. Timeaveraged normalized mean velocity components vs. vertical direction. panels A and B show and at 15°. The second and third rows show the mean horizontal velocities for sets C06 and C46, respectively, at the latitudes indicated in the legends. 
For completeness, we show the vertical component of the normalized mean flow in the bottom row of Fig. 4. All runs show a similar pattern of highfrequency oscillations for regardless of latitude and Ω_{⋆} with amplitudes of about 10^{−4}u_{rms}. These oscillations are identified as longitudinal sound waves, as expected for a compressible system in a vertically confined cavity.
5.4. Reynolds stresses
For zero rotation, it is expected that , see Eq. (8). However, we find that always has a small but nonzero value that also persists in the longest time series of our data. We find that this contribution is present and its magnitude remains unchanged for higher resolutions as well (288^{3} and 576^{3} grids). This issue therefore does not appear to be caused by a numerical convergence issue. We have been unable to identify whether the cause is due to some as yet unidentified physical effect, for example, caused by compressibility, the forcing, inhomogeneities in the system, or a combination thereof. Because this contribution is systematically present, we used the nonrotating run (C0) and subtracted from that run from the results of the runs with rotation.
We show representative results of the offdiagonal stresses at five selected Ω_{⋆} at 15° (left column) and spatially averaged (−0.5 ≤ zk_{1} ≤ 0.5) as function of latitude (right column) in Fig. 6. The Reynolds stresses for all sets are available in Appendix A. The vertical Reynolds stress at all latitudes shows similar profiles as at 15°, see panel A. The stress is nearly constant in the interior of the domain and tends to zero at the boundaries. is always negative, independent of Ω_{⋆} and latitude, as shown in panel B. Thus, the vertical angular momentum transport is inward, in agreement with previous studies (e.g., Pulkkinen et al. 1993; Chan 2001; Käpylä et al. 2004; Kitchatinov & Rüdiger 2005; Käpylä & Brandenburg 2008; Käpylä 2019b). Independent of Ω_{⋆}, the vertical stress vanishes at the pole and has its minimum and maximum amplitude at the equator and 15°, respectively, after which it gradually decreases toward the pole. For a given Ω_{⋆}, its amplitude is about twice larger at 30° latitude than 60°. The latitudinal dependence of Q_{yz} is different from that in previous studies by Pulkkinen et al. (1993) and Käpylä et al. (2004) at Ω_{⋆} ≈ 1, in which they measured Q_{yz} from local box convection simulations. Pulkkinen et al. (1993) reported that the latitudinal dependence was almost constant up to 60° and decreased toward higher latitudes. Q_{yz} has a vshape profile in latitude with the minimum at 45° in Käpylä et al. (2004). The main ingredient that is missing in our forced turbulence simulation in comparison with theirs is density stratification. Moreover, Käpylä et al. (2004) included the overshooting layer below the CZ. This makes it difficult to determine what causes our results to be different from theirs.
Fig. 6. Left column: timeaveraged offdiagonal Reynolds stresses vs. vertical direction at five selected Ω_{⋆} indicated by the legend at a latitude of 15°. Right column: stresses shown in the left panels, further spatially averaged (−0.5 ≤ zk_{1} ≤ 0.5), at different latitudes. The rows from top to bottom show , , and . 
The middle panels C and D in Fig. 6 show horizontal stress Q_{xy}. The signature of turbulent fluctuations at the forcing scale is more clearly visible in this component, and the measurement is quite noisy. The values of Q_{xy} are close to zero up to Ω_{⋆} = 0.46, above which they slowly start to become positive (negative) values in the middle (close to the boundaries). We note here that Q_{xy} obtains notable values at similar Ω_{⋆} as A_{H}. This clearly shows that the generation of horizontal Reynolds stress depends on the presence of the horizontal anisotropy in the flow. However, A_{H} being maximum at certain latitude does not necessarily mean that Q_{xy} also has its maximum value at that latitude. This can be seen by comparing the maximum values in panel D of Fig. 6 with the maximum values of A_{H} shown in panel C of Fig. 3. At a given Ω_{⋆}, the profile of Q_{xy} is similar at all latitudes. Its amplitude is maximum at 30° and gradually decreases toward the equator and the pole, as shown in panel D. This result agrees with the observational measurements of Q_{xy} using sunspot proper motions (Ward 1965; Gilman & Howard 1984; Pulkkinen & Tuominen 1998), but not with the one measured using supergranulation motions, see Fig. 10 in Hanasoge et al. (2016). The horizontal stress has always positive values independent of Ω_{⋆} and latitude, in agreement with previous studies in the slowrotation regime (e.g. Kitchatinov & Rüdiger 2005; Käpylä & Brandenburg 2008; Käpylä 2019b). The latitudinal profile of Q_{xy} measured by Pulkkinen et al. (1993) is very similar to our results, but has negative values because their box is located in the southern hemisphere; see their Fig. 6.
The meridional stress is shown in the last row of Fig. 6. In contrast to the other stresses, shows a complicated profile, in particular, close to the boundaries. Moreover, it has positive or negative values depending on both Ω_{⋆} and θ, as shown in panel E. The latitudinal dependence of the meridional stress is shown in panel F. At Ω_{⋆} < 0.1, is positive at low latitudes and above 45°. By increasing Ω_{⋆}, moves toward negative values and its absolute value increases. For Ω_{⋆} > 0.24, it has its maximum amplitude at about 45° and decreases toward the pole and the equator, similar to . The meridional stress in Pulkkinen et al. (1993) also shows a sign change in agreement with ours, but in midlatitude. However, the sign change occurs at Ω_{⋆} ≈ 1 while our results show only negative values at that Ω_{⋆}. Comparing the absolute amplitude of the stresses in right column of Fig. 6, we see that is always larger than and . For example, at Ω_{⋆} = 0.64, is about 2 to 10 times larger than and 5 to 20 times larger than depending on the latitudes. When we also compare the absolute amplitude of and , we see that for all Ω_{⋆}. These results show that although Q_{xy} increases as a function of Ω_{⋆}, its values are still much lower than the vertical stresses, which agrees with the assumption of KR05 regarding the NSSL.
Although our model is quite simple in comparison to the GDNS, it is of interest to compare the Reynolds stresses with simulations such as those in Käpylä et al. (2011b). These authors modeled turbulent convection in a spherical wedge for a variety of rotation rates. Considering the runs of Käpylä et al. (2011b) with Ω_{⋆} < 1, we find good agreement for the horizontal stress Q_{xy}, which is small and positive for small Ω_{⋆}, and which has appreciable values only for Ω_{⋆} > 0.5. However, we find maximu values at 30° instead of at 10…15° in Käpylä et al. (2011b). We also observe a similar trend for Q_{xz} such that it is positive for small Ω_{⋆} in the northern hemisphere with a sign change after certain Ω_{⋆}. However, this trend depends on latitude in their case; see their Fig. 8. The profile of Q_{yz} in the convection simulations is quite different from ours, such that it has a strong latitudinal dependence and has both positive and negative values depending on Ω_{⋆} and latitude. This is consistent with earlier studies (e.g., Käpylä 2019b), in which a sign change of Q_{yz} occurs at higher Ω_{⋆} than those considered in the present simulations.
5.5. Role of Reynolds stresses in generating the mean flows
As the Reynolds stresses appear in the MF momentum equation, we start by writing the MF equations for and using Eq. (13). We wrote these equations first by considering that our setup is fully compressible and that the forcing we used here is not solenoidal, which might cause density fluctuations that cannot be ignored in the MF equations. These considerations led to three terms in addition to the Reynolds stresses in the momentum equation (e.g., Käpylä et al. 2020). We compared all of them with the Reynolds stresses. They and their gradients are considerably smaller than the Reynolds stresses. Therefore we can ignore the density fluctuations, and the final set of equations read
Omitting terms proportional to the small quantities ν and , and Ω_{y} = 0, yields the final form of the equations,
We verified the validity of the MF equations by considering the steadystate solution, which reads
We show the horizontal mean velocities in comparison with the RHS of Eqs. (27) and (28) from a latitude of 30° in set C46 in the upper panel of Fig. 7. These results are representative of all nonequatorial cases. Although there are fluctuations in the gradient of the Reynolds stresses, the match is satisfactory.
Fig. 7. Timeaveraged mean velocities and and their corresponding balancing terms in Eqs. (27) and (28) at latitude of 30° (upper panel) and Eq. (30) at the equator (lower panel) in the vertical direction for set C46. In the upper panel, the orange and blue lines show and , respectively. The red and black lines show the RHS of Eqs. (28) and (27), respectively. In the lower panel, the solid and dotted lines show the LHS and RHS of Eq. (30), respectively. 
The equator is a particular case, and Eq. (27) cannot be used because Q_{xz} and Ω_{z} are both zero there. Therefore we need to use the third component of the MF momentum equation. Applying similar elimination of the terms as for Eqs. (23) and (24), we have
The pressure gradient appears in this equation due to horizontal averaging. In the steady state, the zonal flow can be written as
We show both sides of Eq. (30) in the lower panel of Fig. 7. The good correspondence indicates that these equations can be used to investigate the role of the stresses in generation of the mean flows.
We emphasize that although in the steady state, for example, at the equator, the two terms on the RHS of Eq. (30) balance, these terms are not the generators of the mean flow. They do, however, determine the final amplitude of the flow. Instead, the mean flows are generated by the gradient of the vertical stress Q_{yz} at the vertical boundaries, as can be seen from Eq. (26). This flow then slowly penetrates to the middle of domain. This behavior is also clearly shown in the first panel of Fig. 4, where we show the time evolution of .
The generation of mean flows is straightforward at the equator because the meridional stress, and hence the meridional flow, vanish there. At other latitudes, the meridional stress and flow have to be included, but it is clear that the Reynolds stresses are the main driver of mean flows in the current setups.
6. Parameterization of Reynolds stresses in terms of meanfield hydrodynamics
Based on the Λeffect theory explained in Sect. 2, the vertical and horizontal Reynolds stresses given in Eqs. (7) and (8), respectively, can be written in the simulation domain as
Measuring the Λeffect coefficients V and H from a single experiment is not possible because the turbulent viscosities ν_{∥} and ν_{⊥} are also unknown. We compensated for this by running another set of otherwise identical simulations, but suppressed the horizontal mean flows artificially at each time step. Therefore the first terms in Eqs. (31) and (32) go to zero. From these simulations, we directly measured Q^{(Λ)}. However, we needed to validate this approach because the velocities can be affected by the nonlinearity of the NavierStokes equations. To do this, we performed yet another set of otherwise identical simulations, but used periodic boundary conditions (PBC) in all directions instead of the stressfree boundary condition (SFBC) in the vertical direction. Then we compared the two sets of stresses obtained with these sets of boundary conditions. This comparison of varying boundary conditions is also important with respect to interpreting the Ω_{⋆} dependence as a depth dependence. This approach is somewhat artificial because we practically enforced unrealistic BCs within the convection zone.
As an example, we show the horizontal mean velocities for the PBC version of C46 at a latitude of 30° in panel A of Fig. 8. Clearly, no notable mean flow is generated in this run. Therefore the first term in both Eqs. (31) and (32) goes to zero, similar to the cases in which the mean flows are suppressed. In panel B of Fig. 8, we show the results of comparing the Reynolds stresses in PBC and SFBC cases. The difference caused by varying boundary conditions is confined to a very narrow layer near the boundary. These results suggest that our method for the separation of different effects and enforcing artificial SFBC at different depths is valid.
Fig. 8. Panel A: timeaveraged normalized mean velocity vs. vertical direction of the PBC run for set C46 at a latitude of 30°. The black and red lines show and , respectively. Panel B: comparison of the timeaveraged normalized stresses obtained from periodic and stressfree boundary condition of the same run. The solid and dashed lines show the measured (red), (blue) and (black) from SFBC and PBC runs, respectively. 
Considering Eq. (31), the subtraction of the Reynolds stresses obtained from these simulations from the total ones gives
Measuring the vertical gradient of , the value of ν_{∥} can be determined by performing an errorweighted linear leastsquares fit to Eq. (33). Inserting the measured values of ν_{∥} back into Q^{(Λ)} of both Eqs. (31) and (32), we can measure V and H provided that ν_{⊥} ≪ ν_{∥}.
6.1. Properties of the diffusive and nondiffusive parts of Reynolds stresses
Similar to Sect. 5.4, we first measured Q_{ij} from a nonrotating run and then subtracted its mean value from the corresponding stress in other sets. We show the different contributions to the Reynolds stresses in Fig. 9. In the left column we show stresses from one or two simulation sets, and in the right column we show the dependence of volume averages, excluding the boundaries (−2 ≤ zk_{1} ≤ 2), of on both latitude and Ω_{⋆}. In panel A we show the vertical stresses for set C24 at the equator and at a latitude of 30°. With these results, we can explain the minimum of Q_{yz} at the equator: the diffusive and nondiffusive parts of the stresses are comparable but of opposite signs, leading to a low negative value for the total. With Eq. (31), we see that ν_{∥} > 0, which in combination with , gives . Moreover, the final negative value of Q_{yz} also shows that generates the zonal flow. The profile of for all other latitudes is similar to the profile at a latitude of 30°, and it shows that the main contribution from the diffusive part occurs close to the boundaries at zk_{1}≳2 with positive values. Furthermore, the amplitude of decreases toward higher latitudes (not shown). In the middle of the domain, it has negative values, which fits well, as is shown in Fig. 5. The nondiffusive part of the vertical stress is always . Its absolute value increases from the pole toward the equator and increases with Ω_{⋆}. We also find that is linearly dependent on Ω_{⋆} in the slowrotation regime, in agreement with previous numerical results (Käpylä 2019a).
Fig. 9. Panels A, C, and E: timeaveraged diffusive and nondiffusive parts of the Reynolds stresses vs. vertical direction. The black and blue lines in panel A show the normalized vertical stresses at the equator and a latitude of 30° for set C24, respectively. In panel C the horizontal stresses are shown at a latitude of 30° for set C64. The blue and black lines in panel E show the meridional stresses for sets C06 and C64 at a latitude of 30°, respectively. The vertical lines denotes the z range used for volume averages. Solid, dotted, and dashdotted lines show , , and , respectively. Panels B, D, and F: Volume averages over −2 ≤ zk_{1} ≤ 2, of vs. Ω_{⋆} at different latitudes as indicated by the legend. 
We show corresponding results for Q_{xy} in panel C of Fig. 9 for Ω_{⋆} = 0.64. has positive values in the whole domain, while is almost zero in the middle of the domain, and its contribution to Q_{xy} is confined to the boundaries at zk_{1}≳2. This also shows that is the main contributor to the negative values of Q_{xy} close to the boundaries shown in Fig. 6 C. The volumeaveraged values of , excluding the boundaries, are shown in Fig. 9 D as a function of both Ω_{⋆} and latitude. Its value is almost zero at the equator and at the pole. It is significantly nonzero above Ω_{⋆} > 0.24 and increases with increasing Ω_{⋆} independent of latitude. Independent of Ω_{⋆}, it has maximum value at a latitude of 30°. We note here that the amplitude of is also significantly smaller than . The measured profile of is almost identical to the profile obtained by Käpylä (2019b).
Our results for Q_{xz} are shown in Fig. 9. At low Ω_{⋆}, contributes almost nothing to the total stresses. For Ω_{⋆} > 0.15, the contribution of disappears in the middle of the domain, but maintains its positive value close to the boundaries. This is shown in panel E, where we plot Q_{xz} for low and high Ω_{⋆} for the sake of comparison. In panel F, we show volume averages of at all Ω_{⋆} and latitudes. The value of is almost zero at the equator and at the pole. At other latitudes, its absolute value increases with increasing Ω_{⋆}. It has always negative values, independent of both Ω_{⋆} and latitude. These results agree with those of Käpylä (2019b) in the slowrotation regime.
6.2. Measuring turbulent viscosity
The diagonal turbulent viscosity ν_{∥}, normalized by ℓu_{rms}, and its dependence on both Ω_{⋆} and latitude is shown in panels A and B of Fig. 10, respectively. Except for the highest latitudes, where measurements are unreliable, the turbulent viscosity decreases monotonically as a function of Ω_{⋆} such that for the largest Ω_{⋆}, corresponding to bottom of the NSSL, its value has decreased by roughly a factor of two. The method used here to measure turbulent viscosity relies on the presence of mean flows. As these diminish toward high latitudes, it is very difficult to obtain reliable estimates of ν_{∥} near the pole. We note that the measurements of also suffer from numerical noise at Ω_{⋆} < 0.1 at low latitudes. In particular, the latitudinal dependence of for θ ≲ 60°, shown in panel B, is probably not reliable. According to the results at lower latitudes, we conclude that the latitude dependence is weaker than the rotational dependence. As ν_{∥} is measured with high confidence, we consider its profile at the equator applicable for other latitudes. We used it to measure V and H at other latitudes. The ratio of turbulent to kinematic viscosity is ν_{∥}/ν ∼ 10–20, as expected for the fluid Reynolds numbers in the current simulations.
Fig. 10. Normalized turbulent viscosity and Λeffect coefficients as a function of Ω_{⋆} and latitude. Panels A, C, and D show , V, and H as a function of Ω_{⋆} from the equator to a latitude of 75°, respectively. Panel B shows as a function of latitude for five selected Ω_{⋆}. 
Käpylä et al. (2020) measured turbulent viscosity from nonrotating isotropically forced turbulence with shear. Using the same normalization as in their study with ν_{t0} = u_{rms}/3k_{f}, we obtain ν_{t}/ν_{t0} ≈ 3.5…3.8 in the slowly rotating simulations in sets C02, C04, and C06. These values are roughly twice higher than those in Käpylä et al. (2020). The cause of the difference is unclear.
We also compared the profile of ν_{∥} with an analytical expression for the rotation dependence of the viscosity obtained under SOCA by Kitchatinov et al. (1994, hereafter KPR94). We considered the first term in Eq. (34) of their work, which is relevant to our simulations in which ν_{∥} = ν_{0}ϕ_{1}(Ω_{⋆}), where ν_{0} = 4/15ℓu_{rms} is the turbulent viscosity obtained for the isotropic nonrotating case, and where ϕ_{1} is a function of Ω_{⋆} given in the appendix of KPR94. We scaled the analytical result by a factor of κ = 0.68 to make it comparable with our numerical result. In Fig. 11 we show the result of this comparison. This result shows that except for the κ factor, the rotation dependence agrees fairly well between the theory and numerical simulations.
Fig. 11. Comparison of the obtained turbulent viscosity with the analytical result of KPR94. The solid and dashed lines show the normalized turbulent viscosity and rescaled analytical expression κϕ_{1}, respectively. 
Considering the offdiagonal turbulent viscosity ν_{⊥}, we failed to measure it because the two terms that constitute it, and , are too small. The measurement error in the former is large.
6.3. Measurements of the vertical Λeffect coefficient
We measured the vertical Λcoefficient by substituting the volume averages of shown in panel A of Fig. 9 and ν_{∥} at the equator using
Our results are shown in panel C of Fig. 10. V is about 0.75 and gradually increases to ≈0.95 for latitudes ≤45°. However, the value of V at the lowest Ω_{⋆} is lower at all latitudes, but it has large error bars. In contrast to low latitudes, values of V at latitudes of 60° and 75° decrease for Ω_{⋆} > 0.3. Considering the large errors in the measurements at low Ω_{⋆}, we might consider V to be roughly constant for Ω_{⋆} ≤ 0.15 independent of latitude, but it shows a strong latitudinal and rotational dependence for Ω_{⋆} > 0.15. This means that considering the first term V^{(0)} in Eq. (9) in the NSSL condition alone is not enough, as is assumed in the theoretical model by KR05 explained in Sect. 2. Moreover, the increase in V toward higher Ω_{⋆} at low latitudes is in contrast with the decrease predicted in KR05 model. The same applies to the results of Käpylä (2019b), who did not consider that ν_{t} = ν_{t}(Ω_{⋆}).
6.4. Measurements of the horizontal Λeffect coefficient
We measured the horizontal Λeffect coefficient similarly to the vertical coefficient using
The results are shown in panel D of Fig. 10. The values of H are always positive, independent of Ω_{⋆} and latitude. Its values are one order of magnitude lower than V up to Ω_{⋆} = 0.6, above which H begins to increase at latitudes < 45°. We also note that its value is zero at the equator and at the pole. H is largest at a latitude of 15° and gradually decreases toward higher latitudes. These results show that close to the surface, H does not play any role in transporting the angular momentum, which validates the assumption applied in the NSSL model by KR05.
7. Conclusions
We applied an alternative approach to the MF and GDNS, that is, we ran direct numerical simulations of forced turbulence in local boxes, to primarily determine whether the assumptions and approximations applied in MF theory to explain the formation of the NSSL are valid. In contrast to the GDNS, we were able to isolate and study the role and contribution of the Reynolds stresses in the rotational regime relevant for the NSSL. Additionally, we were able to measure the turbulent viscosity. Our results show that the three required conditions explained in Sect. 2, which are necessary to generate the NSSL in the RK05 model, are insufficient. In particular, the meridional component of the Reynolds stress cannot be ignored. However, our results are in accordance with Q_{xy} → 0 in the upper part of the NSSL, whereas Q_{xy} obtains small but nonzero values close to the bottom of the NSSL, in agreement with theoretical predictions. The role of the vertical Reynolds stress in transporting the angular momentum radially inward agrees with the theory. However, its profile differs from that predicted by the theory. In particular, Kitchatinov et al. (2013) and Kitchatinov (2016) assumed that only the term V^{(0)} survives in the expansion of V in the NSSL. However, our results indicate that higherorder terms in the expansion of V need to be considered. Moreover, it is also expected from theory that the inward angular momentum flux (V) decreases with increasing Ω_{⋆} at all latitudes, but our results show that this expectation is fulfilled only at high latitudes. We also note here that the rotational quenching of the turbulent viscosity, ν_{∥}, adds another degree of complexity to the problem that was not considered previously in the models of NSSL. From the theoretical MF prediction (Kitchatinov et al. 1994), this behavior agrees qualitatively well with our results, however.
Although these local box simulations have a moderate value of Re ≈ 13 and there is no connection between different latitudes, our results are largely consistent with the stresses and mean flows obtained in the GDNS. On the other hand, the theoretical works used SOCA, which should be valid at Reynolds or Strouhal numbers of up to unity, which is in the vicinity of the parameter regime of the current models. It is expected therefore that we find a relatively good match between the measured turbulent viscosity and the viscosity predicted by SOCA.
Q_{xz} cannot be disregarded in the NSSL. Its role can be further investigated in more realistic setups using spherical geometry where the artifact of the discontinuity between latitudes can be removed. We also note here that we considered only a single modest Reynolds number and one forcing scale, the effects of which need to be explored with wider parameter studies. The other important physics that needs to be investigated are the effects of stratification, compressibility, and magnetic fields, and this needs to be compared with previous studies that have studied them in turbulent convection, namely Pulkkinen et al. (1993), Chan (2001), and Käpylä et al. (2004).
It is worthwhile to note here that a set of companion laboratory experiments is being proposed to test several aspects of our model. In these experiments, a rotating waterfilled apparatus will be used to simulate regions of finite latitudinal extent, including βplane effects, and forcing will be introduced by pumpdriven nozzles at the boundaries (Burin et al. 2019). The relative variation of the system rotation rate and the nozzle exit velocity will allow both the Ω_{⋆} > 1 and Ω_{⋆} < 1 regimes to be explored. The forcing scale length and isotropy will be changed by opening and closing nozzles and by altering the nozzle shapes and orientations. Timeresolved measurements of the components of the flow velocity will allow the mean flows and stresses to be computed and compared with numerical results and theoretical models. Although the details of the forcing and the fluid boundary conditions will be different in the experiment compared to our simulations, it is expected that meaningful results will be obtained as the rotation rate of the system is varied and the experimental data are analyzed to search for signatures of the Λeffect.
Acknowledgments
The simulations have been carried out on supercomputers at GWDG and on the Max Planck supercomputer at RZG in Garching. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project “UniSDyn”, grant agreement No. 818665). A. Barekat and E. Gilson acknowledge funding by the MaxPlanck/Princeton Center for Plasma Physics. P. J. Käpylä acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG) Heisenberg programme (grant No. KA 4825/21).
References
 Barekat, A., Schou, J., & Gizon, L. 2014, A&A, 570, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barekat, A., Schou, J., & Gizon, L. 2016, A&A, 595, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brandenburg, A. 2001, ApJ, 550, 824 [Google Scholar]
 Brandenburg, A. 2005, ApJ, 625, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Johansen, A., Bourdin, P. A., et al. 2021, J. Open Source Softw., 6, 2807 [CrossRef] [Google Scholar]
 Burin, M. J., Caspary, K. J., Edlund, E. M., et al. 2019, Phys. Rev. E, 99, 023108 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. L. 2001, ApJ, 548, 1102 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: Clarendon) [Google Scholar]
 Duvall Jr., T. L. 1979, Sol. Phys., 63, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Foukal, P., & Jokipii, J. R. 1975, ApJ, 199, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A. 1977, Geophys. Astrophys. Fluid Dynam., 8, 93 [CrossRef] [Google Scholar]
 Gilman, P. A. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A., & Howard, R. 1984, Sol. Phys., 93, 171 [NASA ADS] [Google Scholar]
 Glatzmaier, G. A. 1985, ApJ, 291, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Greer, B. J., Hindman, B. W., Featherstone, N. A., & Toomre, J. 2015, ApJ, 803, L17 [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2016, ApJ, 819, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Gunderson, L. M., & Bhattacharjee, A. 2019, ApJ, 870, 47 [CrossRef] [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Ann. Rev. Fluid Mech., 48, 191 [Google Scholar]
 Hathaway, D. H. 1996, ApJ, 460, 1027 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., & Upton, L. 2014, J. Geophys. Res., 119, 3316 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 798, 51 [Google Scholar]
 Käpylä, P. J. 2019a, Astron. Nachr., 340, 744 [CrossRef] [Google Scholar]
 Käpylä, P. J. 2019b, A&A, 622, A195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., & Brandenburg, A. 2008, A&A, 488, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2006, Astron. Nachr., 327, 884 [CrossRef] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2004, A&A, 422, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011a, Astron. Nachr., 332, 883 [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Guerrero, G., Brandenburg, A., & Chatterjee, P. 2011b, A&A, 531, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Rheinhardt, M., Brandenburg, A., & Käpylä, M. J. 2020, A&A, 636, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kichatinov, L. L., & Rüdiger, G. 1993, A&A, 276, 96 [NASA ADS] [Google Scholar]
 Kitchatinov, L. L. 2013, in IAU Symposium, eds. A. G. Kosovichev, E. de Gouveia Dal Pino, & Y. Yan, 294, 399 [Google Scholar]
 Kitchatinov, L. L. 2016, Astron. Lett., 42, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., Pipin, V. V., & Ruediger, G. 1994, Astron. Nachr., 315, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 2005, Astron. Nachr., 326, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, F., & Rädler, K. H. 1980, Meanfield magnetohydro dynamics and dynamo theory (Oxford: Pergamon Press) [Google Scholar]
 Lebedinski, A. I. 1941, Astron. Zh., 18, 10 [Google Scholar]
 Matilsky, L. I., Hindman, B. W., & Toomre, J. 2019, ApJ, 871, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Hindman, B. W. 2011, ApJ, 743, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1955, ApJ, 122, 293 [Google Scholar]
 Pulkkinen, P., & Tuominen, I. 1998, A&A, 332, 755 [NASA ADS] [Google Scholar]
 Pulkkinen, P., Tuominen, I., Brandenburg, A., Nordlund, A., & Stein, R. F. 1993, A&A, 267, 265 [NASA ADS] [Google Scholar]
 Robinson, F. J., & Chan, K. L. 2001, MNRAS, 321, 723 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G. 1980, Geophys. Astrophys. Fluid Dynam., 16, 239 [CrossRef] [Google Scholar]
 Rüdiger, G. 1989, Differential rotation and stellar convection. Sun and the solar stars (Berlin: Akademie Verlag) [CrossRef] [Google Scholar]
 Rüdiger, G., Küker, M., Käpylä, P. J., & Strassmeier, K. G. 2019, A&A, 630, A109 [EDP Sciences] [Google Scholar]
 Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [Google Scholar]
 Stix, M. 2002, The sun: an introduction (Berlin: Springer) [Google Scholar]
 Thompson, M. J., Toomre, J., Anderson, E. R., et al. 1996, Science, 272, 1300 [Google Scholar]
 Ward, F. 1965, ApJ, 141, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2016, A&A, 596, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yoshimura, H. 1975, ApJ, 201, 740 [Google Scholar]
Appendix A: Reynolds stresses
We show the measured , and at all latitudes and Ω_{⋆} in Fig. A.1, Fig. A.2, and Fig. A.3, respectively.
Fig. A.1. vs. vertical direction at all Ω_{⋆} indicated by the legend. Each panel shows at a certain latitude. 
Fig. A.2. vs. vertical direction at all Ω_{⋆} indicated by the legend. Each panel shows at a certain latitude. 
Fig. A.3. vs. vertical direction at all Ω_{⋆} indicated by the legend. Each panel shows at a certain latitude. 
All Tables
All Figures
Fig. 1. Schematic representation of the geometry of the current models and their relation to the NSSL. The depth of the layer is magnified for clarity. The simulation boxes are located at nine depths (not all shown) and seven latitudes. Ω_{⋆} increases gradually from the surface to the bottom of the NSSL. 

In the text 
Fig. 2. Streamlines of the velocity field. The color table shows the amplitude of the azimuthal component of the velocity field normalized by the sound speed. Panels A and B show U_{y}/c_{s} at the equator and at θ = 30° for set C46, respectively. 

In the text 
Fig. 3. Panel A: timeaveraged and normalized diagonal components of the Reynolds stress as functions of z. The dotted (solid) line shows ( and ) of set C0. The vertical dashed lines mark the part of the domain from which A_{V} and A_{H} were measured. Anisotropy parameters A_{V}panel B and A_{H}panel C are shown as functions of Ω_{⋆} at the latitudes indicated in the legend. The diamonds in panel B show the values of A_{V} at the equator from the runs from which the mean flow was removed. 

In the text 
Fig. 4. Normalized mean components of the velocity field vs. time in terms of turnover time in representative runs in set C24. The rows from top to bottom show , , and . The left and right columns show the mean velocities at the equator and at a latitude of 15°, respectively. To make the comparison of the velocity components feasible, we clipped the values of the color table of panel A at 50% of the maximum value. 

In the text 
Fig. 5. Timeaveraged normalized mean velocity components vs. vertical direction. panels A and B show and at 15°. The second and third rows show the mean horizontal velocities for sets C06 and C46, respectively, at the latitudes indicated in the legends. 

In the text 
Fig. 6. Left column: timeaveraged offdiagonal Reynolds stresses vs. vertical direction at five selected Ω_{⋆} indicated by the legend at a latitude of 15°. Right column: stresses shown in the left panels, further spatially averaged (−0.5 ≤ zk_{1} ≤ 0.5), at different latitudes. The rows from top to bottom show , , and . 

In the text 
Fig. 7. Timeaveraged mean velocities and and their corresponding balancing terms in Eqs. (27) and (28) at latitude of 30° (upper panel) and Eq. (30) at the equator (lower panel) in the vertical direction for set C46. In the upper panel, the orange and blue lines show and , respectively. The red and black lines show the RHS of Eqs. (28) and (27), respectively. In the lower panel, the solid and dotted lines show the LHS and RHS of Eq. (30), respectively. 

In the text 
Fig. 8. Panel A: timeaveraged normalized mean velocity vs. vertical direction of the PBC run for set C46 at a latitude of 30°. The black and red lines show and , respectively. Panel B: comparison of the timeaveraged normalized stresses obtained from periodic and stressfree boundary condition of the same run. The solid and dashed lines show the measured (red), (blue) and (black) from SFBC and PBC runs, respectively. 

In the text 
Fig. 9. Panels A, C, and E: timeaveraged diffusive and nondiffusive parts of the Reynolds stresses vs. vertical direction. The black and blue lines in panel A show the normalized vertical stresses at the equator and a latitude of 30° for set C24, respectively. In panel C the horizontal stresses are shown at a latitude of 30° for set C64. The blue and black lines in panel E show the meridional stresses for sets C06 and C64 at a latitude of 30°, respectively. The vertical lines denotes the z range used for volume averages. Solid, dotted, and dashdotted lines show , , and , respectively. Panels B, D, and F: Volume averages over −2 ≤ zk_{1} ≤ 2, of vs. Ω_{⋆} at different latitudes as indicated by the legend. 

In the text 
Fig. 10. Normalized turbulent viscosity and Λeffect coefficients as a function of Ω_{⋆} and latitude. Panels A, C, and D show , V, and H as a function of Ω_{⋆} from the equator to a latitude of 75°, respectively. Panel B shows as a function of latitude for five selected Ω_{⋆}. 

In the text 
Fig. 11. Comparison of the obtained turbulent viscosity with the analytical result of KPR94. The solid and dashed lines show the normalized turbulent viscosity and rescaled analytical expression κϕ_{1}, respectively. 

In the text 
Fig. A.1. vs. vertical direction at all Ω_{⋆} indicated by the legend. Each panel shows at a certain latitude. 

In the text 
Fig. A.2. vs. vertical direction at all Ω_{⋆} indicated by the legend. Each panel shows at a certain latitude. 

In the text 
Fig. A.3. vs. vertical direction at all Ω_{⋆} indicated by the legend. Each panel shows at a certain latitude. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.