On the origin of the double cell meridional circulation in the solar convection zone

Recent advances in helioseismology, numerical simulations and mean-field theory of solar differential rotation have shown that the meridional circulation pattern may consist of two or more cells in each hemisphere of the convection zone. According to the mean-field theory the double-cell circulation pattern can result from the sign inversion of a nondiffusive part of the radial angular momentum transport (the so-called $\Lambda$-effect) in the lower part of the solar convection zone. Here, we show that this phenomenon {can result} from the radial inhomogeneity of the Coriolis number, which depends on the convective turnover time. We demonstrate that if this effect is taken into account then the solar-like differential rotation and the double-cell meridional circulation are both reproduced by the mean-field model. The model is consistent with the distribution of turbulent velocity correlations determined from observations by tracing motions of sunspots and large-scale magnetic fields, indicating that these tracers are rooted just below the shear layer.

reproduced by the mean-field model. The model is consistent with the distribution of turbulent velocity correlations determined from observations by tracing motions of sunspots and large-scale magnetic fields, indicating that these tracers are rooted just below the shear layer.

INTRODUCTION
Angular momentum transport on the Sun is tightly related to the influence of the Coriolis force on motions inside the convection zone. The mean-field hydrodynamic theory predicts that in addition to viscous redistribution of the angular momentum by convective motions there is a nondissipative contribution to turbulent stresses. This contribution is called the Λ-effect (Kitchatinov & Rüdiger 1993).
According to this theory, the differential rotation profile in stellar convection zones is established as a result of a nonlinear balance between the turbulent stresses and meridional circulation (Durney 1999). The meridional circulation itself is driven by perturbation of the Taylor-Proudman balance between the centrifugal and baroclinic forces. The baroclinic forces result from the nonuniform distribution of the mean entropy because of the anisotropic heat transport in a rotating convective zone (Durney 1999).
Along with the solar-like angular velocity profiles, the standard mean-field models predict the meridional circulation structure with one circulation cell in each hemisphere. Contrary, direct numerical simulations, in most cases, predict the meridional circulation structure with multiple cells (see, Miesch & Hindman 2011;Gastine et al. 2013;Guerrero et al. 2013;Käpylä et al. 2014;Featherstone & Miesch 2015;Hotta et al. 2015). The double-cell (or multiple-cell) structure of the meridional circulation has been suggested by recent helioseismology inversions (Zhao et al. 2013;Schad et al. 2013;Kholikov et al. 2014), but is still under debate. For example, Rajaguru & Antia (2015) found that the meridional circulation can be approximated by a single-cell structure with the return flow deeper than 0.77R . However, their results indicate an additional weak cell in the equatorial region, and contradict to the recent results of Böning et al. (2017) who confirmed a shallow return flow at 0.9R .
Recent mean-field hydrodynamic modeling by Pipin & Kosovichev (2016) and Bekki & Yokoyama (2017) (hereafter BY17) showed that the double-cell meridional circulation structure can be reproduced by tuning the nondissipative turbulent stresses, i.e., the Λ-effect. In particular, BY17 argued that the double-cell meridional circulation structure can be explained if the radial transport of angular momentum by the Λ-effect changes sign at some depth of the convection zone. This effect was demonstrated by prescribing ad hoc radial profile of the Λ-tensor, changing sign at the midpoint of the convection zone, at r = 0.825 R . Direct numerical simulations by Käpylä et al. (2011) also showed that the radial Λ-effect can inverse sign in a case of high Coriolis number Ω * = 2Ωτ c , where Ω is the global mean rotation rate, and τ c is the local turnover time of convection. These results motivated us to search for the turbulent mechanism that can explain the sign inversion of the Λ-effect in the solar convection zone.
In this paper, we show that the sign inversion of the radial Λ-effect follows naturally from the standard mean-field hydrodynamics theory (Kitchatinov & Rüdiger 1993;Kitchatinov 2004). By employing the standard solar model and the mixing-length theory of convective energy transport, we show that the sign inversion of the Λ-effect is located in the lower convection zone, at r 0.78 R .
The key point is that in addition to the density stratification considered in the previous mean-field models, the radial gradients of the functions that describe effects of the Coriolis force have to be taken into account. We demonstrate that when these gradients are taken into account then both the solar-like distribution of the angular velocity profile and the double-layer meridional circulation structure are both reproduced by the mean-field model.

The angular momentum balance
We decompose the axisymmetric mean velocity into poloidal and toroidal components: U = U m + r sin θΩφ, whereφ is the unit vector in the azimuthal direction. The mean flow satisfies the stationary continuity equation, Similarly to Rüdiger (1989), the conservation of the angular momentum is expressed as follows: where the turbulent stresses tensor,T, is written in terms of small-scale fluctuations of velocity: where, in following the mean-field hydrodynamic framework (see Rüdiger 1989;Kitchatinov et al. 1994), the turbulent stress tensor is expressed as a sum of two major parts: the first term,T i,j = Λ ijk U k , represents the nondissipative part (the Λ-effect), and the second term,T describes the eddy viscosity tensor contribution. The analytical expressions forT is given in the Appendix.
To determine the meridional circulation, we consider the azimuthal component of the large-scale vorticity, ω = ∇ × U m φ , which is governed by the following equation: where ∂/∂z = cos θ∂/∂r − sin θ/r · ∂/∂θ is the gradient operator along the axis of rotation. Turbulent stresses affect generation and dissipation of large-scale flows, and, in turn, they are affected by global rotation and magnetic field. In this paper, we neglect magnetic field effects. The magnitude of kinetic coefficients in tensorT depends on the convective turnover time, τ c , and on the RMS of the convective velocity, u . The radial profile of τ c is obtained from the standard solar interior model calculated using the MESA code (Paxton et al. 2011(Paxton et al. , 2013. The RMS velocity, u , is determined in the mixing-length approximations from the gradient of the mean entropy, s, where = α MLT H p is the mixing length, α MLT = 2.2 is the mixing-length theory parameter, and H p is the pressure scale height. For a nonrotating star, the u (r) profile corresponds to results of the MESA code. The mean-field equation for the heat transport takes into account effects of rotation: where, ρ and T are the mean density and temperature,and F conv and F rad are the convective and For calculation of the heat eddy-conductivity tensor, χ ij , we take into account effects of global rotation: where functions φ and φ were defined in KPR94.
The eddy conductivity and viscosity are determined from the mixing-length approximation: where Pr T is the turbulent Prandtl number. We found that Pr T = 3 4 gives the magnitude of the latitudinal differential rotation on the surface in agreement with solar observations.

The Λ-effect profiles
In analytical derivations of the Λ effect within the mean-field hydrodynamics framework, it is assumed that the background turbulence can be modeled as a randomly forced quasi-isotropic spatially inhomogeneous turbulent flow. Furthermore, the mean flow generation effect appears in the secondorder terms of the Taylor expansions in terms of the large-scale inhomogeneity parameter, /L, where is a characteristic size of convective vortexes, and L is a spatial scale of mean-field parameters. In is the density scale height. This approximation is different from the assumptions employed in derivations of the heat eddy-conductivity tensor, χ ij (see Eq.7), or the eddy viscosity tensor N ijkl . In this case, it is assumed that the background turbulent flow is spatially homogeneous. The dissipative effects appear when gradients of the large-scale flow are taken into account, i.e., additional terms of the order of /L in Eq(4) (e.g. see Kitchatinov et al. 1994). A detailed discussion of approximations and assumptions of the mean-field hydrodynamics can be found in Rüdiger (1989).
For accurate description of the angular momentum transfer, it is important to take into account the effect of the radial profile of the Coriolis number, Ω * , on the Λ effect. To illustrate this, we consider the case of fast rotation, Ω * 1. The nondiffusive flux of angular momentum can be expressed as follows: where J 1 = − π 4Ω * , and and Λ V , Λ H are the vertical and horizontal components of the Λ effect. We reproduce the original results of Kitchatinov & Rüdiger (1993) if the contribution of the Coriolis number gradient, ∂ log Ω * ∂r , is neglected. We would like to stress that the dependence of the Λ effect coefficients on the Coriolis number is included in models M1 and M2. Those models disregard the effect of the spatial derivatives of the Ω * which seems to be important near the bottom of the convection zone.
However, because the convective overturn time, τ c , in stellar convection zones changes with depth, the radial gradient of the Coriolis number, Ω * = 2Ωτ c , can be essential for the amplitude and direction of the Λ effect. The τ c -profile can be obtained from results of the MESA code using the standard mixing-length approximation, τ c = /u . Note that with the use of the mixing-length theory for the solar interior model, this formula is equivalent to another expression (see Kippenhahn & Weigert 1994): where F con is the amount of heat flux transported by convection, and δ = − ∂ log ρ ∂ log T P (δ = 1 for an ideal mono-atomic gas). The profile of the Coriolis number calculated from the solar interior model is shown in Figure 1a. It is seen that parameter | ∂ log Ω * ∂r | is greater than unity in the lower part of the convection zone, and thus it should be taken into account.
Below, we consider results for the mean-field hydrodynamic models of the solar differential rotation based on three models of the Λ effect, listed in Table 1. We would like to stress that the dependence of the Λ effect coefficients on the Coriolis number is included in models M1 and M2, but these models disregard the radial derivative of the Ω * . Note that models M1 and M2 differ only in the parameter of anisotropy, which is a = 0 in M1 and a = 2 in M2. Model M3 includes both the radial derivative of Ω * and the convective anisotropy.  (9) and (11)) in three models: a) M1; b) M2; c) M3, listed in Table 1.  In all of these cases, the direction of the nondiffusive angular momentum transport is predominantly along the axis of rotation from higher to lower latitudes in the main part of the convection zone.
Model M1 shows the outward transport by the Λ effect in the upper part of the convection zone.
where ζ = r sin θ and we take into account cylinder-like distribution of L. Note, that ∇ ζ L > 0.
The double-cell meridional circulation like that found by Zhao et al. (2014) gives the negative U m ζ at the bottom, positive in the middle, and negative at the top of the convection zone. Then, the LHS of Eq(14) gets the same signs. Therefore, Eq (14) can be satisfied for a radially converging vector fieldT φ . The pattern demonstrated in Figure 2c satisfies this condition as well as the Λ-effect model introduced by BY17 (see Fig. 1 in their paper). As it is seen in Figure 2 that for the case of the moderate and fast rotation regimes, Ω * > 1, the mean-field theory predicts a nearly cylinderlike distribution of the angular momentum fluxes produced by the Λ-effect in the rotating stratified convective media. We note that effects of the convective velocity anisotropy are strongly quenched with the increase of the Coriolis number (Kitchatinov et al. 1994;Kitchatinov 2004).

RESULTS
Numerical solutions of Eq. (1)-(5) for the three models of Table 1  Contrary to results of BY17 signs of the torque from the Λ-effect are not fully balanced with the meridional circulation. This means that the effect of the diffusive part of the angular momentum is also important. Our models take into account the anisotropic eddy viscosity. Its importance for the mean-field models of the solar differential rotations was addressed previously by Küker et al.  Our results suggest that the model of the meridional circulation is sensitive to the choice of the radial Coriolis number profile in the convection zone. In particular, we found that there is no sign inversion of the Λ-effect for the profile of Stix (2002), , where x = r/R , and Ω * 0 is the Corioils number at x 0 = 0.8 (e.g., Ω * 0 ≡ 3). Using this profile and taking into account the inhomogeneity of the Coriolis number we get no inversion of the Λ-effect.
Results of that model are very similar to our model M2. Another need for the future development follows from direct numerical simulations. They often show multiple meridional circulation cells, see, e.g., Guerrero et al. (2016) or Warnecke et al. (2018).
These results are difficult to explain within our model. In our paper, we discussed components of the Λ-effect, which contribute to the generation of the azimuthal large-scale flow. In addition to them, there is a theoretical possibility for the Λ-effect components which generate the meridional circulation. This effect results from the nondissipative part of the off-diagonal turbulent stresses: b θ b r , where u and b are the fluctuating velocity and magnetic field. The origin of this effect can be tightly related to inhomogeneities of the kinetic and magnetic helicities of turbulent flows (Yokoi & Brandenburg 2016). It is interesting that the so-called "anisotropic kinetic alpha -effect" or the "AKA-effect" (Frisch et al. 1987;Pipin et al. 1996;Brandenburg & Rekowski 2001) also results in the nondissipative part of theT rθ (as well as the azimuthal components ofT).
Therefore, the theory of the Λ-effect, which is employed in our paper is rather incomplete.
The model correlation u θ u φ can be compared with motions of sunspots (Sudar et al. 2017) and the large-scale magnetic fields (Latushko 1993). Both of these measurements show the positive correlation in the northern hemisphere of the Sun. The sign of u θ u φ is the main reason for the solar-like differential rotation with the equator rotating faster than the poles (Rüdiger 1989). Model Expression of the turbulent stress tensor results from the mean-field hydrodynamics theory (see, Kitchatinov et al. 1994;Kitchatinov 2004) as follows: where u is fluctuating velocity. Application the mean-field hydrodynamic framework leads to the Taylor expansion given by Eq(4). The viscous part of the azimuthal components of the stress tensor is determined following Kitchatinov et al. 1994 in the following form: where the eddy viscosity, ν T , is determined from the mixing-length theory assuming the turbulent Prandl number P r T = 3 4 : The viscosity quenching functions, Φ ⊥ and Φ , depend nonlinearly on the Coriolis number, Ω * = 2Ωτ c and they are determined by Kitchatinov et al. (1994).
The nondiffusive flux of angular momentum can be expressed as follows (Rüdiger 1989): and H (1) = −V (1) . We employ the parameter of the turbulence anisotropy a = u 2 h − 2u 2 r u 2 r =2, where u h and u r are the horizontal and vertical RMS velocities (Kitchatinov 2004).
The first RHS term of Eq.(5) describes dissipation of the mean vorticity, ω. Similarly to Rempel (2005) we approximate it as follows, where the rotational quenching function φ 1 is given by Kitchatinov et al. (1994). We have tried to apply a more general formalism including all components of the eddy-viscosity tensor for rotating turbulence provided by Kitchatinov et al. (1994), and obtained results that are similar to the approximation given by Eq (28).