EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 655 (November 2021) A&A, 655 (2021) A79 Full HTML
Open Access
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/0004-6361/202040052
Published online 23 November 2021

© A. Barekat et al. 2021

Licence Creative CommonsOpen 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 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 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 near-surface 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

(1)

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 time-averaged 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 mean-field 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 mean-field models, where it can be added or removed by hand. In contrast, global 3D convection simulations typically fail to self-consistently 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 large-scale 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

(2)

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:

(3)

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

(4)

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 near-surface 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 density-stratified 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 |AV| decreases when Ω⋆ increases such that the maximum of |AV| is achieved for Ω⋆ = 0 (e.g., Chan 2001; Käpylä et al. 2004). |AV| 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), AH 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, AH → 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 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 derived an extension to the quasi-linear turbulence model of Kichatinov & Rüdiger (1993) that did not take the slow-rotation regime into account that is required to treat the NSSL. In the new model (KR05), the vertical anisotropy parameter AV 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 |AV|≫AH for Ω⋆ ≲ 1 applies in the surface layers, although the turbulent transport coefficients depend on latitude. They used a hydrodynamic mean-field (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 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 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 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 for a higher-density 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 large-scale 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 NSSL-like 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 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 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 mean-field 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,

(5)

(6)

where Nijkl and Λijk are fourth- and third-rank tensors describing the turbulent viscosity and Λ-effect, respectively. In spherical geometry, Qrϕ, Qθϕ, and Qrθ 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

(7)

(8)

where ν∥ and ν⊥ are the diagonal and off-diagonal components of the turbulent viscosity tensor Nijkl, 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 AV and AH (Rüdiger 1980). These coefficients are typically expanded in latitude in powers of sin2θ as

(9)

(10)

In the NSSL, Ω⋆ ≤ 1 and AH ≈ 0, such that due to the Λ-effect vanishes. The off-diagonal 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 slow-rotation regime, only the first term in the expansion of the vertical coefficient V(0) survives and tends to a constant. Furthermore, when a stress-free boundary condition is applied at the radial boundaries, Qrϕ = Qrθ = 0. Using this in Eq. (7) and equating the diffusive and nondiffusive stresses, we obtain

(11)

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 |AV|≫AH, Ω⋆ ≤ 1, and a stress-free 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 stress-free 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.

thumbnail 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 self-consistently 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 cs 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 = urms/cs = 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 fast-rotating regimes by Käpylä & Brandenburg (2008) and Käpylä (2019b).

The hydrodynamic equations that we solved directly are

(12)

(13)

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

(14)

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

(15)

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 f0 and f1 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, fk is given by

(16)

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 CODE1 (Brandenburg et al. 2021) to run the simulations. We considered a cubic box with size (2π)3 discretized over 1443 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 stress-free conditions are imposed at vertical boundaries with

(17)

where zbot and ztop represent the bottom and top of the domain. The box size is represented by the wave number k1 = 2π/L, and we chose a forcing wave number kf/k1 = 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 f0 = 10−6 and f1 = 0.04 were chosen such that the effects of compressibility are weak with a Mach number Ma = urms/cs ≈ 0.04 in all simulations. Moreover, with f1 ≫ f0, we fulfilled the NSSL condition in which |AV|≫AH; see Table 1. The vigor of turbulence is quantified by the Reynolds number,

(18)

Table 1.

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 τ = ℓ/urms, where ℓ is the size of the eddies. In our simulations the energy-carrying scale of turbulence is the forcing scale ℓ = 2π/kf. The Coriolis number in the simulations is therefore given by

(19)

The corresponding input parameter is the Taylor number,

(20)

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θϕ → Qxy, Qθr → Qxz, and Qrϕ → Qyz. 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 time-averaged 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 |AV|≫AH in Sect. 5.2. In Sect. 5.3, by measuring the mean velocity field at all depths and latitudes, we search for the NSSL-like flow. Next, we focus on the properties of the off-diagonal 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 kf/k1 = 10. The expected large-scale 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.

thumbnail 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 Uy/cs 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

(21)

(22)

We show the time-averaged 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 AV ≈ −0.5. This is shown in panel B of Fig. 3, where we show volume averages of AV excluding the boundaries (−2 ≤ zk1 ≤ 2) to avoid boundary effects. As expected, |AV| has its maximum value in the absence of rotation at Ω⋆ = 0.

thumbnail Fig. 3.

Panel A: time-averaged 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 AV and AH were measured. Anisotropy parameters AVpanel B and AHpanel C are shown as functions of Ω⋆ at the latitudes indicated in the legend. The diamonds in panel B show the values of AV at the equator from the runs from which the mean flow was removed.

In the slow-rotation 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. AH 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 Ω⋆. AH increases significantly by about one order of magnitude from the top to the bottom of the NSSL at low latitudes. In contrast, |AV| decreases by only about 15% at low latitudes and stays constant at the equator. This shows that the latitudinal dependence of |AV| is much weaker than that of AH in Ω⋆ < 1.

Although the latitudinal dependence of AV is weak, it shows a nonmonotonous profile. For a specific Ω⋆, the largest magnitude of AV is obtained at the equator. |AV| 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 AV 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 stress-free 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 AV and AH that the model fulfills the NSSL condition |AV|≫AH. 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 time-averaged 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 urms everywhere except at the boundaries.

thumbnail 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 time-averaged 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 time-averaged 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 NSSL-like 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.

thumbnail Fig. 5.

Time-averaged 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 high-frequency oscillations for regardless of latitude and Ω⋆ with amplitudes of about 10−4urms. 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 (2883 and 5763 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 off-diagonal stresses at five selected Ω⋆ at 15° (left column) and spatially averaged (−0.5 ≤ zk1 ≤ 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 Qyz is different from that in previous studies by Pulkkinen et al. (1993) and Käpylä et al. (2004) at Ω⋆ ≈ 1, in which they measured Qyz 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. Qyz has a v-shape 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.

thumbnail Fig. 6.

Left column: time-averaged off-diagonal 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 ≤ zk1 ≤ 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 Qxy. 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 Qxy 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 Qxy obtains notable values at similar Ω⋆ as AH. This clearly shows that the generation of horizontal Reynolds stress depends on the presence of the horizontal anisotropy in the flow. However, AH being maximum at certain latitude does not necessarily mean that Qxy 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 AH shown in panel C of Fig. 3. At a given Ω⋆, the profile of Qxy 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 Qxy 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 slow-rotation regime (e.g. Kitchatinov & Rüdiger 2005; Käpylä & Brandenburg 2008; Käpylä 2019b). The latitudinal profile of Qxy 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 mid-latitude. 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 Qxy 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 Qxy, 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 Qxz 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 Qyz 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 Qyz 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

(23)

(24)

Omitting terms proportional to the small quantities ν and , and Ωy = 0, yields the final form of the equations,

(25)

(26)

We verified the validity of the MF equations by considering the steady-state solution, which reads

(27)

(28)

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.

thumbnail Fig. 7.

Time-averaged 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 Qxz 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

(29)

The pressure gradient appears in this equation due to horizontal averaging. In the steady state, the zonal flow can be written as

(30)

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 Qyz 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 mean-field 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

(31)

(32)

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 Navier-Stokes 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 stress-free 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.

thumbnail Fig. 8.

Panel A: time-averaged 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 time-averaged normalized stresses obtained from periodic and stress-free 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

(33)

Measuring the vertical gradient of , the value of ν∥ can be determined by performing an error-weighted linear least-squares 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 Qij 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 ≤ zk1 ≤ 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 Qyz 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 Qyz 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 |zk1|≳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 slow-rotation regime, in agreement with previous numerical results (Käpylä 2019a).

thumbnail Fig. 9.

Panels A, C, and E: time-averaged 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 dash-dotted lines show , , and , respectively. Panels B, D, and F: Volume averages over −2 ≤ zk1 ≤ 2, of vs. Ω⋆ at different latitudes as indicated by the legend.

We show corresponding results for Qxy 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 Qxy is confined to the boundaries at |zk1|≳2. This also shows that is the main contributor to the negative values of Qxy close to the boundaries shown in Fig. 6 C. The volume-averaged 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 Qxz 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 Qxz 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 slow-rotation regime.

6.2. Measuring turbulent viscosity

The diagonal turbulent viscosity ν∥, normalized by ℓurms, 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.

thumbnail 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 = urms/3kf, 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ℓurms 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.

thumbnail 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 off-diagonal 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

(34)

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

(35)

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 Qxy → 0 in the upper part of the NSSL, whereas Qxy 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 higher-order 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.

Qxz 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 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). 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. 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. 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 Max-Planck/Princeton Center for Plasma Physics. P. J. Käpylä acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG) Heisenberg programme (grant No. KA 4825/2-1).

References

  1. Barekat, A., Schou, J., & Gizon, L. 2014, A&A, 570, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Barekat, A., Schou, J., & Gizon, L. 2016, A&A, 595, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Brandenburg, A. 2001, ApJ, 550, 824 [Google Scholar]
  4. Brandenburg, A. 2005, ApJ, 625, 539 [NASA ADS] [CrossRef] [Google Scholar]
  5. Brandenburg, A., Johansen, A., Bourdin, P. A., et al. 2021, J. Open Source Softw., 6, 2807 [CrossRef] [Google Scholar]
  6. Burin, M. J., Caspary, K. J., Edlund, E. M., et al. 2019, Phys. Rev. E, 99, 023108 [NASA ADS] [CrossRef] [Google Scholar]
  7. Chan, K. L. 2001, ApJ, 548, 1102 [NASA ADS] [CrossRef] [Google Scholar]
  8. Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: Clarendon) [Google Scholar]
  9. Duvall Jr., T. L. 1979, Sol. Phys., 63, 3 [NASA ADS] [CrossRef] [Google Scholar]
  10. Foukal, P., & Jokipii, J. R. 1975, ApJ, 199, L71 [NASA ADS] [CrossRef] [Google Scholar]
  11. Gilman, P. A. 1977, Geophys. Astrophys. Fluid Dynam., 8, 93 [CrossRef] [Google Scholar]
  12. Gilman, P. A. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
  13. Gilman, P. A., & Howard, R. 1984, Sol. Phys., 93, 171 [NASA ADS] [Google Scholar]
  14. Glatzmaier, G. A. 1985, ApJ, 291, 300 [NASA ADS] [CrossRef] [Google Scholar]
  15. Greer, B. J., Hindman, B. W., Featherstone, N. A., & Toomre, J. 2015, ApJ, 803, L17 [Google Scholar]
  16. 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]
  17. Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
  18. Gunderson, L. M., & Bhattacharjee, A. 2019, ApJ, 870, 47 [CrossRef] [Google Scholar]
  19. Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Ann. Rev. Fluid Mech., 48, 191 [Google Scholar]
  20. Hathaway, D. H. 1996, ApJ, 460, 1027 [NASA ADS] [CrossRef] [Google Scholar]
  21. Hathaway, D. H., & Upton, L. 2014, J. Geophys. Res., 119, 3316 [NASA ADS] [CrossRef] [Google Scholar]
  22. Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 798, 51 [Google Scholar]
  23. Käpylä, P. J. 2019a, Astron. Nachr., 340, 744 [CrossRef] [Google Scholar]
  24. Käpylä, P. J. 2019b, A&A, 622, A195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Käpylä, P. J., & Brandenburg, A. 2008, A&A, 488, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2006, Astron. Nachr., 327, 884 [CrossRef] [Google Scholar]
  27. Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2004, A&A, 422, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011a, Astron. Nachr., 332, 883 [CrossRef] [Google Scholar]
  29. 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]
  30. Käpylä, P. J., Rheinhardt, M., Brandenburg, A., & Käpylä, M. J. 2020, A&A, 636, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Kichatinov, L. L., & Rüdiger, G. 1993, A&A, 276, 96 [NASA ADS] [Google Scholar]
  32. Kitchatinov, L. L. 2013, in IAU Symposium, eds. A. G. Kosovichev, E. de Gouveia Dal Pino, & Y. Yan, 294, 399 [Google Scholar]
  33. Kitchatinov, L. L. 2016, Astron. Lett., 42, 339 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kitchatinov, L. L., Pipin, V. V., & Ruediger, G. 1994, Astron. Nachr., 315, 157 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kitchatinov, L. L., & Rüdiger, G. 2005, Astron. Nachr., 326, 379 [NASA ADS] [CrossRef] [Google Scholar]
  36. Krause, F., & Rädler, K. H. 1980, Mean-field magnetohydro dynamics and dynamo theory (Oxford: Pergamon Press) [Google Scholar]
  37. Lebedinski, A. I. 1941, Astron. Zh., 18, 10 [Google Scholar]
  38. Matilsky, L. I., Hindman, B. W., & Toomre, J. 2019, ApJ, 871, 217 [NASA ADS] [CrossRef] [Google Scholar]
  39. Miesch, M. S., & Hindman, B. W. 2011, ApJ, 743, 79 [NASA ADS] [CrossRef] [Google Scholar]
  40. Parker, E. N. 1955, ApJ, 122, 293 [Google Scholar]
  41. Pulkkinen, P., & Tuominen, I. 1998, A&A, 332, 755 [NASA ADS] [Google Scholar]
  42. Pulkkinen, P., Tuominen, I., Brandenburg, A., Nordlund, A., & Stein, R. F. 1993, A&A, 267, 265 [NASA ADS] [Google Scholar]
  43. Robinson, F. J., & Chan, K. L. 2001, MNRAS, 321, 723 [NASA ADS] [CrossRef] [Google Scholar]
  44. Rüdiger, G. 1980, Geophys. Astrophys. Fluid Dynam., 16, 239 [CrossRef] [Google Scholar]
  45. Rüdiger, G. 1989, Differential rotation and stellar convection. Sun and the solar stars (Berlin: Akademie Verlag) [CrossRef] [Google Scholar]
  46. Rüdiger, G., Küker, M., Käpylä, P. J., & Strassmeier, K. G. 2019, A&A, 630, A109 [EDP Sciences] [Google Scholar]
  47. Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [Google Scholar]
  48. Stix, M. 2002, The sun: an introduction (Berlin: Springer) [Google Scholar]
  49. Thompson, M. J., Toomre, J., Anderson, E. R., et al. 1996, Science, 272, 1300 [Google Scholar]
  50. Ward, F. 1965, ApJ, 141, 534 [NASA ADS] [CrossRef] [Google Scholar]
  51. Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2016, A&A, 596, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. 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.

thumbnail Fig. A.1.

vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows at a certain latitude.

thumbnail Fig. A.2.

vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows at a certain latitude.

thumbnail Fig. A.3.

vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows at a certain latitude.

All Tables

Table 1.

Summary of the runs in which we varied the Taylor number and latitude.

In the text

All Figures

thumbnail 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
thumbnail 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 Uy/cs at the equator and at θ = 30° for set C46, respectively.
In the text
thumbnail Fig. 3.

Panel A: time-averaged 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 AV and AH were measured. Anisotropy parameters AVpanel B and AHpanel C are shown as functions of Ω⋆ at the latitudes indicated in the legend. The diamonds in panel B show the values of AV at the equator from the runs from which the mean flow was removed.
In the text
thumbnail 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
thumbnail Fig. 5.

Time-averaged 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
thumbnail Fig. 6.

Left column: time-averaged off-diagonal 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 ≤ zk1 ≤ 0.5), at different latitudes. The rows from top to bottom show , , and .
In the text
thumbnail Fig. 7.

Time-averaged 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
thumbnail Fig. 8.

Panel A: time-averaged 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 time-averaged normalized stresses obtained from periodic and stress-free 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
thumbnail Fig. 9.

Panels A, C, and E: time-averaged 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 dash-dotted lines show , , and , respectively. Panels B, D, and F: Volume averages over −2 ≤ zk1 ≤ 2, of vs. Ω⋆ at different latitudes as indicated by the legend.
In the text
thumbnail 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
thumbnail 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
thumbnail Fig. A.1.

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

vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows at a certain latitude.
In the text
thumbnail 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.