{'fraction_non_alphanumeric': 0.03207916709617565, 'fraction_numerical': 0.018121843108957838, 'mean_word_length': 4.371055254124681, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generation of mean flows in rotating anisotropic turbulence: The case of solar near-surface shear layer

The radial gradient of the rotation rate in the near-surface shear layer (NSSL) of the Sun is independent of latitude and radius. Theoretical mean-field models have been successful in explaining this property of the solar NSSL, while global direct convection models have been unsuccessful. We investigate reason for this discrepancy by measuring the mean flows, Reynolds stress, and turbulent transport coefficients under NSSL conditions. Simulations have minimal ingredients. These ingredients are inhomogeneity due to boundaries, anisotropic turbulence, and rotation. Parameters of the simulations are chosen such they match 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 is varied by changing the rotation rate such that the resulting Coriolis numbers<1. We measure the turbulent transport coefficient relevant for the non-diffusive and diffusive parts of the Reynolds stress and compare them with predictions of current mean-field theories. A negative radial gradient of mean flow similar to the solar NSSL 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 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. Our local simulations do not validate the explanation for the generation of the NSSL from mean-field theory where meridional stresses are neglected. However, the rotational dependence of turbulent viscosity in our simulations is in good agreement with theoretical prediction. Our results are in qualitative agreement with global convection simulations in that a NSSL obtained near the equator.


Introduction
The convection zone (CZ) of the Sun, despite being highly turbulent, shows a well-organized large-scale 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 large-scale circulation in the meridional plane, known as the meridional flow (MC), is also present. The amplitude of MC is about 15-20 ms −1 which is two orders of magnitude smaller than the rotational velocity (Duvall 1979;Hathaway 1996).
The near-surface shear layer (NSSL) occupies about 17% of the CZ, or roughly 35Mm 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 d ln Ω d ln r ≈ −1 in the upper 13 Mm of the NSSL independent of latitude up to 60 • ( Barekat et al. 2014). Second, the gradient is evolving over time, by an amount between 5-10% of its time-averaged value, following closely the magnetic activity cycle ). 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. Such negative shear exists only in the NSSL in the solar CZ. The effect of the NSSL has been tested numerically in mean-field dynamo models by Käpylä et al. (2006) where it was found to aid equatorward migration. More observational and theoretical arguments for the NSSL strongly shaping the solar dynamo process were presented by Brandenburg (2005). The role of NSSL can be easily investigated in mean-field models where Article number, page 1 of 13 arXiv:2012.06343v1 [astro-ph.SR] 10 Dec 2020 it can be added or removed by hand. On the contrary, global 3D convection simulations typically fail in generating a realistic NSSL self-consistently (e.g. Guerrero et al. 2013;Hotta et al. 2015) and thus its role on the resulting dynamo solutions is unclear. Therefore, understanding the role that the NSSL plays for the dynamo requires that we first understand its formation mechanism and why global simulations do not capture it.
The equations governing the generation of large-scale flows in the solar CZ are the following: First, 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 A and fluctuations around the mean, a, and where averages are taken over the azimuthal direction. Then, we obtain the equation where ρ, U m = (U r , U θ , 0), Ω = U φ /r sin θ, ν, µ 0 , and B 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 S is the mean rate of strain tensor. The Reynolds and Maxwell stresses are the correlations of fluctuating components Q φj = u φ u j and M φj = b φ b j /µ 0 , respectively. Density fluctuations are omitted corresponding to an anelastic approximation.
Second, the azimuthally averaged equation for the azimuthal component of vorticity, describes the time evolution of MC: where w = ∇ × U 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 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 major players generating stellar differential rotation are the first two terms in both Eqs. (2) and (3) (Rüdiger 1989;Kitchatinov 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. It describes the role rotation plays in different parts of the CZ, in particular leading 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 10 −3 Ω NSSL 1 Ω CZ 10. Non-rotating density-stratified convection is dominated by vertical motions in which case the vertical anisotropy parameter A V ∝ u 2 H −u 2 r < 0, where u H and u r are the turbulent horizontal and radial velocities. Rotation tends to suppress convection (e.g. Chandrasekhar 1961) 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). On the other hand, rotation introduces horizontal anisotropy A H ∝ u 2 φ − u 2 θ , where u φ and u θ are the longitudinal and latitudinal velocities. Typically A H is positive and it is increasing with Ω . Furthermore, A H → 0 as Ω → 0. Thus, Ω in the solar CZ reflects also 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 non-zero off-diagonal 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 utilized a hydrodynamic mean-field (MF) model, considering the properties of the turbulent flow explained above and parameterized the Reynolds stresses in the form of turbulent transport coefficients (Rüdiger 1980(Rüdiger , 1989, see also Section 2). They obtain the NSSL by taking the anisotropy of turbulence near the surface into account such that A V A H for Ω 1. This leads to strong inward transport of the angular momentum in the NSSL and ultimately to the generation of the radial shear. The remarkable agreement of recent observed latitudinal independence of the gradient with their model brought motivation to develop the theory further including the effect of the magnetic field in the NSSL (Kitchatinov 2016). This leads to the prediction of the time variation of the angular velocity gradient during the solar cycle, qualitatively agreeing 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 for measuring the sub-surface magnetic field. In their model, however, the Reynolds stresses were computed using second-order 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 utilized since the 1970s (e.g. Gilman 1977Gilman , 1983Glatzmaier 1985), but reproducing the NSSL has turned out to be a serious challenge for these models. With such global convection simulations it is possible to generate a shear layer close to the equator, mostly confined outside the tangent cylinder, where rotation-aligned, elongated large-scale 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 when higher density stratification has been used, a shear layer extending to 60 • latitudes has been found Hotta et al. (2015). In this case, however, the gradient of the radial shear was positive in the range 0 • < θ < 45 • , contrary to the helioseismic inferences of the NSSL. They 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, hence the influence of rotation on convection (Coriolis number) was overestimated, and they also speculated about an unfavourable influence of boundary conditions to their results. Hence, it is unclear whether these results really are applicable to the NSSL. Overall, using global direct numerical simulations (GDNS) to study the origin of the NSSL is cumbersome due to high computational cost, the multitude of effects present, Article number, page 2 of 13 and the difficulty to reliably separate them from each other. For such an approach, a simpler modelling strategy is required which is attempted in this paper.
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 them, Gunderson & Bhattacharjee (2019), considers the formation of the NSSL in a magnetohydrostatic equilibrium model, being driven by a poleward meridional flow near the surface. In addition to the assumption of stationarity, while the magnetic field of the Sun is oscillatory, the model considers only non-turbulent states; nevertheless, a large-scale poloidal flow, when inserted on top of the equilibrium configuration, is seen to reduce the rotational velocity near the surface, hence leading to NSSL-like condition there. In the study of Miesch & Hindman (2011), the Reynolds and Maxwell stresses were accounted for in the governing equations, hence allowing for turbulent effects. They considered a case, where 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 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 it 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, such angular momentum transport is not a sufficient condition to sustain the NSSL, but another necessary ingredient is the meridional force balance in between the turbulent stresses and centrifugally-driven circulation within the NSSL. In the bulk of the convection zone, 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 NSSL, the force caused by the turbulent stresses was balanced by the Coriolis force.
In this paper we adopt an entirely different approach to those reviewed above. We formulate a model with minimal ingredients for the generation of large-scale flows to study the role of rotation-induced 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 unambiguous identification of the drivers of mean flows which can be used to assess the generation mechanisms of the solar NSSL.

The NSSL in terms of mean-field hydrodynamics
In this section we briefly explain the theory of the Λ-effect and its relevance for formation of the NSSL (Kitchatinov 2013(Kitchatinov , 2016. We refer the reader to Rüdiger (1989) for a thorough treatise. In this theory, rotating and anisotropic turbulence contributes to diffusive and non-diffusive transport of angular momentum. The non-diffusive part is known as the Λ-effect (Lebe-dinski 1941). Therefore, the Reynolds stress consists of two parts where N ijkl and Λ ijk are fourth-and third-rank 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 stresses appear only in the vorticity equation and in the model by KR05 they do not play a role in the generation of the NSSL. Ignoring magnetic fields, the vertical and horizontal stresses are given by Q θφ = ν cos θHΩ−sin θ ∂Ω ∂θ +ν ⊥ Ω 2 sin θ 2 cos θr ∂Ω ∂r , where ν and ν ⊥ are the diagonal and off-diagonal components of the turbulent viscosity tensor N ijkl , respectively. The latter component ν ⊥ appears due to 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 are, to the lowest order, 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 Q (Λ) θφ due to the Λeffect vanishes. The off-diagonal viscosity ν ⊥ is non-zero but small such that its influence is negligible (Rüdiger et al. 2019). It has been shown analytically (Kitchatinov & Rüdiger 2005) and numerically (Käpylä 2019b) that in the slow rotation regime only the first term in the expansion of the vertical coefficient V (0) survives and tends to a constant. Furthermore, applying a stressfree boundary condition at the radial boundaries, one realizes that Q rφ = Q rθ = 0. Using this in Eq. (7) and equating the diffusive and non-diffusive stresses we get which shows a reasonable agreement with observational results where the radial rotational gradient is independent of latitude (Barekat et al. 2014).

The model
We use a similar hydrodynamic model in Cartesian domain as in Käpylä & Brandenburg (2008) and Käpylä (2019b). We explain it here briefly and refer the reader to relevant parts of the above-mentioned works for details. In this model gravity is neglected and an external random forcing by non-helical transversal waves with direction-dependent amplitude is applied. The ensuing flow is turbulent and anisotropic. The medium is considered to be isothermal and to obey the ideal gas equation. The governing equations are where D/Dt = ∂/∂t + U · ∇ is the advective derivative, ρ and c s are the density and sound speed, respectively, and Ω = Ω 0 (− cos θ, 0, sin θ) T is the rotation vector. The viscous force is given by where S ij = 1 2 (U i,j + U j,i ) − 1 3 δ ij U k,k is the traceless rate of 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, wavevector, and a random phase, respectively. The desired vertical (z) anisotropy can be enforced using a tensorial normalization factor 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 makes the forcing transversal waves;ê is an arbitrary unit vector. The details of the forcing can be found in Brandenburg (2001). Notes. The grid resolution of all runs is 144 3 , forcing amplitudes f0 = 10 −6 and f1 = 0.04, Re ≈ 13, and ν = 3.3 · 10 −4 (csk 3 1 ) −1 .

Simulation setup
We used the PENCIL CODE 1 ) to run the simulations. We consider a cubic box with size (2π) 3 discretized over 144 3 grid points. z corresponds to vertical, x to latitudinal, and y to azimuthal direction, respectively, the latter two being referred to as the horizontal directions. Horizontal boundaries are periodic and stress-free 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 wavenumber k 1 = 2π/L and we choose a forcing wavenumber k f /k 1 = 10. The units of length, time, and density are k −1 1 , (c s k 1 ) −1 and ρ 0 , respectively, where ρ 0 is the initial uniform value of density. The forcing parameters f 0 = 10 −6 and f 1 = 0.04 are 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 fulfill the NSSL condition in which |A V | |A H |; see Table 1. The vigor of turbulence is quantified by the Reynolds number where u rms = (U 2 −U 2 ) 1/2 , is the root mean square of the fluctuating velocity field. Using a fixed value of the kinematic viscosity, ν = 3.3 · 10 −4 (c s k 3 1 ) −1 , the Reynolds number is about 13 for all simulations. We place 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 is determined by the value of Ω 0 which is varied such that the range of Ω from Eq. (4) is relevant for the NSSL. The turnover time is defined as τ = /u rms , where is the size of the eddies. In our simulations the energy-carrying scale of turbulence is the forcing scale = 2π/k f . Hence, the Coriolis number in the simulations is 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 performed to remove a contribution to the Reynolds stress appearing in the non-rotating case; see Section 5.4. Mean quantities are defined as horizontal (xy) averages. The local Cartesian quantities are related to their counterparts in spherical polar coordinates via (r, θ, φ) → (z, x, y), We normalize quantities such that U i = U i /u rms and Q ij = Q ij /u 2 rms , tilde denoting this operation. Additionally, the error on the measured physical quantities, which are obtained directly from the simulations, is estimated by dividing the time series into three parts and comparing their time averaged values with the one obtained from the whole time series. The maximum deviation from the latter is considered to be the error of the measurement.

Velocity field
A statistically stationary turbulent state appears after about few τ independent of Ω everywhere except at the equator, where the statistically stationary state is reached between few to about 300 τ from the lowest to highest Ω , respectively. As an example, we show snapshots of the zonal flow normalized by the sound speed at about 1000 τ for the set C46 at the equator and at 30 • latitude in panel A and B of Fig. 2, respectively. The other components of the velocity field are very similar to the zonal one shown in panel B. The dominant scale of the turbulence is the forcing scale k f /k 1 = 10. The expected large-scale zonal flow similar to the actual NSSL is generated only at the equator shown in panel A. All other sets show similar behaviour.

Anisotropy of the flow
We start our analysis by measuring the diagonal components of the Reynolds stresses and the anisotropy parameters which are given by Article number, page 5 of 13 We show the volume averaged A V and A H as a function of Ω at different latitudes in the middle and bottom panels of Fig. 3, respectively. We consider data between −2 ≤ zk 1 ≤ 2 for the volume averages to avoid boundary effects. The vertical anisotropy parameter A V is always at least two orders of magnitude greater than A H . Neither shows appreciable variation as a function of Coriolis number for Ω 0.15. The vertical anisotropy parameter is almost independent of Ω at the equator in contrast to other latitudes where its absolute value decreases with increasing Ω . It decreases about 15% at the bottom of the NSSL at 15 • and about 5% at latitudes above 45 • . The horizontal anisotropy parameter shows almost no dependence on latitude above 45 • but it becomes 100 times greater from the top to the bottom of the NSSL below this latitude. The behaviour of both anisotropy parameters is similar to the ones obtained by Käpylä & Brandenburg (2008) in which they have similar set-up as ours but applied fully periodic boundary conditions. This shows that anisotropy of the flow is insensitive to the boundary conditions.

Mean flows
The development of mean flows in rotating cases means that reaching a statistically steady state takes significantly longer than in non-rotating runs. Furthermore, long time averages are needed for statistical convergence of the turbulent quantities. We run 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 the set C24 at the equator and at θ = 15 • in Fig. 4. At the equator, a large zonal flow U y with a negative vertical gradient developed gradually over 100τ as shown in panel (A). All other sets show similar zonal flow profile at the equator, but both the amplitude and steepness of the gradient increase with increasing Ω . Moving away from the equator, the amplitude of the mean zonal flow reduces significantly and the negative gradient disappears as shown in panel (B) of Fig. 4. The dependence of mean zonal flow on rotation can be seen in the panel (A) of Fig. 5 where we show the time-averaged U y at selected Ω at 15 • latitude. By increasing Ω , the gradient of U y changes sign and becomes steeper up to Ω = 0.46, then it becomes shallower and slowly vanishes in the middle at Ω = 1. The latitudinal dependence of U y is shown for sets C06 and C46 in the panels (C) and (E) of Fig. 5, respectively. We find that U y decreases as a function of latitude, vanishes at the poles, and that Article number, page 6 of 13 For completeness, we show the vertical component of the normalized mean flow U z in the bottom row of Fig. 4. All runs show a similar pattern of high-frequency oscillations for U z irrespective of latitude and Ω with amplitudes of the order of 10 −4 u rms . These oscillations are identified as longitudinal sound waves as expected for a compressible system in a confined cavity.

Reynolds stresses
For zero rotation, it is expected that Q xy = 0, see Eq. (8). However, we find that Q xy always has a small but non-zero value which persists also in the longest time series of our data. We find that such contribution is present and its magnitude remains unchanged also for higher resolutions (288 3 and 576 3 grids). Hence, this issue seems not to be caused by a numerical convergence issue, but we have been unable to identify whether the cause is due to some, yet unidentified, physical effect for example due to compressibility, effects due to the forcing, inhomogeneities in the system, or a combination thereof. Given that this contribution is systematically present, we made a non-rotating run (C0) and subtracted Q xy from that run from the results of the runs with rotation.
We show representative results of the off-diagonal stresses at five selected Ω at 15 • (left column) and as function of latitude (right column) in Fig. 6. The data for all sets is available as online electronic material. 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. Q yz 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 decreases gradually towards 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 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. In Pulkkinen et al. (1993) it is almost constant up to 60 • and decreases toward higher lati-tudes. Q yz has a v-shape profile in latitude with the minimum at 45 • in Käpylä et al. (2004). The major ingredient which is missing in our forced turbulence simulation in comparison with theirs is density stratification. Moreover, Käpylä et al. (2004) includes the overshooting layer below the CZ. Therefore, it is difficult to find out what makes our results different from theirs.
The middle panels (C) and (D) in Fig. 6 show horizontal stress Q xy . The signature of turbulent fluctuations at the forcing scale are seen more clearly in this component and the measurement is quite noisy. The values of Q xy are close to zero up to Ω = 0.46, above which it slowly starts to get positive (negative) values in the middle (close to the boundaries). This is the same behavior as with A H seen in Fig. 3. At a given Ω , the profile of Q xy is similar at all latitudes. Its amplitude is maximum at 30 • and decreases gradually towards the equator and the pole as shown in panel (D). This result is in agreement 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 Figure 10 in Hanasoge et al. (2016). The horizontal stress has always positive values independent of Ω and latitude in agreement with previous studies in slow rotation 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, albeit with negative values as their box is located at the southern hemisphere; see their Figure 6.
The meridional stress is shown in the last row of Fig. 6. In contrast to the other stresses, Q xz shows a complicated profile, in particular close to the boundaries. Moreover, it has positive or negative values depending on both Ω and θ shown in panel (E). The latitudinal dependency of the meridional stress is shown in panel (F). At Ω < 0.1, Q xz is positive at low latitudes and Q xz → 0 above 45 • . By increasing Ω , Q xz moves toward negative values and its absolute value increases. For Ω > 0.24, it has its maximum amplitude at about 45 • and it decreases toward the pole and the equator similar to Q xy . The meridional stress in Pulkkinen et al. (1993) also shows a sign change in agreement with ours albeit in mid latitude. However, the sign change occurs at Ω ≈ 1 while our shows only negative values at that Ω .
Comparing the absolute amplitude of the stresses in right column of Fig. 6 we see that Q yz always larger than Q xz and Q xy . For example at Ω = 0.64, Q yz is about two to ten times larger than Q xz and five to twenty times larger than Q yz depending on latitudes. Comparing also the absolute amplitude of Q xy and Q xz , we see that Q xz > Q xy for all Ω . These results show that, in spite of the fact that Q xy is increasing as a function of Ω , its values are still much smaller than vertical stresses which is in agreement with the assumption of KR05 regarding the NSSL.
Although our model is quite simple in comparison to GDNS, it is of interest to compare the Reynolds stresses with simulations such as those in Käpylä et al. (2011b). These authors modelled 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 maximal 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 Ω on the northern hemisphere with a sign change after certain Ω . However, this trend depends on latitude in their case; see their Figure 8. The profile of Q yz in the convection simulations is quite different from ours such that it has a strong latitudinal dependency and gets both positive and negative values depending on Ω and latitude. This is consistent with earlier studies (e.g. Käpylä 2019b) where a sign change of Q yz occurs at higher Ω than those considered in the present simulations.

The role of Reynolds stresses in the generation of the mean flows
As the Reynolds stresses appear in the MF momentum equation, we start by writing the MF equations for U x and U y using Eq. (13). We wrote these equations first by considering the facts that our setup is fully compressible and the forcing we used here is not solenoidal which might cause density fluctuations which cannot be ignored in the MF equations. These considerations lead to the presence of three additional terms to the Reynolds stresses in the momentum equation (e.g. Käpylä et al. 2020). We compared all of them with the Reynolds stresses, and it turns out that 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 U z , and Ω y = 0, yields the final form of the equations: Article number, page 8 of 13 We double-checked the validity of the MF equations by considering the steady-state solution which reads We show the horizontal mean velocities in comparison with the RHS of Eqs. (27) and (28) 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 emphasise, 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. Such behavior can also be clearly seen in the first panel of Fig. 4, where we show the time evolution of U y .
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 has to be included, but it is clear that the Reynolds stresses are the main driver of mean flows in the current setups.

Parameterization of Reynolds stresses in terms of mean-field hydrodynamics
Based on the Λ-effect theory explained in Section 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 also the turbulent viscosities ν and ν ⊥ are unknown. Our strategy around this is to run another set of otherwise identical simulations, but where the horizontal mean flows are artificially suppressed at each time step. Therefore, the first terms in Eqs. (31) and (32) go to zero. Then from these simulations we can directly measure Q (Λ) . However, we need to validate this approach because the velocities can be affected through the non-linearity of the Navier-Stokes equations. Therefore, we perform yet another set of otherwise identical simulations, but use periodic boundary conditions (PBC) in all directions instead of stress-free boundary condition (SFBC) in the vertical direction. Then we compare the two sets of stresses obtained with these sets of boundary conditions. Such a comparison of varying boundary conditions is important also in the respect of interpreting the Ω dependence as depth dependencethis approach is somewhat artificial, as we practically enforce unrealistic BCs within the convection zone.
As an example, we show the horizontal mean velocities for the PBC version of C46 at 30 • latitude 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 where the mean flows are suppressed. In panel (B) of Fig. 8, we show the results of the comparison of the Reynolds stresses between 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 depth is valid.
Considering Eq. (31), the subtraction of the Reynolds stresses obtained from these simulations from the total ones gives Measuring the vertical gradient of U y , the value of ν can be determined by performing an error-weighted linear least-squares fit to Eq. (33). Putting the measured values of ν back into Q (Λ) of both Eqs. (31) and (32), we can measure V and H provided that ν ⊥ ν .

Properties of the diffusive and non-diffusive parts of Reynolds stresses
Similar to Section 5.4, we first measure Q ij from a non-rotating run and then subtract its mean value from 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 of Q ij on both latitude and Ω . In panel (A), we show the vertical stresses for set C24 at the equator and at 30 • latitude. With these results we can explain the minimum of Q yz at the equator: the diffusive and non-diffusive parts of the stresses are comparable but of opposite signs, leading to a small negative value for the total. With Eq. (31), we see that ν > 0 which in combination with ∂ z U y < 0, gives Q (ν) yz > 0. Moreover, the final negative value of Q yz also shows that Q (Λ) yz is responsible for the generation of the zonal flow. The profile of Q (ν) yz for all other latitudes is similar to the one at 30 • , and shows that the major contribution from the diffusive part is happening close to the boundaries at |zk 1 | 2 with positive values. Furthermore, the amplitude of Q (ν) yz decreases towards higher latitudes (not shown). In the middle of the domain it has negative values, which fits well with ∂ z U y > 0 that can be seen in Fig. 5. The non-diffusive part of the vertical stress is always Q (Λ) yz < 0. Its absolute value increases from the pole towards the equator and increases with Ω . We also find that Q (Λ) yz is linearly dependent on Ω in slow rotation regime Ω 1 in agreement with previous numerical results (Käpylä 2019a).
We show corresponding results for Q xy in panel (C) of Fig. 9 for Ω = 0.64. Q This also shows that Q (ν) xy is the main contributor to the negative values of Q xy close to the boundaries shown in Fig. 6 (C). The volume averaged values of Q (Λ) xy , excluding the boundaries, are shown in Fig. 9 (D) as a function of both Ω and latitude. Its value is almost zero both at the equator and at the pole. It is significantly non-zero above Ω > 0.24 and increases with increasing Ω independent of latitude. Independent of Ω , it has maximum value at 30 • latitude. We note here that the amplitude of Q (Λ) xy is also significantly smaller than Q (Λ) yz . The measured profile of Q (Λ) xy is almost identical to the one obtained by Käpylä (2019b).
Our results for Q xz are shown in Fig. 9. At low Ω , there is almost no contribution of Q is almost zero at the equator and at the pole. In other latitudes, its absolute value increases with increasing Ω . It has always negative values independent of both Ω and latitude. These results are in agreement with those of Käpylä (2019b) in the slow rotation regime.

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. Apart from 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, we think that the latitudinal dependence of ν for θ 60 • , shown in panel (B), is not reliable. According to the results at lower latitudes, we conclude that the latitude dependence is weak in comparison to the rotational dependence. Hence, we consider the profile of ν at the equator applicable for other latitudes which is measured with high confidence and use it for measuring 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. Käpylä et al. (2020) measured turbulent viscosity by imposing a weak sinusoidal shear flow in a non-rotating isotropically forced turbulent medium in a Cartesian geometry and measured the response of the system. The response is an off-diagonal Reynolds stress that is assumed to be proportional to the imposed shear flow according to the Boussinesq ansatz. They defined a shear number Sh as where U 0 is the amplitude of the flow and k U is the wavenumber of the imposed sinusoidal shear. To obtain Sh for the shear flows generated at the equator in the present simulations, we set k U = k 1 /2 and U 0 /u rms = max( U y ). We consider only the slow rotation regime where Ω < 0.1, corresponding to sets C02, C04 and C06, where Sh = 0.002, 0.004 and 0.006 which is within the range of Sh values used in Käpylä et al. (2020). We normalized our ν with the same normalization factor as in Käpylä et al. (2020), that is ν t0 = u rms /3k f ; see their Section 2.2. This differs from the currently used normalization by a factor of 6π such that the current normalized values, for example in Fig. 10, are smaller than theirs by this factor. Using their normalization, we obtain values of ν t /ν t0 ≈ 3.5 . . . 3.8 which are roughly twice larger than in Käpylä et al. (2020). The difference is likely caused by the presence of strong anisotropy of turbulence in our simulations due to the forcing and the rotation that were absent in the study of Käpylä et al. (2020). We also compare 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 consider the first term in Equation (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 non-rotating case, and where φ 1 is a function of Ω given in the Appendix of KPR94. We scale 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 apart from κ factor the rotation dependence is in fair agreement between the theory and numerical simulations.
Considering the off-diagonal turbulent viscosity ν ⊥ , we failed measuring it as both terms constituting it, Q (ν) xy and Ω 2 sin θ cos θ∂U y /∂z, are too small, and the measurement error in the former is large.

Measurements of the vertical Λ-effect coefficient
We measure the vertical Λ-coefficient by substituting the volume averages of Q Our results are shown in panel (C) of Fig. 10. The absolute value of V is about 0.75 and gradually increases to ≈ 0.95 for latitudes ≤ 45 • . However, the value of V at the lowest Ω are smaller at all latitudes, but they have large error bars. In contrast to low latitudes, the absolute values of V at 60 • and 75 • decrease for Ω > 0.3. Considering the large errors in the measurements at low Ω , we might consider V being roughly constant for Ω ≤ 0.15 independent of latitude, but it shows strong latitudinal and rotational dependency for Ω > 0.15. This means that considering only the first term V (0) in Eq. (9) in the NSSL condition is not enough as it is assumed by the theoretical model by KR05 explained in Section 2. Moreover, the increase of V towards higher Ω at low latitudes is in contrast with the decrease predicted by KR05 model. The same applies to the results of Käpylä (2019b) who did not consider that ν t = ν t (Ω ).

Measurements of the horizontal Λ-effect coefficient
We measure the horizontal Λ-effect coefficient similarly to the vertical one 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 smaller 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 the largest at 15 • and decreases gradually towards higher latitudes. These results show that H does not play any role close to the surface in transporting the angular momentum which validates the assumption applied in the NSSL model by KR05.

Conclusions
We applied an alternative approach to MF and GDNS, namely running direct numerical simulations of forced turbulence in lo- cal boxes, to primarily find out if the assumptions and approximations applied in MF theory to explain the formation of the NSSL are valid. In contrast to GDNS, we could 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 applying the three required conditions, explained in Section 2, that 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 non-zero values close to the bottom of the NSSL in agreement with the theoretical predictions. Regarding the vertical Reynolds stress, its role in transporting the angular momentum radially inward is in agreement with theory. However, its profile differs from that predicted by theory. In particular, it was assumed in Kitchatinov (2013) and Kitchatinov (2016) that only the term V (0) survives in the expansion of V in the NSSL. However, our results indicate that higher order terms in the expansion of V need to be considered. Moreover, it is also expected from theory that the vertical transport of angular momentum decreases by increasing Ω independent of latitude, 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 which was not considered previously in the models of NSSL. This behavior, however, from a theoretical MF prediction is in good qualitative agreement with our results (Kitchatinov et al. 1994). 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 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. Hence, it is not surprising that we find a relatively good match in between the measured turbulent viscosity and the one predicted by SOCA.
Concerning the fact that Q xz cannot be disregarded in the NSSL, its role can be further investigated in more realistic setup using spherical geometry where the artifact of discontinuity between latitudes can be removed. We also note here that in this work we consider 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 need to be investigated are the effect of stratification, compressibility and magnetic fields, and comparing with previous studies that have studied these 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 water-filled apparatus will be used to simulate regions of finite latitudinal extent, including β-plane effects, and forcing will be introduced by pump-driven nozzles at the boundaries (Burin et al. 2019). 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/closing nozzles and altering the nozzle shapes and orientations. Time resolved 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. Despite that the details of the forcing and the fluid boundary conditions will be different in the experiment than in the present simulations, it is expected that meaningful results will be obtained as the rotation rate of the system is varied and the experimental data is analyzed to look for signatures of the Λ-effect.