Anisotropic turbulent transport with horizontal shear in stellar radiative zones

Context. Turbulent transport in stellar radiative zones is a key ingredient of stellar evolution theory, but the anisotropy of the transport due to the stable stratification and the rotation of these regions is poorly understood. The assumption of shellular rotation, which is a cornerstone of the so-called rotational mixing, relies on an efficient horizontal transport. However, this transport is included in many stellar evolution codes through phenomenological models that have never been tested. Aims. We investigate the impact of horizontal shear on the anisotropy of turbulent transport. Methods. We used a relaxation approximation (also known as {\tau} approximation) to describe the anisotropising effect of stratification, rotation, and shear on a background turbulent flow by computing velocity correlations. Results. We obtain new theoretical scalings for velocity correlations that include the effect of horizontal shear. These scalings show an enhancement of turbulent motions, which would lead to a more efficient transport of chemicals and angular momentum, in better agreement with helio- and asteroseismic observations of rotation in the whole Hertzsprung-Russell diagram. Moreover, we propose a new choice for the non-linear time used in the relaxation approximation, which characterises the source of the turbulence. Conclusions. For the first time, we describe the effect of stratification, rotation, and vertical and horizontal shear on the anisotropy of turbulent transport in stellar radiative zones. The new prescriptions need to be implemented in stellar evolution calculations. To do so, it may be necessary to implement non-diffusive transport.


Introduction
The transport of angular momentum and chemical species plays a key role in stellar evolution. In particular, rotation induces transport both through large-scale motions such as meridional circulation and through turbulent, small-scale motions generated by hydrodynamical instabilities (Zahn 1983(Zahn , 1992Maeder & Zahn 1998;Maeder 2009). These multi-dimensional small-scale transport processes cannot be resolved in one-dimensional (1D) stellar evolution codes, so their secular effects have to be modelled in those codes.
In the last decades, helio-and asteroseismology have provided us with many constraints on the internal rotation of the Sun and distant stars. In the case of the Sun, helioseismic data bring to light an almost flat rotation profile in the solar radiative zone down to 0.2 R ⊙ , and a possibly faster core (Brown et al. 1989; Thompson et al. 2003;García et al. 2007;Fossat et al. 2017). Benomar et al. (2015) detected only weak differential rotation on a few solar-type stars. For a large number of subgiant and red giant stars, asteroseismic data show a weak core-surface rotation contrast (Beck et al. 2012;Mosser et al. 2012;Deheuvels et al. 2012Deheuvels et al. , 2014Deheuvels et al. , 2015Triana et al. 2017;Gehan et al. 2018). Weak differential rotation has also been detected in intermediate-mass and massive stars (Kurtz et al. 2014;Saio et al. 2015;Triana et al. 2015;Murphy et al. 2016;Van Reeth et al. 2016;Aerts et al. 2017;Van Reeth et al. 2018;Ouazzani et al. 2019). Finally, a strong transport of angular momentum is needed to explain data obtained for white dwarves (Suijs et al. 2008;Hermes et al. 2017) and neutron stars (Heger et al. 2005;Hirschi & Maeder 2010). Zahn (1992) introduced a formalism for turbulent transport accounting for both advection by the meridional circulation and the anisotropic turbulent transport induced by differential rotation, which is modelled thanks to vertical and horizontal turbulent diffusion coefficients. The anisotropy of the transport is due to the combined effect of rotation and stable stratification, which limit horizontal and vertical motions through the buoyancy force and the Coriolis acceleration, respectively. Those coefficients are based on phenomenological arguments and thus introduce large uncertainties in stellar evolution models. A lot of work has been done to improve the physics of the vertical diffusion coefficient by adding physical ingredients such as gradients of mean molecular weight (Maeder & Meynet 1996;Talon & Zahn 1997), and local numerical simulations have been performed to test existing models (Prat & Lignières 2013Garaud & Kulenthirarajah 2016) and investigate the effect of viscosity (Prat et al. 2016;Garaud et al. 2017;Gagnier & Garaud 2018). Current models predict less transport than what is needed to explain observations of internal rotation (Eggenberger et al. 2012;Ceillier et al. 2013;Marques et al. 2013;Cantiello et al. 2014;Ouazzani et al. 2019).
Concerning the horizontal diffusion coefficient, since the inital model of Zahn (1992) that was based on phenomenological arguments, other models based on energetic arguments (Maeder 2003) or results of an unstratified Taylor-Couette experiment ) have been proposed, but this coefficient is still poorly constrained, and its impact on stellar evolution is clearly overlooked. The first numerical simulations of the turbulent transport triggered by the instabilities of a horizontal shear in stellar conditions have also been performed (Cope et al. 2019;Garaud 2020). Finally, Park et al. (2020a,b) studied the linear in-A&A proofs: manuscript no. article stabilities of the horizontal differential rotation in stellar radiation zones as a function of stratification, rotation, and thermal diffusivity.
However, the widely used formalism of Zahn (1992) assumes that the turbulent transport of angular momentum is viscous, which is not necessarily the case. Indeed, anisotropic turbulent motions can generate non-zero turbulent fluxes of angular momentum even in the absence of differential rotation, which is a characteristic of non-diffusive transport. The existence of such non-diffusive turbulent transport associated with the Coriolis force, known as the Λ effect (Rüdiger 1989), has been confirmed by direct numerical simulations (see e.g. Käpylä 2019). Originally, this effect was introduced in the context of solar and stellar convective zones to explain the persistance of differential rotation despite the presence of a magnetic field, which was supposed to erase it. Kitchatinov & Brandenburg (2012, hereafter referred to as KB12) later investigated the effect of the stable stratification of stellar radiative zones on the Reynolds stress.
Based on this work, Mathis et al. (2018, hereafter referred to as M+18) proposed a new model of the horizontal transport induced by the vertical shear instability. In particular, M+18 added a vertical shear (due to radial differential rotation) to the formalism of KB12 and show that it has no significant effect on the anisotropy of turbulence. Besides, they used estimates of mean velocity correlations to propose new prescriptions for the horizontal turbulent diffusion coefficient of chemical elements. The implementation of these prescriptions in stellar evolution computations shows a slightly enhanced transport of angular momentum throughout the main sequence, but not enough to fit helioand asteroseismic observations. Furthermore, a major caveat of this implementation is that it assumed that the transport of angular momentum had the same (diffusive) nature as the transport of chemical elements, and thus completely ignored the Λ effect.
In the present work, we perform a new generalisation of the formalism of KB12 in the presence of a general (both vertical and horizontal) shear. It allows us to propose new prescriptions for the transport of both angular momentum and chemical elements that include the Λ effect.
In Sect 2, we derive spectral properties of the flow anisotropised by stratification, rotation, and shear. Then, in Sect. 3, we deduce scalings for velocity correlations as a function of stratification, rotation, and shear. We propose possible interpretations of the new scalings in Sect. 4. We explain how non-viscous turbulent transport of angular momentum should be implemented in stellar evolution codes to account for the Λ effect in Sect. 5. Finally, we conclude on the impact of these results on the modelling of turbulence in stellar evolution codes in Sect. 6.

Turbulent spectrum
In the Boussinesq approximation, which is justified here because turbulent motions are at length scales much smaller than the pressure scale height and velocity scales much smaller than the sound speed, the equations governing the flow are ∂s ∂t where V is the velocity, P ′ and ρ ′ are pressure and density fluctuations, respectively, ρ is the background density, g is the gravity vector, ν is the kinematic viscosity, f is a forcing term, s is the specific entropy, and κ is the thermal diffusivity. If one neglects the large-scale meridional circulation, velocity can be split into a mean rotation part and a fluctuating part: where r, θ, and ϕ are the standard spherical coordinates, (e r , e θ , e ϕ ) is the associated orthonormal basis, Ω is the rotation rate, and u is the velocity fluctuation. The continuity equation (1) implies that The left-hand side of the momentum equation (2) can be developed into where s is the distance to the rotation axis and e s the associated unit vector. The second term, which is the centrifugal acceleration, will be neglected in the following. Then, in the frame that locally rotates at the same rate as the fluid, Eq.
(2) becomes where S = r sin θ∇Ω is the shear rate. Similarly, the entropy can be split into a mean part s and a fluctuating one s ′ : Equation (3) yields Non-linear and dissipative terms are usually difficult to model analytically in turbulent flows. To simplify the problem, we use here a relaxation approximation, also known as τ approximation, which assumes that in a stationary steady state, the main effect of these terms is to tend the flow to relax with a given time constant τ. Equations (7) and (9) can then be approximated by For a complete treatment of thermal diffusion, we refer the reader to Park et al. (2020a,b). These two equations can be combined into a single one by expressing entropy fluctuations as a function of density fluctuations: where c P is the specific heat capacity at constant pressure. This yields where the Brunt-Väisälä frequency N is defined by Article number, page 2 of 7 V. Prat and S. Mathis: Anisotropic turbulent transport with horizontal shear in stellar radiative zones Using the Fourier decomposition the continuity equation (5) becomes and pressure fluctuations can be eliminated from Eq. (13): whereS = τS,k = k/k, σ =k·e z ,Ω = 2τΩ,Ñ = τN, µ =k r and u (0) = τ f . This way of writing the forcing term means that we assume that it comes from a pre-existing background turbulent flow. The previous equation can be written as a matrix relation where M is the matrix andS =S r e r +S θ e θ . When M is inversible, Eq. (18) implieŝ Using the fact thatû (0) also verifies the continuity equation, the previous equation reduces tô where is the determinant of M, η =k ·S/S , D is the matrix given by and δ i j and ε i jk are the Kronecker and Levi-Civita tensors, respectively. The spectrum of the background turbulence is assumed to be characterised by the correlation tensor where denotes a statistical average. In the continuity of Mathis et al. (2018), we assume a purely horizontal background turbulence, for which the Fourier transform of Q (0) i j iŝ where E(k) is the kinetic energy spectrum (see e.g. Davidson 2004) such that The turbulent spectrum of the sheared turbulence is deduced from that of the background turbulence using the relation mn .
After some algebra, this leads tô where

Mean correlations
Mean velocity correlations can be computed using the relation Stellar radiative zones are usually strongly stratified. It implies thatÑ ≫Ω (for example, the typical value of the ratio N/(2Ω) is around 500 for the Sun; see also André et al. 2018André et al. , 2019, for the evolution of this ratio along stellar evolution for low-and intermediate-mass stars),Ñ ≫S , andÑ ≫ 1. As a first approximation, where µ = cos α,k θ = sin α cos β andk ϕ = sin α sin β (with these notations, σ = cos θ cos α − sin θ sin α cos β). For the same reason, Using parity arguments and integrating over k (see Appendix A), this leads, when |S θ | < 2, to This term corresponds to the Reynolds stress in the meridional plane. Similarly, The first two terms correspond to the vertical and horizontal transport of angular momentum, respectively. The other two describe the horizontal turbulent velocity scales. We note that whenS θ tends towards zero, the scalings for u r u θ and u θ u ϕ (in Eqs. (34) and (36), respectively) vanish, while the others tend towards finite values. These values, as well as the zero limit of u θ u ϕ , are consistent with the results of KB12. In contrast, KB12 found a non-zero value for u r u θ : which is of lower order inÑ than the new scaling of Eq. (34).
Article number, page 3 of 7 A&A proofs: manuscript no. article The case of u 2 r is more complicated, since it contains three terms that come from and a priori have the same order of magnitude. The first term leads to which is consistent with the scaling derived in KB12 and M+18. The second term vanishes for parity reasons. The third term leads to When |S θ | ≪ 1, this reduces to which is small compared to u 2 r 1 . In total,

Interpretation
Interestingly, none of the terms computed in Sect. 3 explicitely depend on the radial shearS r , although it was included in the governing equations. In contrast, all share a common feature: they have a denominator in (1 −S 2 θ /4) 3/2 . This implies that when |S θ | gets close to 2, the turbulence is enhanced. Physically, this might be related to the limit of the Rayleigh-Taylor instability. The Rayleigh-Taylor stability criterion is Assuming that Ω > 0, this criterion can be rewritten cos θS θ > −(2Ω + sin θS r ).
When 2Ω + sin θS r > 0, the last inequality is verified for all positive values of cos θS θ , but only for negative values verifying By comparison with |S θ | < 2, this suggests that This expression is similar to τ = (2Ω+S r ) −1 , proposed by M+18. However, the presence of cos θ at the numerator means that τ would vanish at the equator, in which case the scalings derived in Sect. 3 are not valid there. Besides, the comparison with the Rayleigh-Taylor instability stands only in the special case where 2Ω + sin θS r > 0 and cos S θ < 0. Thus, the criterion |S θ | < 2 could be linked to some other hydrodynamical instability instead, such as the inertial or the inflectional instability studied by Park et al. (2020a,b). The form of the expression of u θ u ϕ in Eq. (36), which describes the horizontal transport of angular momentum, suggests that this flux term is viscous, with a turbulent viscosity coefficient equal to We note that this coefficient has no explicit dependence on the Brunt-Väisälä frequency or on the vertical shear. However, τ most likely depends on S r , and not necessarily on N. Thus, the viscosity coefficient probably depends on both S r and S θ . The dependence on u 2 (0) is problematic for the use in stellar evolution codes, since there is no easy way to reliably prescribe this quantity. One way to roughly estimate it is to assume that horizontal motions generated by the background turbulence are global. This can be written where R is the radial extent of the radiative zone. Introducing this expression in Eq. (48) yields In the case where the dynamics is dominated by the horizontal shear, assuming τ ∼ 1/S θ further leads to a horizontal viscosity coefficient that is proportional to the horizontal shear rate, which is consistent with the viscosity coefficient derived by Park et al. (2020b). We note that, in contrast, Garaud (2001) found that the horizontal shear instability can lead to an anti-diffusive horizontal transport of angular momentum. Although the expression in Eq. (48) is not directly usable in stellar evolution codes, it can be used to express other flux terms as a function of ν h . Thus, In the present work, we consider the case of a horizontally isotropic turbulence (i.e. with identical statistical properties in the latitudinal and azimuthal directions). This is rigourously justified for the radiative zone of slowly rotating stars, in which turbulent structures would have the shape of a horizontal pancake (see e.g. Fig. 1 of M+18). Because of that, there is no horizontal Λ effect. In rapidly rotating stars, turbulent structures are elongated along the rotation axis, so the horizontal isotropy is broken (see e.g. Rüdiger 1980) and a non-zero horizontal flux of angular momentum due to the Λ effect is present. However, in stellar radiative zones, provided that the Brunt-Väisälä frequency is still much larger than the rotation rate, the horizontal flux due to the Λ effect remains negligle compared to the leading diffusive term, which reduces to Eq. (48) in the horizontally symmetric case.
In contrast with the horizontal transport, the vertical transport may not be viscous, since the term u r u ϕ given in Eq. (35) is not explicitly proportional to S r . Instead, it is proportional to the rotation rate Ω, which is a characteristic of the Λ effect. However, when the dynamics is dominated by the vertical shear, one can assume that τ ∼ 1/S r , in which case Eq. (35) becomes Article number, page 4 of 7 V. Prat and S. Mathis: Anisotropic turbulent transport with horizontal shear in stellar radiative zones This then leads to a vertical turbulent viscosity coefficient which is negative. This corresponds to an anti-diffusive transport of angular momentum. Cope et al. (2019) and Garaud (2020) estimated the vertical turbulent diffusion coefficient associated with the horizontal shear instability using scaling laws but they did not provide any proof that the transport of chemical elements is indeed diffusive. In addition, they did not discuss the transport of angular momentum. Using the approximation of Eq. (49) leads to The same approximation can be used in the non-viscous case to express the flux from Eq. (35): Following KB12 and M+18, we estimate the turbulent diffusion coefficients using a mixing-length approximation D v ≃ τ u 2 r and D h ≃ τ u 2 h , where u 2 h = u 2 θ + u 2 ϕ . This allows us to estimate the ratio between vertical and horizontal diffusion coefficients: which has the same scaling with N as found by M+18. Besides, in the limit whereS θ ≪ 1, Eq. (57) reduces to which also has the same scaling with τ as found by M+18. Moreover, in contrast with M+18, where u θ u ϕ was zero and no horizontal turbulent viscosity coefficient could be directly derived, we can now compare horizontal viscosity and diffusion coefficients. We thus find that D h ≃ 4ν h . This suggests that the common assumption that they are equal may be qualitatively valid, but given the uncertainties introduced by the approximation we used, we cannot really conclude on the quantitative validity. In the vertical direction, we cannot directly compute the turbulent viscosity coefficient, so there is no evidence of such an equality.
We note that the term of horizontal transport is no longer present since it has a zero average. We showed in Sect. 4 that the horizontal transport can be seen as viscous. Therefore, the equation for the fluctuationsΩ is the same as in , for example: whereΩ(r, θ) = Ω 2 (r)[P 2 (cos θ) + 1/5] and

Conclusion
In the current paper, we investigated the impact of horizontal shear on the anisotropy of the transport in stably stratified, rotating stellar radiative zones. We were able to derive new scalings for mean velocity correlations that confirmed that the main effect is to significantly enhance turbulent motions. Thus, horizontal shear can be considered as an additional source of transport of angular momentum and chemical elements both in the horizontal and vertical directions.
In the regime where the stratification has a much stronger effect than rotation, we found that the horizontal transport of angular momentum is diffusive, while the vertical transport is dominated by the Λ effect. Additionally, the horizontal turbulent viscosity and diffusion coefficients are proportional to each other, which partly justifies the approximation that they are equal, wellspread in stellar evolution codes. The expressions of all the mean Reynolds stresses computed in the current paper depend on a property of the background turbulent flow, which is unknown a priori. Nevertheless, if one of the components is known by another way (for example using a phenomenological model for the horizontal turbulent viscosity coefficient), this allows one to compute fluxes of angular momentum. We finally show how the non-viscous transport of angular momentum predicted by our model should be implemented in stellar evolution codes to go beyond the current diffusive and viscous formalism.
Although the relaxation approximation has been validated in special cases by numerical simulations (see e.g. KB12), its validity in the configuration considered in the present work is uncertain. In particular, the dependence of the non-linear time τ on the A&A proofs: manuscript no. article physical parameters is undetermined. In M+18, three different physically motivated but phenomenological expressions were proposed. Here, we derived from more rigourous arguments a new expression for the non-linear turbulent time scale which is close to one of those proposed in M+18, τ = (2Ω + S r ) −1 , which corresponds to a turbulence dominated by the effects of rotation and vertical shear. Ultimately, direct numerical simulations should be performed to test the formalism and the new prescriptions in the presence of a general shear.
The results presented in this work have been obtained in the Boussinesq approximation, which neglects density fluctuations except in the buoyancy term. This treatment is valid as long as turbulent structures are small compared to the density scale height. When this condition is not satisfied, for example in regions with strong density gradients, similar calculations should be performed in the more general framework of the anelastic approximation, which accounts for density gradients, but filters pressure waves out.