© The Authors 2026

1. Introduction

Observations of the Sun and other stars with convective envelopes show that differential rotation is ubiquitous (e.g. Rüdiger et al. 2013; Reinhold & Gizon 2015). The main driver of differential rotation has long been thought to be non-vanishing Reynolds stress due to the interaction of rotation and convective turbulence (e.g. Rüdiger 1989):

$$\begin{array}{c}\hfill Q_{\mathit{ij}}=\overline{u_i{u}_j},\end{array}$$(1)

where $u_i=U_i-{\overline{U}}_i$ is the fluctuating velocity and where the overbar denotes a suitably defined average. More specifically, the off-diagonal components of Q_ij have non-diffusive contributions in anisotropic turbulence that are proportional to the angular velocity, Ω, which lead to the generation of differential rotation. This is more commonly known as the Λ effect (e.g. Rüdiger 1980).

The Λ effect has been studied with simulations of anisotropically forced turbulence (Käpylä & Brandenburg 2008; Käpylä 2019b,a; Barekat et al. 2021). These results have confirmed the existence of the Λ effect and are in broad agreement with mean-field hydrodynamics (Rüdiger 1989). Convection simulations show qualitatively similar results at slow rotation (e.g. Pulkkinen et al. 1993; Rüdiger et al. 2005, 2019; Hupfer et al. 2005; Rüdiger & Küker 2021), but differences arise at rapid rotation where, for example, the radial angular momentum flux changes sign (Chan 2001; Käpylä et al. 2004). The latter arises due to large-scale wave-like convective modes also known as Busse columns, banana cells, or thermal Rossby waves that appear as tilted columns in an equatorial latitude belt (e.g. Busse 2002; Aurnou et al. 2007); see Käpylä (2023) and Mori & Hotta (2023) for recent analyses of the associated angular momentum transport. Mean-field theories of stellar angular momentum transport (e.g. Kitchatinov & Rüdiger 1995, 2005; Kleeorin & Rogachevskii 2006) do not include such large-scale convective modes, and therefore the accuracy of their predictions in the case of convection is uncertain. The main goal of the current study is to compute the Λ effect from convection using a comprehensive set of simulations at varying rotation rates for the first time.

The paper is organised as follows. Section 2 introduces the Λ effect in the framework of mean-field hydrodynamics. Section 3 describes the model. In Sect. 4 the results of the study are discussed, and in Sect. 5 the conclusions are presented.

2. Mean-field theory and the Λ effect

The Reynolds stress in rotating turbulence can be written as

$$\begin{array}{c}\hfill Q_{\mathit{ij}}=Q_{\mathit{ij}}^{(0)}+Q_{\mathit{ij}}^{(\mathrm{\Omega })}+Q_{\mathit{ij}}^{(S)},\end{array}$$(2)

where $Q_{\mathit{ij}}^{(0)}$ is present in the absence of rotation or shear, $Q_{\mathit{ij}}^{(\mathrm{\Omega })}$ is due to rotation, and $Q_{\mathit{ij}}^{(S)}$ is due to shear. Assuming the large-scale fields vary slowly in space and time, we write an ansatz,

$$\begin{array}{c}\hfill Q_{\mathit{ij}}=Q_{\mathit{ij}}^{(0)}+{\mathrm{\Lambda }}_{\mathit{ijk}}{\mathrm{\Omega }}_k-{\mathcal{N}}_{\mathit{ijkl}}\frac{\partial {\overline{U}}_k}{\partial x_l},\end{array}$$(3)

where Λ_ijk and 𝒩_ijkl are third and fourth rank tensors, respectively, and Ω is the mean angular velocity. The appearance of cross-correlations requires two preferred directions (e.g. Rüdiger 1989). The gravity (g) and angular velocity (Ω) vectors fulfill this requirement in rotating convection. However, in such cases large-scale flows are produced due to the Λ effect, and disambiguation of $Q_{\mathit{ij}}^{(\mathrm{\Omega })}$ and $Q_{\mathit{ij}}^{(S)}$ is impossible using a single simulation. To circumvent this complication the mean flows can be removed artificially (e.g. Rüdiger et al. 2019; Barekat et al. 2021). Thus $Q_{\mathit{ij}}^{(S)}$ vanishes, and the Reynolds stress reads

$$\begin{array}{c}\hfill Q_{\mathit{ij}}=Q_{\mathit{ij}}^{(0)}+Q_{\mathit{ij}}^{(\mathrm{\Omega })},\end{array}$$(4)

where the rotation-generated stress is

$$\begin{array}{c}\hfill Q_{\mathit{ij}}^{(\mathrm{\Omega })}=Q_{\mathit{ij}}-Q_{\mathit{ij}}^{(0)}.\end{array}$$(5)

The non-diffusive stress corresponding to the Λ effect is then

$$\begin{array}{c}\hfill Q_{\mathit{ij}}^{(\mathrm{\Omega })}={\mathrm{\Lambda }}_{\mathit{ijk}}{\mathrm{\Omega }}_k,\end{array}$$(6)

which can be extracted from a single experiment because off-diagonal components of $Q_{\mathit{ij}}^{(0)}$ vanish. The horizontal, meridional, and vertical components of the Λ effect are given by

$$\begin{array}{cc}\hfill Q_{\theta \varphi }^{(\mathrm{\Omega })}& ={\nu }_{\mathrm{t}0}\mathrm{\Omega }Hcos\theta ,\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill Q_{r\theta }^{(\mathrm{\Omega })}& ={\nu }_{\mathrm{t}0}\mathrm{\Omega }Msin\theta cos\theta ,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill Q_{r\varphi }^{(\mathrm{\Omega })}& ={\nu }_{\mathrm{t}0}\mathrm{\Omega }Vsin\theta ,\hfill \end{array}$$(9)

where ν_t0 is a reference value of turbulent viscosity. The factor ν_t0Ω is used for normalization, whereas V, H, and M are dimensionless position-dependent functions.

It is important to bear in mind that Eqs. (7)–(9) are parametrisations derived from symmetry properties but do not specify the physical origin of the stress. Furthermore, Eqs. (7)–(9) are not limited to low Reynolds numbers as the results derived under the first-order smoothing approximation (FOSA) are.

3. The model

The model is the same as in several earlier studies (e.g. Käpylä 2024, and references therein). The equations for compressible hydrodynamics,

$$\begin{array}{cc}\hfill \frac{Dln\rho }{\mathit{Dt}}& =-\mathbf{\nabla }\mathbf{·}\mathit{U},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill \frac{D\mathit{U}}{\mathit{Dt}}& =\mathit{g}-\frac{1}{\rho }(\mathbf{\nabla }p-\mathbf{\nabla }\mathbf{·}2\nu \rho \mathsf{S})-2\mathbf{\Omega }\times \mathit{U},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill T\frac{\mathit{Ds}}{\mathit{Dt}}& =-\frac{1}{\rho }[\mathbf{\nabla }\mathbf{·}({\mathit{F}}_{\mathrm{rad}}+{\mathit{F}}_{\mathrm{SGS}})+\mathcal{C}]+2\nu {\mathsf{S}}^2,\hfill \end{array}$$(12)

were solved, where D/Dt = ∂/∂t + U⋅∇ is the advective derivative, ρ is the density, U is the velocity, $\mathit{g}=-g{\widehat{\mathit{e}}}_z$ is the acceleration due to gravity with g > 0, p is the pressure, ν is the constant kinematic viscosity, S is the traceless rate-of-strain tensor with ${\mathsf{S}}_{\mathit{ij}}=\frac{1}{2}(U_{i,j}+U_{j,i})-\frac{1}{3}{\delta }_{\mathit{ij}}\mathbf{\nabla }\mathbf{·}\mathit{U},$ and Ω = Ω₀(− sin θ, 0, cos θ)^T is the rotation vector, where θ is the angle between Ω and g. T is the temperature, s is the specific entropy, F_rad and F_SGS are the radiative and turbulent sub-grid scale (SGS) fluxes, respectively, and 𝒞 describes surface cooling. The gas is assumed to be optically thick and fully ionised such that radiation is modelled via diffusion approximation. The ideal gas equation of state, $p=(c_{\mathrm{P}}-c_{\mathrm{V}})\rho T=\mathcal{R}\rho T,$ applies, where ℛ is the gas constant, and c_P and c_V are the specific heats at constant pressure and volume, respectively. The radiative flux is F_rad = −K∇T, where K is the radiative heat conductivity given by K = 16σ_SBT³/(3κρ); here, σ_SB is the Stefan-Boltzmann constant and κ is the opacity. Kramers opacity law, κ = κ₀(ρ/ρ₀)^a(T/T₀)^b, with a = 1 and b = −7/2 was used (Weiss et al. 2004); see also Edwards (1990), Brandenburg et al. (2000), and Dorch & Nordlund (2001) for early applications to convection simulations.

Turbulent SGS diffusivity is applied for the entropy fluctuations with F_SGS = −ρTχ_SGS∇s′, where $s\prime (\mathit{x})=s(\mathit{x})-\overline{s}(z)$ and the overbar indicates horizontal averaging. χ_SGS is constant throughout the domain and ${\overline{\mathit{F}}}_{\mathrm{SGS}}\approx 0$ because the SGS diffusivity operates on the entropy fluctuations. The surface cooling is given by $\mathcal{C}=f(z)c_{\mathrm{P}}\rho (T-T_{\mathrm{cool}})/{\tau }_{\mathrm{cool}}$, where ${\tau }_{\mathrm{cool}}=0.4\sqrt{d/g}$ is a cooling timescale; T = e/c_V is the temperature where e is the internal energy; and T_cool = T_top corresponds to the fixed value at the upper boundary. Horizontal averages of horizontal flows ${\overline{U}}_z(z)$ and ${\overline{U}}_y(z)$ were calculated at every time step and removed from the velocity field, so at all times we had U = u. The effects of mean flow removal are discussed in Appendix A.

3.1. Geometry, initial, and boundary conditions

The computational domain is a rectangular box where z_bot ≤ z ≤ z_top is the vertical coordinate, with z_bot/d = −0.45, z_top/d = 1.05, and where d is the depth of the initially isentropic layer (see below). The horizontal coordinates x and y are given by −2d ≤ (x, y)≤2d. The initial stratification consists of three layers. The two lower layers are polytropic, with the polytropic indices n₁ = 3.25 (z_bot/d ≤ z/d < 0) and n₂ = 1.5 (0 ≤ z/d ≤ 1). The uppermost layer above z/d = 1 is initially isothermal. The initial velocity follows a Gaussian-noise distribution with amplitude of the order of ${10}^{-4}\sqrt{\mathit{dg}}$.

The horizontal boundaries are periodic. The vertical boundaries are impenetrable and stress-free for the flow, such that U_x, z = U_y, z = U_z = 0. The temperature gradient at the bottom boundary is given by ∂_zT = −F_bot/K_bot, where F_bot is the fixed input flux and K_bot(x, y, z_bot) is the heat conductivity at z = z_bot. A constant temperature of T = T_top was assumed on the top boundary.

3.2. Units, control parameters, and simulation strategy

The units of length, time, density, and entropy are given by [x]=d, $[t]=\sqrt{d/g}$, [ρ]=ρ₀, [s]=c_P, where ρ₀ is the initial value of density at z = z_top. The profile f(z) = 1 above z/d = 1 and f(z) = 0 below z/d = 1, connecting smoothly over a width of 0.025d, and ${\xi }_0=H_{\mathrm{p}}^{\mathrm{top}}/d=\mathcal{R}{T}_{\mathrm{top}}/gd$ sets the initial pressure scale height at the surface. The Prandtl number based on the radiative heat conductivity is Pr(x, t) = ν/χ(x, t), where χ(x, t) = K(x, t)/c_Pρ(x, t). SGS Prandtl number is given by Pr_SGS = ν/χ_SGS.

The Rayleigh number based on the energy flux is given by

$$\begin{array}{c}\hfill {\mathrm{Ra}}_{\mathrm{F}}=\frac{g{H}^4{F}_{\mathrm{bot}}}{c_{\mathrm{P}}\rho T\nu {\chi }^2}·\end{array}$$(13)

A diffusion-free Rayleigh number can be constructed using Re_cs = c_sH/ν, via

$$\begin{array}{c}\hfill {\mathrm{Ra}}_{\mathrm{F}}^{\mathrm{c}}=\frac{{\mathrm{Ra}}_{\mathrm{F}}{\mathrm{Pr}}^2}{{\mathrm{Re}}_{{\mathrm{c}}_{\mathrm{s}}}^3}=\frac{gH{F}_{\mathrm{bot}}}{c_{\mathrm{P}}\rho T{c}_{\mathrm{s}}^3}·\end{array}$$(14)

Assuming H = H_T ≡ c_PT/g, ${\mathrm{Ra}}_{\mathrm{F}}^{\mathrm{c}}$ reduces to the dimensionless normalised flux (e.g. Brandenburg et al. 2005):

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{n}}=F_{\mathrm{bot}}/\rho c_{\mathrm{s}}^3.\end{array}$$(15)

In rotating simulations, a flux-based, diffusion-free modified Rayleigh number is given by (e.g. Christensen 2002; Christensen & Aubert 2006; Käpylä 2024)

$$\begin{array}{c}\hfill {\mathrm{Ra}}_{\mathrm{F}}^{★}=\frac{{\mathrm{Ra}}_{\mathrm{F}}}{{\mathrm{Pr}}^2{\mathrm{Ta}}^{3/2}},\end{array}$$(16)

where Ta = 4Ω₀²d⁴/ν² is the Taylor number. This can be recast into a flux-based Coriolis number (Käpylä 2024; Bekki 2025):

$$\begin{array}{c}\hfill {\mathrm{Co}}_{\mathrm{F}}=\frac{2{\mathrm{\Omega }}_{0}H}{u_{★}}=2{\mathrm{\Omega }}_{0}H{\left(\frac{\rho }{F_{\mathrm{bot}}}\right)}^{1/3}={({\mathrm{Ra}}_{\mathrm{F}}^{★})}^{-1/3},\end{array}$$(17)

with u_★ ≡ (F_bot/ρ)^1/3 and H = H_T.

In addition to explicit diffusion terms, the advective terms in Eqs. (10)–(12) are written in terms of a fifth-order upwinding derivative with a hyperdiffusive sixth-order correction with a flow-dependent diffusion coefficient; see Appendix B of Dobler et al. (2006). The PENCIL CODE¹ was used to produce the simulations (Pencil Code Collaboration 2021).

3.3. Diagnostics quantities

The global Reynolds, SGS Péclet, and Coriolis numbers,

$$\begin{array}{c}\hfill \mathrm{Re}=\frac{u_{\mathrm{rms}}}{\nu k_1},{\mathrm{Pe}}_{\mathrm{SGS}}={\mathrm{Pr}}_{\mathrm{SGS}}\mathrm{Re},\mathrm{Co}=\frac{2{\mathrm{\Omega }}_0}{u_{\mathrm{rms}}{k}_1},\end{array}$$(18)

describe the importance of viscosity, SGS diffusion, and rotation relative to advection, respectively. Here, u_rms is the volume averaged rms velocity and k₁ = 2π/d is an estimate of the largest eddies in the system. Typically, we find χ ≪ χ_SGS in the convection zone in the current simulations.

The overall anisotropy of the flow is characterised by the parameters

$$\begin{array}{c}\hfill A_{\mathrm{V}}=\frac{Q_{\mathit{xx}}+Q_{\mathit{yy}}-2Q_{\mathit{zz}}}{Q},A_{\mathrm{H}}=\frac{Q_{\mathit{yy}}-Q_{\mathit{xx}}}{Q},\end{array}$$(19)

where Q = TrQ_ij. The parameters A_V and A_H do not provide any information about the scale dependence of anisotropy. This is revealed by spectral anisotropy parameters, which are based on the power spectra of velocity components (e.g. Käpylä 2019a):

$$\begin{array}{cc}\hfill A_{\mathrm{V}}(k)& =\frac{E_x(k)+E_y(k)-2E_z(k)}{E_{\mathrm{K}}(k)},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill A_{\mathrm{H}}(k)& =\frac{E_y(k)-E_x(k)}{E_{\mathrm{K}}(k)},\hfill \end{array}$$(21)

where u_i² = ∫E_i(k)dk, and E_K(k) = ∑_iE_i(k). Diagnostics are time-averaged over the statistically steady part of the simulations, and horizontal or volume averages are additionally applied. Error estimates are given as the mean error of the mean where the number of turnover times is taken as the number of realisations.

4. Results

Some of the current slowly rotating runs use saturated snapshots of the progenitor Run K3h from Käpylä (2019c) as initial conditions, whereas several sets of runs were run from scratch. In total, eight sets of runs with a different Co_F and seven runs corresponding to different values of θ were made; see Table 1.

4.1. Description of the flow fields

Figure 1 shows representative flow fields from runs with slow (Set B), intermediate (Set E), and rapid rotation (Set H), highlighting the transition from cellular convection at slow rotation to columnar convection at rapid rotation. The effects of rotation are not apparent in the most slowly rotating cases at any latitude; see the left column of Fig. 1. For intermediate rotation the average size of convection cells somewhat reduced, and alignment of the convection cells with the rotation vector is discernible at the equator; see the middle column of Fig. 1. For rapid rotation the convection cells are significantly smaller and strongly aligned with the rotation vector; see right column of Fig. 1. This is particularly apparent at the equator (θ = 90°, Run H7), where structures are aligned with the x direction. The dominant size of convective structures at the pole (θ = 0°) was shown to follow the Coriolis-inertial-Archimedean (CIA) scaling, ℓ ∝ Co^−1/2, in an earlier study (Käpylä 2024). A similar result was found for the equator (θ = 90°), corresponding here to Runs B7, E7, and H7 – shown in the bottom row of Fig. 1 – for spectra taken in the longitudinal (y) direction in Bekki (2025). The evident horizontal anisotropy is a key ingredient leading to the non-zero Λ effect.

Fig. 1. Flow fields from Runs B1, B4, and B7 with CoF ≈ 0.55 (left column), Runs E1, E4, and E7 with CoF ≈ 4.5…4.9 (middle), and from Runs H1, H4, and H7 with CoF ≈ 12…17 (right) at latitudes θ = 0° (top row), θ = 45° (middle), and θ = 90° (bottom).

4.2. Anisotropy of turbulence

The vertical and horizontal Λ effects are, to the lowest order, proportional to the turbulence anisotropies given by A_V and A_H, respectively (e.g. Rüdiger 1980). Representative results for A_V and A_H from runs with slow (Set B), intermediate (Set E), and rapid rotation (Set H) are shown in Fig. 2. In the slowly rotating cases, A_V is negative throughout the convection zone, similarly to non-rotating cases (not shown). For Co_F ≈ 0.55 (Set B) the vertical anisotropy, A_V, remains practically constant in the upper part of the convection zone (z/d ≳ 0.65) as a function of latitude. The absolute value of A_V decreases slightly in the lower part of the convection zone from the pole towards the equator. For intermediate (Set E) and rapid rotation (Set H) the latitudinal variation of A_V increases; see panels (b) and (c) of Fig. 2. A_V generally decreases from the pole to the equator. In the most rapidly rotating cases, A_V at the equator is again somewhat higher than in the low-latitude cases; see panel (c) of Fig. 2. Furthermore, while A_V at the pole decreases for intermediate rotation, it is enhanced in the rapidly rotating cases.

Fig. 2. Anisotropy parameters AV (top row) and AH (bottom) from runs in Set B with CoF ≈ 0.55 (left), Set E with CoF ≈ 4.5…4.9 (middle), and from Set H with CoF ≈ 12…17 (right).

The horizontal anisotropy A_H is much weaker than A_V for slow rotation because it is induced by rotation, while the latter is partly an inherent property of convection. A_H vanishes at the pole due to symmetry and increases monotonically towards the equator for slow and intermediate rotation; see panels (d) and (e) of Fig. 2. A_H is negative for rapid rotation (Set H; see panel (f) of Fig. 2), which to the lowest order would imply poleward angular momentum transport that promotes anti-solar differential rotation (e.g. Rüdiger 1980). However, the latter result is valid for slow rotation corresponding to Co ≪ 1, which is not satisfied in the current rapidly rotating runs. Furthermore, this argument relies only on the magnitude of A_V, while spatial scales of the horizontal flows are also very different in the latitudinal (x) and longitudinal (y) directions; see Fig. 1.

To study the anisotropy in more detail, a spectral decomposition was made; see Eqs. (20) and (21). Representative results for A_V(k) and A_H(k) are shown in Figs. 3 and 4. Large scales are dominated by horizontal flows such that A_V(k) > 0, similarly to the non-rotating case (not shown). The zero-crossing of A_V(k) occurs at a progressively higher normalised wave number of $\stackrel{\sim }{k}=k/k_1$ as the rotation rate increases. The negative A_V(k) at intermediate and large wave numbers are responsible for A_V = ∫A_V(k)dk < 0 in the bulk of the convection zone in all cases. For slow rotation (Runs E[1,4,7]) A_V(k) does not change appreciably as a function of θ, whereas for the intermediate and rapid rotation runs the crossover from positive to negative A_V(k) occurs at a lower wave number; see the middle and right panels of Fig. 3.

Fig. 3. Anisotropy parameter AV(k) from runs with slow (Runs B1, B4, and B7), intermediate (Run E1, E4, and E7), and rapid rotation (Run H1, H4, and H7) at θ = 0° (left panel), θ = 45° (middle), and θ = 90° (right), respectively, near the middle of the convection zone at z/d = 0.49.

Fig. 4. Anisotropy parameter AH(k) for the same runs as in Fig. 3.

As discussed above, the horizontal anisotropy becomes prominent only at rapid rotation, and it vanishes for θ = 0°; see Fig. 4. For slow and intermediate rotation the main contribution to the overall positive A_H comes from large and intermediate scales; see the middle and right panels of Fig. 4. This can be understood as a consequence of rotational alignment of the largest convection cells that have the greatest Coriolis numbers. In the rapid rotation regime a similar trend is seen on large scales, but the behaviour for $\stackrel{\sim }{k}\gtrsim 10$ of the latitudinal (x) component dominates, leading to A_H = ∫A_Hkdk < 0.

Turbulence becomes quasi-two-dimensional for sufficiently rapid rotation leading to diminishing velocity along Ω (e.g. Müller & Thiele 2007). Strong density stratification has a similar effect along the direction of the gravitational acceleration (e.g. Boffetta 2023). The latter has been adopted as one of the building blocks for theories of angular momentum transport in stellar convection zones (e.g. Kitchatinov & Rüdiger 1993). This approach was generalised by Kitchatinov & Rüdiger (2005), the authors of which introduced an anisotropy parameter, a, that characterises the difference between vertical and horizontal motions in the unperturbed case; see their Eq. A14. Assuming for simplicity that the scale ratio ℓ_corr²/L² = 1, the vertical anisotropy parameter, A_V, can be written in terms of a as follows:

$$\begin{array}{c}\hfill A_{\mathrm{V}}=\frac{1}{2}-\frac{1}{5}a,\end{array}$$(22)

such that $a=\frac{5}{2}$ corresponds to isotropy. The range of A_V in Fig. 3 inside the convection zone is roughly −0.2 to −1.2, corresponding to a = 3.5…8.5, which is significantly higher than the values typically used in mean-field approaches (e.g. Kitchatinov & Rüdiger 2005; Pipin & Kosovichev 2018).

4.3. Reynolds stress and Λ effect

Horizontally and temporally averaged off-diagonal Reynolds stresses are shown in Fig. 5. The stresses are assumed to be due to non-diffusive origin, and are attributed here to the Λ effect since horizontally averaged mean flows have been suppressed in the present simulations. The horizontal stress, Q_xy, is responsible for generating horizontal or latitude-dependent differential rotation in spherical coordinates. Q_xy is on average positive almost everywhere corresponding to equatorward flux of angular momentum; see the left column of Fig. 5. For the slowest rotation rates the data are not sufficiently converged to rule out negative values, which, however, are theoretically expected only in cases where shear is retained (Rüdiger et al. 2019). As the rotation rate is increased, the horizontal stress is increasingly concentrated near the equator; see the data for Sets E to H in the four bottom left panels of Fig. 5. This behaviour is commonly seen in f-plane simulations (e.g. Chan 2001; Käpylä et al. 2004; Hupfer et al. 2005) of convection but not in corresponding forced turbulence calculations (e.g. Käpylä & Brandenburg 2008; Käpylä 2019b). The main difference between the forced turbulence and convection set-ups is the lack of stratification and thermodynamics in the former such that, for example, thermal Rossby waves are not excited. Convective structures corresponding to the latter are the main contributor to Q_xy near the equator (see also Käpylä et al. 2011a). On the other hand, global convection simulations do not show a similar concentration of horizontal stress (Q_θϕ) near the equator (e.g. Käpylä et al. 2011a). This is likely explained by the differences in geometry; in the current set-up the horizontal stress cannot generate a mean flow in the horizontally averaged sense because of the horizontal periodicity of the domain, whereas in global simulations the divergence of Q_θϕ leads to a mean flow. However, it is possible to generate mean flows that depend on x or y but vanish under a horizontal average even in the present geometry, although this has not been observed in the current simulations. A possible reason is that the horizontal size of the domain is too small or that the horizontal aspect ratio of the domain needs to be unequal to unity, which effectively introduces a preferred direction (e.g. Guervilly & Hughes 2017; Currie et al. 2020). Furthermore, the formation of large-scale vortices (e.g. Chan 2007; Käpylä et al. 2011b; Guervilly et al. 2014) was not observed in the current simulations, most likely because the fluid Reynolds number is too low; see the discussion in Käpylä (2024).

Fig. 5. Off-diagonal Reynolds stresses ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ (left column), ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ (middle), and ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ (right) from all of the runs. The set of runs and the corresponding CoF are denoted in the left panel of each row.

The meridional stress, Q_xz, does not directly participate in the angular momentum transport, but it has an indirect influence via the meridional flow. This component is predominantly negative for slow rotation (Co ≲ 0.25, Sets A to C), with a magnitude similar to Q_xy. For Co ≳ 0.6, Q_xz first changes sign for θ = 75° and, at even more rapid rotation, for θ = 60°. This is in contrast to earlier forced turbulence simulations, which also yielded Q_xz < 0 for rapid rotation (e.g. Käpylä 2019b). It is possible that the difference is also due to the thermal Rossby waves that were absent in the isothermal forced turbulence models. The magnitude of Q_xz is less than the other off-diagonal components for Co ≳ 0.6. The impact of this stress component for stellar differential rotation has not been studied in detail.

The vertical stress Q_yz generates radial differential rotation. For slow rotation (Co ≲ 0.6) Q_yz is negative almost everywhere and roughly proportional to Ω. This agrees with analytic results and forced turbulence simulations (e.g. Käpylä & Brandenburg 2008; Käpylä 2019b). However, at sufficiently rapid rotation the sign of Q_yz changes near the equator where large positive values are observed; see the data for Sets F to H in the middle column of Fig. 5. This feature is absent in isothermal forced turbulence models and is due to thermal Rossby waves that are excited in rotating convection (e.g. Busse 2002). Analyses of the radial stress Q_rϕ – corresponding to Q_yz in the current Cartesian geometry – in spherical shells and wedges have shown that the large-scale convective modes associated with the thermal Rossby waves produce the dominant contribution to Q_rϕ and produce equatorial acceleration (e.g. Käpylä 2023). This is in contrast to some hydrodynamic mean-field models, where the radial stress vanishes for rapid rotation and equatorial acceleration is achieved via equatorward horizontal flux of angular momentum (e.g. Kitchatinov 2013).

The non-diffusive Reynolds stresses from Eqs. (7) to (9) can be represented in terms of the Λ-effect coefficients as

$$\begin{array}{cc}\hfill \mathcal{H}& \equiv Hcos\theta =\frac{2}{15}{\mathrm{Co}}^{-1}{\stackrel{\sim }{Q}}_{\mathit{xy}},\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill \mathcal{M}& \equiv Msin\theta cos\theta =\frac{2}{15}{\mathrm{Co}}^{-1}{\stackrel{\sim }{Q}}_{\mathit{xz}},\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill \mathcal{V}& \equiv Vsin\theta =\frac{2}{15}{\mathrm{Co}}^{-1}{\stackrel{\sim }{Q}}_{\mathit{yz}},\hfill \end{array}$$(25)

where ${\nu }_{\mathrm{t}0}=\frac{4}{15}{u}_{\mathrm{rms}}{k}_1$ was used, and where the tilde refers to normalisation by u_rms². The functional form of ν_t0 is taken for a non-rotating case from the analytic study of Kitchatinov et al. (1994). Results from numerical simulations of isotropically forced homogeneous turbulence tend to be somewhat higher; see Käpylä et al. (2020). Therefore, the choice of ν_t0 in the normalisation leads to some uncertainty in the absolute values of the coefficients. We followed the same procedure as in Käpylä (2019b) and expanded the coefficients in powers of sin²θ

$$\begin{array}{c}\hfill H={\displaystyle \sum _1^{n_{\max }}}{H}^{(i)}{sin}^{2n}\theta ,\end{array}$$(26)

$$\begin{array}{c}\hfill M={\displaystyle \sum _0^{n_{\max }}}{M}^{(i)}{sin}^{2n}\theta ,\end{array}$$(27)

$$\begin{array}{c}\hfill V={\displaystyle \sum _0^{n_{\max }}}{V}^{(i)}{sin}^{2n}\theta ,\end{array}$$(28)

where H⁽⁰⁾ is assumed to vanish due to symmetry. The Reynolds stresses vary as a function of height, but here we applied volume averaging, denoted by angle brackets ⟨.⟩ over the 0 ≤ z/d ≤ 1 range to simplify the analysis. The coefficients V, H, and M were fitted with up to n_max = 3 using Eqs. (23)–(25) and Eqs. (26)–(28) with volume-averaged numerical data for the Reynolds stresses. The fit was deemed accurate enough when the mean error of the mean between the numerical data and the reconstruction with the fitted coefficients did not decrease by more then ten per cent when an additional coefficient was taken into account. However, an n_max higher than three would be needed to fit the vertical stress $⟨{\stackrel{\sim }{Q}}_{\mathit{yz}}⟩$ for the most rapidly rotating cases (Sets G and H), but this was not pursued. The numerical data as well as the fits are shown in Fig. 6 for all runs. Furthermore, Table 2 summarises the fitted coefficients.

Fig. 6. Volume and time-averaged Λ coefficients ⟨ℋ⟩ (left panels), ⟨ℳ⟩ (middle), and ⟨𝒱⟩ (right) from all sets of runs along with best fits to Eqs. (28)–(27). The thick lines show an adequate fit according to the criteria in the text, and the thin lines show fits with one fewer coefficient for comparison. The colours denote the maximum number of coefficients, nmax, taken into account in the fits such that nmax = 0, 1, 2, and 3 correspond to yellow, orange, red, and purple, respectively.

The horizontal Λ effect is positive for all rotation rates, although in the slowest rotation cases studied (Sets A, B, and C) the results are not statistically significant. These results are in accordance with quasi-linear theory (e.g. Rüdiger 1989), where negative values of ℋ are expected only when large-scale flows (Rüdiger et al. 2019) or magnetic fields (e.g. Kitchatinov & Rüdiger 2004; Käpylä 2019b) are present, both of which are absent in the current simulations. The concentration of the horizontal stress near the equator is reflected in the horizontal Λ effect coefficient, ℋ, although this behaviour is subdued in Fig. 6 due to the volume averaging. Nevertheless, ⟨ℋ⟩ tends to be at its maximum at the nearest latitude away from the equator for Co > 0.6 (Sets E to H). This is reflected by the need to use coefficients up to H⁽³⁾ to already fit the data in Set E; see the second to fourth columns of Table 2. More data points at lower latitudes would be needed to determine the full latitude profile, but this is out of the scope of the present study. In Käpylä (2019b) all of the H⁽ⁱ⁾ coefficients were found to be positive, which is not the case for H⁽¹⁾ in the present study.

The volume-averaged meridional Λ effect, ⟨ℳ⟩, shown in the middle column of Fig. 6, is negative for slow rotation up to around Co = 0.25 (Set C). For more rapid rotation, ⟨ℳ⟩ changes sign near the equator which was not observed in the forced turbulence simulations of Käpylä (2019b), but was present in earlier simulations by Käpylä & Brandenburg (2008). Pulkkinen et al. (1993)² reported M⁽⁰⁾ = −0.03 and M⁽¹⁾ = 0.06, although the quality of their fit is rather modest; see their Fig. 11. The signs coincide with the current Set C, although the values are greater than in the current runs. This is possibly because mean flows were retained in the runs of Pulkkinen et al. (1993) that also contribute to the Reynolds stress; (see e.g. Rüdiger et al. 2019). For more rapid rotation (Sets D onwards), at least three coefficients are needed to fit the data. The sign of M⁽⁰⁾ coincides with that in forced turbulence simulations of Käpylä (2019b), but the behaviour of the higher order coefficients is more complex.

The vertical Λ effect, ⟨𝒱⟩, is negative for for Co_F ≲ 4.56 (Sets A to E). In the two slowest rotating sets A and B, only the fundamental mode of the Λ effect, described by V⁽⁰⁾ is needed to fit the numerical data; see Table 2. The absolute value of V⁽⁰⁾ is of the order of 0.05 in the slowest rotation cases. The sign and functional form agree with forced turbulence simulations (e.g. Käpylä 2019b; Barekat et al. 2021), but the absolute value is an order of magnitude smaller in the current results although the vertical anisotropy of the flow as measured by Λ_V is similar in both cases; compare Fig. 2 with Fig. 3 of Barekat et al. (2021). For Co ≳ 2.5 (Sets F to H), the sign of ⟨𝒱⟩ changes near the equator. In the current results V⁽⁰⁾ < 0 and V⁽³⁾ > 0 everywhere, whereas V⁽¹⁾ > 0 and V⁽²⁾ < 0, apart from in Set E; see the three last columns in Table 2. In Käpylä (2019b), V⁽⁰⁾ is also predominantly negative, whereas V⁽¹⁾ < 0 and V⁽²⁾ > 0 in contrast to the current results. Analytic studies do not typically consider coefficients higher than V⁽¹⁾; see, for example, Kitchatinov & Rüdiger (2005), Pipin & Kosovichev (2018).

The qualitative differences between the analytic theories and isothermal forced turbulence simulations in comparison to the current results for the vertical Λ effect are likely due to the dominating influence of thermal Rossby waves that appear as tilted columnar convection cells near the equator (e.g. Käpylä 2023). Thermal Rossby waves are typically not captured by analytics that rely on simplified turbulence models, and therefore a positive ⟨𝒱⟩ for rapid rotation is absent in such theories. This often also applies to analytic theories where the convective heat flux is taken into account (e.g. Kleeorin & Rogachevskii 2006); see, however Rogachevskii & Kleeorin (2018). Yet, Pipin & Kosovichev (2018) found that the standard mean-field theory does yield a sign change of the vertical Λ effect for rapid rotation when the spatial non-uniformity of the convective turnover time is taken into account.

5. Conclusions

The current simulations show that the radial, non-diffusive angular momentum transport, parameterised by the vertical Λ effect, is negative for slow rotation (Co ≲ 1.3) in rotating density-stratified convection. The horizontal Λ effect is positive, corresponding to equatorward angular momentum flux, for all rotation rates. These findings are in qualitative agreement with isothermal forced turbulence simulations (e.g. Käpylä & Brandenburg 2008; Käpylä 2019b; Barekat et al. 2021) and analytic mean-field theories (e.g. Kitchatinov & Rüdiger 1995; Kleeorin & Rogachevskii 2006). For more rapid rotation the radial Λ effect changes sign, and the horizontal flux is increasingly concentrated at low latitudes. The change of sign of of the vertical Λ effect is due to the emergence of prograde propagating thermal Rossby waves, which are not accounted for in analytic theories of angular momentum transport (e.g. Rüdiger 1989; Kitchatinov & Rüdiger 2005). In global convection simulations the thermal Rossby waves are the cause of solar-like differential rotation. On the other hand, thermal Rossby waves have not been observed in the Sun (e.g. Birch et al. 2024, and references therein), and if they are present, they must be much weaker than in simulations. Therefore, the origin of solar and stellar differential rotation is related to the ongoing debate regarding the nature of convection in the Sun (e.g. Schumacher & Sreenivasan 2020; Hotta et al. 2023, and references therein), often referred to as the convective conundrum (e.g. O’Mara et al. 2016).

The current study probes the hydrodynamic regime with simulations where the density stratification is much weaker than in stars and where the transition to optically thin photosphere is modelled rather crudely. Magnetic fields have been shown to strongly affect the angular momentum transport and large-scale flows in global simulations (e.g. Hotta et al. 2022; Hotta 2025; Käpylä 2023; Soderlund et al. 2025). Magnetic fields have also been shown to significantly influence the Λ effect in theoretical and idealised numerical studies (e.g. Kitchatinov & Rüdiger 2004; Käpylä 2019b,a). The inability to incorporate a realistic stratification or proper radiative surface possibly inhibits non-local entropy rain (e.g. Spruit 1997; Brandenburg 2016; Käpylä 2025) in current studies. The effects on convective angular momentum transport of these aspects need to be addressed in future studies.

Acknowledgments

I thank the anonymous referee, Valery Pipin, and Igor Rogachevskii for their comments on the manuscript. The simulations were performed using the resources granted by the Gauss Center for Supercomputing for the Large-Scale computing project “Cracking the Convective Conundrum” in the Leibniz Supercomputing Centre’s SuperMUC-NG supercomputer in Garching, Germany. This work was supported in part by the Deutsche Forschungsgemeinschaft Heisenberg programme (grant No. KA 4825/4-1).

References

Appendix A: Effects of mean flow removal

According to Eq. (2) the Reynolds stress can be decomposed to contributions from the background turbulence, rotation, and shear. In practice this is complicated by the inherent non-linearity of the Navier-Stokes equations. For example, sufficiently strong shear leads to production of turbulence; see e.g. Käpylä et al. (2020). Removing shear flows is therefore also likely to affect the solutions. Here we consider a few otherwise identical simulations where the mean flows are either a part of the solution or where they are removed. The latter correspond to the runs listed in Table 1.

We use two diagnostics to study the effect of mean flow removal, which are the overall velocity amplitude u_rms and the horizontally averaged off-diagonal Reynolds stress component Q_yz. These are shown for four representative runs in Fig. A.1. The change of the volume-averaged rms velocity reflects the strength of the removed mean flows. The mean flows in the slowly rotating Run D6 are relatively strong due to which the rms-velocity is decreased by roughly 5 per cent. The clearest effect is seen at the equator (Run F7), where strong mean flows are generated and their removal leads to a decrease of the flow amplitude by about 20 per cent. Mean flows in runs at higher latitudes and faster rotation tend to be weaker; see the data for Runs F4 and H5. These results regarding mean flows are in accordance with earlier studies (e.g. Chan 2001; Käpylä et al. 2004). Similarly the Reynolds stresses are modified by the absence of mean flows because the diffusive contribution to the stress is absent. This effect is again strongest in cases with larger mean flows as indicated by the lower panel of Fig. A.1. Assuming that the turbulent viscosity ν_t is positive and because ${\partial }_z{\overline{U}}_y<0$ in Run D6, we would expect the diffusive contribution to the stress to be positive, that is $Q_{\mathit{yz}}^{(S)}>0$. This agrees with the results of the current simulations; see the black lines in the lower panel of Fig. A.1. In Runs F4 and H5 the differences between the runs with and without shear are small and not always suggestive of a positive ν_t. This can be due to insufficient statistical convergence. The results of Run F7, where Q_yz should vanish in the case where the mean flows are retained (e.g. Bekki 2025) are suggestive of this. The non-zero residual (red solid line) in the lower panel of Fig. A.1 can be used as an estimate of the statistical uncertainty of the results.

Fig. A.1. Volume-averaged rms velocity urms/(dg)1/2 as a function of time (upper panel) and Reynolds stress component Qyz/dg in units of 10−4 as a function of z (lower panel) from Runs D6 (dotted black line), F4 (dotted blue), F7 (dotted red), and H5 (dotted orange). The solid lines show the same quantities from corresponding runs where mean flows are retained.

All Tables

All Figures

	Fig. 1. Flow fields from Runs B1, B4, and B7 with CoF ≈ 0.55 (left column), Runs E1, E4, and E7 with CoF ≈ 4.5…4.9 (middle), and from Runs H1, H4, and H7 with CoF ≈ 12…17 (right) at latitudes θ = 0° (top row), θ = 45° (middle), and θ = 90° (bottom).
In the text

	Fig. 2. Anisotropy parameters AV (top row) and AH (bottom) from runs in Set B with CoF ≈ 0.55 (left), Set E with CoF ≈ 4.5…4.9 (middle), and from Set H with CoF ≈ 12…17 (right).
In the text

	Fig. 3. Anisotropy parameter AV(k) from runs with slow (Runs B1, B4, and B7), intermediate (Run E1, E4, and E7), and rapid rotation (Run H1, H4, and H7) at θ = 0° (left panel), θ = 45° (middle), and θ = 90° (right), respectively, near the middle of the convection zone at z/d = 0.49.
In the text

	Fig. 4. Anisotropy parameter AH(k) for the same runs as in Fig. 3.
In the text

	Fig. 5. Off-diagonal Reynolds stresses ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ (left column), ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ (middle), and ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ (right) from all of the runs. The set of runs and the corresponding CoF are denoted in the left panel of each row.
In the text

Fig. 6. Volume and time-averaged Λ coefficients ⟨ℋ⟩ (left panels), ⟨ℳ⟩ (middle), and ⟨𝒱⟩ (right) from all sets of runs along with best fits to Eqs. (28)–(27). The thick lines show an adequate fit according to the criteria in the text, and the thin lines show fits with one fewer coefficient for comparison. The colours denote the maximum number of coefficients, nmax, taken into account in the fits such that nmax = 0, 1, 2, and 3 correspond to yellow, orange, red, and purple, respectively.

In the text

	Fig. A.1. Volume-averaged rms velocity urms/(dg)1/2 as a function of time (upper panel) and Reynolds stress component Qyz/dg in units of 10−4 as a function of z (lower panel) from Runs D6 (dotted black line), F4 (dotted blue), F7 (dotted red), and H5 (dotted orange). The solid lines show the same quantities from corresponding runs where mean flows are retained.
In the text