A&A 425, 229-242 (2004)
DOI: 10.1051/0004-6361:20040278

Transport and mixing in the radiation zones of rotating stars

I. Hydrodynamical processes

S. Mathis - J.-P. Zahn

LUTH, Observatoire de Paris, 92195 Meudon, France

Received 16 February 2004 / Accepted 3 June 2004

Abstract
The purpose of this paper is to improve the modelization of the rotational mixing which occurs in stellar radiation zones, through the combined action of the thermally driven meridional circulation and of the turbulence generated by the shear of differential rotation. The turbulence is assumed to be anisotropic, due to the stratification, with stronger transport in the horizontal directions than in the vertical. The main difference with the former treatments by Zahn (1992, A&A, 265, 115) and Maeder & Zahn (1998, A&A, 334, 1000) is that we expand here the departures from spherical symmetry to higher order, and include explicitly the differential rotation in latitude, to first order. This allows us to treat simultaneously the bulk of a radiation zone and its tachocline(s). Moreover, we take fully into account the non-stationarity of the problem, which will enable us to tackle the rapid phases of evolution. The system of partial differential equations, which govern the transport of angular momentum, heat and chemical elements, is written in a form which makes it ready to implement in a stellar evolution code. Here the effect of a magnetic field is deliberately ignored; it will be included in forthcoming papers.

Key words: turbulence - stars: evolution - stars: rotation

1 Introduction - Motivation

In rotating stars, radiation zones are generally not in radiative equilibrium, as was shown by Von Zeipel (1924). Eddington (1925) and Vogt (1925) concluded that this thermal imbalance drives a large-scale meridional circulation, which was described in detail by Sweet (1950) for a given rotation law. Mestel (1953) discussed the role of composition inhomogeneities, which build up through the nuclear burning of hydrogen; he showed that they would slow down the circulation, and possibly prevent any mixing. Partly for this reason, rotational mixing was largely ignored in further work on stellar evolution. The next significant step was made by Busse (1982) who showed that angular momentum would be advected in such a way as to achieve a state of zero circulation, at least in an non-evolving, inviscid star with no stress applied to its surface.

The subject has been re-examined in Zahn (1992), who assumed that hydrodynamical instabilities due to the shear of differential rotation would enforce a rotation law with much larger variation in the vertical than in the horizontal direction. With such "shellular'' rotation, the problem reduces to one dimension (or rather 1.5 dimension, since departures from spherical symmetry are taken in account), and this renders it much more tractable. In particular one can then treat consistently the redistribution of angular momentum through the meridional circulation, and thus its feedback on the rotation profile, while preserving its advective nature. This makes it possible to prove that the boundary condition applied at the surface on the angular momentum flux plays a crucial role. In a star which loses angular momentum through a wind, the meridional circulation is enforced by the requirement of transporting this momentum to the surface. On the contrary, when the star loses no angular momentum, it relaxes to a state where the advection of momentum is balanced by its transport through turbulent motions. Thereafter this treatment was extended to inhomogeneous stars, with a general equation of state, and partly adapted to phases of fast stellar evolution (Maeder & Zahn 1998).

The first evolutionary sequence calculated according to these prescriptions was that of a 9 $M_\odot$ star (Talon et al. 1997); it was followed by extensive grids of models built by Meynet & Maeder (2000), who refined the prescription for the loss of angular momentum by taking into account the latitudinal variation of the wind. For massive stars, these models are in better agreement with observational data than models without rotational mixing: they account for the observed surface composition, and they predict the correct ratio between red and blue supergiants observed in open clusters and in the Magellanic clouds (Maeder & Meynet 2001). They also explain the absence of light elements depletion on the blue side of the Li dip, as shown by Talon & Charbonnel (1998). But they fail to correctly reproduce the flat rotation profile in the radiative interior of solar-type stars (Mathias & Zahn 1997). In such stars, which are slow rotators, other physical processes must therefore contribute to the transport of angular momentum; the most plausible candidates are either magnetic stresses, as advocated already by Mestel (1953), or internal gravity waves emitted at the interface with convection zones (Talon et al. 2002; Talon & Charbonnel 2003).

The purpose of the present paper is to improve the modeling of rotational mixing, in particular during the phases of fast stellar evolution, when steep vertical gradients develop, both of angular velocity and of chemical composition. To achieve this, our aim is threefold:

• we extend the linear expansion of all variables to higher order spherical harmonics, in order to capture the tachocline circulation which is induced in its vicinity by the differential rotation of a convection zone;
• for the same reason, we treat explicitly the departure from shellular rotation, in the linear approximation;
• we retain all time derivatives, except those describing the relaxation to hydrostatic and geostrophic balance, in order to better resolve the phases of fast stellar evolution.

In forthcoming papers, we shall introduce a magnetic field, and examine its impact on the thermally driven circulation.

2 Main assumptions

Our treatment is based on the assumption that the rotating star is only weakly two-dimensional, and this for two reasons. The first is that the rotation rate is sufficiently moderate to allow the centrifugal force to be considered as a perturbation, compared to the gravitational field. The second reason rests on the conjecture that the differential rotation induced by the meridional circulation gives rise to turbulent motions which are strongly anisotropic due to the stable stratification, with much stronger transport in the horizontal directions than in the vertical: i.e.  $\nu_{\rm v} \ll \nu_{\rm h}$ and  $D_{\rm v} \ll D_{\rm h}$ respectively for the turbulent viscosity and diffusivity. Such anisotropic turbulence is observed in the Earth's atmosphere and oceans; in a star we expect it to smooth the horizontal variations of angular velocity and of chemical composition, a property we shall invoke to discard certain non-linear terms. The prescriptions to be used for these turbulent diffusivities are discussed and updated in a companion paper (Mathis et al. 2004).

We thus consider an axisymmetric star, and assume that the horizontal variations of all quantities are small and smooth enough to allow their linearization and their expansion in a modest number of spherical harmonics  $P_{l}(\cos\theta)$. As reference surface, we chose either the sphere or the isobar, and write all scalar quantities either as

$$X(r,\theta)=X_{0}(r)+\sum_{l}\widehat{X}_{l}(r)P_{l}\left(\cos\theta\right)$$ (1)

or

$$X(P,\theta)=\overline{X}(P)+\sum_{l>0}\widetilde{X}_{l}(P)P_{l}\left(\cos\theta\right).$$ (2)

Let us establish the relation between those two modal expansions. We expand the pressure around the sphere as:

$$P(r,\theta)=P_{0}(r)+\sum_{l>0}\widehat{P}_{l}(r)P_{l}\left(\cos\theta\right),$$ (3)

and introduce the radial coordinate of the isobar

$$r_P(r, \theta) = r + \sum_{l>0}\xi_{l}(r) P_{l}(\cos\theta),$$ (4)

where r is the mean value of the radius of an isobar. Taking the Taylor expansion of P to first order, we have:

$$P\left(r+\sum_{l>0}\xi_{l}(r)P_{l}(\cos\theta),\theta\right)=... ...}{ {{\rm d} r}}\right) \sum_{l>0}\xi_{l}(r)P_{l}(\cos\theta),$$ (5)

and since by definition the pressure is constant on the isobar, we conclude that

$$\xi_{l}(r)=-\frac{\widehat{P}_{l}(r)} {{{\rm d} P_{0}/{\rm d}r}}\cdot$$ (6)

If we apply the same procedure to any other variable X, we get

$$X\left(r+\sum_{l>0} \xi_{l}(r)P_{l}(\cos\theta), \theta\right... ...\rm d}P_{0}}\right)\widehat{P}_{l}(r)\right]P_{l}(\cos\theta)$$ (7)

and therefore

$$\tilde{X}_{l}({r})=\widehat{X}_{l}(r)- \left( \frac{{\rm d} X_{0}} {{\rm d} P_{0}} \right) \widehat{P}_{l}(r).$$ (8)

This relation will serve below in Sect. 5.1 to calculate the effective gravity.

3 Transport of angular momentum

We start from the momentum equation

$$\rho\left[ \partial_{t}\vec V+\left( \vec V \cdot \vec \nabla... ...- \rho\vec\nabla\phi+\vec\nabla\cdot \vert\vert\tau\vert\vert$$ (9)

where $\rho$ is the density, $\phi$ the gravitational potential, and $\vert\vert\tau\vert\vert$ represents the turbulent stresses. The macroscopic velocity field  $\vec V$ is the sum of a zonal flow with angular velocity  $\Omega(r,\theta)$ and a meridional flow  $\vec{\mathcal{U}}(r, \theta)$:

$$\vec V=r\sin\theta ~ \Omega(r,\theta) ~ {\vec e}_{\varphi}+\vec{\mathcal{U}}(r, \theta).$$ (10)

The latter can be split into a spherically symmetric part, which represents the contraction or dilation of the star during its evolution, plus the meridional circulation

$$\vec{\mathcal{U}}=\dot{r} ~ {\vec e}_{r}+\vec{\mathcal{U}}_{M}(r, \theta),$$ (11)

which we expand in spherical functions

$$\vec{\mathcal{U}}_{M}=\sum_{l>0}\left[U_{l}(r)P_{l}\left(\cos... ... d}P_{l}(\cos\theta)}{{\rm d}\theta}{\vec e}_{\theta}\right].$$ (12)

Using the continuity equation in the anelastic approximation, i.e.  $\vec{\nabla} \cdot (\rho~\vec {\mathcal{U}}_{M})=0$, we obtain the following relation between U_l(r) and V_l(r):

$$V_{l}=\frac{1}{l(l+1)\rho r}{{\rm d} \over {\rm d}r} \left(\rho r^2U_{l}\right).$$ (13)

Making again use of the continuity equation, the azimuthal component of (9) can be written as an advection/diffusion equation for the angular momentum:

$$\rho {{\rm d} \over {\rm d}t} \left(r^{2} \sin^{2} \theta ~ \... ...nu_{\rm h} \sin^{3} \theta ~ \partial_{\theta} \Omega \right),$$ (14)

where as in Zahn (1992) we assume that the effect of the turbulent stresses on the large scale flow are adequately described by an anisotropic eddy viscosity, whose components are  $\nu_{\rm v}$ and  $\nu_{\rm h}$ respectively in the vertical and horizontal directions. As discussed in Mathis et al. (2004), they act to reduce the cause of turbulence, as observed in laboratory experiments of rotating flows, namely here the vertical and horizontal gradients of angular velocity. This explains why the respective fluxes of angular momentum contain only these gradients, in contrast with the treatment of Kippenhahn (1963), who considered the effect of an imposed anisotropic turbulence, due to thermal convection; in that case the rotation becomes non-uniform, as is actually observed in the solar convection zone. Note the introduction of the Lagrangian time-derivative  ${{\rm d} / {\rm d}t}$, meaning that the radial coordinate r is the mean radius of the layer (either sphere or isobar) which encloses the mass M_r, with  ${\rm d}M_{r}=4\pi\rho r^{2}{\rm d}r$.

The form of Eq. (14) incites to expand the angular velocity as

$$\Omega(r,\theta)=\overline{\Omega}(r)+ \widehat{\Omega}(r,\th... ... \overline{\Omega}(r)+ \sum_{l>0}\Omega_{l}(r) ~ Q_l(\theta),$$ (15)

with the horizontal average being defined as

$$\overline{\Omega}(r)=\frac{\int_{0}^{\pi}\Omega(r,\theta)\sin... ...rm d} \theta}{\int_{0}^{\pi}\sin^{3}\theta ~ {\rm d} \theta},$$ (16)

and where the horizontal functions  $Q_l(\theta)$ satisfy the orthogonality condition

$$\frac{\int_{0}^{\pi} Q_l(\theta) \sin^{3}\theta ~ {\rm d} \theta} {\int_{0}^{\pi}\sin^{3}\theta ~ {\rm d} \theta} =0.$$ (17)

These horizontal functions are readily identified:

$$Q_l(\theta)= P_{l}\left(\cos\theta\right) - I_l \quad \hbox{w... ...{\rm d} \theta} = \delta_{l,0} - {1 \over 5} \delta_{l,2} ~ ;$$ (18)

except for l=2, these functions Q_l reduce to the Legendre polynomials.

3.1 Vertical transport of angular momentum

Taking the horizontal average of Eq. (14), and using the assumption that  $\overline{\Omega} \gg \Omega_{l}$, we obtain the following vertical advection/diffusion equation for the mean angular velocity  $\overline{\Omega}$:

$$\rho {{\rm d} \over {\rm d}t} (r^2\overline{\Omega})=\frac{1}... ...\left(\rho\nu_{\rm v}r^4\partial_{r}\overline{\Omega}\right).$$ (19)

Note that only the l=2 component of the circulation is able to advect a net amount of angular momentum; the higher order components of  $\vec{\mathcal{U}}_{M}$ (for instance those induced in its tachocline by a differentially rotating convection zone) do not contribute to the vertical transport of angular momentum, as was explained in Spiegel & Zahn (1992).

3.2 Horizontal transport of angular momentum

We establish the equation governing the horizontal transport of angular momentum by multiplying Eq. (19) through  $\sin^{2}\theta$ and subtracting it from the original form (14):

$$\rho {{\rm d} \over {\rm d}t} \left(r^2\sin^2\theta ~ \wideha... ...m h}\sin^3\theta ~ \partial_{\theta}\widehat{\Omega}\right) ;$$ (20)

in the advection term we have again neglected the fluctuation  $\widehat{\Omega}$ compared to the mean  $\overline{\Omega}$.

The next step is to replace  $\widehat{\Omega}(r, \theta)$ by its expansion (15) in the horizontal functions  $Q_l(\theta)$. For l=2, the equation separates neatly into

$$\rho {{\rm d} \over {\rm d}t} \left(r^2 \Omega_{2}\right)-2\r... ...r^4\partial_{r}\Omega_{2}\right)-10\rho\nu_{\rm h}\Omega_{2},$$ (21)

$$\hbox{with} \quad V_{2}=\frac{1}{6\rho r}\frac{{\rm d}}{{\rm ... ...{{\rm d}\ln\left(r^2\overline{\Omega}\right)}{{\rm d}\ln r},$$ (22)

which we can simplify by assuming that the turbulent transport is much more efficient in the horizontal than in the vertical direction, i.e.   $\nu_{\rm v} \ll \nu_{\rm h}$:

$$\rho {{\rm d} \over {\rm d}t} \left(r^2 \Omega_{2}\right)-2\r... ...left[2V_{2}-\alpha U_{2}\right]=-10\rho\nu_{\rm h}\Omega_{2}.$$ (23)

In the asymptotic regime   $t \gg r^2/\nu_{\rm h}$, we retrieve the equation for  $\Omega_{2}$ given by Zahn (1992):

$$\nu_{\rm h}\Omega_{2}(r)=\frac{1}{5}r\left[2V_{2}(r)-\alpha(r)U_{2}(r)\right]\overline{\Omega}(r).$$ (24)

For l>2 the situation is more intricate, with couplings between terms of different l. However, when making the three following assumptions: (i) stationarity ( $t \gg r^2/\nu_{\rm h}$); (ii) strong anisotropy of the turbulent diffusion ( $\nu_{\rm v} \ll \nu_{\rm h}$); (iii) vertical advection negligible compared to the horizontal one (thin layer approximation), the l-components separate, and one recovers the result given by Spiegel & Zahn (1992) in their treatment of the solar tachocline:

$$\nu_{\rm h} \Omega_{l}(r)= \overline{\Omega} rV_{l}(r) ,$$ (25)

namely that horizontal diffusion balances horizontal advection.

4 Transport of chemicals

The transport of chemicals species is governed by an advection/diffusion equation similar to that for the transport of angular momentum:

$$\rho \partial_{t}c_{i}+ \rho~ \vec{\mathcal{U}} \cdot \vec\na... ...left(\rho D_{\rm h}\sin\theta~ \partial_{\theta}c_{i}\right),$$ (26)

where c_i is the concentration of a given element, and where we assume again that the turbulent diffusivity is anisotropic. For simplicity, we have not included here the term representing the production or destruction of the species due to nuclear reactions, which can be easily added. We follow the same procedure as for the angular momentum: we split the concentration in its horizontal average and its fluctuation

$$c_{i}(r, \theta)=\overline{c_{i}}(r)+c_{i}'(r,\theta),$$ (27)

and explicit the velocity field in

$$\vec{\mathcal{U}} = {\vec e}_r \dot{r} + \vec{\mathcal{U}}_M + {\vec e}_r U_i^{\rm diff}$$ (28)

where, to the radial displacement accompanying the evolution of the star and to the meridional circulation, we add the microscopic diffusion velocity  $U_i^{\rm diff}$, which depends on the considered species. We insert these expressions into Eq. (26), which leads us to

$$\rho {{\rm d} \over {\rm d}t} \left( \overline{c_{i}}+{c_{i}'... ...t(\rho D_{\rm h}\sin\theta ~\partial_{\theta}{c}_{i}'\right).$$ (29)

Then we take its horizontal average, making use of the anelastic continuity equation   $\vec{\nabla} \cdot (\rho~\vec {\mathcal{U}}_{M})=0$; we thus obtain the equation governing the evolution of the mean concentration  $\overline{c_{i}}$:

$$\rho {{\rm d} \over {\rm d} t } \overline{c_{i}}+\frac{1}{r^2... ...} \theta \over \int_0^\pi \sin \theta ~ {\rm d} \theta} \cdot$$ (30)

The subtraction of Eq. (30) from Eq. (29) yields the following equation for the fluctuations of the concentration:

$$\rho {{\rm d} \over {\rm d}t} {c_{i}'}+\rho~ \mathcal{U}_{M, ... ...eft(\rho D_{\rm h}\sin\theta\partial_{\theta}{c_{i}'}\right).$$ (31)

Expanding c'_i and  $\mathcal{U}_{M, r}$ in Legendre polynomials:

$$c_{i}'=\sum_{l>0}\widetilde{c_{i}}_{,l}(r)P_{l}(\cos\theta) \... ...quad \mathcal{U}_{M, r}=\sum_{l>0}U_{l}(r)P_{l}(\cos\theta),$$ (32)

and assuming that the gradient of the concentration fluctuations is negligible compared to that of the mean concentration, i.e.  $\vert\vec\nabla c'_i\vert \ll \partial_r \overline{c_i}$, this equation can be recast into

$$\rho\ {{\rm d} \over {\rm d}t} \widetilde{c_{i}}_{,l}+\rho U_... ...ight)-\frac{l(l+1)}{r^2}\rho D_{\rm h}\widetilde{c_{i}}_{,l}.$$ (33)

As in Chaboyer & Zahn (1992), we further assume that the horizontal diffusivity is much stronger than the vertical one, more precisely that  $D_{\rm h} \gg D_{\rm v} ({l_{\rm h}}/{l_{\rm v}})^2$, where  $l_{\rm h} (l_{\rm v})$ is the distance over which c_i' changes significantly in the horizontal (vertical) direction, and we finally obtain:

$${{\rm d} \over {\rm d}t} \widetilde{c_{i}}_{,l}+U_{l}\partial... ...line{c_i}=-\frac{l(l+1)}{r^2}D_{\rm h}\widetilde{c_{i}}_{,l}.$$ (34)

For  $t \gg r^2/D_{\rm h}$ there is an asymptotic state where

$$\widetilde{c_{i}}_{,l}=-\frac{r^2}{l(l+1)D_{\rm h}}U_{l}\partial_{r}\overline{c_i}$$ (35)

and then

$$\rho {{\rm d} \over {\rm d}t} \overline{c_i} + \frac{1}{r^2}\... ...}\frac{r^2\left(U_{l}\right)^{2}}{l(l+1)(2l+1)D_{\rm h}}\cdot$$ (36)

However, in order to better resolve the fast evolution phases, one may prefer to retain the time derivative in (34), and solve the equation for the mean concentration (30) in the form:

$$\rho {{\rm d} \over {\rm d} t } \overline{c_{i}}+\frac{1}{r^2... ...r}\left[r^2\rho D_{\rm v}\partial_{r}\overline{c_{i}}\right].$$ (37)

Finally, let us recall how the molecular weight $\mu$ is related to the concentrations c_i:

$$\frac{1}{\mu}=\sum_{i}\frac{1+Z_{i}}{A_{i}}c_{i}$$ (38)

where A_i is the mass number (protons+neutrons) and Z_i the number of electrons. From Eq. (34) we can thus derive the following advection/diffusion equation for the horizontal fluctuation of the molecular weight  $\Lambda_{l}=\widetilde{\mu_{l}}/\overline{\mu}$:

$${{\rm d} \over {\rm d}t} \Lambda_{l} - {{\rm d} \ln \overline... ..._{\rm p}}\nabla_{\mu}-\frac{l(l+1)}{r^2}D_{\rm h}\Lambda_{l},$$ (39)

where we have introduced the pressure scale-height  $H_{\rm p}=\vert{\rm d} r/{\rm d} \ln P \vert$ and the logarithmic gradient  $\nabla_{\mu}={\rm d} \ln\overline{\mu}/ {\rm d}\ln P$. This equation will play a key role in the derivation of the meridional circulation, which we shall discuss in Sect. 6.

5 Structural properties of the differentially rotating star

5.1 Baroclinic relation

In our case where we take a non-uniform and a non-cylindrical rotation law, the centrifugal force does not derive from a potential. We are in the baroclinic configuration, where the isobars and the surfaces of constant density do not coincide anymore. To find how the density varies on an isobar, we start from the hydrostatic equation:

$$\frac{1}{\rho}\vec\nabla P=\vec g= - \vec\nabla \phi+\vec{\ma... ...{1}{2}{\Omega}^{2}\vec\nabla\left(r^{2}\sin^{2}\theta\right)$$ (40)

where $\phi$ is the gravitational potential, and where the local effective gravity $\vec g$ includes the centrifugal force  $\vec{\mathcal F}_{\mathcal{C}}$. Taking the curl of this equation, we get

$$-\frac{1}{{\rho}^{2}}\vec\nabla\rho\wedge\vec\nabla P=-\frac{... ...}{2}\vec\nabla (\Omega)^{2}\wedge\vec\nabla(r\sin\theta)^{2},$$ (41)

which to first order takes the form

$$- \frac{\vec\nabla {\rho}^{'}\wedge\vec g}{\overline{\rho}}=\... ...l_{\theta}(\Omega^{2})\sin^{2}\theta\right]\vec{e}_{\varphi},$$ (42)

with $\rho'$ being the variation of the density on the isobar, which we write as

$$\rho^{'}(r, \theta)=\sum_{l>0} \widetilde{\rho_{l}}(r)P_{l}(\cos\theta).$$ (43)

Likewise we expand $\Omega^2$ in spherical harmonics; since (cf. Eqs. (15), (18))

$$\Omega(r,\theta)=\overline{\Omega}(r)+\sum_{l>0}\Omega_{l}(r)\left(P_{l}(\cos\theta)-I_{l}\right)$$ (44)

we have

$$\Omega^{2}(r,\theta)=\left[{\overline{\Omega}}^{2}-2\overline... ...ight]+2\overline{\Omega}\sum_{l>0}\Omega_{l}P_{l}(\cos\theta)$$ (45)

where we kept only the terms linear in $\Omega_l$. We expect the largest departures from shellular rotation to occur in the tachoclines, where they are enforced by the differential rotation of the adjacent convection zone; judging by the solar case, such departures should be rather small, of the order of 1/10, thus justifying the linear approximation.

Inserting these expansions in Eq. (42), we reach after some algebra the following expression for the modal amplitudes of the relative density fluctuation on an isobar

$$\frac{\widetilde{\rho}_{l}(r)}{\overline{\rho}}=\frac{r}{\overline{g}}\mathcal{D}_{l}(r),$$ (46)

where  $\overline{g}$ is the horizontal average, on the isobar, of the modulus of $\vec g$, and

 

$$\mathcal{D}_{l}(r) \mathcal{N}_{l}^{0}\left\{r\partial_{r}\left[\overline{\Omega}^{2... ...(r)\Omega_{2}(r)I_{2} \right] \frac{1}{3\mathcal{N}_{2}^{0}}\delta_{l,2}\right.$$ =

$$+2r\sum_{s>0}\partial_{r}\left(\overline{\Omega}(r)\Omega_{s}(r)\... ..._{s}^{0} \left(-C_{s+1}^{0}\delta_{l,s}+D_{s+1}^{0}\delta_{l,s+2}\right)\right]$$

$$-{\left. 2\overline{\Omega}(r)\sum_{s>0}\frac{\Omega_{s}(r)}{\mat... ...left(-C_{s-1}^{0}\delta_{l,s-2}+D_{s-1}^{0}\delta_{l,s}\right)\right]\right\}}.$$ (47)

All numerical coefficients involved ( $A^0_l, B^0_l, C^0_l, D^0_l, G^0_l, H^0_l, \mathcal{N}^0_l$) are given in Appendix A. We shall display here the result for l=2, 4, where we keep only the first term  $\Omega_{2}$ of the expansion of $\Omega$:

$$\mathcal{D}_{2} = \frac{1}{3}\left[r\partial_{r}\overline{\Om... ...mega_{2}\right)\right]+\frac{8}{7}\overline{\Omega}\Omega_{2}$$ (48)

$$\mathcal{D}_{4} = \frac{6}{35} \left[r\partial_{r} \left(\ove... ...a}\Omega_{2}\right) -2 \overline{\Omega} \Omega_{2}\right].$$ (49)

For  $\Omega_2=0$ we recover the expression given in Zahn (1992):  $\widetilde{\rho}_{2}/\overline\rho = (r^2/ 3\overline g) \partial_r \overline{\Omega}^{2}$.

This baroclinic Eqs. (46)-(47) plays a key role in linking the density fluctuation on an isobar with the rotation profile. It will allow us to close the system formed by the equation for the transport of angular momentum, that for the transport of the chemical species and that for the transport of heat, which we shall establish in Sect. 6.

5.2 Effective gravity

Next we examine the redistribution of masses by the centrifugal force and its effect on the gravity inside the star. As all other scalars, we expand the gravity as:

$$g(P,\theta)=\overline{g}(P)+\sum_{l>0}\widetilde{g}_{l}(P)P_{l}(\cos\theta).$$ (50)

The goal of this section is to determine the amplitude of the fluctuation on an isobar,  $\widetilde{g}_{l}$. The first step will be to calculate the perturbation of the gravitational potential $\phi$ on the sphere of radius r, and thus to calculate the functions  $\widehat \phi_{l}(r)$, where we have expanded $\phi$ as:

$$\phi(r, \theta)=\phi^{(0)}(r)+\phi^{(1)}(r,\theta)=\phi_{0}(r... ...a) \quad \hbox{with} \quad \phi_{0}(r)=-\frac{GM(r)}{r}\cdot$$ (51)

We follow the method of linearization developed in Sweet (1950), and expand the pressure and the density around the sphere in the same way as $\phi$:

$$P(r, \theta)=P^{(0)}(r)+P^{(1)}(r,\theta)=P_{0}(r)+\sum_{l>0}\widehat{P}_{l}(r)P_{l}(\cos\theta)$$ (52)

$$\rho(r, \theta)=\rho^{(0)}(r)+\rho^{(1)}(r,\theta)=\rho_{0}(r)+\sum_{l>0}\widehat{\rho}_{l}(r)P_{l}(\cos\theta).$$ (53)

Then, we take the hydrostatic equation (using the classical definition of a potential, contrary to Zahn 1992):

$$\frac{\vec\nabla P}{\rho}= - \vec\nabla\phi+\vec{\mathcal F}_... ...athcal{C}} =\frac{1}{2}\Omega^{2}\vec\nabla(r^2\sin^2\theta),$$ (54)

which we expand to first order:

$$\vec\nabla P^{(0)}= - \rho^{(0)}\vec\nabla\phi^{(0)}$$ (55)

$$\vec\nabla P^{(1)}= - \rho^{(0)}\vec\nabla\phi^{(1)} - \rho^{... ...c\nabla\phi^{(0)}+ \rho^{(0)}\vec{\mathcal F}_{\mathcal{C}}.$$ (56)

We eliminate the pressure fluctuation by taking the curl of the latter equation, which gives us:

$$\partial_{\theta}\rho^{(1)} = \frac{\partial_{r}\rho^{(0)}}{g... ...(r\rho^{(0)}{\mathcal F}_{{\mathcal C},\theta}\right)\right].$$ (57)

Next we insert the modal expansion of  $\rho^{(1)}$ and those of the components of the centrifugal force:

$${\mathcal F}_{{\mathcal C}, r}(r,\theta)=\sum_{l}a_{l}(r)P_{l... ...,\theta)=-\sum_{l}b_{l}(r)\partial_\theta P_{l}(\cos\theta) ;$$ (58)

after integration in $\theta$, this yields the modal amplitude of the density fluctuation over the sphere:

$$\widehat{\rho}_{l}(r) = \frac{1}{g_{0}} \left[ \frac{{\rm d}... ...}+\frac{{\rm d}}{{\rm d}r}\left(r\rho_{0}b_{l}\right)\right].$$ (59)

It remains to insert this expression in the perturbed Poisson equation  $\nabla^2 \widehat{\phi}_{l} = 4 \pi G \widehat{\rho}_{l}$ to retrieve Sweet's result:

$$\frac{1}{r}\frac{{\rm d}^2}{{\rm d}r^2}\left(r\widehat{\phi}_... ...}+\frac{{\rm d}}{{\rm d}r}\left(r\rho_{0}b_{l}\right)\right].$$ (60)

The functions a_l(r) and b_l(r) are given by

 

$$\frac{2}{3}r\left[\overline{\Omega}^{2}(r)-2\overline{\Omega}(r)\Omega_{2}(r)I_{2}\right]\left(\delta_{l,0}-\delta_{l,2}\right)$$                              a_l(r) = (61)

$$+ 2r\overline{\Omega}(r)\sum_{s>0}\Omega_{s} (r)\left[\frac{1}{\m... ...0}\delta_{l,s+2}-H_{s+1}^{0}\mathcal{N}_{s}^{0}\delta_{l,s}\right)\Big\}\right]$$

$$\frac{1}{3}r\left[\overline{\Omega}^{2}(r)-2\overline{\Omega}(r)\Omega_{2}(r)I_{2}\right]\delta_{l,2}$$                              b_l(r) = (62)

$$+2r\overline{\Omega}(r)\sum_{s>0}\frac{\Omega_{s}ñ(r)}{\mathcal{N... ...s}^{0}\delta_{l,s}+D_{s+1}^{0}\mathcal{N}_{s+2}^{0}\delta_{l,s+2}\right)\Big\},$$

where we have used the expansion (45) of $\Omega^2$. The numerical coefficients are given in Appendix A; here we shall quote only the results for l=2, 4, with the expansion of $\Omega$ limited to the term $\Omega_{2}$:

$$a_{0}=\frac{2}{3}r \overline{\Omega}^{2} \qquad \qquad \qquad \qquad b_{0}$$ = 0 (63)

$$a_{2}=-\frac{2}{3} r\overline{\Omega}^{2} + \frac{24}{35}r\overline{\Omega}\Omega_{2} \qquad\;~ b_{2} \frac{1}{3}r \overline{\Omega}^{2}+\frac{8}{35}r\overline{\Omega}\Omega_{2}$$ = (64)

$$a_{4}=-\frac{24}{35}r\overline{\Omega}\Omega_{2} \qquad \qquad\qquad b_{4} \frac{6}{35}r\overline{\Omega}\Omega_{2} \cdot$$ = (65)

With this multipolar expansion of the gravitational potential, we are ready to calculate the horizontal fluctuation of the effective gravity:

$$\vec g= - \vec\nabla\phi+ \vec{\mathcal F}_{\mathcal{C}}= - \... ...{l}(\cos\theta)\right)\right]+\vec{\mathcal F}_{\mathcal{C}}.$$ (66)

Taking the scalar product of $\vec g$ by itself, and retaining only the first order terms, we get the following expansion for g² over the sphere of radius r:

$$g^{2}(r,\theta)=\left[\left(\frac{{\rm d}\phi_{0}(r)}{{\rm d}... ...t{\phi}_{l}(r)}{{\rm d}r} - a_{l}(r)\right)P_{l}(\cos\theta).$$ (67)

Comparing with

$$g^2(r,\theta)=\left[g_{0}(r)+\sum_{l\neq 0}\widehat{g}_{l}(r)... ...um_{l\neq 0}\widehat{g}_{l}(r)P_{l}\left(\cos\theta\right) ~,$$ (68)

we identify:

$$\widehat{g}_{l}(r)= \frac{{\rm d}\widehat{\phi}_{l}(r)}{{\rm d} r} - a_{l}(r).$$ (69)

Making use of (8) we get likewise for the variation along an isobar:

$$\widetilde{g}_{l}(P)=\frac{{\rm d} g_{0}}{{\rm d} r}\frac{1}{... ...\right) + \frac{{\rm d}\widehat{\phi}_{l}}{{\rm d}r} - a_{l}.$$ (70)

It remains to replace  $\widehat{P}_{l}$ by its expression drawn from the $\theta$-component of the hydrostatic equation

$$\frac{\widehat{P}_{l}}{\rho_{0}}= - \widehat{\phi}_{l}-r b_{l}$$ (71)

to finally obtain the variation of the effective gravity along an isobar:

$$\frac{\widetilde{g}_{l}}{\overline{g}}=-\left[\frac{{\rm d} g... ...m d}}{{\rm d}r}\left(\frac{\widehat{\phi}_{l}}{g_{0}}\right).$$ (72)

Keeping only the $\Omega_{2}$ term and carrying the expansion to l=4, we have:

$$\frac{\widetilde{g}_{2}}{\overline{g}}=-\frac{{\rm d}g_{0}}{{... ...rm d}}{{\rm d}r}\left(\frac{\widehat{\phi}_{2}}{g_{0}}\right)$$ (73)

$$\frac{\widetilde{g}_{4}}{\overline{g}}=-\frac{{\rm d}g_{0}}{{... ...}{{\rm d}r}\left(\frac{\widehat{\phi}_{4}}{g_{0}}\right)\cdot$$ (74)

For strictly shellular rotation ( $\Omega_2=0$) we retrieve the expression of Zahn (1992) (after changing the sign of  $\widehat{\phi}_{2}$):

$$\frac{\widetilde{g}_{2}}{\overline{g}}= {1 \over 3} \overline... ...}{{\rm d}r}\left(\frac{\widehat{\phi}_{2}}{g_{0}}\right)\cdot$$ (75)

The two terms represent respectively the contribution of the centrifugal force and that of the modified mass distribution. The latter is obtained by integrating the perturbed Poisson Eq. (60) with the appropriate boundary conditions (see Sect. 7).

6 Thermal imbalance and transport of heat

Following closely the procedure described in Zahn (1992), we now undertake to link the meridional circulation with the rotation profile and the inhomogeneities of chemical composition. We start from the equation stating the conservation of thermal energy:

$$\rho T \left[ {\partial S \over \partial t} + \vec{\mathcal{U... ...\nabla T\right)+\rho\epsilon-\vec\nabla\cdot{\vec{F}}_{\rm h}$$ (76)

where S is the entropy per unit mass, $\chi$ the thermal conductivity and $\epsilon$ the nuclear energy production rate per unit mass. As Maeder & Zahn (1998), we include the flux  ${\vec F}_{\rm h}$ carried by horizontal turbulence.

We project this equation on the base of spherical harmonics, starting with the left hand side:

$$\rho T \left[ {\partial S \over \partial t} + \vec{\mathcal{U... ...ine{S} \over {\rm d}r} \right] P_{l}\left(\cos\theta\right) ,$$ (77)

where we have used again the expression of the velocity field (28)  $\vec{\mathcal{U}} = \widehat{\vec e}_r \dot{r} + \vec{\mathcal{U}}_M~$, together with the modal expansion of the meridional flow (12). Recognizing in the horizontal average the production of energy related with the contraction or expansion of the star during its evolution, i.e.  $\overline{T} {{\rm d}\overline{S}}/{{\rm d}t}=-\epsilon_{\rm grav}$, we rewrite Eq. (76) as

$$\overline{\rho}\overline{T}\sum_{l>0} \left[\frac{{\rm d}\wid... ...verline{\epsilon}_{\rm grav}-\vec\nabla\cdot{\vec F}_{\rm h},$$ (78)

thus introducing the Lagrangian time derivative. Note that by construction the right hand side has now a zero horizontal average. Then we expand the temperature in

$$T(P,\theta)=\overline{T}(P)+\sum_{l>0}\widetilde{T}_{l}(P)P_{l}(\cos\theta),$$ (79)

from which we deduce the gradient

$$\vec\nabla T=\rho\left[\frac{{\rm d}\overline{T}}{{\rm d}P} +... ...o}+ \sum_{l>0}\widetilde{T}_{l}\vec\nabla P_{l}(\cos\theta).$$ (80)

Multiplying by the conductivity $\chi$ and taking the divergence, we obtain

$$\vec\nabla\cdot\left(\chi\vec\nabla T\right) \rho\chi\left[\frac{{\rm d}\overline{T}}{{\rm d}P}+\sum_{l>0}\fra... ...{T}_{l}}{{\rm d}P}P_{l}(\cos\theta)\right]\right)\cdot\frac{\vec\nabla P}{\rho}$$ =

$$+\sum_{l>0}\vec\nabla\left(\chi\widetilde{T}_{l}\right)\cdot\vec\... ...a P_{l}(\cos\theta)+\sum_{l>0}\chi\widetilde{T}_{l}\nabla^{2}P_{l}(\cos\theta).$$ (81)

Next we replace

$$\frac{\vec\nabla P}{\rho}=\vec g= - \vec \nabla \phi +\vec{\m... ...al{C}} = {1 \over 2} \Omega^2 \vec{\nabla} (r \sin \theta)^2 ,$$ (82)

and expand

$$\rho \chi = \overline{\rho \chi} +\sum_{l>0}\widetilde{\rho \... ...C}}+\sum_{l>0}\widetilde{f}_{\mathcal{C},l}P_{l}(\cos\theta) .$$ (83)

The functions  $\overline{f}_{\mathcal{C}}$ and  $\widetilde{f}_{\mathcal{C},l}$ are easily derived. Expressing the components of the centrifugal force in terms of a_l and b_l introduced in (58), we get

$$\vec\nabla\cdot\vec {\mathcal F}_{\mathcal{C}}=\sum_{l=0}^{\i... ...{l}(r)\right)+\frac{l(l+1)}{r}b_{l}(r)\right]P_{l}(\cos\theta)$$ (84)

and therefore:

$$\overline{f}_{\mathcal{C}}=\frac{1}{r^2}\partial_{r}\left(r^2... ...2}\partial_{r}\left(r^2a_{l}\right)+l(l+1)\frac{b_{l}}{r}\cdot$$ (85)

Then we explicit the turbulent flux; with the assumption that the turbulent diffusivity is much stronger in the horizontal than in the vertical direction ( $D_{\rm h} \gg D_{\rm v}$), its divergence reduces to

$${\vec\nabla}\cdot {\vec F_{\rm h}}=-\frac{1}{r^2\sin\theta} ... ...\theta D_{\rm h}\rho T ~ \partial_{\theta} S(r,\theta)\right]$$ (86)

which we expand in

$$\vec\nabla\cdot\vec F_{\rm h}=\sum_{l>0}\frac{l(l+1)}{r^2}\overline{\rho T}D_{\rm h}\widetilde{S}_{l}P_{l}(\cos\theta).$$ (87)

After some arrangements, and keeping only first order terms, this leads us to

 

$$\vec\nabla\cdot\left(\chi\vec\nabla T\right)+\rho\epsilon+\overli... ...\overline{\rho} \left(\overline{\epsilon}+\overline{\epsilon}_{g}\right)\rangle$$ (88)

$$+\sum_{l>0}\Bigg\{\left(\overline{\rho\chi}\frac{{\rm d}\widetild... ...\chi} \frac{{\rm d}\overline{T}}{{\rm d}P}\right)2\overline{g}\widetilde{g}_{l}$$

$$+{ \overline{\rho}\frac{{\rm d}}{{\rm d}P}\left(\widetilde{\rho\c... ...ver r^2} \overline{\rho T} D_{\rm h} \widetilde{S}_l \Bigg\}}P_{l}(\cos\theta).$$

As mentioned above, the horizontal average <...> is zero on a level surface, which means that the star is in average radiative equilibrium, including the production of nuclear and gravitational energy.

It remains to replace all tilded quantities in terms of two only (for each l-component). Here we choose these two variables to be the horizontal variations of temperature and molecular weight, instead of Zahn (1992) and Maeder & Zahn (1998) who took the density variation. We define

$$\Psi_l = {\widetilde{T}_l \over \overline{T}} ~, \qquad \Lamb... ...qquad \Theta_l= {\widetilde{\rho}_l \over \overline{\rho}} ~,$$ (89)

and use the equation of state in differential form

$${{\rm d}\rho \over \rho} = \alpha {{\rm d}P \over P} - \delta {{\rm d}T \over T} + \varphi {{\rm d} \mu \over \mu}$$ (90)

with the straightforward definitions $\delta = - ({\partial \ln \rho / \partial \ln T})_{P, \mu}$ etc., to express

$$\Theta_l = \varphi \Lambda_l - \delta \Psi_l.$$ (91)

We operate likewise with the other variables:

$$\frac{\widetilde{\chi}_{l}}{\chi}=\chi_{T}\Psi_{l}+\chi_{\mu}... ...l}}{\epsilon}=\epsilon_{T}\Psi_{l}+\epsilon_{\mu}\Lambda_{l},$$ (92)

and introduce

$$K = {\overline{\chi} \over {\overline{\rho} C_P}}, \qquad \rh... ...verline{\epsilon}+\overline{\epsilon}_{\rm grav}\right)}\cdot$$ (93)

Switching to r as independent variable, with  ${\rm d}P = - \overline{\rho g} ~ {\rm d}r$, we cast (88) in its final form:

$$\frac{\vec\nabla\cdot\left(\chi\vec\nabla T\right)+\rho\epsil... ...{\rho}}=\sum_{l>0}\frac{L}{M}\mathcal{T}_{l}P_{l}(\cos\theta)$$ (94)

with

 

$$\mathcal{T}_{l} 2\left[1-\frac{\overline{f}_{\mathcal{C}}}{4\pi G\overline{\rho}}... ...hcal{C}}}{4\pi G\overline{\rho}}\left(-\delta\Psi_{l}+\varphi\Lambda_{l}\right)$$ =

$$-\frac{\rho_{m}}{\overline{\rho}}{r \over 3} \partial_r \left(-H_... ...al_r\Psi_{l}+(1-\delta+\chi_{T})\Psi_{l}+(\varphi+\chi_{\mu})\Lambda_{l}\right)$$

$$+\frac{\left(\overline{\epsilon}+\overline{\epsilon}_{\rm grav}\r... ...l} + \left(\frac{D_{\rm h}}{K}\right){\widetilde{S_{l} } \over C_P}\right]\cdot$$ (95)

Here we have used the property that

$$\overline{\rho}\overline{\chi}\frac{{\rm d}\overline{T}}{{\rm... ...}}\right]{\rm d}V}=\frac{L(r)}{4\pi G M(r) - r^2 a_0(r)}\cdot$$ (96)

To lowest order  $r^2 a_0= (2/3) r^3 \overline \Omega^2$, according to Eq. (61), which justifies the approximation

$$\overline{\rho}\overline{\chi}\frac{{\rm d}\overline{T}}{ {\rm d}P} = \frac{L(r)}{4\pi G M(r)}\cdot$$ (97)

The last step is to put in final form the left hand side of the heat Eq. (78). In a medium of varying composition the entropy of mixing must be taken into account (Maeder & Zahn 1998). In the simplest case, where the stellar material can be approximated by a simple mixture of hydrogen and helium with a fixed abundance of metals, the entropy of mixing can be expressed in terms of the molecular weight only; then we have

$${\rm d}S=C_{\rm p}\left[\frac{{\rm d}T}{T}-\nabla_{\rm ad}\frac{{\rm d}P}{P}+\Phi(P,T,\mu)\frac{{\rm d}\mu}{\mu}\right] ~,$$ (98)

where  $\nabla_{\rm ad}$ is the adiabatic gradient and $\Phi$ is a function of the metal mass fraction and of $\mu$, the mean molecular weight (see Maeder & Zahn 1998). Applying this relation to the fluctuation of entropy around an isobar leads to:

$$\widetilde{S}_{l}=C_{\rm p}\left[\Psi_{l}+\Phi\Lambda_{l}\right].$$ (99)

In the same way, we get for the mean entropy gradient:

$$\partial_{r}\overline{S}=\frac{C_{\rm p}}{H_{\rm p}}\left(\nabla_{\rm ad}-\nabla-\Phi\nabla_{\mu}\right)$$ (100)

with:  $\nabla={\rm d} \ln T/ {\rm d} \ln P$ and  $\nabla_{\mu}={\rm d} \ln\mu /{\rm d} \ln P$.

We refer back to Eq. (78) to express the left hand side of the thermal energy equation, and take into account that $\Lambda_l$ obeys a similar advection/diffusion Eq. (39); this permits the elimination of  ${\rm d} \Lambda_l$/${\rm d}t$ on the l.h.s. and leads us to

$$\overline{T} C_{\rm p}\left[{{\rm d} \Psi_{l} \over {\rm d} t... ...{\rm ad}-\nabla\right)\right]=\frac{L}{M}{\mathcal T}_{l}(r),$$ (101)

where we cast  ${\mathcal T}_{l}(r)$ (95) in its final form, having replaced  $\widetilde{S}_l$ by (99) on the right hand side:

 

$$\mathcal{T}_{l} 2\left[1-\frac{\overline{f}_{\mathcal{C}}}{4\pi G\overline{\rho}}... ...hcal{C}}}{4\pi G\overline{\rho}}\left(-\delta\Psi_{l}+\varphi\Lambda_{l}\right)$$ = (102)

$$+\frac{\rho_{m}}{\overline{\rho}}\left[ \frac{r}{3}\partial_{r}\l... ...t) - \frac{l(l+1)H_{T}}{3r} \left(1 + \frac{D_{\rm h}}{K}\right){\Psi_l}\right]$$

$$+\frac{\left(\overline{\epsilon}+\overline{\epsilon}_{\rm grav}\r... ...l}+(f_{\epsilon}\epsilon_{\mu}+f_{\epsilon}\varphi - \varphi)\Lambda_{l}\Big\}.$$

We recall that  $\overline{f}_{\mathcal{C}} + \sum {\widetilde{f}_{\mathcal{C},l}}P_l(\cos \theta)$ is the divergence of the centrifugal force (cf. Eq. (82)); the functions  $\overline{f}_{\mathcal{C}}~$ and  $\widetilde{f}_{\mathcal{C},l}$ are given in (85).

This expression is similar to that derived in Maeder & Zahn (1998), where  $\Theta_2 = \widetilde{\rho}_2 / \overline {\rho}$ was used as dependent variable instead of  $\Psi_l = \widetilde{T}_l / \overline {T}$ here. However we note an important difference, namely that the l.h.s. involves here the only the thermal part of the subadiabatic gradient  $\left(\nabla_{\rm ad}-\nabla\right)$, whereas in Maeder & Zahn (1998) it was   $\left(\nabla_{\rm ad}-\nabla + (\varphi/ \delta) \nabla_\mu \right)$. Also, here the highest order derivative at the r.h.s. operates only on $\Psi_l$, whereas in Maeder & Zahn (1998) it applied both on $\Theta_2$ and on $\Lambda_2$. We expect therefore the present form to be less sensitive, numerically, to steep composition gradients.

The temperature fluctuation on an isobar is thus governed by an advection/diffusion equation, from which we can derive the radial component of the meridional circulation. It allows us to link the circulation with its causes, namely both the rotation profile and that of chemical composition. The fact that this equation has been derived for a general law of rotation allows us to treat simultaneously the bulk of the radiation zones and the tachoclines.

The system of equations is now completed: we have an advection/diffusion equation for the transport of the angular momentum (mean and fluctuating), another for the transport of chemical elements (mean and fluctuating), another for the temperature (mean and fluctuating), and we have the baroclinic relation which allows us to close the system. It remains to specify the boundary conditions of this system, to be applied on the limits of the radiation zone.

7 Boundary conditions

We shall now review the differential equations we have established for the various quantities which enter in the problem, and state for each of them the appropriate boundary conditions, for a given radiation zone. To be specific, we consider a star with a radiation zone located between a convective core and an upper convection zone, and designate by r_b and r_t the radius respectively of the base and of the top of that radiative zone. Of course, in a solar-type main-sequence star we have r_b=0, whereas for a massive main sequence star r_t=R, R being the radius of the star.

The equation for the mean angular velocity  $\overline{\Omega}$ (19) is of second order in r, and therefore it requires two boundary conditions. They are obtained by evaluating the budget of angular momentum in each adjacent convective zone:

$${ {\rm d} \over {\rm d}t} \left[\int_{0}^{r_{b}}r^{4}\rho \Om... ...\Omega}U_{2}+\mathcal{F}_{\Omega} \quad \hbox{at} \; r=r_{t},$$ (103)

where  $\mathcal{F}_{\Omega}$ is the angular momentum flux which is lost at the surface by the stellar wind. (The perturbation $\Omega_{2}$ obeys an evolution equation which does not include any derivative in r, at least with the approximation we made; therefore it needs no boundary condition.)

Likewise, the equation for the mean concentration (30) is of second order in r, and its two boundary conditions are obtained by calculating the evolution in time of that mean concentration  $\overline{c_{i}}$ in the two adjacent convection zones (cf. Palacios et al. 2001). We thus have respectively for r_t and r_b:

$${ {\rm d} \over {\rm d}t}\left[\overline{c}_{i}\int_{r_{t}}^{R}r^2\rho {\rm d} r\right] r^{2}\rho\left(U^{\rm diff}_i \overline{c}_{i}\right)-r^{2}\rho\l... ...\partial_{r}\overline{c}_{i}-\dot{M}\overline{c}_{i} \quad \hbox{at}\; r=r_{t},$$ =

$${ {\rm d} \over {\rm d}t}\left[\overline{c}_{i}\int_{0}^{r_{b}}r^{2}\rho {\rm d}r\right] -r^{2}\rho\left(U^{\rm diff}_i \overline{c}_{i}\right)+r^{2}\rho\... ...rm v}+D_{\rm eff}\right)\partial_{r}\overline{c}_{i}\quad \hbox{at} \; r=r_{b},$$ = (104)

where $\dot{M}$ is the rate of mass loss through the wind.

Let us examine now the heat Eqs. (101)-(102), which is of second order in r for the variable $\Psi_l$. The variable $\Lambda_l$ need not to be considered here, since it is determined through its evolution Eq. (39). The boundary conditions on $\Psi_l$ are deduced from the baroclinic Eqs. (46)-(47), which involves  $\overline{\Omega}$ and $\Omega_{2}$, and their first order derivatives. We must therefore specify all these functions at the boundaries of the radiation zone. For  $\overline{\Omega}$ this is already done (Eq. (103)), and for  $\Omega_{l}(r)$ we require the fluctuations to be continuous at the boundaries:

$$\Omega_{l}(r)=\Omega_{l,b} \quad \hbox{at} \; r=r_{b}, \qquad \Omega_{l}(r)=\Omega_{l,t} \quad \hbox{at} \; r=r_{t},$$ (105)

thus assuming that there are no stresses between the radiative and the convective zone. From the baroclinic relation (42) we deduce that the gradient of the angular velocity must also be continuous, and therefore that

$$\partial_{r} \overline{\Omega}(r)= \partial_{r}\overline{\Omega}_{b} \quad \hbox{at} \; r=r_{b}, \qquad \;~ \partial_{r} \overline{\Omega}(r) \partial_{r}\overline{\Omega}_{t} \quad \hbox{at} \; r=r_{t},$$ =

$$\partial_{r} \Omega_{l}(r)= \partial_{r}\Omega_{l,b} \quad \hbox{at} \; r=r_{b}, \qquad \partial_{r} \Omega_{l}(r) \partial_{r}\Omega_{l,t}\quad \hbox{at} \; r=r_{t}.$$ = (106)

Since we do not solve here for the rotation law in the convective regions (a formidable task!), the best we can afford is to apply a horizontal profile which is inspired from helioseismology, for the base of the outer convection zone. Since the rotation rate  $\Omega(r,\theta)$ seems to be independent of r in the solar convection zone, we can take

$$\partial_{r} \overline{\Omega}(r)= 0 \quad \hbox{at} \; r=r_{... ...d \partial_{r} \Omega_{l}(r)= 0 \; \hbox{at} \quad r=r_{t}.$$ (107)

For the condition at the boundary of a convective core, we will have to rely on numerical simulations (see the recent results of Browning et al. 2004). In the case where we stop the expansion of rotation at the $\Omega_{2}$ term, and where we assume that the vertical gradient of $\Omega$ is zero, we obtain:

$$\varphi\Lambda_{2}-\delta\Psi_{2} \frac{8}{7}\frac{r}{\overline{g}}\overline{\Omega}\Omega_{2},$$ =

$$\varphi\Lambda_{4}-\delta\Psi_{4} -\frac{12}{35}\frac{r}{\overline{g}}\overline{\Omega}\Omega_{2}.$$ = (108)

Finally, let us state the boundary conditions for the perturbed Poisson equation. Unlike the precedent equations above, this second order equation must be integrated over the whole star. We require regularity at origin, and continuity with a potential field at the surface (r=R):

$$\widehat{\phi}_l = 0 \; \hbox{at} \; r=0, \qquad (l+1) \wideh... ...d} \over {\rm d}r} \widehat{\phi}_l = 0 \; \hbox{at} \; r=R.$$ (109)

8 Conclusion

The work presented here marks another step in the description of rotational mixing in stellar radiation zones, after the papers by Zahn (1992) and Maeder & Zahn (1998). When applied to massive stars, which are fast rotators, such mixing yields much better agreement between models and observations, as was discussed by Meynet & Maeder (2000). However numerical problems arise during the advanced stages of evolution, when steep gradients of composition and rotation develop. For this reason, our main purpose here was to increase the accuracy of the modelization in latitude by adding higher order terms to the expansion in spherical harmonics, and also to better resolve the phases of rapid evolution by retaining all time derivatives, except those related to the dynamical relaxation.

This allows us to include the transport occurring in the tachoclines, where the differential rotation imposed by the adjacent convection zone generates an octupolar circulation, with no net advection of angular momentum. Until now, these layers required an ad-hoc treatment, such as performed by Brun et al. (1999, 2002), who found that the mixing in the solar tachocline contributes indeed to the depletion of lithium, and that it modifies sufficiently the structure of the star to be detected through helioseismology.

From a more technical point of view, we have written here the heat Eqs. (101)-(102), which determines the meridional circulation, as a relaxation equation for the temperature fluctuation $\Psi_l$. This breaks the apparent symmetry between the circulation driven by the differential rotation profile (hence by  $\Psi_l$) and the counteracting flow (Mestel's $\mu$-current) which is induced by the inhomogeneities in chemical composition ($\Lambda_l$), since now the expression giving the circulation velocity involves only the first derivative in $\Lambda_l$, while keeping the second derivative of $\Psi_l$. A benefit of that choice is that it renders this differential equation less sensitive, numerically, to the steep composition gradients which arise in the course of stellar evolution.

One weakness of the modelization of rotational mixing remains the description of the turbulence which is generated by the differential rotation. Such turbulence is certainly anisotropic, due to the stable stratification, and it tends therefore to smooth the angular velocity and the chemical composition on horizontal surfaces. This problem is addressed in the companion paper (Mathis et al. 2004), where we discuss the prescriptions for the turbulent transport which are presently available.

A major shortcoming of the present modelization of rotational mixing is that it predicts a rapidly spinning solar core, contrary to the findings of helioseismology. This proves that other physical processes contribute to the transport of angular momentum in solar-type stars, which are presently slow rotators because they have been spun down by their wind. Internal gravity waves, which are generated by turbulent convection at the interface with a convective region, are one possible cause. Their effect has been studied by Talon et al. (2002), and is now being introduced in stellar evolution codes (Talon & Charbonnel 2003).

Magnetic stresses are another serious candidate for the angular momentum transport, as was pointed out already by Mestel (1953). The alternating dynamo field does not penetrate enough into the solar radiation zone to produce any impact on the rotation law, as was shown by Garaud (1999). But even a weak fossil field could enforce nearly uniform rotation, although the outcome depends sensitively on the poloidal field which is assumed, and on whether it penetrates into the differentially rotating convection zone, as illustrated by the calculations performed by Ruediger & Kitchanitov (1996, 1997) and by MacGregor & Charbonneau (1999). Therefore the problem must be treated consistently, taking into account the advection of the field by the meridional circulation, itself being modified by the field, and imposing proper boundary conditions. We shall address this problem in two forthcoming papers, one dealing with an axisymmetric field and the second with a non-axisymmetric field, as observed in oblique rotators.

Acknowledgements

We are grateful to C. Charbonnel, A. Maeder, P. Morel, A. Palacios and S. Talon for their careful reading of the manuscript and their helpful remarks. This work was supported in part by the Centre National de la Recherche Scientifique (Programme National de Physique Stellaire).

References

• Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [NASA ADS] [CrossRef] (In the text)
• Brun, A. S., Turck-Chièze, S., & Zahn, J.-P. 1999, ApJ, 525, 1032 [NASA ADS] [CrossRef] (In the text)
• Brun, A. S., Antia, H. M., Chitre, S. M., & Zahn, J.-P. 2002, A&A, 391, 725 [EDP Sciences] [NASA ADS] (In the text)
• Busse, F. H. 1982, ApJ, 259, 759 [NASA ADS] [CrossRef] (In the text)
• Chaboyer, B., & Zahn, J.-P. 1992, A&A, 253, 173 [NASA ADS] (In the text)
• Eddington, A. S. 1925, Observatory, 48, 78 [NASA ADS] (In the text)
• Garaud, P. 1999, MNRAS, 304, 583 [NASA ADS] [CrossRef] (In the text)
• Kippenhahn, R. 1963, ApJ, 137, 664 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• MacGregor, K. B., & Charbonneau, P. 1999, ApJ, 519, 911 [NASA ADS] [CrossRef] (In the text)
• Maeder, A. 2003, A&A, 399, 263 [EDP Sciences] [NASA ADS]
• Maeder, A., & Zahn, J.-P. 1998, A&A, 334, 1000 [NASA ADS] (In the text)
• Maeder, A., & Meynet, G. 2001, A&A, 373, 555 [EDP Sciences] [NASA ADS] (In the text)
• Meynet, G., & Maeder, A. 2000, A&A, 361, 101 [NASA ADS] (In the text)
• Mathias, J., & Zahn, J.-P. 1997, IAU Symp. 181, ed. J. Provost & F.-X. Schmider; poster volume, 103 (In the text)
• Mathis, S., Palacios, A., & Zahn, J.-P. 2004, A&A, 425, 243 [EDP Sciences] [NASA ADS] (In the text)
• Mestel, L. 1953, MNRAS, 113, 716 [NASA ADS] (In the text)
• Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603 [EDP Sciences] [NASA ADS]
• Ruediger, G., & Kitchatinov, L. L. 1996, ApJ, 466, 1078 [NASA ADS] [CrossRef] (In the text)
• Ruediger, G., & Kitchatinov, L. L. 1997, Astron. Nachr., 318, 273 [NASA ADS] (In the text)
• Spiegel, E., & Zahn, J.-P. 1992, A&A, 265, 106 [NASA ADS] (In the text)
• Sweet, P. A. 1950, MNRAS, 110, 548 [NASA ADS] [MathSciNet] (In the text)
• Talon, S., & Charbonnel, C. 1998, A&A, 335, 959 [NASA ADS] (In the text)
• Talon, S., & Charbonnel, C. 2003, A&A, 405, 1025 [EDP Sciences] [NASA ADS] (In the text)
• Talon, S., Kumar, P., & Zahn, J.-P. 2002, ApJ, 574, L175 [NASA ADS] [CrossRef] (In the text)
• Talon, S., Zahn, J.-P., Maeder, A., & Meynet, G. 1997, A&A, 322, 209 [NASA ADS] (In the text)
• Vogt, H. 1925, Astron. Nachr., 223, 229 [NASA ADS] (In the text)
• Von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] (In the text)
• Zahn, J.-P. 1992, A&A, 265, 115 [NASA ADS] (In the text)

 

  
Online Material

Appendix A: Algebra related to the spherical harmonics

The spherical harmonics are defined by:

$$Y_{l}^{m}(\theta,\varphi)=\mathcal{N}_{l}^{m}P_{l}^{\vert m\vert}\left(\cos\theta\right){\rm e}^{{\rm i}m\varphi}$$ (A.1)

with the normalization:

$$\mathcal{N}_{l}^{m}=(-1)^{\frac{\left(m+\vert m\vert\right)}{... ...l-\vert m\vert)!}{(l+\vert m\vert)!}\right]^{\frac{1}{2}}\cdot$$ (A.2)

We recall that they obey the well-known differential equation:

$$\frac{1}{\sin\theta}\partial_{\theta}\left(\sin\theta\partial... ...right)=-l\left(l+1\right)Y_{l}^{m}\left(\theta,\varphi\right).$$ (A.3)

In deriving certain equations in Sect. 5 we used the following recursion relations for m=0:

$$\cos\theta Y_{l}^{0}(\theta)=A_{l}^{0}Y_{l-1}^{0}(\theta)+B_{... ...uad \hbox{and} \; B_{l}^{0}=\frac{(l+1)}{\sqrt{(2l+3)(2l+1)}},$$ (A.4)

$$\sin\theta Y_{l}^{0}(\theta)=C_{l}^{0}\partial_{\theta}Y_{l-1... ...} \quad \hbox{and} \; D_{l}^{0}=\frac{1}{\sqrt{(2l+3)(2l+1)}},$$ (A.5)

$$\cos\theta\partial_{\theta}Y_{l}^{0}(\theta)=E_{l}^{0}\partia... ...} \quad \hbox{and} \; F_{l}^{0}=\frac{l}{\sqrt{(2l+3)(2l+1)}},$$ (A.6)

$$\sin\theta\partial_{\theta}Y_{l}^{0}(\theta)=G_{l}^{0}Y_{l+1}... ...hbox{and} \; H_{l}^{0}=\frac{l(l+1)}{\sqrt{(2l+1)(2l-1)}}\cdot$$ (A.7)

The identities (A.5) and (A.6) have been deduced from the two others (A.4)-(A.7), with the help of (A.3).

Appendix B: System of transport equations for l = 2, 4

In this appendix we give the explicit form of the transport equations for l=2, 4, with only the first term $\Omega_{2}$ being kept in the expansion of the rotation rate $\Omega$.

B.1 Mean equations

B.1.1 mean rotation rate

$$% latex2html id marker 2525 \rho\frac{{\rm d}}{{\rm d}t}\left... ...\overline{\Omega}\right) \qquad \hbox{(Eq.~(\ref{mean-AM}))} ;$$ (B.1)

B.1.2 mean chemical composition

$$% latex2html id marker 2546 \rho {{\rm d} \over {\rm d} t } \... ...verline{c_{i}}\right] \qquad \hbox{(Eq.~(\ref{mean-c-fin}))}.$$ (B.2)

B.2 System for l = 2

B.2.1 meridional circulation

We rewrite (101)-(102) in two first order equations as

$$U_{2}=\frac{L}{M \overline{g}}\left(\frac{P}{\overline{\rho}C... ...line{T}}\right)\frac{1}{\nabla_{\rm ad}-\nabla}\mathcal{B}_{2}$$ (B.3)

where we have defined:

$$\mathcal{B}_{2} -2\left[1-{\partial_r (r^3 \overline{\Omega}^2) \over 6 \pi G \ov... ...t) - {{\rm d} \over {\rm d}r} \left({\widehat \phi_2 \over g_0} \right) \right]$$ =

$$+\frac{1}{\pi G\overline{\rho}} \left[\frac{6}{7}\overline{\Omega... ...athcal{A}_{2}-\frac{2H_{T}}{r}\left(1+\frac{D_{\rm h}}{K}\right)\Psi_{2}\right]$$ (B.4)

$$+\frac{\overline{\epsilon}+\overline{\epsilon}_{\rm grav}}{\epsil... ...r {\rm d}t} + \Phi {{\rm d} \ln \overline \mu \over {\rm d}t} \Lambda_2 \right)$$

and:

$$\mathcal{A}_{2}=H_{T}\partial_{r}\Psi_{2} -\left(1-\delta + \... ...}\right)\Psi_{2} -\left(\varphi+ \chi_{\mu}\right)\Lambda_{2},$$ (B.5)

and where we have replaced the gravity fluctuation by (74),  $\overline{f}_{\mathcal{C}}$ by (84) and  $\widetilde{f}_{{\mathcal C}ñ,2}$ by (85).

B.2.2 baroclinic relation

From Eqs. (46)-(48) we have:

$$\varphi\Lambda_{2}-\delta\Psi_{2}=\frac{r}{\overline{g}}\left... ...}+\frac{8}{35}r\overline{\Omega}\partial_{r}\Omega_{2}\right].$$ (B.6)

B.2.3 horizontal fluctuation of the molecular weight

We apply (39) to l=2:

$$\frac{{\rm d}\Lambda_{2}}{{\rm d}t}-\frac{{\rm d}\ln\overline... ...{2}}{H_{\rm p}}\nabla_{\mu}-\frac{6}{r^2}D_{\rm h}\Lambda_{2}.$$ (B.7)

B.2.4 Poisson equation

Likewise, for (60):

$$\frac{1}{r}\frac{{\rm d}^2}{{\rm d}r^2}\left(r\widehat{\phi}_... ...}+\frac{{\rm d}}{{\rm d}r}\left(r\rho_{0}b_{2}\right)\right].$$ (B.8)

B.3 System for l = 4

B.3.1 horizontal shear

$$% latex2html id marker 2765 \rho \frac{{\rm d}\left(r^{2}\Ome... ...rho\nu_{\rm h}\Omega_{2} \qquad \hbox{(Eq.~(\ref{omega2}))}.$$ (B.9)

B.3.2 meridional circulation

We proceed as for l=2:

$$U_{4}=\frac{L}{M\overline{g}}\left(\frac{P}{\overline{\rho}C_... ...line{T}}\right)\frac{1}{\nabla_{\rm ad}-\nabla}\mathcal{B}_{4}$$ (B.10)

where we have defined:

$$\mathcal{B}_{4} -2\left[1- {\partial_r (r^3 \overline{\Omega}^2) \over 6 \pi G \o... ...ght)- {{\rm d} \over {\rm d}r} \left({\widehat \phi_4 \over g_0} \right)\right]$$ =

$$+\frac{1}{\pi G \overline{\rho}}\left[\frac{6}{35}\left(2\overlin... ...hcal{A}_{4}-\frac{20H_{T}}{3r}\left(1+\frac{D_{\rm h}}{K}\right)\Psi_{4}\right]$$ (B.11)

$$+\frac{\overline{\epsilon}+\overline{\epsilon}_{\rm grav}}{\epsil... ...r {\rm d}t} + \Phi {{\rm d} \ln \overline \mu \over {\rm d}t} \Lambda_4 \right)$$

and:

$$\mathcal{A}_{4}=H_{T}\partial_{r}\Psi_{4}-\left(1-\delta + \c... ...}\right)\Psi_{4} -\left(\varphi+ \chi_{\mu}\right)\Lambda_{4}.$$ (B.12)

B.3.3 baroclinic relation

We restate (46)-(48) as

$$\varphi\Lambda_{4}-\delta\Psi_{4}=\frac{6}{35}\frac{r}{\overl... ...ga}\partial_{r}\Omega_{2}-2\overline{\Omega}\Omega_{2}\right].$$ (B.13)

B.3.4 horizontal fluctuation of the molecular weight

$$% latex2html id marker 2902 \frac{{\rm d}\Lambda_{4}}{{\rm d}... ...}D_{\rm h}\Lambda_{4} \qquad \hbox{(Eq.~(\ref{lamba-evol}))}.$$ (B.14)

B.3.5 Poisson equation

$$% latex2html id marker 2920 \frac{1}{r}\frac{{\rm d}^2}{{\rm ... ...{4}\right)\right] \qquad \hbox{(Eq.~(\ref{poisson-pert}))}.$$ (B.15)

These equations are ready to be implemented in a stellar evolution code, together with the boundary conditions discussed in Sect. 7.