A&A 440, 653-666 (2005)
DOI: 10.1051/0004-6361:20052640

Transport and mixing in the radiation zones of rotating stars

II. Axisymmetric magnetic field

S. Mathis - J.-P. Zahn

LUTH, Observatoire de Paris, 92195 Meudon, France

Received 4 January 2005 / Accepted 11 May 2005

Abstract
The purpose of this paper is to improve the modeling of the mixing of chemical elements that occurs in stellar radiation zones. In addition to the classical rotational mixing considered in our previous paper, which results of the combined action of the thermally-driven meridional circulation and of the turbulence generated by the shear of differential rotation, we include here the effect of an axisymmetric magnetic field in a self-consistent way. We treat the advection of the field by the meridional circulation, its Ohmic diffusion, and the production of its toroidal component through the shear of differential rotation. The Lorentz force is assumed not to exceed the centrifugal force; it acts on the baroclinic balance and therefore on the meridional flow, and it has a strong impact on the transport of angular momentum. All variables and governing equations are expanded in spherical or spherical vectorial functions, to arbitrary order: this yields a system of partial differential equations in time and in the radial coordinate, which is ready to be implemented in a stellar structure code.

Key words: magnetohydrodynamics (MHD) - turbulence - stars: rotation - stars: magnetic fields - stars: evolution

1 Motivation

The need for improved stellar structure models, which go beyond the so-called standard model, is now well recognized. The radiation zones can no longer be treated as motionless regions, where no mixing occurs, but one has to include physical processes which lead to the transport of matter and angular momentum. So far the focus has been on the effect of rotation, which generates a thermally driven meridional circulation; by advecting angular momentum, that circulation renders the rotation non-uniform and prone to shear turbulence. This rotational mixing has been formulated by Zahn (1992) and Maeder & Zahn (1998), assuming that the angular velocity varies much less in the horizontal than in the vertical direction (shellular rotation). For rapidly rotating stars, models built along these lines are in much better agreement than standard models; this was first shown by Talon et al. (1997), and confirmed by Maeder & Meynet (2000, 2001), who in addition refined the description for the mass loss by taking into account the latitudinal variation of the wind.

However, these models fail to correctly reproduce the flat rotation profile observed in the solar radiative interior (Matias & Zahn 1997). In solar-type stars, at least once they have been spun down through their magnetized wind, other physical processes are therefore responsible for the transport of angular momentum, and the most plausible candidates are now being investigated. One is the transport by internal gravity waves emitted at the base of the convection zones, associated with turbulence which smoothes the differential rotation; this mechanism was described by Kumar et al. (1999), Talon et al. (2002) and applied by Talon & Charbonnel (2003, 2004) to stars of various types. The other possibility is magnetic torquing, as advocated already by Mestel (1953): he showed that even a very modest field could enforce rigid rotation. In the present paper we shall address the latter effect, and examine how rotational mixing is affected by an axisymmetric field.

Moss (1974) was the first to map the thermally-driven meridional circulation in a star containing an axisymmetric poloidal field; when rotation was added, it was assumed constant and uniform. Various generalizations of this work followed, by Mestel & Moss (1977) and by Moss (1977, 1984, 1985, 1987). The a priori assumption of rigid rotation was relaxed by Tassoul & Tassoul (1986), who claimed that no magnetic field could render the rotation nearly uniform in presence of meridional circulation; however their treatment failed to correctly represent the generation of toroidal field through the shearing of poloidal field, which provides the major feed-back in this problem. This was pointed out by Mestel et al. (1988), who carried out numerical simulations that illustrated how a weak magnetic field could enforce nearly uniform rotation within a few Alfvén times; for sake of simplification, both the poloidal field and the meridional circulation were taken as given, and constant in time.

Recently Garaud (2002) treated the problem in a fully consistent way, and she applied it to the description of the solar tachocline: she let the field adjust through Ohmic diffusion and be advected by the meridional circulation. However, due to numerical limitations, her circulation was driven mainly through Ekman pumping, rather than through thermal diffusion; also, only stationary solutions were constructed.

The purpose of the present paper is to model the time-dependent rotational mixing in presence of an axisymmetric magnetic field, thus extending the results obtained by Mathis & Zahn (2004; from now on referred to as Paper I) in the purely hydrodynamical case. To achieve this goal:

• we allow for a general axisymmetric magnetic field, including all its effects;
• we expand all variables in spherical harmonics, formally to unlimited order;
• we treat explicitely the departures from shellular rotation, in the linear approximation, in order to capture the tachocline(s);
• we filter out the short times, and keep only the time derivative in the heat equation and in that describing the transport of angular momentum.

In a forthcoming paper, we shall generalize the present work to a non-axisymmetric magnetic field.

2 Assumptions and expansions

We focus here our attention to magnetic fields of moderate strength, whose Alfvén speed $V_{\rm A}$ does not exceed the rotation velocity $R \Omega$, but is larger than the meridional circulation velocity $V_{\rm M}$. Thus the field will not affect much the hydrostatic balance in the meridional plane, but it may play a major role in the transport of angular momentum, by tending to render the angular velocity $\Omega$ constant along the field lines of the poloidal field (Ferraro 1937). We also assume that the field does not vary on a timescale which is shorter than the Alfvén time  $R/V_{\rm A}$; thus we shall deal with fields of primordial origin, and exclude possible dynamo fields produced in the adjacent convection zones, whose penetration would be damped anyway by skin effect (Garaud 1999).

Our treatment is not able to address the question whether the field is stable or not, because it is axisymmetric, and because we do not resolve the Alfvén waves. We will have to rely on our intuition, guided by 3D simulations such as recently performed by Braithwaite & Spruit (2004), to choose the magnetic configurations we implement in our code. As for instabilities which may lead to a turbulent state, for instance the magneto-rotational instability when the equator rotates slower than the poles (Balbus & Hawley 1994), their transport properties will have to be accounted for by a suitable parametrization, as was done previously in the hydrodynamical case.

As in Paper I, we assume that the rotating star is only weakly two-dimensional, and this for two reasons. The first is that the rotation rate and the magnetic field, which is taken axisymmetric, are sufficiently moderate to allow the centrifugal and the Lorentz forces to be considered as perturbations compared to gravity. The second reason rests on the less justified hypothesis that the shear instabilities due to the differential rotation give 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. In the radiation zones, we expect it to smooth out the horizontal variations of angular velocity and of chemical composition, a property we shall use to discard certain non-linear terms. Let us emphasize also that the influence of the magnetic field on such turbulence is not taken into account, as seems to be allowed by the condition $V_{\rm A} < R \Omega$.

Hence, we 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. As reference surface, we choose 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)

where $\tilde{X}_{l}({r})$ and $\widehat{X}_{l}(r)$ are related by (cf. Paper I, Sect. 2)

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

with r being here the mean radius of an isobar.

Likewise, we expand all axisymmetric vector fields using the method outlined by Rieutord (1987):

$$\vec{u}(r,\theta)=\sum_{l=0}^{\infty}\left\{u_{0}^{l}(r)\vec ... ... S_{l}^{0}(\theta)+w_{0}^{l}(r)\vec T_{l}^{0}(\theta)\right\}.$$ (4)

The orthogonal axisymmetric vectorial spherical harmonics $\vec R_{l}^{0}\left(\theta\right)$, $\vec S_{l}^{0}\left(\theta\right)$, $\vec T_{l}^{0}\left(\theta\right)$ are defined by:

$$\vec R_{l}^{0}(\theta)=Y_{l}^{0}(\theta)~\vec{\widehat{e}}_{r... ...}(\theta)=\vec\nabla_{\mathcal{S}}\wedge\vec R_{l}^{0}(\theta)$$ (5)

where $Y_{l}^{0}\left(\theta\right)$ is the classical spherical harmonic with m=0 (cf. Edmonds 1974) and $\vec\nabla_{\mathcal{S}}$ is the horizontal gradient

$$\vec\nabla_{\mathcal{S}}=\vec{\widehat{e}}_{\theta}\partial_{... ...{\widehat{e}}_{\varphi}\frac{1}{\sin\theta}\partial_{\varphi}.$$ (6)

Their detailed properties are given in Appendix A.

Using Eqs. (A.1) and (A.6), $\vec u$ may be written explicitly on the usual vector basis in spherical coordinates:

$$\vec u \sum_{l=0}^{\infty}\mathcal{N}_{l}^{0}\left[u_{0}^{l}\left(r\righ... ...)\partial_{\theta}P_{l}\left(\cos\theta\right)\vec{\widehat e}_{\varphi}\right]$$ =

$$\mathcal{N}_{0}^{0}u_{0}^{0}\left(r\right)\vec{\widehat e}_{r}+\s... ...left(r\right)P_{l}^{1}\left(\cos\theta\right)\vec{\widehat e}_{\varphi}\right];$$ = (7)

$\vec{\widehat e}_{r}$, $\vec{\widehat e}_{\theta}$ and $\vec{\widehat e}_{\varphi}$ are the unit-vectors repectively in the r, $\theta$ and $\varphi$ directions.

This expansion of the vector fields allows us to separate explicitly the spatial coordinates $(r, \theta)$ in the vectorial partial differential equations which govern the problem: using this decomposition, we cast the problem into partial differential equations in t and r only. This point is very important numerically, because the existing stellar evolution codes in which we have to introduce the transport equations are all 1-dimensional, and also because we need to achieve much higher accuracy in the radial direction than in the horizontal.

3 Axisymmetric magnetic field

We start with expansions and equations involving the magnetic field.

3.1 Definitions

The magnetic field, being divergenceless, is conveniently split in its poloidal and toroidal parts:

$$\vec B\left(r,\theta\right)=\vec B_{\rm P}\left(r,\theta\righ... ...\xi_{\rm T}\left(r,\theta\right)\vec{\widehat{e}}_{r}\right) ,$$ (8)

where $\xi_{\rm P}$ and $\xi_{\rm T}$ are respectively the poloidal and the toroidal magnetic stream-functions. Using the classical method introduced by Bullard & Gellman (1954) for the spectral treatment of the geodynamo problem (see also James 1973, 1974 and Serebrianaya 1988), we get the following expansion in radial functions and in spherical harmonics which are defined in Sect. A.1.1:

$$\xi_{\rm P}\left(r,\theta\right)=\sum_{l=1}^{\infty}\xi_{0}^{... ...\infty}\chi_{0}^{l}\left(r\right)Y_{l}^{0}\left(\theta\right).$$ (9)

Since $\vec{R}_{l}^{0}\left(\theta\right)=Y_{l}^{0}\left(\theta\right)\vec{\widehat{e}}_{r}$, as we have seen in the previous section where we have defined $\vec{R}_{l}^{0}\left(\theta\right), \vec{S}_{l}^{0}\left(\theta\right), \vec{T}_{l}^{0}\left(\theta\right)$, the field may be written as

$$\vec B_{\rm P}=\sum_{l=1}^{\infty} \vec\nabla\wedge\vec\nabla... ...dge\left(\chi_{0}^{l}\vec R_{l}^{0}\left(\theta\right)\right).$$ (10)

Next, we develop the curl operator as explained in the appendix (Eqs. (A.32) and (A.33)) to obtain the following expansions for the magnetic field:

$$\vec B\left(r,\theta\right)= \vec B_{\rm P}\left(r,\theta\rig... ...}^{l}}{r}\right]\vec T_{l}^{0}\left(\theta\right)\right\}\cdot$$ (11)

We proceed likewise for the current density $\vec j$, which in the framework of magnetohydrodynamics is related with $\vec B$ through the Maxwell-Ampère equation, neglecting the displacement current:

$$\vec j\left(r,\theta\right)=\frac{1}{\mu_{0}}\vec\nabla\wedge\vec B\left(r,\theta\right),$$ (12)

$\mu_{0}$ being the magnetic permeability of vacuum.

Therefore, using expression (8) for the field, we get:

$$\vec j=\frac{1}{\mu_{0}}\vec\nabla\wedge\vec B=\frac{1}{\mu_{... ...bla\wedge\left(\xi_{\rm T}\vec{\widehat{e}}_{r}\right)\right],$$ (13)

which, using the relations (A.33) and (A.36), we can project again on $\vec R_{l}^{0}\left(\theta\right)$, $\vec S_{l}^{0}\left(\theta\right)$ and $\vec T_{l}^{0}\left(\theta\right)$, and split into its poloidal and toroidal parts:

$$\vec j_{\rm P}\left(r,\theta\right)=\frac{1}{\mu_{0}}\vec\nab... ...r}\chi_{0}^{l}\right]\vec S_{l}^{0}\left(\theta\right)\right\}$$ (14)

$$\vec j_{\rm T}\left(r,\theta\right)=\frac{1}{\mu_{0}}\vec\nab... ...r}\right)\right]\vec T_{l}^{0}\left(\theta\right)\right\}\cdot$$ (15)

with $\Delta_{l}$ being the Laplacian operator

$$\Delta_{l}=\partial_{r,r}+\frac{2}{r}\partial_{r}-\frac{l(l+1)}{r^2}\cdot$$ (16)

We are now ready to examine the aspects of transport and mixing in the radiation zones of rotating stars that are related to the presence of an axisymmetric magnetic field. To achieve this, we shall proceed in three steps

• first, we introduce the Lorentz force;
• next, we derive the transport equation for $\vec B$, which is the classical induction equation;
• finally, we take into account the energy losses due to Ohmic heating.

3.2 Lorentz force

The Lorentz force $\vec{\mathcal F}\!\!_{\mathcal L}$ plays a crucial role in a rotating star, since it acts to render the angular velocity constant on the field lines of the poloidal field $\vec B_{\rm P}$. In this section, we establish the formal expression of $\vec{\mathcal F}\!\!_{\mathcal L}$ in terms of vectorial spherical harmonics.

Starting from the definition of the Lorentz force:

$$\vec{\mathcal{F}}\!\!_{\mathcal{L}}\left(r,\theta\right)=\vec... ...[\frac{1}{\mu_{0}}(\vec\nabla\wedge\vec B)\right]\wedge\vec B,$$ (17)

we replace $\vec B$ by its expansions (11), and use the algebra ruling the vector product of two general axisymmetric vectors in Appendix A, whereby we obtain the following formal projection:

 

$$\vec{\mathcal{F}}\!\!_{\mathcal{L}}\left(r,\theta\right) \mathcal{X}_{{\vec {\mathcal F}}_{\mathcal{L}};0}(r)\vec{\widehat... ...F}}_{\mathcal{L}};l}(r)P_{l}^{1}(\cos\theta)\vec{\widehat{e}}_{\varphi}\right\}$$ =

$$\left[\frac{\mathcal{X}_{{\vec {\mathcal F}}_{\mathcal{L}};0}(r)}... ...thcal{L}};l}(r)}{\mathcal{N}_{l}^{0}}\right]\vec T_{l}^{0}(\theta)\right\}\cdot$$ = (18)

In Appendix B, we give the explicit expressions for the radial functions $\mathcal{X}_{\vec{\mathcal F}\!\!_{\mathcal L};l}\left(r\right)$, $\mathcal{Y}_{\vec{\mathcal F}\!\!_{\mathcal L};l}\left(r\right)$, $\mathcal{Z}_{\vec{\mathcal F}\!\!_{\mathcal L};l}\left(r\right)$ in terms of the magnetic stream-functions $\xi_{0}^{l}$ and $\chi_{0}^{l}$, the normalization coefficient $\mathcal{N}_{l}^{0}$ being given in (A.2).

As we have done previously for the magnetic field $\vec B$ and for the current $\vec j$, we split the magnetic force into its poloidal and toroidal parts:

$$\vec{\mathcal{F}}\!\!_{\mathcal{L},\rm P} \vec j_{\rm T}\wedge\vec B_{\rm P}+\vec j_{\rm P}\wedge\vec B_{\r... ...cal F}}_{\mathcal{L}};l}P_{l}^{1}(\cos\theta)\vec{\widehat{e}}_{\theta}\right\}$$ =

$$\left[\frac{\mathcal{X}_{{\vec {\mathcal F}}_{\mathcal{L}};0}}{\m... ...F}}_{\mathcal{L}};l}}{\mathcal{N}_{l}^{0}}\right]\vec S_{l}^{0}(\theta)\right\}$$ = (19)

$$\vec{\mathcal{F}}\!\!_{\mathcal{L},\rm T}=\vec j_{\rm P}\wedg... ...\mathcal{N}_{l}^{0}}\right]\vec T_{l}^{0}(\theta)\right\}\cdot$$ (20)

As it was underlined by Mestel et al. (1988) and as we shall see in Sects. 4-6, this decomposition is very useful: the poloidal Lorentz force $\vec{\mathcal F}\!\!_{\mathcal L,\rm P}$ operates on the hydrostatic balance and thus contributes to the thermal imbalance, while the toroidal component $\vec{\mathcal F}\!\!_{\mathcal L,\rm T}$ acts on the transport of angular momentum through its torque.

The expression of $\vec{\mathcal F}\!\!_{\mathcal L}$ in the case where only the $l\!=\!\left\{1,2,3\right\}$ terms are kept in the expansion of $\vec B$ can be found in Appendix E. One should note that if the expansion of $\vec B$ is terminated at mode $l_{\rm max}=N_{\rm m}$ in Eq. (11), then the expansion of $\vec{\mathcal F}\!\!_{\mathcal L}$ involves terms up to $l_{\rm max}=2N_{\rm m}$, due to selection rules (see Appendix B).

3.3 Induction equation

We now turn to the evolution of the magnetic field in stellar radiation zones, which is governed by the induction equation:

$$\partial_{\it t}\vec B=\vec\nabla\wedge\left(\vec V\wedge\vec... ...\vert\vert\eta\vert\vert\otimes\vec\nabla\wedge\vec B\right) ,$$ (21)

where we have allowed for an anisotropic eddy diffusivity $\vert\vert\eta\vert\vert$. We recall (cf. Paper I) that the macroscopic velocity field $\vec V$ is the sum of a zonal flow $\vec{\mathcal U}_{\varphi} = r\sin\theta ~\Omega\left(r,\theta\right) \vec{\widehat{e}}_{\varphi}$ and of a meridional flow $\vec{\mathcal U}\left(r,\theta\right)$:

$$\vec V\left(r,\theta\right)=r\sin\theta ~ \Omega\left(r,\thet... ...\widehat{e}}_{\varphi}+\vec{\mathcal{U}}\left(r,\theta\right).$$ (22)

The latter can be split into a spherically symmetric part, which represents the contractions and dilatations of the star during its evolution, plus the thermally-driven circulation

$$\vec{\mathcal{U}}\left(r,\theta\right)=\dot{r}~\vec{\widehat{e}}_{r}+\vec{\mathcal{U}}_{\rm M}\left(r,\theta\right) ,$$ (23)

the meridional flow being expanded in spherical functions:

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

The radial functions U_l(r) and V_l(r) are related by the continuity equation, i.e. $\vec{\nabla} \cdot (\rho~\vec {\mathcal{U}}_{\rm M})=0$ in the anelastic approximation:

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

We introduce the expanded form of the velocity field in the induction Eq. (21):

$$\partial_{\it t}\vec B-\vec\nabla\wedge\left(\dot{r}~ \vec{\w... ...\vert\vert\eta\vert\vert\otimes\vec\nabla\wedge\vec B\right) ,$$ (26)

and as before we project this equation on the vectorial spherical harmonics $\vec R_{l}^{0}\left(\theta\right)$, $\vec S_{l}^{0}\left(\theta\right)$, $\vec T_{l}^{0}\left(\theta\right)$. The time derivative on the left-hand side of Eq. (26) is readily derived (cf. Eq. (11)):

$$\partial_{t}\vec B=\sum_{l=1}^{\infty}\left\{\left[l(l+1)\fra... ..._{t}\chi_{0}^{l}}{r}\right]\vec T_{l}^{0}(\theta)\right\}\cdot$$ (27)

Then, using the algebra related to the $\vec R_{l}^{0}\left(\theta\right)$, $\vec S_{l}^{0}\left(\theta\right)$, $\vec T_{l}^{0}\left(\theta\right)$ and the identity (A.32), the following expansion of the contractions/dilatations term is obtained:

$$\vec\nabla\wedge\left(\dot{r}\hbox{ }\vec{\widehat e}_{r}\wed... ...l}\right)\right]\vec T_{l}^{0}\left(\theta\right)\right\}\cdot$$ (28)

In the same way, we get for the dissipation term:

 

$$\vec\nabla\wedge\left(\vert\vert\eta\vert\vert\otimes\vec\nabla\wedge\vec B\right) -\sum_{l=1}^{\infty}\left\{\left[\eta_{\rm h}\frac{1}{r}\Delta_{l... ...ac{\xi_{0}^{l}}{r}\right)\right)\right]\vec S_{l}^{0}\left(\theta\right)\right.$$ =

$$+{\left.\left[\frac{1}{r}\partial_{r}\left(\eta_{\rm h}\partial_{... ...ight)\frac{\chi_{0}^{l}}{r^3}\right]\vec T_{l}^{0}\left(\theta\right)\right\}},$$ (29)

where we recall that $\Delta_{l}$ is the Laplacian operator: $\Delta_{l}=\partial_{r,r}+({2}/{r})\partial_{r}-{l(l+1)}/{r^2}$.

The task is more difficult with the non-linear expressions describing the stretching and advection of the magnetic field: $\vec\nabla\wedge\left[\vec{\mathcal U}_{\varphi}\wedge\vec B\right]+\vec\nabla\wedge\left[\vec{\mathcal U}_{\rm M}\wedge\vec B\right]$. The first term corresponds to the generation of toroidal field trough the shearing of the poloidal component by the differential rotation, while the second one represents the advection of the field by the meridional circulation.

The first step is to expand ( $\vec{\mathcal U}_{\varphi} + \vec{\mathcal U}_{\rm M})$ on the $\vec R_{l}^{0}\left(\theta\right)$, $\vec S_{l}^{0}\left(\theta\right)$ and $\vec T_{l}^{0}\left(\theta\right)$:

$$\vec{\mathcal U}_{\varphi}+\vec{\mathcal U}_{\rm M}=\sum_{l=0... ...0}^{l}\left(r\right)\vec T_{l}^{0}\left(\theta\right)\right\}.$$ (30)

From (24), we get immediately for the meridional part:

$$u_{0}^{l}\left(r\right)=\frac{U_{l}\left(r\right)}{\mathcal{N... ...(r\right)=\frac{V_{l}\left(r\right)}{\mathcal{N}_{l}^{0}}\cdot$$ (31)

Next, we expand the rotation law as in Paper I:

$$\Omega\left(r,\theta\right)=\overline{\Omega}\left(r\right)+\... ...nfty}\Omega_{l}^{*}\left(r\right)P_{l}\left(\cos\theta\right),$$ (32)

where $\overline{\Omega}\left(r\right)\!\!=\!\!\int_{0}^{\pi}\Omega\left(r,\theta\right)\sin^{3}\theta{\rm d}\theta/\int_{0}^{\pi}\sin^{3}\theta {\rm d}\theta$ the horizontal functions are given by $Q_{l}\left(\theta\right)=P_{l}\left(\cos\theta\right)-I_{l}$, with
$I_{l}={\int_{0}^{\pi}P_{l}\left(\cos\theta\right)\sin^{3}\theta ~ {\rm d}\theta... ...nt_{0}^{\pi}\sin^{3}\theta ~{\rm d}\theta}=\delta_{l,0}-\frac{1}{5}\delta_{l,2}$. Thus, we get:

$$\left\{ \begin{array}{l@{\quad}l} \Omega_{0}^{*}\left(r\right... ...t)=\Omega_{l}\left(r\right)\hbox{ for }l>0. \end{array}\right.$$ (33)

Using the identity (A.9), we obtain the zonal flow $\vec{\mathcal U}_{\varphi}=\sum_{l>0}w_{0}^{l}\left(r\right){\vec T}_{l}^{0}\left(\theta\right)$, where $w_{0}^{l}\left(r\right)$ is given by

$$w_{0}^{l}(r)=r\left[\frac{D_{l-1}^{0}}{\mathcal{N}_{l-1}^{0}}... ...+1}^{0}}{\mathcal{N}_{l+1}^{0}}\Omega_{l+1}^{*}(r)\right]\cdot$$ (34)

The expression of the w₀^l when we keep only the two first terms of the expansion (32), $\overline{\Omega}$ and $\Omega_{2}$ is given in (E.3).

It remains to project the vector product of the two axisymmetric vectors, ( $\vec{\mathcal U}_{\varphi}+\vec{\mathcal U}_{\rm M}$) and $\vec B$, on the vector harmonics $\vec R_{l}^{0}\left(\theta\right)$, $\vec S_{l}^{0}\left(\theta\right)$ and $\vec T_{l}^{0}\left(\theta\right)$. This task is accomplished in Sect. A.2.3; applying it to the advection term, we get:

$$\left(\vec{\mathcal U}_{\varphi}+\vec{\mathcal U}_{\rm M}\right)\wedge\vec B \sum_{l=1}^{\infty}\left\{\mathcal{X}_{{\rm Ad};l}(r)P_{l}(\cos\t... ...hcal{Z}_{{\rm Ad};l}(r)P_{l}^{1}(\cos\theta)\vec{\widehat{e}}_{\varphi}\right\}$$ =

$$\sum_{l=1}^{\infty}\left\{\left[\frac{\mathcal{X}_{{\rm Ad};l}(r)... ...Z}_{{\rm Ad};l}(r)}{\mathcal{N}_{l}^{0}}\right]\vec T_{l}^{0}(\theta)\right\} ,$$ = (35)

where the radial functions $\mathcal{X}_{{\rm Ad};l}\left(r\right)$, $\mathcal{Y}_{{\rm Ad};l}\left(r\right)$ and $\mathcal{Z}_{{\rm Ad};l}\left(r\right)$, which are explicit functions of $\overline{\Omega}$, $\Omega_{l}$, U_l, $\xi_{0}^{l}$ and $\chi_{0}^{l}$, are given in Appendix C. Finally, using once again Eq. (A.32), we obtain:

 

$$\vec\nabla\wedge\left[\left(\vec{\mathcal U}_{\varphi}+\vec{\mathcal U}_{\rm M}\right)\wedge\vec B\right] \sum_{l=1}^{\infty}\left\{\left[\frac{l(l+1)}{\mathcal{N}_{l}^{0}... ...{\rm Ad};l}\right)\right)\right]\vec T_{l}^{0}\left(\theta\right)\right\}}\cdot$$ = (36)

We are now ready to put in their final form the one-dimensional advection/diffusion equations respectively for the poloidal $\xi_{0}^{l}$ and for the toroidal magnetic stream-functions $\chi_{0}^{l}$, by collecting the $\vec R_{l}^{0}\left(\theta\right)$ and the $\vec T_{l}^{0}\left(\theta\right)$ components of Eqs. (27)-(29) and (36):

$$\frac{{\rm d}\xi_{0}^{l}}{{\rm d}t}=\frac{1}{\mathcal{N}_{l}^... ...};l}+\eta_{\rm h}r\Delta_{l}\left(\frac{\xi_{0}^{l}}{r}\right)$$ (37)

$$\frac{{\rm d}\chi_{0}^{l}}{{\rm d}t}+\partial_{r}\left(\dot{r... ...\right)-\eta_{\rm v}l(l+1)\frac{\chi_{0}^{l}}{r^2}\right]\cdot$$ (38)

Note that these equations involve the Lagrangian time-derivative, which makes them suitable for their implementation in stellar evolution codes. They are the equivalent in our formalism of the two classical equations for $\vec{B}_{\rm P}\!=\!\vec\nabla\wedge\vec A$, $\vec A\!=\!A\vec{\widehat e}_{\varphi}$ being the potential vector, and $\vec{B}_{\rm T}\!=\!B_{\rm T}\vec{\widehat e}_{\varphi}$:

$$\partial_{t}A+\frac{1}{s}\vec{\mathcal{U}}_{\rm M}\cdot\vec\n... ...ga+\eta\left(\nabla^2 B_{\rm T}-\frac{B_{\rm T}}{s^2}\right) ,$$ (39)

where $s=r\sin\theta$ and $\eta$ is taken isotropic and uniform (Campbell 1997; Mestel 1999).

Finally, we note that in the ideal case, in the absence of shear, meridional circulation and Ohmic diffusion, Eqs. (37) and (38) reduce to:

$$\frac{{\rm d}\xi_{0}^{l}}{{\rm d}t}=0\hbox{ and }\frac{{\rm d... ...ac{{\rm d}}{{\rm d}t}\left(\frac{B_{\varphi}}{r\rho}\right)=0.$$ (40)

These equations express the Lagrangian flux conservation of respectively

• the magnetic field through the sphere of radius r, as the star expands or contracts;
• the toroidal field in a homothetic contraction or expansion where the density varies as r^-3 (cf. Cowling 1957).

In Appendix E.1 we give the explicit expressions of Eqs. (37) and (38) for the dipole, the quadrupole and the octupole ( $l=\left\{1,2,3\right\}$), retaining the lowest order terms in the rotation law: $\Omega(r,\theta)=\overline{\Omega}\left(r\right)+\Omega_2 \left(r\right)Q_{2}\left(\theta\right)$, and the associated meridional circulation.

3.4 Ohmic heating

The last point we shall consider for completeness is Ohmic heating, which contributes somewhat to the thermal balance/imbalance. In the case of an anisotropic eddy-magnetic diffusivity, its input rate is given by

$$\mathcal{J}\left(r,\theta\right)=\frac{1}{\mu_{0}}\left[\vert... ...\vec B\right)\right]\cdot\left(\vec\nabla\wedge\vec B\right) ,$$ (41)

which reduces to the classical result $\mathcal{J}={{\vec j}^{2}}/{\sigma}$, with $\sigma={1}/\left({\mu_{0}\eta}\right)$ being the conductivity, when the diffusivity is isotropic.

The scalar product of two axisymmetric vectors is treated in Sect. A.2.3; applying the result to $\mathcal{J}$, we get:

$$\mathcal J\left(r,\theta\right)=\sum_{l>0}\mathcal{J}_{l}\left(r\right)P_{l}\left(\cos\theta\right),$$ (42)

where the radial functions $\mathcal{J}_{l}(r)$, which are functions of the magnetic stream-functions $\xi_{0}^{l}$ and $\chi_{0}^{l}$, are given in Appendix D, and in Appendix E for the special case where only the $l=\left\{1,2,3\right\}$ modes are kept. Note again here that if the expansion of $\vec B$ is cut at mode $l_{\rm max}=N_{\rm m}$ in (11), then the expansion of $\mathcal{J}$ involves terms up to $l_{\rm max}=2N_{\rm m}$.

We are now ready to introduce the Lorentz force in the equations governing the transport of momentum and of heat.

4 Transport of angular momentum

We start from the momentum equation

$$\rho\left[ \partial_{t}\vec V+\left( \vec V \cdot \vec \nabla... ...ert\vert\tau\vert\vert+\vec{\mathcal{F}}\!\!_{\mathcal{L}} ~ ,$$ (43)

where $\rho$ is the density, $\phi$ the gravitational potential, $\vert\vert\tau\vert\vert$ the turbulent stresses and $\vec{\mathcal F}\!\!_{\mathcal L}$ the Lorentz force. We insert in it the expression for the macroscopic velocity $\vec V$ given above in Eqs. (22), (23), and take its azimuthal component, which yields an advection/diffusion equation for the angular momentum density:

$$\rho {{\rm d} \over {\rm d}t} \left( r^{2} \sin^{2} \theta ~ ... ...amma_{\vec{\mathcal F}\!\!_{\mathcal L}}\left(r,\theta\right).$$ (44)

Note that here again we have introduced the Lagrangian time derivative, making use of the anelastic continuity equation. This equation is of course identical to that given in Paper I (Eq. (14)), except for the magnetic torque $\Gamma_{\vec{\mathcal F}\!\!_{\mathcal L}}\left(r,\theta\right)= r \sin \theta ~ \vec e _\varphi \cdot \vec{\mathcal F}\!\!_{\mathcal L}\left(r,\theta\right)$. As in Zahn (1992), we assume that the effect of the turbulent stresses on the large scale flow is adequately described by an anisotropic eddy viscosity, whose components are $\nu_{\rm v}$ and $\nu_{\rm h}$ respectively in the vertical and horizontal directions. In the absence of circulation, turbulence and magnetic field, we retrieve the Lagrangian conservation of angular momentum.

It remains to project this equation on spherical harmonics; the expansions of the meridional circulation and the angular velocity were established in Paper I and they have been recalled above in Eqs. (24), (32). We proceed likewise for the magnetic torque:

$$\Gamma_{\vec{\mathcal F}\!\!_{\mathcal L}}\left(r,\theta\righ... ...thcal{F}}_{\mathcal L};l}(r)P_{l}^{1}\left(\cos\theta\right) ,$$ (45)

where $\mathcal{F}_{\mathcal{L},\varphi}$ is the azimuthal component of $\vec{\mathcal F}\!\!_{\mathcal L}$ obtained from (18). Then, using the following property of the associated Legendre functions:

$$P_{l}^{1}\left(\cos\theta\right)=-\frac{{\rm d}}{{\rm d}\thet... ...}{{\rm d}\mu}P_{l}\left(\mu\right)\hbox{ where }\mu=\cos\theta$$ (46)

and the expression for ${{\rm d}P_{l}\left(\mu\right)}/{\rm d\mu}$ which is given in Eq. (A.45), we put $\Gamma_{\vec{\mathcal F}\!\!_{\mathcal L}}$ in its final form:

$$\Gamma_{\vec{\mathcal F}\!\!_{\mathcal L}}\left(r,\theta\righ... ...c{k-1}{2}\right]}\left(2k-4j-1\right)\delta_{l,k-2j-1}\right],$$ (47)

$E\left[x\right]$ is the integer part of x, and $\mathcal{Z}_{\vec{\mathcal F}\!\!_{\mathcal L};k}$ is defined in Eq. (18) and spelled out in explicit form in Appendix B.

4.1 Evolution equation for the angular mean angular velocity

Taking the horizontal average of Eq. (44) over an isobar and using the assumption that $\overline{\Omega}(r) \gg \Omega_{l}(r)$, 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}... ...erline{\Gamma}_{\vec{\mathcal{F}}_{\mathcal L}}\left(r\right),$$ (48)

$$\hbox{where the mean Lorentz torque is} \quad \overline{\Gamm... ...\left(r\right)-\frac{1}{5}\Gamma_{2}\left(r\right)\right]\cdot$$ (49)

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}}_{\rm 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 pointed out in Spiegel & Zahn (1992). The explicit form of Eq. (48), keeping only the $l=\left\{1,2,3\right\}$ terms in the expansion of $\vec B$ (cf. Eq. (11)) is derived in Appendix E.2.1.

4.2 Evolution equation for the differential rotation in latitude

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

$$\rho {{{\rm d} \over {\rm d}t} \left(r^2\sin^2\theta ~ \widehat{\... ...}{r^2}\partial_{r}\left(\rho\nu_{\rm v}r^4\partial_{r}\widehat{\Omega}\right) }$$

$$+\frac{1}{\sin\theta}\partial_{\theta}\left(\rho\nu_{\rm h}\sin^3... ...thcal L}}-\sin^2\theta ~ \overline{\Gamma}_{\vec{\mathcal F}\!\!_{\mathcal L}}.$$ (50)

In the advection term we have again neglected the fluctuation $\widehat{\Omega}(r,\theta)$ compared to the mean $\overline{\Omega}$.

The next step is to replace $\widehat{\Omega}(r,\theta)$ by its expansion (32) 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\rh... ...al_{r}\Omega_{2}\right)-10\rho\nu_{\rm h}\Omega_{2}+\Gamma_{2}$$ (51)

$$\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}\cdot$$ (52)

It can be simplified by assuming that the turbulent transport is 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\rh... ...}-\alpha U_{2}\right]=-10\rho\nu_{\rm h}\Omega_{2}+\Gamma_{2}.$$ (53)

In the asymptotic regime $t \gg {r^2}/{\nu_{\rm h}}$, a stationary state is reached:

$$\nu_{\rm h}\Omega_{2}=\frac{1}{5}r\left[2V_{2}-\alpha U_{2}\right]\overline{\Omega}+\frac{\Gamma_{2}}{10\rho} ,$$ (54)

where horizontal diffusion balances horizontal advection and the torque due to the Lorentz force.

For $l\!\!>\!\!2$ the situation is more intricate, because there are couplings between terms of different l, which prevent a clean separation for each l. This is mainly due to the magnetic torque. Indeed, we recall that in the hydrodynamical case the hypothesis of Spiegel & Zahn (1992) allows such separation (cf. Paper I, Sect. 3.2). Therefore, we choose here to stop the expansion of the rotation law at $\Omega_{2}$. The explicit form of Eq. (53), keeping only the $l=\left\{1,2,3\right\}$ terms in the expansion of $\vec B$ (cf. Eq. (11)) is derived in Appendix E.2.1.

5 Structural properties of the differentially rotating magnetic star

5.1 Baroclinic relation

As in Paper I, we consider a non-uniform and a non-cylindrical rotation law, and add here a general axisymmetric magnetic field. Then neither the centrifugal force nor the Lorentz force derive from a potential, and therefore the isobars and the surfaces of constant density in general do not coincide. In order to determine 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... ...cal{\rm C}}+\frac{\vec{\mathcal F}_{\mathcal{L},\rm P}}{\rho},$$ (55)

where $\phi$ is the gravitational potential and $\vec g$ the local effective gravity, which includes both the centrifugal force
$\vec{\mathcal F}_{\mathcal{\rm C}}= \frac{1}{2}{\Omega}^{2}\vec\nabla\left(r^{2}\sin^{2}\theta\right)$ and the meridional magnetic force $\vec{\mathcal{F}}_{\mathcal L,\rm P}$. Taking the curl of this equation, we get

$$-\frac{1}{{\rho}^{2}}\vec\nabla\rho\wedge\vec\nabla P=-\frac{... ...left(\frac{\vec{\mathcal F}_{\mathcal{L},\rm P}}{\rho}\right),$$ (56)

which to first order reduces to

$$-\frac{\vec\nabla {\rho}{'}\wedge\vec g}{\overline{\rho}}=\le... ...vec{{\mathcal F}}_{\mathcal L,\rm P}}{\overline{\rho}}\right),$$ (57)

with $\rho' (r, \theta)$ being the variation of the density on the isobar.

We now expand all terms in spherical harmonics. For the density fluctuation this is readily done:

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

Then, differentiating the equation of state:

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

we can write the modal amplitude of the density fluctuation as

$${ \widetilde{\rho_{l}} \over \overline\rho} = \Theta_l = \varphi \Lambda_l - \delta \Psi_l .$$ (60)

For the centrifugal force, we recall that

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

which leads us to

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

keeping only the terms linear in $\Omega_{l}$ (beside $\overline{\Omega}^{2}$). Finally, we draw the expansion of the poloidal Lorentz force from (18):

$$\vec{\mathcal{F}}\!\!_{\mathcal{L},\rm P}\left(r,\theta\right... ...l}(r)P_{l}^{1}(\cos\theta)\vec{\widehat{e}}_{\theta}\right\} .$$ (63)

It remains to insert these expansions in Eq. (57), noting that it projects itself only on the azimuthal vectorial spherical harmonics $\vec T_{l}^{0}\left(\theta\right)$. Using the algebra related with these harmonics, we reach the following expression for the modal amplitudes of the relative fluctuation of density on an isobar:

$$\frac{\widetilde{\rho_{l}}\left(r\right)}{\overline\rho}=\var... ...ft(r\right)}{\overline{\rho}\left(r\right)}\right)\right]\cdot$$ (64)

$\overline{g}$ is the horizontal average of the modulus of $\vec g$ and $\mathcal{D}_{l}\left(r\right)$ has been derived in Paper I (Eq. (47)):

 

$$\mathcal{D}_{l}(r) \mathcal{N}_{l}^{0}\left\{r\partial_{r}\left[\overline{\Omega}^{2... ...A_{s}^{0}\left(-C_{s-1}^{0}\delta_{l,s-2}+D_{s-1}^{0}\delta_{l,s}\right)\right.$$ =

$$\left.+B_{s}^{0} \left(-C_{s+1}^{0}\delta_{l,s}+D_{s+1}^{0}\delta... ...left(-C_{s-1}^{0}\delta_{l,s-2}+D_{s-1}^{0}\delta_{l,s}\right)\right]\right\}},$$ (65)

where all the numerical coefficients ( A⁰_l, B⁰_l, C⁰_l, D⁰_l, G⁰_l, H⁰_l) are given in Appendix A.

This baroclinic Eqs. (64)-(65) plays a key role in linking the density fluctuation on an isobar with the differential rotation and the magnetic field. It allows us to close the system formed by the induction equation, that 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 and perturbed potential

It remains to examine how the redistribution of mass we just described modifies the local gravity. In the absence of magnetic field, the perturbing force is just the centrifugal force $\vec{\mathcal F}_{\mathcal{\rm C}} =\frac{1}{2}\Omega^{2}\vec\nabla(r^2\sin^2\theta)$, whose r and $\theta$ components are expanded as

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

where 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) =

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

- 1pc

 

$$\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) =

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

To this we have here to add the meridional components of the Lorentz force

$${\mathcal F}_{{\mathcal L},{\rm P};r}\left(r,\theta\right)=\s... ...}\left(r\right)\partial_{\theta}P_{l}\left(\cos\theta\right).$$ (69)

In the expressions that we have derived in Paper I, it suffices then replace $\rho_0 a_l$ by $\rho_0 a_l +\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}$ and $\rho_0 b_l$ by $\rho_0 b_l + \mathcal{Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}$ to establish the version including the effect of the magnetic field. Thus for the gravity perturbation along an isobar we now get (cf. Paper I, Eq. (72))

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

The last term involves the fluctuation of the gravity field along the sphere, $\widehat{\phi}_{l}$, which is obtained by integrating the Poisson equation (cf. Paper I, Eq. (60)), modified along the same rules:

$$\frac{1}{r}\frac{{\rm d}^2}{{\rm d}r^2}\left(r\widehat{\phi}_... ...Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}\right)\right] \cdot$$ (71)

Let us recall that Sweet (1950) was the first to establish this result for the most general perturbing force.

6 Thermal imbalance and transport of heat

It remains to implement the magnetic field in the heat equation:

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

where S is the entropy per unit mass, $\chi$ the thermal conductivity, $\epsilon$ the nuclear energy production rate per unit mass and $\vec F_{\rm h}$ the flux carried in the horizontal direction by the anisotropic turbulence. Note that in a medium of varying composition, we have to take into account the entropy of mixing (cf. Maeder & Zahn 1998). In the simplest case, applicable to main-sequence stars, where the stellar material can be approximated by a mixture of hydrogen and helium with a fixed abundance of metals, it can be expressed in terms of the mean molecular weight only; then we have

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

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.

The magnetic field manifests itself in the heat Eq. (72) through the Ohmic heating term $\mathcal{J}$, but also in the divergence of the thermal flux, because it involves the divergence of the perturbing force. This force includes both the centrifugal force $\vec{\mathcal{F}}_{\mathcal{\rm C}}$ and the Lorentz force per unit volume, ${\vec{\mathcal F}_{\mathcal{L},\rm P}}/{\rho_{0}}$, and their divergence will again be expanded in spherical harmonics as

$$\vec\nabla\cdot\vec{\mathcal{F}}_{\mathcal{\rm C}}=\overline{... ...l>0}\widetilde{f}_{\mathcal{L},l}P_{l}\left(\cos\theta\right).$$ (74)

In Paper I (Eq. (85)) we gave the expressions of $\overline{f}_{\mathcal{\rm C}}$ and $\widetilde{f}_{\mathcal{\rm C},l}$ in terms of a_l and b_l:

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

As above in Sect. 5.2, it suffices to replace a_l by $\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}/\rho_0$ and b_l by $\mathcal{Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}/\rho_0$ to obtain the equivalent expressions for the divergence of the Lorentz force:

$$\overline{f}_{\mathcal{L}}=\frac{1}{r^2}\partial_{r}\left(r^2... ...mathcal{Y}_{\vec{\mathcal F}_{\mathcal{L};l}}}{r\rho_{0}}\cdot$$ (76)

This leads us to the modal form of the heat equation, which can be implemented directly into a stellar structure code:

$$\overline{T} C_{\rm p}\left[{{\rm d} \Psi_{l} \over {\rm d} t... ..._{l}(r)+\frac{\mathcal{J}_{l}\left(r\right)}{\overline{\rho}},$$ (77)

where $\mathcal{T}_{l}$ is given by:

 

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

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

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

We recall that L is the luminosity, M the mass, $\overline{T}$ the horizontal average of the temperature, $C_{\rm p}$ the specific heat at constant pressure, and $\nabla$ the radiative gradient. We have also introduced the temperature scale-height $H_{\rm T}=\left\vert{{\rm d}r}/{{\rm d}\ln \overline{T}}\right\vert$, the thermal diffusivity $K={\overline{\chi}}/{\overline{\rho}C_{\rm p}}$, the horizontal eddy-diffusivity $D_{\rm h}$ and $f_{\epsilon}={\overline{\epsilon}}/{\left(\overline{\epsilon}+\overline{\epsilon}_{\rm grav}\right)}$, with $\overline{\epsilon}$ and $\overline{\epsilon}_{\rm grav}$ being respectively the mean nuclear and gravitational energy release rates, whereas $\epsilon_{\mu}$ and $\chi_{\mu}$ are the logaritmic derivatives of $\epsilon$ and of the radiative conductivity $\chi$ with respect to $\mu$, their derivatives with respect to T being noted as $\epsilon_{\rm T}$ and $\chi_{\rm T}$. Moreover, we have $\epsilon_{\rm m}=L\left(r\right)/M\left(r\right)$ and $\rho_{\rm m}$ is the mean density inside the considered level surface.

Three remarks before we conclude this section.

• First, we have noted previously that if the expansion of $\vec B$ is stopped at mode $l_{\rm max}=N_{\rm m}$ in (11), then the expansions of $\vec{\mathcal F}\!\!_{\mathcal L}$ and of $\mathcal{J}$ must include all terms up to $l_{\rm max}=2N_{\rm m}$, due to selection rules. Therefore, keeping $\Omega\left(r,\theta\right)=\overline{\Omega}\left(r\right)+\Omega_{2}\left(r\right)Q_{2}\left(\theta\right)$, the expansion of the meridional flow will extend to $l_{\vec{\mathcal U}_{\rm M};{\rm max}}={\rm max}\left\{4,2N_{\rm m}\right\}$.
• Next, the relative importance of the terms contributed respectively by the magnetic force and by the centrifugal force depends on the ratio $\left({V_{\rm A}}/{V_{\Omega}}\right)^{2}$ where $V_{\rm A}={B}/{\sqrt{\mu_{0}\rho_{0}}}$ is the Alfvén speed and $V_{\Omega}=R\Omega$, R being the radius of the star. If the magnetic field is weak enough, this ratio will be small except just below the surface because of the decay of the density. In this case, for the implementation in existing stellar evolution codes of the equations related to the thermal imbalance and to the structural properties of the star (cf. Sect. 5), we keep only the modes $l=\left\{2,4\right\}$, which allows to describe the tachocline(s) circulation to first order (see Appendix E). However, in the case of a strong field higher order modes of the Lorentz force must be taken into account.
• Finally, the ratio between the term representing the Ohmic heating, ${\mathcal{J}_{l}\left(r\right)}/{\overline{\rho}}$, and the highest-order derivative term of $\vec\nabla\cdot\left(\chi\vec\nabla T\right)$ in Eqs. (77) and (78), namely $({L}/{M}) ({\rho_{\rm m}}/{\overline{\rho}}) ({r}/{3}) \partial_{r}\left(H_{\rm T}\partial_{r}\Psi_{l}\right)$, is given by $({\eta}/{K}) \left({{V_{\rm A}}}/{V_{\Omega}}\right)^{2}$. Therefore the Ohmic term can be neglected since $V_{\rm A} < R \Omega$ and $\eta \ll K$.

In conclusion, the link between the circulation and its cause, namely the thermal imbalance due to the rotation, the magnetic field and the chemical composition, is established through Eqs. (77) and (78) where the temperature fluctuation on an isobar is governed by an advection/diffusion equation from which we can derive the radial component of $\vec{\mathcal{U}}_{\rm M}$. These equations have been established for a rotation law which depends on r and $\theta$, and that allows us to treat simultaneously the bulk of the radiation zones and the tachoclines in presence of a general axisymmetric magnetic field.

7 Boundary conditions

The system of equations is now complete: we have the induction equation, an advection/diffusion equation for the transport of angular momentum (mean and fluctuating) including the action of the magnetic torque, another for the temperature (mean and fluctuating), and the baroclinic relation which allows us to close the system. The equation for the transport of chemical species (or alternatively for the transport of molecular weight) is unchanged. It thus remains to specify the boundaries conditions of this system, to be applied on the limits of the radiation zone. To be specific, we consider a star with a radiation zone located between a convective core and an upper convection zone. We designate by $r_{\rm b}$ and $r_{\rm t}$ the respective radius of the base and of the top of that radiative zone. Of course, in a solar-type main-sequence star we have $r_{\rm b}=0$, whereas for a massive main sequence star $r_{\rm t}=R$.

The boundary conditions for the equations already present in the hydrodynamic case, which was examined in Paper I, are unchanged, except those to be applied to the equation of angular momentum transport, since it includes the Lorentz torque. The novelty here is the induction equation, which we have split in two equations, respectively for the poloidal and toroidal magnetic stream-functions $\xi_{0}^{l}$ (Eq. (37)) and $\chi_{0}^{l}$ (Eq. (38)). These equations are of second order in r, and therefore each of them requires two boundary conditions. In principle they can be obtained by solving the induction equation, together with the equation of momentum, etc. in each adjacent convective zone, but this represents a formidable task which is well beyond our scope here. Instead, we make the following - simplifying but reasonable - hypothesis, valid only in the convection zone(s):

(i)

we do not take into account the meridional velocity field. The only motion we retain is the zonal flow associated with the differential rotation: $\vec{\mathcal U}_{\varphi}=r\sin\theta ~ \Omega_{\rm CZ}\left(r,\theta\right)\vec{\widehat{e}}_{\varphi}$, with $\Omega_{\rm CZ}(r, \theta)$ being the angular velocity in the convection zone;

(ii)

we put $\eta_{\rm v}=\eta_{\rm h}=\eta_{\rm {turb}}$ where $\eta_{\rm {turb}}$ is the magnetic eddy-diffusivity, which we assume constant inside the considered convective region;

(iii)

this eddy-diffusivity is high enough to allow for a stationary state (compared to the slow evolution in the radiation zone).

With these assumptions, the poloidal field is potential, and Eq. (37) reduces to:

$$\frac{{\rm d^{2}}\xi_{0}^{l}}{{\rm d}r^{2}}-\frac{l\left(l+1\right)}{r^2}\xi_{0}^{l}=0;$$ (79)

it admits the general solution:

$$\xi_{0}^{l}=Ar^{l+1}+\frac{B}{r^{l}}\cdot$$ (80)

When applying it to a convective core where we must take the non-singular solution at the origin:

$$r\frac{{\rm d}\xi_{0}^{l}}{{\rm d}r}-\left(l+1\right)\xi_{0}^{l}=0\quad \hbox{at} \; r=r_{\rm b} ,$$ (81)

which becomes $\xi_{0}^{l}=0$ if $r_{\rm b}=0$. In the case of a convective envelope, surrounded by vacuum, we have

$$r\frac{{\rm d}\xi_{0}^{l}}{{\rm d}r}+l\xi_{0}^{l}=0\quad \hbox{at} \; r=r_{\rm t}.$$ (82)

With the same assumptions, the toroidal field obeys the simplified form of Eq. (38):

$$\frac{{\rm d^{2}}\chi_{0}^{l}}{{\rm d}r^{2}}-\frac{l\left(l+1... ...l{Y}_{{\rm Ad};l}\left(\Omega_{\rm {CZ}}\right)\right)\right].$$ (83)

$\mathcal{X}_{{\rm Ad};l}\left(\Omega_{\rm {CZ}}\right)$ and $\mathcal{Y}_{{\rm Ad};l}\left(\Omega_{\rm {CZ}}\right)$ are derived using the previous multipolar solutions for $\xi_{0}^{l}$ (cf. Eq. (80)), expanding $\Omega_{\rm {CZ}}$ in spherical functions like $\Omega_{\rm {CZ}}\left(r,\theta\right)=\sum_{l}\Omega_{l}\left(r\right)P_{l}\left(\cos\theta\right)$ where the $\Omega_{l}$ are taken from the results of the inversion of helioseismological data (Corbard et al. 2002), using (34) and Appendix C. Then, the solution of (83) inside a convection zone is derived using Green's functions:

$$\chi_{0}^{l}=\frac{1}{2l+1} {1 \over \eta_{\rm {turb}}} \left... ...mathcal{S}_{l}\left(\Omega_{\rm {CZ}}\right) ~ {\rm d}x\right]$$ (84)

where $R_{{\rm int}\rm CZ}$ and $R_{\rm {extCZ}}$ are respectively the radius at its base and at its top. Applying Eq. (84) to a convective core where $R_{{\rm int}\rm CZ}=0$ and $R_{\rm {extCZ}}=r_{\rm b}$, we get:

$$\chi_{0}^{l}=\frac{1}{2l+1} {1 \over \eta_{\rm {turb}}} r_{\r... ...a_{\rm {CZ}}\right) ~ {\rm d}x \quad \hbox{at} \; r=r_{\rm b}.$$ (85)

Finally, for a convective envelope where $R_{{\rm int}CZ}=r_{\rm t}$ and $R_{\rm {extCZ}}=R$, we obtain:

$$\chi_{0}^{l}=\frac{1}{2l+1}{1 \over \eta_{\rm {turb}}} r_{\rm... ...ga_{\rm {CZ}}\right) ~ {\rm d}x\quad \hbox{at} \; r=r_{\rm t}.$$ (86)

In the case where $\mathcal{S}_{l}\left(\Omega_{\rm {CZ}}\right)$ is negligible, the previous results leads to:

$$\chi_{0}^{l}=0\hbox{ }\hbox{at}\hbox{ }r=r_{\rm b} \quad \hbo... ...\rm t}, \quad \hbox{and thus to} \quad \vec B_{\rm T}=\vec 0 .$$ (87)

It remains to write down the boundary conditions to be applied on the transport of angular momentum. Since the equation for the mean angular velocity $\overline{\Omega}$ (48) is of second order in r, 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_{\rm b}}r^{4}\rho... ...athcal{F}_{\rm B}(r_{\rm t}) \quad \hbox{at} \; r=r_{\rm t} .$$ (88)

$\mathcal{F}_{\Omega}$ is the (signed) flux of angular momentum which is lost at the surface by the stellar wind, and $\mathcal{F}_{\rm B}(r)$ the flux carried by the magnetic field through the considered surface:

$$\mathcal{F}_{\rm B}(r) = -{1 \over \mu_0} r^3 \int_0^\pi B_r ... ...l-1\right)\left(2l+1\right)}\chi_{0}^{l-1}\right]\right\}\cdot$$ (89)

Here again the perturbation $\Omega_{2}$ obeys an evolution equation which does not include any derivative in r, and therefore it needs no boundary condition.

8 Conclusion

The work presented here is the continuation of Paper I (Mathis & Zahn 2004), where the star was assumed without magnetic field. Here we allow for a magnetic field of moderate strength in the radiation zone(s), with an Alfvén speed not exceeding the rotational velocity: $(V_{\rm A})^2 < (R\Omega)^2$. Such a field has little impact on the hydrostatic balance, since we also assume, as in Paper I, that the centrifugal force is small compared to gravity: $R \Omega^2 \ll g$. But it will compete with the centrifugal force in the baroclinic balance (56), and therefore it will participate in governing the meridional flow. Moreover, such a field tends to enforce uniform rotation along the field lines of the poloidal field (Ferraro's law). The question then arises whether or not the poloidal field threads into the convection zone(s), where such differential rotation is probably maintained through the turbulent motions, as observed in the Sun. To answer that question one has to treat consistently the evolution of both angular velocity and magnetic field, with special care for the dynamics in the tachocline(s), and for this reason we expanded them in spherical harmonics to arbitrary order.

The overall problem of rotational mixing in a magnetized star is highly non-linear, with multiple feed-backs, as illustrated by the diagram displayed in Fig. 1.

	Figure 1: Rotational mixing in stellar interiors with magnetic field: a highly non-linear problem (the action of the convection through penetration and overshoot which enhance the mixing has been added).
Open with DEXTER

Mixing is achieved through the meridional circulation, which is due to the thermal imbalance caused by the centrifugal force and by the Lorentz force (mainly through the baroclinic balance), and also through the turbulence generated by the shear of differential rotation and possibly by the magnetic field. These motions modify both the rotation profile and the magnetic field through large-scale advection and turbulent diffusion. As the star evolves, molecular gradients build up through microscopic diffusion (including radiative levitation and gravitational settling), and through nuclear burning. Large-scale advection turns them into horizontal gradients, which react on the meridional circulation and are smoothed by the turbulence.

The boundary conditions play a crucial role: in a non-magnetic star, the circulation is driven by the loss (or gain) of angular momentum, as explained in Zahn (1992); in the absence of such loss, the circulation tends to vanish, as was shown by Busse (1982). It remains to be seen whether this is still the case when angular momentum is transported mainly through magnetic stresses, with the poloidal field anchored in the convection zone(s).

The main weakness of our modeling remains the description of the turbulence. As in Paper I, it is assumed to be anisotropic, due to the stable stratification, and that it tends to smooth angular velocity and chemical composition on horizontal surfaces. The action of magnetic field on such turbulence is not taken in account, which seems reasonable as long as the Alfvén speed does not exceed the rotational velocity. Moreover, we have deliberately ignored here the possibility of MHD instabilities generating their own turbulence (Balbus & Hawley 1994; Spruit 1999).

From the technical point of view, we adopted here the formalism developed in Paper I and introduced the spherical vectorial harmonics, which obey the Racah-Wigner angular momentum algebra widely used in quantum mechanics. This allowed us to separate the colatitude $\theta$ in the vectorial partial differential equations which govern the problem, and to reduce the problem to solve coupled partial differential equations in $t\hbox{ and }r$ only. It is then straightforward to introduce these equations in existing stellar structure codes, a task we have now undertaken.

Among the first applications we plan to model a solar-type star, to check whether a fossil magnetic field, which is able to extract angular momentum from the radiative interior as the star spins down, yields a profile of angular velocity that is compatible with the helioseismic data. Next we will undertake the modeling of magnetized massive stars, since they are progenitors of neutron stars which are the seat of strong magnetic fields (Maeder & Meynet 2003, 2004; Heger et al. 2004).

Most of the appendix is dedicated to the algebra involving the spherical and vectorial spherical harmonics used in this paper. The numerical values of the coupling coefficients given in the tables have been computed with Mathematica. In Appendix E, we state the equations in a form ready to be implemented in a stellar structure code, when the magnetic field is expanded up to the octupole (l=3) and only the first term $\Omega_{2}$ of the differential rotation is retained.

Acknowledgements

This work was supported by the Centre National de la Recherche Scientifique (Programme National de Physique Stellaire and GdR Dynamo).

References

• Balbus, S. A., & Hawley, J. F. 1994, MNRAS, 266, 769 [NASA ADS] (In the text)
• Braithwaite, J., & Spruit, H. C. 2004, Nature, 431, 819 [NASA ADS] [CrossRef] (In the text)
• Bullard, E. C., & Gellman, H. 1954, Phil. Trans. R. Soc. London A, 247, 213 [NASA ADS] (In the text)
• Busse, F. H. 1982, ApJ, 259, 759 [NASA ADS] [CrossRef] (In the text)
• Campbell, C. G. 1997, Magnetohydrodynamics in binary stars, Astrophysics & Space Science Library (Kluwer Academic Publishers), 216 (In the text)
• Corbard, T., & Thompson, M. J. 2002, Sol. Phys., 205, 211 [NASA ADS] [CrossRef] (In the text)
• Cowling, T. G. 1957, Magnetohydrodynamics, Interscience tracts on physics and astronomy (In the text)
• Edmonds, A. R. 1974, Angular Momentum in Quantum Mechanics (Third printing with corrections) (Princeton University Press) (In the text)
• Ferraro, V. C. A. 1937, MNRAS, 97, 458 [NASA ADS] (In the text)
• Garaud, P. 1999, MNRAS, 304, 583 [NASA ADS] [CrossRef] (In the text)
• Garaud, P. 2002, MNRAS, 329, 1 [NASA ADS] [CrossRef] (In the text)
• Gradshteyn, I. S., & Ryzhik, I. M. 1980, Table of Integrals, Series and Products - Corrected and Enlarged Edition (Academic Press)
• Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] (In the text)
• James, R. W. 1973, Proc. R. Soc. London A, 331, 469 [NASA ADS] (In the text)
• James, R. W. 1974, Proc. R. Soc. London A, 340, 287 [NASA ADS] (In the text)
• Kumar, P., Talon, S., & Zahn, J.-P. 1999, ApJ, 520, 859 [NASA ADS] [CrossRef] (In the text)
• Maeder, A., & Zahn, J.-P. 1998, A&A, 334, 1000 [NASA ADS] (In the text)
• Maeder, A., & Meynet, G. 2000, ARA&A, 38, 143 [NASA ADS] (In the text)
• Maeder, A., & Meynet, G. 2001, A&A, 373, 555 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Maeder, A., & Meynet, G. 2003, A&A, 411, 543 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Maeder, A., & Meynet, G. 2004, A&A, 422, 225 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Matias, J., & Zahn, J.-P. 1997, ed. J. Provost, & F.-X. Schmider; poster volume, IAU Symp., 181, 103 (In the text)
• Mathis, S., & Zahn, J.-P. 2004, A&A, 425, 229 [EDP Sciences] [NASA ADS] [CrossRef] (Paper I) (In the text)
• Mestel, L. 1953, MNRAS, 113, 716 [NASA ADS] (In the text)
• Mestel, L., & Moss, D. 1977, MNRAS, 178, 27 [NASA ADS] (In the text)
• Mestel, L., Tayler, R. J., & Moss, D. 1988, MNRAS, 231, 873 [NASA ADS] (In the text)
• Mestel, L. 1999, in International series of monographs on Physics, 99, Stellar Magnetism, ed. J. Birman, S. F. Edwards, R. H. Friend, C. H. Llewellyn Smith, M. Rees, D. Sherrington, & G. Veneziano (Oxford: Clarendon Press) (In the text)
• Moss, D. 1974, MNRAS, 168, 61 [NASA ADS] (In the text)
• Moss, D. 1977, MNRAS, 178, 51 [NASA ADS] (In the text)
• Moss, D. 1984, MNRAS, 209, 607 [NASA ADS] (In the text)
• Moss, D. 1985, MNRAS, 213, 575 [NASA ADS] (In the text)
• Moss, D. 1987, MNRAS, 226, 281 [NASA ADS] (In the text)
• Rieutord, M. 1987, Geophys. Astrophys. Fluid Dynamics, 39, 163 (In the text)
• Serebrianaya, P. M. 1988, Geophys. Astrophys. Fluid Dynamics, 44, 141 (In the text)
• Spiegel, E., & Zahn, J.-P. 1992, A&A, 265, 106 [NASA ADS] (In the text)
• Spruit, H. C. 1999, A&A, 349, 189 [NASA ADS] (In the text)
• Sweet, P. A. 1950, MNRAS, 110, 548 [NASA ADS] (In the text)
• Talon, S., Zahn, J.-P., Maeder, A., & Meynet, G. 1997, A&A, 322, 209 [NASA ADS] (In the text)
• Talon, S., Kumar, P., & Zahn, J.-P. 2002, ApJ, 574, L175 [NASA ADS] [CrossRef] (In the text)
• Talon, S., & Charbonnel, C. 2003, A&A, 405, 1025 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Talon, S., & Charbonnel, C. 2004, A&A, 418, 105 [NASA ADS] (In the text)
• Tassoul, J.-L., & Tassoul, M. 1986, ApJ, 310, 786 [NASA ADS] [CrossRef] (In the text)
• Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1975, Quantum Theory of Angular Momentum (World Scientific) (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

A.1 Scalar quantities

A.1.1 Definition and basic properties

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)

where $P_{l}^{\vert m\vert}\left(\cos\theta\right){\rm e}^{{\rm i}m\varphi}$ is the associated Legendre function, and the normalization coefficient being

$$\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)

They obey the orthogonality relation:

$$\int_{\Omega}\left(Y_{l_{1}}^{m_{1}}\left(\theta,\varphi\righ... ...i\right){\rm d}\Omega=\delta_{l_{1},l_{2}}\delta_{m_{1},m_{2}}$$ (A.3)

where ${\rm d}\Omega=\sin\theta ~ {\rm d}\theta ~ {\rm d}\varphi$ and where the complex conjugate spherical harmonic is given by:

$$\left(Y_{l}^{m}\left(\theta,\varphi\right)\right)^{*}=\left(-1\right)^{m}Y_{l}^{-m}\left(\theta,\varphi\right).$$ (A.4)

Using these properties, every function $f(\theta,\varphi)$ can be expanded as:

$$f\left(\theta,\varphi\right)=\sum_{l=0}^{\infty}\sum_{m=-l}^{... ...(Y_{l}^{m}\left(\theta,\varphi\right)\right)^{*}{\rm d}\Omega.$$ (A.5)

A.1.2 Special case of axisymmetric spherical harmonics

In that case we have:

$$Y_{l}^{0}\left(\theta,\varphi\right)=\mathcal{N}_{l}^{0}P_{l}... ...^{1}\left(\cos\theta\right)\hbox{ for }l>0. \end{array}\right.$$ (A.6)

A.1.3 Linear differential and recursion relations

We recall that the spherical harmonics obey to the differential equation:

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

In deriving the functions which are related to the centrifugal force in Sect. 5 (cf. Eqs. (65), (67), (68)) we have used the following recursion relations for m=0:

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

$$\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.9)

$$\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.10)

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

The identities (A.9) and (A.10) have been deduced from the two others (A.8)-(A.11) with the help of (A.7).

A.1.4 Expansion of products of the spherical harmonics

Using the normalization and the orthogonality of spherical harmonics (cf. Eqs. (A.3)-(A.5)) and their complex conjugate (cf. Eq. (A.4)), we can write:

$$Y_{l_{1}}^{m_{1}}\left(\theta,\varphi\right)Y_{l_{2}}^{m_{2}}... ...1}+m_{2}\right)}Y_{l}^{m_{1}+m_{2}}\left(\theta,\varphi\right)$$ (A.12)

where we define the integral $\mathcal{I}_{l_{1},l_{2},l}^{m_{1},m_{2},m}$ like in Edmonds (1968) or Varshalovich et al. (1975):

$$\mathcal{I}_{l_{1},l_{2},l}^{m_{1},m_{2},m}=\int_{\Omega}Y_{l... ...{array}{ccc} l_{1} & l_{2} & l\\ 0 & 0 & 0 \end{array}\right)$$ (A.13)

with the 3j-Wigner coefficients which are related to the classical Clebsch-Gordan coefficients by:

$$\left(\begin{array}{ccc} l_{1} & l_{2} & l \\ m_{1} & m_{2} ... ..._{1}-l_{2}-m}}{\sqrt{2l+1}}C_{l_{1},m_{1},l_{2},m_{2}}^{l,-m}.$$ (A.14)

So finally, we get the following expansion for the product of two spherical harmonics:

$$Y_{l_{1}}^{m_{1}}\left(\theta,\varphi\right)Y_{l_{2}}^{m_{2}}... ...l_{2},m_{2}}^{l}Y_{l}^{m_{1}+m_{2}}\left(\theta,\varphi\right)$$ (A.15)

where the coefficient c_l1,m1,l2,m2^l is given by:

$$c_{l_{1},m_{1},l_{2},m_{2}}^{l}=\left(-1\right)^{\left(m_{1}+... ...array}{ccc} l_{1} & l_{2} & l\\ 0 & 0 & 0 \end{array}\right).$$ (A.16)

Then, using the initial definition of spherical harmonics (cf. Eqs. (A.1) and (A.2)), we deduce the expansion for the product of two associated Legendre functions:

$$P_{l_{1}}^{m_{1}}\left(\cos\theta\right)P_{l_{2}}^{m_{2}}\lef... ...{1},l_{2},m_{2}}^{l}P_{l}^{m_{1}+m_{2}}\left(\cos\theta\right)$$ (A.17)

where

$$d_{l_{1},m_{1},l_{2},m_{2}}^{l}=(-1)^{\left(m_{1}+m_{2}\right... ...array}{ccc} l_{1} & l_{2} & l\\ 0 & 0 & 0 \end{array}\right).$$ (A.18)

Note that the previous expression is symmetric; therefore

d_l1,m1,l2,m2^l=d_l2,m2,l1,m1^l. (A.19)

A.2 Vector fields

A.2.1 Definitions and basic properties

Following Rieutord (1987), we expand any vector field $\vec{u}(r,\theta,\phi)$ in vectorial spherical harmonics as

$$\vec{u}(r,\theta,\phi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\lef... ...a,\varphi)+w_{m}^{l}(r)\vec T_{l}^{m}(\theta,\varphi)\right\},$$ (A.20)

where the vectorial spherical harmonics $\vec R_{l}^{m}\left(\theta,\varphi\right)$, $\vec S_{l}^{m}\left(\theta,\varphi\right)$, $\vec T_{l}^{m}\left(\theta,\varphi\right)$ are defined as:

$$\vec R_{l}^{m}(\theta,\varphi)=Y_{l}^{m}(\theta,\varphi)\vec{... ...\vec\nabla_{\mathcal{S}}\wedge\vec R_{l}^{m}(\theta,\varphi) ,$$ (A.21)

$$\hbox{with the horizontal gradient} \; \vec\nabla_{\mathcal{S... ...c{\widehat e}_{\varphi}\frac{1}{\sin\theta}\partial_{\varphi}.$$ (A.22)

These vector functions obey the following orthogonality relations:

$$\int_{\Omega}\vec R_{l_{1}}^{m_{1}}\cdot\vec S_{l_{2}}^{m_{2}... ...c S_{l_{1}}^{m_{1}}\cdot\vec T_{l_{2}}^{m_{2}}{\rm d}\Omega=0,$$ (A.23)

$$\int_{\Omega}\vec R_{l_{1}}^{m_{1}}\cdot\left(\vec R_{l_{2}}^... ...\Omega=l_{1}(l_{1}+1)\delta_{l_{1},l_{2}}\delta_{m_{1},m_{2}}.$$ (A.24)

The vector function $\vec u$ may also be projected on the classical spherical vectorial basis:

$$\vec u=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left\{u_{m}^{l}(r)Y... ...{m}(\theta,\varphi)\right]\vec{\widehat e}_{\varphi}\right\} ,$$ (A.25)

$\vec{\widehat e}_{r}$, $\vec{\widehat e}_{\theta}$ and $\vec{\widehat e}_{\varphi}$ beeing the unit-vectors respectively in the r, $\theta$ and $\varphi$ directions. These expansions of vector fields allow us to separate explicitly the angular variables $\theta$ and $\varphi$ in the vectorial partial differential equations which govern the problem. We thus reduce the problem to solve partial differential equations in $t\hbox{ and }r$ only, which are easy to implement in existing stellar structure codes. Note that $\vec R_{l}^{m}\left(\theta,\varphi\right)$ and $\vec S_{l}^{m}\left(\theta,\varphi\right)$ represent the poloidal part, and $\vec T_{l}^{m}\left(\theta,\varphi\right)$ the toroidal part of $\vec u$.

A.2.2 Expansions of differential operators

As stated in the previous section, the expansion of the vector fields in $\vec R_{l}^{m}\left(\theta,\varphi\right)$, $\vec S_{l}^{m}\left(\theta,\varphi\right)$ and $\vec T_{l}^{m}\left(\theta,\varphi\right)$ allows us to separate the variables in the vectorial partial differential equations which govern the problem. This will prove particularly useful when dealing with magnetic field, and with the non-linear expressions where it is involved. We start by expanding the classical linear vectorial operators: gradient, divergence, curl and Laplacian (scalar or vectorial).

Gradient:

Taking a scalar function $f(r,\theta,\varphi)$ expanded in the $Y_{l}^{m}\left(\theta,\varphi\right)$:

$$f(r,\theta,\varphi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}f_{m}^{l}(r)Y_{l}^{m}\left(\theta,\varphi\right),$$ (A.26)

it is straightforward to derive its gradient:

$$\vec\nabla f(r,\theta,\varphi)=\sum_{l=0}^{\infty}\sum_{m=-l}... ...hi)+\frac{f_{m}^{l}}{r}\vec S_{l}^{m}(\theta,\varphi)\right\}.$$ (A.27)

Divergence:

Taking a vector field $\vec u$ expanded as in (A.20), its divergence is given by:

$$\vec\nabla\cdot\vec u(r,\theta,\varphi)=\sum_{l=0}^{\infty}\s... ...ht)-l(l+1)\frac{v_{m}^{l}}{r}\right]Y_{l}^{m}(\theta,\varphi),$$ (A.28)

where whe note that, if $\vec u$ is divergence-free, such as the magnetic field $\vec B$ or the momentum density $\rho\vec{\mathcal{U}_{\rm M}}$ in the anelastic approximation, we have the following relation between u_m^l and v_m^l:

$$v_{m}^{l}=\frac{1}{l(l+1)}\frac{1}{r}\partial_{r}\left(r^2 u_{m}^{l}\right).$$ (A.29)

Laplacian of a scalar function:

With these expressions for the gradient of a scalar quantity and for the divergence of a vector field, one can easily derive the Laplacian of a scalar function. Using the well-known property $\nabla^{2}f=\vec\nabla\cdot\left(\vec\nabla f\right)$, we get:

$$\nabla^{2}f=\vec\nabla\cdot(\vec\nabla f)=\sum_{l=0}^{\infty}... ...ty}\sum_{m=-l}^{l}\Delta_{l}f_{m}^{l}Y_{l}^{m}(\theta,\varphi)$$ (A.30)

where $\Delta_{l}$ is the Laplacian operator:

$$\Delta_{l}=\partial_{r,r}+\frac{2}{r}\partial_{r}-\frac{l\left(l+1\right)}{r^2}\cdot$$ (A.31)

Curl:

This operator is found in the expression of the magnetic field in terms of a stream function, in that of the current density, etc. If we take $\vec u$ expanded as in Eq. (A.20), we retrieve Rieutord's (1987) result:

$$\vec\nabla\wedge\vec u=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\lef... ...v_{m}^{l}\right)\right]\vec T_{l}^{m}(\theta,\varphi)\right\}.$$ (A.32)

Taking this result, the expansion for the Laplacian of a vector field is derived. First, we have:

$$\vec\nabla\wedge\left(\vec\nabla\wedge\vec u\right)=\sum_{l=0... ...elta_{l}w_{m}^{l}\right]\vec T_{l}^{m}(\theta,\varphi)\right\}$$ (A.33)

and therefore

 

$$\nabla^{2}\vec u \vec\nabla\left(\vec\nabla\cdot\vec u\right)-\vec\nabla\wedge\left(\vec\nabla\wedge\vec u\right)$$ =

$$\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left\{\left[\Delta_{l}u_{m}^{l... ...varphi)+\left[\Delta_{l}w_{m}^{l}\right]\vec T_{l}^{m}(\theta,\varphi)\right\},$$ = (A.34)

result which becomes in the case where $\vec u$ is divergence-free (cf. Eq. (A.29) ):

 

$$\nabla^{2}\vec u -\vec\nabla\wedge\left(\vec\nabla\wedge\vec u\right)$$ =

$$\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left\{\left[\frac{1}{r}\Delta_... ...varphi)+\left[\Delta_{l}w_{m}^{l}\right]\vec T_{l}^{m}(\theta,\varphi)\right\}.$$ = (A.35)

Finally, using Eq. (A.32) once again, we get the "triple curl'' operator:

$$\left(\vec\nabla\wedge\right)^{3}\vec u=\sum_{l=0}^{\infty}\s... ...elta_{l}z_{m}^{l}\right]\vec T_{l}^{m}(\theta,\varphi)\right\}$$ (A.36)

where:

$$z_{m}^{l}=\frac{1}{r}\partial_{r}(r v_{m}^{l})-\frac{u_{m}^{l}}{r},$$ (A.37)

which becomes in the case where $\vec u$ is divergence-free (cf. Eq. (A.29)):

 

$$\left(\vec\nabla\wedge\right)^{3}\vec u=\sum_{l=0}^{\infty}\sum_{... ...\Delta_{l}\Delta_{l}(r u_{m}^{l})\right]\vec T_{l}^{m}(\theta,\varphi)\right\}.$$ (A.38)

   
A.2.3 Products of axisymmetric vectorial spherical harmonics

Before we turn to the non-linear terms involving the magnetic field: the advection term in the induction equation and the Lorentz force, we shall first perform the projections of the scalar product of two general axisymmetric vectors and of their vector product, and derive the associated coupling coefficients.

Scalar product:

We take two general axisymmetric vectors $\vec X_{1}\left(r,\theta\right)$ and $\vec X_{2}\left(r,\theta\right)$, which we expand as in Eq. (A.20):

$$\left\{ \begin{array}{l@{\quad}l} \vec X_{1}(r,\theta)=\sum_{... ...c T_{l_{2}}^{0}\left(\theta\right)\right\}, \end{array}\right.$$ (A.39)

and perform their scalar product. After some algebra involving the expansion of products of spherical harmonics (cf. Sect. A.1.4), we obtain

$$\vec X_{1}\left(r,\theta\right)\cdot\vec X_{2}\left(r,\theta\... ...\vec X_{2}\right);l}\left(r\right)P_{l}\left(\cos\theta\right)$$ (A.40)

with the following expression for $\mathcal{P}_{\left(\vec X_{1}\cdot\vec X_{2}\right);l}(r)$:

 

$$\mathcal{P}_{\left(\vec X_{1}\cdot\vec X_{2}\right);l}(r)=\sum_{l... ...right)\mathcal{F}_{0}^{l_{2}}\left(r\right)\right]X_{l_{1},l_{2}}^{l}\right\} .$$ (A.41)

Here X_l1,l2^lis defined in terms of the coupling coefficients d_l1,m1,l2,m2^l (A.18):

$$X_{l_{1},l_{2}}^{l}=\sum_{j=I(l_{1},1,l_{2},1)}^{l_{1}+l_{2}}... ...j-1}{2}\right]}\left[\left(2j-4p-1\right)\cdot X\right]\right)$$ (A.42)

with

$$X=\sum_{q=0}^{E\left[\frac{\left(j-2p-1\right)-1}{2}\right]}\... ...j-2p-1)-2q-1\right],0}^{r}\delta_{l,r}\right]\right)\right\} .$$ (A.43)

We have used the classical notation $\delta_{i,j}$ for the usual Kronecker symbol, $E\left[x\right]$ is the integer part of x and $I(l_{1},m_{1},l_{2},m_{2})=\max\left(\vert l_{1}-l_{2}\vert,m_{1}+m_{2}\right)$.

Vector product:

We operate likewise for the vector product of two general axisymmetric vectors $\vec X_{1}\left(r,\theta\right)$ and $\vec X_{2}\left(r,\theta\right)$, again expanded as in Eq. (A.20). We reach the following result:

 

$$\vec X_{1}(r,\theta)\wedge\vec X_{2}(r,\theta) \mathcal{X}_{\left(\vec{X}_{1}\wedge\vec{X}_{2}\right);0}(r)\wide... ...}_{2}\right);l}(r)P_{l}^{1}\left(\cos\theta\right)\widehat{e}_{\varphi}\right\}$$ =

$$\sum_{l=0}^{\infty}\left\{\left[\frac{\mathcal{X}_{\left(\vec{X}_... ...;l}(r)}{\mathcal{N}_{l}^{0}}\right]\vec{T}_{l}^{0}\left(\theta\right)\right\} ,$$ = (A.44)

where we have also used the following property of the Legendre function (Gradshteyn & Ryzhik 1967):

$$\frac{{\rm d}P_{l}(\mu)}{{\rm d}\mu}=\sum_{k=0}^{E\left(\frac{l-1}{2}\right)}\left(2l-4k-1\right)P_{l-2k-1}(\mu).$$ (A.45)

The radial function $\mathcal{X}_{\left(\vec{X}_{1}\wedge\vec{X}_{2}\right);l}\left(r\right)$ is given by:

$$\mathcal{X}_{\left(\vec{X}_{1}\wedge\vec{X}_{2}\right);l}(r)=... ...\mathcal{F}_{0}^{l_{2}}(r)\right]X_{l_{1},l_{2}}^{l}\right\} ,$$ (A.46)

with the same coefficient X_l1,l2^l derived above in (A.42). The two other radial functions $\mathcal{Y}_{\left(\vec{X}_{1}\wedge\vec{X}_{2}\right);l}\left(r\right)$ and $\mathcal{Z}_{\left(\vec{X}_{1}\wedge\vec{X}_{2}\right);l}\left(r\right)$ are given by

$$\mathcal{Y}_{\left(\vec{X}_{1}\wedge\vec{X}_{2}\right);l}(r)=... ...},1)}^{l_{1}+l_{2}}d_{l_{1},0,l_{2},1}^{j}\delta_{l,j}\right\}$$ (A.47)

and

$$\mathcal{Z}_{\left(\vec{X}_{1}\wedge\vec{X}_{2}\right);l}(r)=... ...,1)}^{l_{1}+l_{2}}d_{l_{1},0,l_{2},1}^{j}\delta_{l,j}\right\}.$$ (A.48)

Appendix B: Lorentz force

In this section we shall project the Lorentz force,

$$\vec{\mathcal{F}}\!\!_{\mathcal{L}}=\vec j\wedge\vec B=\left[\frac{1}{\mu_{0}}(\vec\nabla\wedge\vec B)\right]\wedge\vec B ,$$ (B.1)

on spherical vectorial harmonics, and write the result in terms of the poloidal and the toroidal magnetic stream-functions, respectively $\xi_{0}^{l}$ and $\chi_{0}^{l}$. We have just seen in Sect. A.2.3 how to expand the vector product of two axisymmetric vectors in the $\vec R_{l}^{0}(\theta)$, $\vec S_{l}^{0}(\theta)$, $\vec T_{l}^{0}(\theta)$. We apply the method to (B.1), with (cf. Eqs. (14) and (15))

$$\vec X_{1}\left(r,\theta\right) \vec j\left(r,\theta\right)=\frac{1}{\mu_{0}}\left(\vec\nabla\wed... ...l_{1}}^{0}(\theta)+\mathcal{C}_{0}^{l_{1}}(r)\vec T_{l_{1}}^{0}(\theta)\right\}$$ =

$$\frac{1}{\mu_{0}}\sum_{l_{1}=1}^{\infty}\left\{\left[l_{1}(l_{1}+... ...\left(\frac{\xi_{0}^{l_{1}}}{r}\right)\right]\vec T_{l_{1}}^{0}(\theta)\right\}$$ = (B.2)

and (cf. Eq. (11))

$$\vec X_{2}\left(r,\theta\right) \vec B\left(r,\theta\right)=\sum_{l_{2}=1}^{\infty}\left\{\mathca... ..._{2}}^{0}(\theta)+ \mathcal{F}_{0}^{l_{2}}(r)\vec T_{l_{2}}^{0}(\theta)\right\}$$ =

$$\sum_{l_{2}=1}^{\infty}\left\{\left[l_{2}(l_{2}+1)\frac{\xi_{0}^{... ...eta)+\left[\frac{\chi_{0}^{l_{2}}}{r}\right]\vec T_{l_{2}}^{0}(\theta)\right\}.$$ = (B.3)

So, we have, using the definition of Sect. A.2.3:

$$\left\{ \begin{array}{l@{\quad}l} \mathcal{A}_{0}^{l_{1}}=\fr... ...}}=\left[\frac{\chi_{0}^{l_{2}}}{r}\right]. \end{array}\right.$$ (B.4)

Therefore, we get using Eqs. (A.44)-(A.46)-(A.48):

$$\vec{\mathcal{F}}\!\!_{\mathcal{L}}(r,\theta) \vec X_{1}(r,\theta)\wedge\vec X_{2}(r,\theta)$$ =

$$\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};0}(r)\widehat{e}... ...\mathcal{L}};l}(r)P_{l}^{1}\left(\cos\theta\right)\widehat{e}_{\varphi}\right\}$$ =

$$\left[\frac{\mathcal{X}_{{\vec {\mathcal F}}_{\mathcal{L}};0}}{\m... ...};l}(r)}{\mathcal{N}_{l}^{0}}\right]\vec{T}_{l}^{0}\left(\theta\right)\right\},$$ = (B.5)

where the radial functions $\mathcal{X}_{\vec{\mathcal F}\!\!_{\mathcal L};l}\left(r\right)$, $\mathcal{Y}_{\vec{\mathcal F}\!\!_{\mathcal L};l}\left(r\right)$ and $\mathcal{Z}_{\vec{\mathcal F}\!\!_{\mathcal L};l}\left(r\right)$ are:

$$\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}=\frac{2}{... ...{\chi_{0}^{l_{2}}}{r}\right]\right)X_{l_{1},l_{2}}^{l}\right\}$$ (B.6)

$$\mathcal{Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}=\frac{1}{... ...},1)}^{l_{1}+l_{2}}d_{l_{1},0,l_{2},1}^{j}\delta_{l,j}\right\}$$ (B.7)

$$\mathcal{Z}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};l}=\frac{1}{... ...,1)}^{l_{1}+l_{2}}d_{l_{1},0,l_{2},1}^{j}\delta_{l,j}\right\}.$$ (B.8)

Explicit values for X_l1,l2^l (A.42) are given in Table B.1, and in Table B.2 for the coupling coefficients d_l1,1,l2,0^l.

Table B.1: Values of the X_l1,l2^l coefficient involved in the r-component of the Lorentz force and in the Ohmic heating.

Table B.2: Values of the d_l1,1,l2,0^l coefficients involved in the $\theta$- and in the $\varphi$-components of the Lorentz force and in the Ohmic heating.

Appendix C: Advection term of the induction equation

We deal likewise with the advection term in the induction equation:

$$\left(\vec{\mathcal{U}}_{\varphi}+\vec{\mathcal{U}}_{\rm M}\right)\wedge\vec B,$$ (C.1)

where $\vec{\mathcal{U}}_{\varphi}(r,\theta)=r\sin\theta~\Omega(r,\theta)\vec{\widehat e}_{\varphi}$ and $\vec{\mathcal{U}}_{\rm M}$ represents the meridional flow, whose expansions are found in (30). Therefore

$$\vec X_{1}\left(r,\theta\right) \vec{\mathcal{U}}_{\varphi}(r,\theta)+\vec{\mathcal{U}}_{\rm M}(r... ...l_{1}}^{0}(\theta)+\mathcal{C}_{0}^{l_{1}}(r)\vec T_{l_{1}}^{0}(\theta)\right\}$$ =

$$\sum_{l_{1}=0}^{\infty}\left\{u_{0}^{l_{1}}\left(r\right)\vec R_{... ..._{1}}^{0}(\theta)+w_{0}^{l_{1}}\left(r\right)\vec T_{l_{1}}^{0}(\theta)\right\}$$ = (C.2)

and (cf. Eq. (11))

$$\vec X_{2}\left(r,\theta\right) \vec B\left(r,\theta\right)=\sum_{l_{2}=1}^{\infty}\left\{\mathca... ..._{2}}^{0}(\theta)+ \mathcal{F}_{0}^{l_{2}}(r)\vec T_{l_{2}}^{0}(\theta)\right\}$$ =

$$\sum_{l_{2}=1}^{\infty}\left\{\left[l_{2}(l_{2}+1)\frac{\xi_{0}^{... ...eta)+\left[\frac{\chi_{0}^{l_{2}}}{r}\right]\vec T_{l_{2}}^{0}(\theta)\right\}.$$ = (C.3)

So, we have, expliciting the velocity field:

$$\left\{ \begin{array}{l@{\quad}l} \mathcal{A}_{0}^{l_{1}}=u_{... ...}}=\left[\frac{\chi_{0}^{l_{2}}}{r}\right]. \end{array}\right.$$ (C.4)

Hence we get, using Eqs. (A.44)-(A.46)-(A.48):

$$\left(\vec{\mathcal{U}}_{\varphi}(r,\theta)+\vec{\mathcal{U}}_{\rm M}(r,\theta)\right)\wedge\vec B(r,\theta) \vec X_{1}(r,\theta)\wedge\vec X_{2}(r,\theta)$$ =

$$\sum_{l=1}^{\infty}\left\{\mathcal{X}_{\rm {\rm Ad};l}(r)P_{l}\le... ...\rm {\rm Ad};l}(r)P_{l}^{1}\left(\cos\theta\right)\widehat{e}_{\varphi}\right\}$$ =

$$\sum_{l=1}^{\infty}\left\{\left[\frac{\mathcal{X}_{\rm {\rm Ad};l... ...};l}(r)}{\mathcal{N}_{l}^{0}}\right]\vec{T}_{l}^{0}\left(\theta\right)\right\},$$ = (C.5)

where $\mathcal{X}_{\rm {\rm Ad};l}\left(r\right)$, $\mathcal{Y}_{\rm {\rm Ad};l}\left(r\right)$ and $\mathcal{Z}_{\rm {\rm Ad};l}\left(r\right)$:

$$\mathcal{X}_{\rm {\rm Ad};l}=\frac{2}{3}\sum_{l_{1}=0}^{\inft... ...\chi_{0}^{l_{2}}}{r}\right]\right)X_{l_{1},l_{2}}^{l}\right\},$$ (C.6)

$$\mathcal{Y}_{{\rm {\rm Ad}};l}=\sum_{l_{1}=0}^{\infty}\sum_{l... ...,1)}^{l_{1}+l_{2}}d_{l_{1},0,l_{2},1}^{j}\delta_{l,j}\right\},$$ (C.7)

$$\mathcal{Z}_{{\rm {\rm Ad}};l}=\sum_{l_{1}=0}^{\infty}\sum_{l_{2}... ...j=I(l_{1},0,l_{2},1)}^{l_{1}+l_{2}}d_{l_{1},0,l_{2},1}^{j}\delta_{l,j}\right\};$$ (C.8)

Explicit values for X_l1,l2^l (A.42) are given in Table C.1, and in Tables C.2-C.4 for the coupling coefficients d_l1,1,l2,0^l.

Table C.1: Values of the X_l1,l2^l coefficients which are involved in the r-component of the induction equation.

Table C.2: Values of the coupling coefficients which are involved in the $\theta$- and in the $\varphi$- components of the induction equation for l=1.

Table C.3: Values of the coefficients involved in the $\theta$- and in the $\varphi$- components of the induction equation for l=2.

Table C.4: Same as in Table C.3, but for l=3.

Appendix D: Ohmic heating

For sake of completeness, we also calculate the Ohmic heating rate, which we shall express in terms of the magnetic stream functions $\xi_{0}^{l}$ and $\chi_{0}^{l}$. From Sect. 2.4, we get:

$$\mathcal{J}=\frac{1}{\mu_{0}}\left[\vert\vert\eta\vert\vert\o... ...\vec B\right)\right]\cdot\left(\vec\nabla\wedge\vec B\right) ,$$ (D.1)

where we allow for different eddy-diffusivities ( $\eta_{\rm h}, \eta_{\rm v}$) respectively in the horizontal and vertical direction. We apply again the method of Sect. A.2.3, and identify the two vector functions which enter the scalar product:

$$\vec X_{1}\left(r,\theta\right) \frac{1}{\mu_{0}}\left[\vert\vert\eta\vert\vert\otimes\left(\vec\... ...thcal{C}_{0}^{l_{1}}\left(r\right)\vec T_{l_{1}}^{0}\left(\theta\right)\right\}$$ =

$$\sum_{l_{1}=1}^{\infty}\left\{\frac{\eta_{\rm v}}{\mu_{0}}\left[l... ...{\xi_{0}^{l_{1}}}{r}\right)\right]\vec T_{l_{1}}^{0}\left(\theta\right)\right\}$$ = (D.2)

and

$$\vec X_{2}\left(r,\theta\right) \vec\nabla\wedge\vec B\left(r,\theta\right)=\sum_{l_{2}=1}^{\inft... ...thcal{F}_{0}^{l_{2}}\left(r\right)\vec T_{l_{2}}^{0}\left(\theta\right)\right\}$$ =

$$\sum_{l_{2}=1}^{\infty}\left\{\left[l_{2}(l_{2}+1)\frac{\chi_{0}^... ...\xi_{0}^{l_{2}}}{r}\right)\right]\vec T_{l_{2}}^{0}\left(\theta\right)\right\}.$$ = (D.3)

We thus obtain, using Eqs. (A.40) and (A.41):

$$\mathcal{J}\left(r,\theta\right)=\vec X_{1}\left(r,\theta\rig... ...nfty}\mathcal{J}_{l}\left(r\right)P_{l}\left(\cos\theta\right)$$ (D.4)

with the following expression for $\mathcal{J}_{l}\left(r\right)$:

 

$$\mathcal{J}_{l} \frac{1}{\mu_{0}}\sum_{l_{1}=1}^{\infty}\sum_{l_{2}=1}^{\infty}\m... ...l_{1},0,l_{2},0\right)}^{l_{1}+l_{2}}d_{l_{1},0,l_{2},0}^{j}\delta_{l,j}\right.$$ =

$$+{\left. \frac{2}{3}\frac{\eta_{\rm h}}{\mu_{0}}\left(\left[\frac... ...ft(\frac{\xi_{0}^{l_{2}}}{r}\right)\right]\right)X_{l_{1},l_{2}}^{l}\right\}} ,$$ (D.5)

recalling that X_l1,l2^l is given in (A.42), and in explicit form in Table B.1.

Appendix E: Equations to be implemented in stellar evolution codes

In this section we shall give the equations ready to be implemented in a stellar structure code, with all coupling coefficients being replaced by their explicit value, in the special case where only the dipole (l=1), quadrupole (l=2) and octupole (l=3) are kept in the magnetic field, and where the differental rotation is reduced to its first term (l=2).

E.1 Induction equation

The linear terms of this equation readily separate into their poloidal and toroidal components, and they project on a single multipole, whereas the advection term is the result of various couplings.

E.1.1 Equations for the dipole

The induction equation translates into two evolution equations for the magnetic stream functions:

$$\left\{ \begin{array}{l@{\quad}l} \frac{{\rm d}\xi_{0}^{1}}{{... ...lta_{1}=\partial_{r,r}+\frac{2}{r}\partial_{r}-\frac{2}{r^2} ,$$ (E.1)

and where the advective terms are given by:

$$\left\{ \begin{array}{l@{\quad}l} \mathcal{X}_{{\rm Ad};1}=\f... ...A}_{0}^{4}\mathcal{E}_{0}^{3}\right)\right] \end{array}\right.$$ (E.2)

with

$$\left\{ \begin{array}{l@{\quad}l} \mathcal{A}_{0}^{2}=2\sqrt{... ...0}^{3}=\left[\frac{\chi_{0}^{3}}{r}\right]. \end{array}\right.$$ (E.3)

E.1.2 Equations for the quadrupole

Likewise

$$\left\{ \begin{array}{l@{\quad}l} \displaystyle\frac{{\rm d}\... ...Delta_{2}=\partial_{r,r}+\frac{2}{r}\partial_{r}-\frac{6}{r^2}$$ (E.4)

and

$$\left\{ \begin{array}{l@{\quad}l} \mathcal{X}_{{\rm Ad};2}= ... ...}_{0}^{4}\mathcal{E}_{0}^{2}\right)\right]. \end{array}\right.$$ (E.5)

All functions $\mathcal{A}_{0}^{2}$ . . . $\mathcal{F}_{0}^{3}$ have the same meaning as in (E.3).

E.1.3 Equations for the octupole

The result is similar for the octupole:

$$\left\{ \begin{array}{l@{\quad}l} \displaystyle\frac{{\rm d}\... ...elta_{3}=\partial_{r,r}+\frac{2}{r}\partial_{r}-\frac{12}{r^2}$$ (E.6)

and

$$\left\{ \begin{array}{l@{\quad}l} \displaystyle\mathcal{X}_{{... ...}_{0}^{4}\mathcal{E}_{0}^{3}\right)\right]. \end{array}\right.$$ (E.7)

E.2 Mean equations

E.2.1 Mean rotation rate

We recall the transport equation for the mean angular momentum (48):

$$\rho {{\rm d} \over {\rm d}t} (r^2\overline{\Omega})=\frac{1}... ...ega}\right)+\overline{\Gamma}_{\vec{\mathcal{F}}_{\mathcal L}}$$ (E.8)

where the mean magnetic torque is given by:

$$\overline{\Gamma}_{\vec{\mathcal{F}}_{\mathcal L}}=\Gamma_{0}... ...ec{\mathcal{F}}\!\!_{\mathcal{L}};5}\right) \end{array}\right.$$ (E.9)

with

$$\left\{ \begin{array}{l@{\quad}l} \displaystyle\mathcal{Z}_{\... ...}_{0}^{3}\mathcal{E}_{0}^{2}\right)\right]. \end{array}\right.$$ (E.10)

The $\mathcal{A}_{0}^{l}$, $\mathcal{B}_{0}^{l}$, $\mathcal{C}_{0}^{l}$, $\mathcal{D}_{0}^{l}$, $\mathcal{E}_{0}^{l}$ and $\mathcal{F}_{0}^{l}$ with $l=\left\{1,2,3\right\}$ are given by:

$$\left\{ \begin{array}{l@{\quad}l} \displaystyle\mathcal{A}_{0... ...0}^{3}=\left[\frac{\chi_{0}^{3}}{r}\right]. \end{array}\right.$$ (E.11)

From (E.9), the final form of (E.8) is immediately derived:

$$\rho {{\rm d} \over {\rm d}t} (r^2\overline{\Omega})=\frac{1}... ...}\right)+r\mathcal{Z}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};1}.$$ (E.12)

E.2.2 Mean chemical composition

We restate the mean transport equation for the concentration of chemical species, which remains unchanged by the introduction of the magnetic field (Paper I, Eq. (B.2)):

$$\rho {{\rm d} \over {\rm d} t } \overline{c_{i}} +\frac{1}{r^... ...o (D_{\rm v}+D_{\rm eff})\partial_{r}\overline{c_{i}}\right] .$$ (E.13)

E.3 System for l = 2

E.3.1 Meridional circulation

We split (77)-(78) in two first order equations as:

$$U_{2}=\frac{L}{M\overline{g}}\left(\frac{P}{\overline{\rho}C_... ...e{T}}\right)\frac{1}{\nabla_{\hbox{ad}}-\nabla}\mathcal{B}_{2}$$ (E.14)

where

$$\mathcal{B}_{2} 2\left[1-\frac{\overline{f}_{\mathcal{\rm C}}+\overline{f}_{\math... ...{A}_{2}-\frac{2H_{\rm T}}{r}\left(1 + \frac{D_{\rm h}}{K}\right){\Psi_2}\right]$$ =

$$+\frac{\left(\overline{\epsilon}+\overline{\epsilon}_{\rm grav}\r... ...\rm d}t}+\Phi\frac{{\rm d}\ln\overline{\mu}}{{\rm d}t}\Lambda_{2}\right)\right]$$ (E.15)

and

$$\mathcal{A}_{2}=H_{\rm T}\partial_{r}\Psi_{2}-(1-\delta+\chi_{\rm T})\Psi_{2}-(\varphi+\chi_{\mu})\Lambda_{2}.$$ (E.16)

The coefficients related respectively to the centrifugal force and to the Lorentz force are given by:

$$\overline{f}_{\mathcal{\rm C}}=\frac{1}{r^2}\partial_{r}\left... ...r^2}\partial_{r}\left(r^2a_{2}\right)+6\frac{b_{2}}{r}\hbox{,}$$ (E.17)

$$\overline{f}_{\mathcal{L}}=\frac{1}{r^2}\partial_{r}\left(r^2... ...rac{\mathcal{Y}_{\vec{\mathcal F}_{\mathcal{L};2}}}{r\rho_{0}}$$ (E.18)

where

$$\left\{ \begin{array}{l@{\quad}l} \displaystyle a_{0}=\frac{2... ...2}+\frac{8}{35}r\overline{\Omega}\Omega_{2} \end{array}\right.$$ (E.19)

and

$$\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};0}=\frac{2}{... ..._{0}^{3}-\mathcal{B}_{0}^{3}\mathcal{F}_{0}^{3}\right)\right\}$$ (E.20)

 

$$\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};2} \frac{2}{3}\frac{1}{2\pi}\left\{13\sqrt{\frac{3}{7}}\left[\left(\... ...{0}^{2}\mathcal{E}_{0}^{2}-\mathcal{B}_{0}^{2}\mathcal{F}_{0}^{2}\right)\right.$$ =

$$+{\left. \frac{239}{7}\left(\mathcal{C}_{0}^{3}\mathcal{E}_{0}^{3}-\mathcal{B}_{0}^{3}\mathcal{F}_{0}^{3}\right)\right\}}$$ (E.21)

$$\mathcal{Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};2}=\frac{1}{... ...{0}^{3}-\mathcal{A}_{0}^{3}\mathcal{F}_{0}^{3}\right)\right\}.$$ (E.22)

The $\mathcal{A}_{0}^{l}$, $\mathcal{B}_{0}^{l}$, $\mathcal{C}_{0}^{l}$, $\mathcal{D}_{0}^{l}$, $\mathcal{E}_{0}^{l}$ and $\mathcal{F}_{0}^{l}$ with $l=\left\{1,2,3\right\}$ are given in (E.11).

The relative fluctuation of the effective gravity is given in (70); we apply it here to l=2:

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

where $\widehat{\phi}_{2}$ is solution of the Poisson equation (cf. Eq. (71)):

$$\frac{1}{r}\frac{{\rm d}^2}{{\rm d}r^2}\left(r\widehat{\phi}_... ...hcal{Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};2}\right)\right].$$ (E.24)

The Ohmic heating term is derived from (D.5):

$$\mathcal{J}_{2} \frac{1}{4\pi}\left\{\left[\mathcal{A}_{0}^{1}{D}_{0}^{1}-\left(\... ...\mathcal{E}_{0}^{3}+\mathcal{C}_{0}^{1}\mathcal{F}_{0}^{3}\right)\right]\right.$$ =

$$+{\left. \frac{2}{7}\left[5\mathcal{A}_{0}^{2}\mathcal{D}_{0}^{2}... ...athcal{E}_{0}^{3}+\mathcal{C}_{0}^{3}\mathcal{F}_{0}^{3}\right)\right]\right\}}$$ (E.25)

where the $\mathcal{A}_{0}^{l}$, $\mathcal{B}_{0}^{l}$, $\mathcal{C}_{0}^{l}$, $\mathcal{D}_{0}^{l}$, $\mathcal{E}_{0}^{l}$ and $\mathcal{F}_{0}^{l}$ with $l=\left\{1,2,3\right\}$ are given by:

$$\left\{ \begin{array}{l@{\quad}l} \displaystyle\mathcal{A}_{0... ...}\left(\frac{\xi_{0}^{3}}{r}\right)\right]. \end{array}\right.$$ (E.26)

E.3.2 Baroclinic relation

We apply Eq. (64) to l=2:

$$\varphi\Lambda_{2}-\delta\Psi_{2}=\frac{r}{\overline{g}}\left... ...thcal{F}}\!\!_{\mathcal{L}};2}}{\overline{\rho}}\right)\right]$$ (E.27)

where $\mathcal{D}_{2}=\frac{1}{3}\left[r\partial_{r}\overline{\Omega}^{2}\right]+\fra... ...overline{\Omega}\Omega_{2}\right)\right]+\frac{8}{7}\overline{\Omega}\Omega_{2}$ and where $\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};2}$ and $\mathcal{Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};2}$ have been given in (E.21)-(E.22).

E.3.3 Horizontal fluctuation of the molecular weight

The evolution equation for the mean molecular weight (Paper I, Eq. (39)) is unchanged by the introduction of the magnetic field:

$$\frac{{\rm d}\Lambda_{2}}{{\rm d}t}-\frac{{\rm d}\ln \overlin... ...U_2}{H_{\rm p}}\nabla_{\mu}-\frac{6}{r^2}D_{\rm h}\Lambda_{2}.$$ (E.28)

E.4 System for l = 4

E.4.1 Horizontal shear

We have:

$$\rho{{\rm d} \over {\rm d}t} \left(r^2 \Omega_{2}\right)-2\rh... ...)-\alpha U_{2}\right]=-10\rho\nu_{\rm h}\Omega_{2}+\Gamma_{2}$$ (E.29)

where the expression for $\Gamma_{2}$ has been given in Eq. (E.9).

E.4.2 Meridional circulation

In the same way as for l=2, we recast (77) and (78) in two first order equations as

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

where

$$\mathcal{T}_{4} 2\left[1-\frac{\overline{f}_{\mathcal{\rm C}}+\overline{f}_{\math... ..._{4}-\frac{20 H_{\rm T}}{3r}\left(1 + \frac{D_{\rm h}}{K}\right){\Psi_4}\right]$$ =

$$+\frac{\left(\overline{\epsilon}+\overline{\epsilon}_{\rm grav}\r... ...\rm d}t}+\Phi\frac{{\rm d}\ln\overline{\mu}}{{\rm d}t}\Lambda_{4}\right)\right]$$

and

$$\mathcal{A}_{4}=H_{\rm T}\partial_{r}\Psi_{4}-(1-\delta+\chi_{\rm T})\Psi_{4}-(\varphi+\chi_{\mu})\Lambda_{4}.$$ (E.31)

The coefficients related respectively to the centrifugal force and to the Lorentz force are given by:

$$\widetilde{f}_{\mathcal{\rm C},4}=\frac{1}{r^2}\partial_{r}\l... ...rac{\mathcal{Y}_{\vec{\mathcal F}_{\mathcal{L}};4}}{r\rho_{0}}$$ (E.32)

where:

$$\left\{ \begin{array}{l@{\quad}l} a_{4}=-\displaystyle\frac{2... ...4}=\frac{6}{35}r\overline{\Omega}\Omega_{2} \end{array}\right.$$ (E.33)

and

$$\left\{ \begin{array}{l@{\quad}l}\displaystyle \mathcal{X}_{\... ..._{0}^{3}\mathcal{F}_{0}^{3}\right)\right\}. \end{array}\right.$$ (E.34)

We get the relative fluctuation of the effective gravity from (70):

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

where $\widehat{\phi}_{4}$ obeys the Poisson equation:

$$% latex2html id marker 7587 \frac{1}{r}\frac{{\rm d}^2}{{\rm ... ...\hbox{ }\hbox{ }\hbox{ }\hbox{(cf. Eq.~(\ref{poisson-pert}))}.$$ (E.36)

We proceed as for l=2 to calculate the Ohmic heating:

$$\mathcal{J}_{4} \frac{1}{4\pi}\left\{4\sqrt{\frac{3}{7}}\left[\mathcal{A}_{0}^{3}... ...\mathcal{E}_{0}^{2}+\mathcal{C}_{0}^{2}\mathcal{F}_{0}^{2}\right)\right]\right.$$ =

$$+{\left. \frac{2}{11}\left[9\mathcal{A}_{0}^{3}\mathcal{D}_{0}^{3... ...athcal{E}_{0}^{3}+\mathcal{C}_{0}^{3}\mathcal{F}_{0}^{3}\right)\right]\right\}}$$ (E.37)

where the $\mathcal{A}_{0}^{l}$, $\mathcal{B}_{0}^{l}$, $\mathcal{C}_{0}^{l}$, $\mathcal{D}_{0}^{l}$, $\mathcal{E}_{0}^{l}$ and $\mathcal{F}_{0}^{l}$ with $l=\left\{1,2,3\right\}$ have been given in the previous section in (E.26).

E.4.3 Baroclinic relation

We apply (64) to l=4. We get:

$$\varphi\Lambda_{4}-\delta\Psi_{4}=\frac{r}{\overline{g}}\left... ...thcal{F}}\!\!_{\mathcal{L}};4}}{\overline{\rho}}\right)\right]$$ (E.38)

where $\mathcal{D}_{4}=\displaystyle\frac{6}{35}\left[r\partial_{r}\left(\overline{\Omega}\Omega_{2}\right)-2\overline{\Omega}\Omega_{2}\right]$ and where $\mathcal{X}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};4}$ and $\mathcal{Y}_{\vec{\mathcal{F}}\!\!_{\mathcal{L}};4}$ have been given in (E.34).

E.4.4 Horizontal fluctuation of the molecular weight

This equation is not affected by the introduction of the magnetic field:

$$\frac{{\rm d}\Lambda_{4}}{{\rm d}t}-\frac{{\rm d}\ln \overlin... ..._4}{H_{\rm p}}\nabla_{\mu}-\frac{20}{r^2}D_{\rm h}\Lambda_{4}.$$ (E.39)

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