A&A 471, 1011-1022 (2007)
DOI: 10.1051/0004-6361:20077418

Angular momentum conservation and torsional oscillations in the Sun and solar-like stars

A. F. Lanza

INAF - Osservatorio Astrofisico di Catania, via S. Sofia 78, 95123 Catania, Italy

Received 6 March 2007 / Accepted 7 June 2007

Abstract
Context. The solar torsional oscillations, i.e., the perturbations of the angular velocity of rotation associated with the eleven-year activity cycle, are a manifestation of the interaction among the interior magnetic fields, amplified and modulated by the solar dynamo, and rotation, meridional flow and turbulent thermal transport. Therefore, they can be used, at least in principle, to put constraints on this interaction. Similar phenomena are expected to be observed in solar-like stars and can be modelled to shed light on analogous interactions in different environments.
Aims. The source of torsional oscillations is investigated by means of a model for the angular momentum transport within the convection zone.
Methods. A description of the torsional oscillations is introduced, based on an analytical solution of the angular momentum equation in the mean-field approach. It provides information on the intensity and location of the torques producing the redistribution of the angular momentum within the convection zone of the Sun along the activity cycle. The method can be extended to solar-like stars for which some information on the time-dependence of the differential rotation is becoming available.
Results. Illustrative applications to the Sun and solar-like stars are presented. Under the hypothesis that the solar torsional oscillations are due to the mean-field Lorentz force, an amplitude of the Maxwell stresses $\vert B_{\rm r}B_{\phi}\vert \ga 8$ $\times$ 10³ G² at a depth of $\sim$0.85 $R_{\odot }$ at low latitude is estimated. Moreover, the phase relationship between $B_{\rm r}$ and $B_{\phi }$ can be estimated, suggesting that $B_{\rm r} B_{\phi} > 0$ below $\sim$0.85 $R_{\odot }$ and $B_{\rm r} B_{\phi} < 0$ above.
Conclusions. Such preliminary results show the capability of the proposed approach to constrain the amplitude, phase and location of the perturbations leading to the observed torsional oscillations.

Key words: Sun: rotation - Sun: activity - Sun: magnetic fields - Sun: interior - stars: rotation - stars: activity

1 Introduction

Doppler measurements of the surface rotation of the Sun show bands of faster and slower zonal flows that appear at midlatitudes and migrate toward the equator with the period of the eleven-year cycle, accompanying the bands of sunspot activity. The amplitude of such velocity perturbations, called torsional oscillations, is $\sim$5 m s^-1 and faster rotation is observed on the side equatorward of the sunspot belt (Howard & LaBonte 1980). Helioseismology has revealed that the torsional oscillations are not at all a superficial phenomenon but involve much of the convection zone, as shown, for example, by Howe et al. (2000), Vorontsov et al. (2002), Basu & Antia (2003) and more recently by Howe et al. (2006,2005). The amplitude of the angular velocity variation is $\delta \Omega/2\pi \sim 0.5{-}1$ nHz at least down to 10%-15% of the solar radius, although the precise depth of penetration of the oscillations is difficult to establish given the present uncertainties of the inversion methods in the lower half of the solar convection zone (e.g., Howe et al. 2006). In addition to such a low-latitude branch of the torsional oscillations, helioseismic studies have detected the presence of a high-latitude branch (above $\sim$ $60^{\circ}$ latitude) that propagates poleward and the amplitude of which is about $\delta \Omega /2\pi \sim 1{-}2$ nHz (e.g., Toomre et al. 2000; Basu & Antia 2001). Such a branch seems to propagate almost all the way down to the base of the convection zone.

A general description of the perturbation of the angular velocity of the torsional oscillations, given the present accuracy of the observations, is provided by the simple formula (e.g., Howe et al. 2005; Vorontsov et al. 2002):

 

$$% \omega (r, \theta) A^{(c)} (r, \theta) \cos (\sigma t) + A^{(s)} (r, \theta) \sin (\sigma t)$$ =

$$A(r, \theta) \sin[\sigma t + \phi(r, \theta)],$$ = (1)

where r is the distance from the centre of the Sun, $\theta$ the colatitude measured from the North pole, $\sigma$ the frequency of the eleven-year cycle, t the time, and the amplitude functions $A^{(c)} \equiv A \sin \phi$, $A^{(s)} \equiv A \cos \phi$depend on the amplitude A and the initial phase $\phi$. Moreover, the velocity perturbation is symmetric with respect to the equator:

$$% \omega( r, \theta) = \omega (r, \pi - \theta).$$ (2)

Several models have been proposed to interpret the torsional oscillations beginning with the pioneering work by Schüssler (1981) and Yoshimura (1981) who considered the Lorentz force associated with the magnetic fields in the activity belts as the cause of the velocity perturbations observed in the solar photosphere. Later models, based on the effects of the Lorentz force on the turbulent Reynolds stresses, were proposed, by e.g., Küker et al. (1996) and Kichatinov et al. (1999), following an original suggestion by Rüdiger & Kichatinov (1990). Spruit (2003) proposed that the low-latitude branch of the torsional oscillations is a geostrophic flow driven by temperature variations due to the enhanced radiative losses in the active region belts.

More recent studies by Covas et al. (2004,2005) present models based on the simultaneous solution of non-linear mean-field dynamo equations and the azimuthal component of the Navier-Stokes equation with a uniform turbulent viscosity. They reproduce the gross features of the torsional oscillations and of the solar activity cycle with an appropriate tuning of the free parameters. Rempel (2006,2007) considers the role of the meridional component of the Navier-Stokes equation in mean-field models and finds that the perturbation of the meridional flow cannot be neglected in the interpretation of the torsional oscillations. His models suggest that the low-latitude branch of the torsional oscillations cannot be explained solely by the effect of the mean-field Lorentz force, but that thermal perturbations in the active region belt and in the bulk of the convection zone do play an active role, as proposed by Spruit (2003).

In the present study, the angular momentum conservation is considered and the relevant equation in the mean-field approximation is solved analytically for the case of a turbulent viscosity that depends on the radial co-ordinate. A general solution is derived independently of any specific dynamo model, allowing us to put constraints on the localization of the torques producing the torsional oscillations. An illustrative application of the proposed methods is presented using the available data.

The observations of young solar-like stars by means of tomographic techniques based on high-resolution spectroscopy have provided recent evidence for time variation of their surface differential rotation (e.g., Donati et al. 2003; Jeffers et al. 2007). Lanza (2006a) has recently shown how such variations can be related to the intensity of the magnetic torque produced by a non-linear dynamo in their convective envelopes, in the case of rapidly rotating stars for which the Taylor-Proudman theorem applies. In the near future, the possibility of measuring the time variation of the rotational splittings of p-mode oscillations in solar-like stars may provide us with information on the changes of their internal rotation, although with limited spatial resolution. In the present study, we extend the considerations of Lanza (2006a) to the case of a generic internal rotation profile, not necessarily verifying the Taylor-Proudman theorem, to obtain hints on the amplitude of the torque leading to the rotation change.

2 The model

2.1 Hypotheses and basic equations

We consider an inertial reference frame with the origin in the barycentre of the Sun and the z-axis in the direction of the rotation axis. A spherical polar co-ordinate system  $(r, \theta, \varphi)$ is adopted, where r is the distance from the origin, $\theta$ the co-latitude measured from the North pole and $\varphi$ the azimuthal angle. We assume that all variables are independent of $\varphi$ and that the solar density stratification is spherically symmetric.

The equation for the angular momentum conservation in the mean-field approach reads (e.g., Rüdiger & Hollerbach 2004; Rüdiger 1989):

$$% \frac{\partial}{\partial t} (\rho r^{2} \sin^{2} \theta ~ \Omega) + \nabla \cdot {\vec \Theta} = 0,$$ (3)

where $\rho(r)$ is the density, $\Omega (r, \theta, t)$ the angular velocity and ${\vec \Theta}$ the angular momentum flux vector given by:

$$% {\vec \Theta} = (\rho r^{2} \sin^{2} \theta ~ \Omega) {\vec u}_... ...{\vec B} B_{\varphi} + \langle {\vec B}^{\prime} B_{\varphi}^{\prime} \rangle),$$ (4)

where ${\vec u}_{\rm (m)}$ is the meridional circulation, ${\vec u}^{\prime}$ the fluctuating velocity field, $\tilde{\mu}$ the magnetic permeability, ${\vec B}$ the mean magnetic field and ${\vec B}^{\prime}$ the fluctuating magnetic field; angular brackets indicate the Reynolds average defining the mean-field quantities. The Reynolds stresses can be written as:

$$% \langle \rho u^{\prime}_{i} u^{\prime}_{j} \rangle = -\eta_... ...\frac{\partial u_{j}}{\partial x_{i}} \right) + \Lambda_{ij},$$ (5)

where ${\vec u}$ is the mean flow field, $\eta_{\rm t}(r)$ is the turbulent viscosity, assumed to be a scalar function of r only, and $\Lambda_{ij}$ indicates the non-diffusive part of the Reynolds stresses due to the velocity correlations in a rotating star (see Rüdiger & Hollerbach 2004; Rüdiger 1989, for details). The conservation of the total angular momentum of the convection zone implies:

$$% \Theta_{r} = 0 \;{\rm for} \; r=r_{\rm b}, \; R_{\odot},$$ (6)

where $r_{\rm b}$ is the radius at the lower boundary of the convection zone and $R_{\odot }$ is the radius of the Sun.

The equation for the conservation of the angular momentum can be recast in the form:

 

$$% \frac{\partial \Omega}{\partial t}\! -\! \frac{1}{\rho r^{4}} \... ...} \left[ (1 - \mu^{2})^{2} \frac{\partial \Omega}{\partial \mu}\right] \!=\! S,$$ (7)

where $\mu \equiv \cos \theta$ and the source term S is given by:

$$% S = - \frac{\nabla \cdot {\vec \tau}}{\rho r^{2} (1- \mu^{2})},$$ (8)

and ${\vec \tau}$ is a vector whose components are:

$$% \tau_{i} = r \sin \theta \left[ \Lambda_{i \varphi} - \frac... ...ht) \right] + \rho r^{2} \sin^{2} \theta \Omega u_{\rm (m) i}.$$ (9)

The boundary conditions given by Eq. (6) can be written as:

$$% \frac{\partial \Omega}{\partial r} = 0 \;{\rm for} \; r=r_{\rm b}, \; R_{\odot},$$ (10)

when we assume $\tau_{\rm r} = 0$ at the surface. Note that helioseismic measurements indicate the presence of a subsurface shear layer with $\frac{\partial \Omega}{\partial r} < 0$ at low latitudes (Corbard & Thompson 2005), but we prefer to adopt the stress-free boundary condition (10) at the surface to ensure the conservation of the total angular momentum of the convection zone in our model.

The solar angular velocity can be split into a time-independent component  $\Omega_{0}$ and a time-dependent component $\omega$, i.e., the torsional oscillations:

$$% \Omega (r, \mu, t) = \Omega_{0} (r, \mu) + \omega (r, \mu, t).$$ (11)

The equation for the torsional oscillations thus becomes:

 

$$% \frac{\partial \omega}{\partial t} \!-\! \frac{1}{\rho r^{4}} \... ... (1 \!-\! \mu^{2})^{2} \frac{\partial \omega}{\partial \mu}\right] \!=\! S_{1},$$ (12)

where the perturbation of the source term is:

$$% S_{1} = - \frac{\nabla \cdot {\vec \tau}_{1}}{\rho r^{2} (1- \mu^{2})},$$ (13)

with

$$% \tau_{1 i} \!=\! r \sin \theta \left[ \Lambda_{i \varphi}^{... ...^{2} \sin^{2} \theta \bar{\Omega}_{0} u^{\rm (p)}_{\rm (m) i},$$ (14)

where $\Lambda_{i \varphi}^{\rm (p)}$ and $u^{\rm (p)}_{\rm (m) i}$ are the time-dependent perturbations of the non-diffusive Reynolds stresses and of the meridional circulation, respectively, and $\bar{\Omega}_{0}$ is the average of the solar angular velocity over the convection zone. Note that the Maxwell stresses appear in Eq. (14), but not in the corresponding equation for $\Omega_{0}$ because the solar magnetic field has no time-independent component. Moreover, in deriving Eq. (14), the variation of the angular velocity over the convection zone has been neglected in the term containing the perturbation of the meridional circulation since $\vert\omega\vert \ll \Omega_{0}$ (e.g., Rempel 2007; Rüdiger 1989). Equation (12) must be solved together with the boundary conditions:

$$% \frac{\partial \omega}{\partial r} = 0 \;{\rm for}\; r=r_{\rm b}, \; R_{\odot}.$$ (15)

   
2.2 Solution of the angular momentum equation

The general solution of Eq. (12) with the boundary conditions (15) can be obtained by the method of separation of the variables and expressed as a series of the form (e.g., Lanza 2006b, and references therein):

$$% \omega (r, \mu , t) = \sum_{n = 0, 2, 4, ...}^{\infty} \sum... ...}^{\infty} \alpha_{nk} (t) \zeta_{nk} (r) P_{n}^{(1,1)} (\mu),$$ (16)

where $\alpha_{nk}(t)$ and $\zeta_{nk}(r)$ are functions that will be specified below and $P_{n}^{(1, 1)}(\mu)$ are Jacobian polynomials, i.e., the finite solutions of the equation:

$$% \frac{{\rm d}}{{\rm d} \mu} \left[ (1- \mu^{2})^{2} \frac{{... ...\rm d} \mu }\right] + n(n+3) (1 - \mu^{2}) P_{n}^{(1, 1)} = 0,$$ (17)

in the interval $-1 \leq \mu \leq 1$ including its ends (e.g., Smirnov 1964a). The Jacobian polynomials form a complete and orthogonal set in the interval [-1, 1] with respect to the weight function  $(1-\mu^{2})$. Only the polynomials of even degree appear in Eq. (16) because the angular velocity perturbation is symmetric with respect to the equator (see Eq. (2)). For $n \gg 1$ the asymptotic expression of the Jacobian polynomials is (e.g., Gradshteyn & Ryzhik 1994):

$$% P_{n}^{(1,1)} (\cos \theta) = \frac{\cos \left[ (n + \frac... ...\frac{\theta}{2}\right)\right]^{3/2}} + O( n^{-\frac{3}{2}}).$$ (18)

The functions $\zeta_{nk}$ are the solutions of the Sturm-Liouville problem defined in the interval $r_{\rm b} \leq r \leq R_{\odot}$ by the equation:

$$% \frac{{\rm d}}{{\rm d}r} \left( r^{4} \eta_{\rm t} \frac{{\... ...ta_{\rm t} \zeta_{nk} + \lambda_{nk} \rho r^{4} \zeta_{nk} = 0$$ (19)

with the boundary conditions (following from Eq. (15)):

$$% \frac{{\rm d} \zeta_{nk}}{{\rm d}r} = 0 \;{\rm at} \; r=r_{\rm b}, \; R_{\odot}.$$ (20)

We shall consider normalized eigenfunctions, i.e.: $\int_{r_{\rm b}}^{R_{\odot}} \rho r^{4} \zeta_{nk}^{2} {\rm d}r$ = 1. The eigenfunctions  $\zeta_{nk}$ for a fixed n, form a complete and orthonormal set in the interval  $[r_{\rm b}, R_{\odot}]$ with respect to the weight function  $\rho r^{4}$ that does not depend on n. We recall from the theory of the Sturm-Liouville problem that the eigenvalues verify the inequality: $\lambda_{n 0} < \lambda_{n 1} <...\lambda_{n k}< \lambda_{n k+1}< ...$ and that the eigenfunction  $\zeta_{nk}$ has k nodes in the interval  $[r_{\rm b}, R_{\odot}]$ for each n. For n=0, the first eigenvalue corresponding to the eigenfunction  $\zeta_{00}$ is zero and the eigenfunction vanishes at all points in $[r_{\rm b}, R_{\odot}]$, as it is evident by integrating both sides of Eq. (19) in the same interval, applying the boundary conditions (20) and considering that $\zeta_{00}$ has no nodes. For n >0, all the eigenvalues  $\lambda_{n0}$ are positive, as can be derived by integrating both sides of Eq. (19) in the interval  $[r_{\rm b}, R_{\odot}]$, applying the boundary conditions (20) and considering that $\zeta_{n0}$ has no nodes. In view of the inequality given above, all the eigenvalues  $\lambda_{nk}$ are then positive for $n\geq 0$. Moreover, it is possible to prove that $\lambda_{n^{\prime} k} \geq \lambda_{n k}$ if $n^{\prime} > n$ and that (see Smirnov 1964b):

$$% \frac{\left(\frac{k^{2} \pi^{2}}{l^{2}}\right) p + n(n+3) q... ...ac{\left(\frac{k^{2} \pi^{2}}{l^{2}} \right) P + n(n+3) Q}{m},$$ (21)

where P, Q and M are the maximum values of the functions  $r^{4}\eta_{\rm t}$, $r^{2}\eta_{\rm t}$ and $\rho r^{4}$ in the interval  $[r_{\rm b}, R_{\odot}]$, respectively, and p, q and m their minimum values in the same interval, respectively; and $l \equiv \int_{r_{\rm b}}^{R_{\odot}} \sqrt{\rho / \eta_{\rm t}} {\rm d}r$. For $k \gg 1$ the asymptotic expression for the eigenfunction  $\zeta_{nk}$ is (e.g., Morse & Feshbach 1953):

$$% \zeta_{nk}(r) \simeq \left( \rho \eta_{\rm t}\right)^{-\fra... ...m b}}^{r} \sqrt{\rho/\eta_{\rm t}} {\rm d} r^{\prime} \right).$$ (22)

The time dependence of the solution in Eq. (16) is specified by the functions  $\alpha_{nk}$ that are given by:

$$% \frac{{\rm d} \alpha_{nk}}{{\rm d}t} + \lambda_{nk} \alpha_{nk}(t) = \beta_{nk} (t),$$ (23)

where the functions $\beta_{nk}$ appear in the development of the perturbation term S₁:

$$% S_{1}(r, \mu, t) = \sum_{n} \sum_{k} \beta_{nk} (t) \zeta_{nk} (r) P_{n}^{(1,1)}(\mu),$$ (24)

and are given by (Lanza 2006b):

 

$$% \beta_{nk} = \frac{(2n+3) (n+2)}{8(n+1)} \int_{r_{\rm b}}^{R_{\... ...mu^{2}) S_{1}(r, \mu, t) \zeta_{nk}(r) P_{n}^{(1,1)} (\mu) {\rm d}r {\rm d}\mu.$$ (25)

The solution of Eq. (23) is:

$$% \alpha_{nk}(t) = \alpha_{nk} (0) + \exp(-\lambda_{nk} t) \i... ...(t^{\prime}) \exp(\lambda_{nk} t^{\prime}) {\rm d} t^{\prime},$$ (26)

which allows us to specify the general solution of Eq. (12) with the boundary conditions (15) when the perturbation term  $S_{1}(r, \mu, t)$ and the initial conditions are given.

   
2.3 Solution for the solar torsional oscillations

To find the solution appropriate to the solar torsional oscillations as specified by Eqs. (1) and (2), it is useful to derive an alternative expression for the functions  $\beta_{nk}$ as follows. Substituting Eq. (13) into Eq. (25) and taking into account that the element of volume is ${\rm d}V = - 2 \pi r^{2} {\rm d}r {\rm d}\mu$, the right hand side of Eq. (25) can be recast in the form of a volume integral extended to the solar convection zone:

$$% \beta_{nk} = F_{n} \int_{V} (\nabla \cdot {\vec \tau}_{1}) \zeta_{nk} P_{n}^{(1,1)} {\rm d}V,$$ (27)

where:

$$% F_{n} \equiv \frac{(2n+3) (n+2)}{ 16\pi (n+1)}$$ (28)

is a factor coming from the normalization of the Jacobian polynomials. It is possible to further simplify Eq. (27) by considering the identity:

$$% \nabla \cdot [\zeta_{nk} P_{n}^{(1,1)} {\vec \tau}_{1}] = \... ...{1}) + {\vec \tau}_{1} \cdot \nabla[\zeta_{nk} P_{n}^{(1,1)}].$$ (29)

Integrating both sides of Eq. (29) over the volume of the convection zone and considering that the integral of the left hand side vanishes thanks to the Gauss's theorem and the condition that the radial component of the stresses  $\tau_{1 r}$ is zero on the boundaries, we find:

$$% \beta_{nk} = - F_{n} \int_{V} {\vec \tau}_{1} \cdot \nabla[\zeta_{nk} P_{n}^{(1,1)}] {\rm d}V.$$ (30)

The time dependence of the solar torsional oscillations specified in Eq. (1) suggests that a similar dependence for the perturbation term should be considered:

$$% {\vec \tau}_{1} (r, \mu, t) = {\vec \tau}_{1}^{(c)}(r, \mu) \cos (\sigma t) + {\vec \tau}_{1}^{(s)}(r, \mu) \sin (\sigma t),$$ (31)

from which a similar expression for the $\beta_{nk}$ follows by substitution into Eq. (27). If we put such an expression for $\beta_{nk}$ into Eq. (26) and perform the integrations with respect to the time, we find the stationary solution:

 

$$% \alpha_{nk} (t) \frac{F_{n}}{\lambda_{nk}^{2} \!+\! \sigma^{2}}$$ =

$$\times \left\{ \left[ \int_{V} \zeta_{nk} P_{n}^{(1,1)} \left( \l... ...la \cdot {\vec \tau}_{1}^{(s)} \right) {\rm d}V \right] \cos (\sigma t) \right.$$

$$\left. +\! \left[ \int_{V} \zeta_{nk} P_{n}^{(1,1)} \left( \lambd... ...cdot\! {\vec \tau}_{1}^{(c)} \right) {\rm d}V \right] \sin (\sigma t) \right\}.$$ (32)

This expression can be substituted into Eq. (16) to give the angular velocity perturbation. It can be written in the form of Eq. (1) with:

 

$$A^{(c)} (r, \mu) \!=\! \! \int_{V} \left[ G_{1}(r, \mu , r^{\prim... ...mu^{\prime}) \nabla \!\cdot\! {\vec \tau}_{1}^{(s)} \right] {\rm d} V^{\prime},$$

$$A^{(s)} (r, \mu) \!=\!\! \int_{V} \!\left[ G_{1}(\!r, \mu , r^{\p... ...mu^{\prime}) \nabla \!\cdot\! {\vec \tau}_{1}^{(c)} \right] {\rm d} V^{\prime},$$ (33)

where the symbol ${\rm d} V^{\prime}$ means that the volume integration is to be performed with respect to the variables  $r^{\prime}$ and $\mu^{\prime}$, and the functions G₁ and G₂ are Green functions defined as:

 

$$% G_{1}(r, \mu , r^{\prime}, \mu^{\prime}) = \sum_{n}^{\infty} F_... ...a_{nk}(r) P_{n}^{(1,1)}(\mu) \zeta_{nk}(r^{\prime}) P_{n}^{(1,1)}(\mu^{\prime})$$

$$G_{2}(r, \mu , r^{\prime}, \mu^{\prime}) = \sum_{n}^{\infty} F_{n... ..._{nk}(r) P_{n}^{(1,1)}(\mu) \zeta_{nk}(r^{\prime}) P_{n}^{(1,1)}(\mu^{\prime}).$$ (34)

The Green functions are continuous with respect to the arguments r, $\mu$, $r^{\prime}$, $\mu^{\prime}$, but their partial derivatives with respect to $r^{\prime}$, $\mu^{\prime}$ have discontinuities of the first kind in the points where $r=r^{\prime}$ or $\mu=\mu^{\prime}$. The convergence of the series in Eqs. (34), here used to represent the Green functions, is assured by the general theory of the Green function (e.g., Smirnov 1964b) and is also proven in Appendix A.

Equations (33) can be used to compute the angular velocity perturbation when ${\vec \tau}_{1}$ is known. Note that in the case in which the Lorentz force due to the mean field and the meridional flow are the only sources of angular momentum redistribution, Eq. (14) gives:

$$% \nabla \cdot {\vec \tau}_{1} = - \frac{1}{\tilde{\mu}}{\vec... ...hi}) + 2 \rho \bar{\Omega}_{0} (r \sin \theta) u_{\rm (m) s},$$ (35)

where ${\vec B}_{\rm p}$ is the mean poloidal magnetic field and $u_{\rm (m) s}$ the component of the meridional flow in the direction orthogonal to the rotation axis. To obtain Eq. (35) we made use of the solenoidal nature of the mean poloidal field and of the continuity equation for the meridional flow.

   
2.4 Localization of the source of the torsional oscillations in the solar convection zone

The results derived above allow us to introduce methods to localize the torques producing the torsional oscillations in the convection zone. Suppose that the observations provide us with the functions  $A^{(c)}(r, \mu)$ and $A^{(s)}(r, \mu)$ appearing in Eq. (1). The functions  $\alpha_{nk}(t)$ can be written as:

$$% \alpha_{nk} = a_{nk}^{(c)} \cos (\sigma t) + a_{nk}^{(s)} \sin (\sigma t),$$ (36)

where the constants a_nk^(c,s) are given by:

$$% a_{nk}^{(c,s)} = 2\pi F_{n} \int_{r_{\rm b}}^{R_{\odot}} \i... ...^{4} (1-\mu^{2}) \zeta_{nk} P_{n}^{(1,1)} {\rm d}r {\rm d}\mu.$$ (37)

Similarly, we can write:

$$% \beta_{nk}(t) = b_{nk}^{(c)} \cos (\sigma t) + b_{nk}^{(s)} \sin (\sigma t),$$ (38)

with the relationships:

 

$$b_{nk}^{(c)} = \lambda_{nk} a_{nk}^{(c)} + \sigma a_{nk}^{(s)},$$

$$b_{nk}^{(s)} = \lambda_{nk} a_{nk}^{(s)} - \sigma a_{nk}^{(c)},$$ (39)

that follow from Eq. (23).

The divergence of ${\vec \tau}_{1}^{(c,s)}$ can be obtained from Eqs. (13) and (24) as:

$$% \nabla \cdot {\vec \tau}_{1}^{(c,s)} = - \rho r^{2} (1-\mu^... ... \sum_{k} b_{nk}^{(c, s)} \zeta_{nk} (r) P_{n}^{(1,1)} (\mu).$$ (40)

Moreover, it is possible to construct a localized estimate of the perturbation term  ${\vec \tau}_{1}$ by considering a function $f(r, \mu)$ the gradient of which is different from zero only within a given volume V_f. It can be developed in series of the eigenfunctions in the form:

$$% f(r, \mu) = \sum_{n} \sum_{k} c_{nk} \zeta_{nk} (r) P_{n}^{(1,1)} (\mu),$$ (41)

where the coefficient c_nk are given by:

$$% c_{nk} = 2\pi F_{n} \int_{r_{\rm b}}^{R_{\odot}} \int_{-1}^... ...^{4} (1-\mu^{2}) \zeta_{nk} P_{n}^{(1,1)} {\rm d}r {\rm d}\mu.$$ (42)

Let us consider the equation:

 

$$\int_{V_{\rm f}} \left( {\vec \tau}_{1}^{(c,s)} \cdot \nabla f \r... ...} \sum_{k} \nabla [c_{nk} \zeta_{nk} (r) P_{n}^{(1,1)} (\mu)] \right) {\rm d}V,$$ (43)

obtained by means of Eq. (41). Considering Eqs. (30) and (38), Eq. (43) can be recast as:

 

$$% \int_{V_{\rm f}} \left( {\vec \tau}_{1}^{(c,s)} \cdot \nabla f \right) {\rm d}V = - \sum_{n} \frac{1}{F_{n}} \sum_{k} c_{nk} b_{nk}^{(c,s)}.$$ (44)

Moreover, if we introduce the volume average of the modulus of the perturbation  $\vert{\vec \tau}_{1}\vert$ with respect to the weight function  $\vert\nabla f\vert$, i.e.:

$$% \langle \vert {\vec \tau}_{1}^{(c,s)} \vert \rangle_{V_{\rm... ...rt {\rm d}V}{ \int_{V_{\rm f}} \vert\nabla f\vert {\rm d}V},$$ (45)

then Eq. (44) gives a lower limit for it in the form:

$$% \langle \vert {\vec \tau}_{1}^{(c,s)} \vert \rangle_{V_{\rm... ...right\vert}{\int_{V_{\rm f}} \vert\nabla f\vert {\rm d}V}\cdot$$ (46)

The minimum dimensions of the volume $V_{\rm f}$ are set by the spatial resolution of the measurements of the angular velocity variations. They depend on the accuracy of the rotational splitting coefficients, the inversion technique and the position within the convection zone (e.g., Schou et al. 1998; Howe et al. 2005). The minimum order of the Jacobian polynomials $N_{\rm m}$ needed to reproduce an angular velocity variation with a latitudinal resolution  $\Delta \theta$ is $N_{\rm m} \sim \frac{2\pi}{\Delta \theta}$. Similarly, the minimum order of the radial eigenfunctions $K_{\rm m}$ is set by the radial resolution $\Delta r$ as $K_{\rm m} \sim \frac{2(R-r_{\rm b})}{\Delta r}$. Therefore, it is possible to truncate the series in Eqs. (44) and (46) to those upper limits for n and k because the coefficients c_nk will decrease rapidly for $n> N_{\rm m}$ and $k > K_{\rm m}$ giving a negligible contribution to the sum.

The statistical errors in the measurements of the angular velocity variations can be easily propagated through the linear Eqs. (37), (39), (40) and (44) to find the errors on the estimates of $\nabla \cdot {\vec \tau}_{1}^{(c,s)}$ or the average of ${\vec \tau}_{1}$. For instance, if we consider the standard deviations  $\sigma_{i}$ of the data d_i, i.e., the rotational splittings or the splitting coefficients from which the internal rotation is derived, the standard deviation  $\sigma_{\rm I}$ of the integral in Eq. (44) is:

$$% \sigma_{\rm I}^{2} \simeq \sum_{n} \frac{1}{F_{n}^{2}} \sum... ...{nk}^{2} \lambda_{nk}^{2} \sum_{i} e_{ink}^{2} \sigma_{i}^{2},$$ (47)

where

$$% e_{ink} \equiv 2 \pi F_{n} \int_{r_{\rm b}}^{R_{\odot}} \in... ... c_{i} (r, \mu) \zeta_{nk} P_{n}^{(1,1)} {\rm d}r {\rm d}\mu,$$ (48)

and the functions $c_{i} (r, \mu)$ are the rotational inversion coefficients defined in Eq. (8) of Schou et al. (1998).

Note that a constant relative error $\epsilon = \Delta \omega/ \omega$ in the measurements of A^(c,s) leads to the same relative error in Eq. (40) and in Eqs. (44) and (46), given the linear equations that relate the corresponding quantities. Indeed, there is also a systematic error in our inversion method related to the poor knowledge of the turbulent viscosity  $\eta_{\rm t}(r)$ that determines the form of the radial eigenfunctions  $\zeta_{nk}$.

2.5 Kinetic energy variation and dissipation

The variation of the kinetic energy of rotation associated with the torsional oscillations, averaged over the eleven-year cycle, can be computed according to Lanza (2006b) and is:

$$% \langle \Delta {\cal T} \rangle_{\rm c}= \sum_{n} \sum_{k} \langle \Delta {\cal T}_{nk} \rangle_{\rm c},$$ (49)

where

$$% \langle \Delta {\cal T}_{nk} \rangle_{\rm c} = \frac{1}{4F_{n}} \left\{ [a_{nk}^{(c)}]^{2} + [a_{nk}^{(s)}]^{2} \right\}.$$ (50)

The average dissipation rate of the kinetic energy of the torsional oscillations due to the turbulent viscosity is:

$$% \left\langle \frac{{\rm d} {\cal T}}{{\rm d}t} \right\rangl... ...{k} \lambda_{nk} \langle \Delta {\cal T}_{nk} \rangle_{\rm c}.$$ (51)

   
2.6 Torsional oscillations due to mean-field Lorentz force

Most models of the torsional oscillations assume that they are due to the Lorentz force produced by the mean field as derived from dynamo models. Therefore, let us consider the case in which only the mean-field Maxwell stresses contribute to the perturbation, i.e.:

$$% {\vec \tau}_{1} = - \frac{1}{\tilde{\mu}} (r \sin \theta) B_{\phi} ~ {\vec B}_{\rm p}.$$ (52)

If the mean radial $B_{\rm r}$ and toroidal fields $B_{\phi }$ are given by:

$$B_{\rm r} = B_{\rm0 r} \cos \left(\frac{1}{2} \sigma t\right),$$

$$B_{\phi} = B_{0 \phi} \cos \left(\frac{1}{2} \sigma t + \Pi_{\rm r}\right),$$ (53)

where $\Pi _{\rm r}$ is the phase lag between the two field components, the components of $\tau_{1 r}$ in Eq. (31) are:

 

$$\tau_{1r}^{(c)} = -\frac{1}{2} \cos \Pi_{\rm r} ~ \frac{r \sin \theta}{\tilde{\mu}} ~ B_{\rm0 r} B_{0 \phi},$$

$$\tau_{1r}^{(s)} = \frac{1}{2} \sin \Pi_{\rm r} ~ \frac{r \sin \theta}{\tilde{\mu}} ~ B_{\rm0 r} B_{0 \phi}.$$ (54)

An estimate of $\tau_{1r}^{(c)}$ and $\tau_{1r}^{(s)}$ can be obtained from the method outlined in Sect. 2.4, considering a localization function f(r) that depends only on the radial co-ordinate. Specifically, Eq. (44) can be used to compute a volume average of $\tau_{1r}^{(c)}$ and $\tau_{1r}^{(s)}$ from which the average stress amplitude  $\vert B_{\rm r} B_{\phi} \vert$ and phase lag  $\Pi _{\rm r}$ can be determined. Analogous considerations are valid for the meridional component of the mean field  $B_{\theta }$ and $B_{\phi }$. Adopting a localization function $f(\theta)$ depending only on $\theta$, it is possible to estimate $\vert B_{\theta} B_{\phi}\vert$ and the phase lag  $\Pi _{\theta }$ between $B_{\theta }$ and $B_{\phi }$. Such results are important to constrain mean-field dynamo models of the solar cycle, as discussed by, e.g., Schlichenmaier & Stix (1995) (see also Sect. 3.3).

   
2.7 Application to solar-like stars

Sequences of Doppler images can be used to measure the surface differential rotation of solar-like stars and its time variability, as done by, e.g., Donati et al. (2003) and Jeffers et al. (2007) in the cases of AB Dor and LQ Hya. Lanza (2006a) discussed the implications of the observed changes of the surface differential rotation on the internal dynamics of their convection zones, assuming that the angular velocity is constant over cylindrical surfaces co-axial with the rotation axis. The present model allows us to relax the Taylor-Proudman constraint on the internal angular velocity, but some different assumptions must be introduced to obtain the internal torques in the active stars. Here we assume that the variation of the surface differential rotation is entirely due to the Maxwell stresses of the internal magnetic fields, localized in the overshoot layer below the convection zone, as in the interface dynamo model by, e.g., Parker (1993). The interior model of the Sun can also be applied to AB Dor and LQ Hya because they have a similar relative depth of the convection zone. Therefore, the basic quantities can be scaled according to the stellar parameters, as explained in Lanza (2006a). Following Donati et al. (2003), we consider a surface differential rotation of the form:

$$% \Omega(\mu, t) = \Omega_{\rm eq} (t) - {\rm d}\Omega(t) \mu^{2},$$ (55)

where $\Omega_{\rm eq}$ is the equatorial angular velocity and ${\rm d}\Omega$ the latitudinal shear, both regarded as time-dependent. Since P₀^(1,1) = 1 and $\mu^{2} = \frac{1}{5} +\frac{4}{15} P_{2}^{(1,1)}$, it can be recast in the form of Eq. (16) leading to:

 

$$\Omega_{\rm eq} - \frac{1}{5} {\rm d}\Omega = \sum_{k} \alpha_{0k}(t) \zeta_{0k}(R),$$

$$- \frac{4}{5} {\rm d}\Omega = \sum_{k} \alpha_{2k}(t) \zeta_{2k}(R),$$ (56)

where R is the radius of the star. The functions  $\beta_{nk}(t)$ can be computed by assuming that the timescale of variation of the differential rotation  $t_{\rm DR}$ is significantly shorter than the timescales for angular momentum transport, as given by $\lambda_{nk}^{-1}$. This assumption is justified in the case of rapidly rotating stars by the observed variation timescales of the order of a few years together with the expected rotational quenching of the viscosity (e.g., Kichatinov et al. 1994). Therefore, Eq. (23) leads to $\alpha_{nk} (t) \simeq t_{\rm DR} \beta_{nk}(t)$. The angular velocity can be written in a form similar to Eq. (33) by introducing an appropriate Green function. For instance, considering the first of Eqs. (56), we find:

$$% \Omega_{\rm eq} - \frac{1}{5} {\rm d}\Omega = \frac{3}{8\pi... ... \tau}_{1} \right) G_{\rm s}(R,r^{\prime}) {\rm d} V^{\prime},$$ (57)

where:

$$% G_{\rm s} (r, r^{\prime}) = \sum_{k} \zeta_{0k} (r) \zeta_{0k} (r^{\prime}),$$ (58)

and the integration is extended over the stellar convection zone. Assuming that $\nabla \cdot {\vec \tau}_{1}$ is different from zero only in an overshoot layer of volume $V_{\rm o}$ and applying considerations similar to those of Sects. 2.3 and 2.4, we find a lower limit for the variation of the averaged perturbation term:

$$% \langle \Delta (\nabla \cdot {\vec \tau}_{1}) \rangle_{V_{\... ...} - \frac{1}{5} \Delta ({\rm d}\Omega) }{V_{\rm o} M} \right],$$ (59)

where $\Delta \Omega_{\rm eq}$ and $\Delta ({\rm d}\Omega)$ are the amplitudes of variation of the differential rotation parameters, and M is the maximum of $G_{\rm s} (R, r^{\prime})$ in the overshoot layer.

This result made use of the limited information we can get from surface differential rotation. However, in the near future, the observations of the rotational splittings of stellar oscillations promise to give information on the internal rotation and its possible time variations. Since only the modes of low degrees ( $\ell \leq 3$) are detectable in disk-integrated measurements, the spatial resolution of the derived internal angular velocity profile is very low. Lochard et al. (2005) considered the case in which only the mean radial profile of $\Omega$ is measurable.

In view of such an information accessible through asteroseismology, let us consider a more general case in which some average of the internal angular velocity of the star, say $\omega_{\rm m}(t)$, can be measured as a function of the time:

$$% \omega_{\rm m}(t) = \int_{r_{\rm b}}^{R} \int_{-1}^{1} \omega(r, \mu ,t) w(r, \mu) {\rm d}r {\rm d}\mu,$$ (60)

where w is an appropriate weight function that takes into account the averaging effects of the limited spatial resolution, and $r_{\rm b}$ is the radius at the base of the stellar convective envelope. We introduce an auxiliary weight function:

$$% w_{1}(r, \mu) \equiv \frac{w(r, \mu)}{(1-\mu^{2}) \rho r^{4}},$$ (61)

which will not diverge toward the poles ( $\mu = \pm 1$) and the surface ($\rho = 0$) because the weight function w is localized in the stellar interior and goes to zero rapidly enough toward the rotation axis and the stellar surface. We can develop the function w₁ as:

$$% w_{1} (r, \mu) = \sum_{n}^{N_{\rm w}} \sum_{k}^{K_{\rm w}} w_{nk} \zeta_{nk}(r) P_{n}^{(1,1)}(\mu),$$ (62)

where the summation can be truncated at some low orders, say, $N_{\rm w}$ and $K_{\rm w}$, because the function w₁ is not sharply localized in the stellar interior. Considering Eq. (60), we find:

$$% \omega_{\rm m} (t) = \sum_{n}^{N_{\rm w}} \sum_{k}^{K_{\rm w}} \frac{1}{2 \pi F_{n}} w_{nk} \alpha_{nk}.$$ (63)

For the sake of simplicity, we assume now that the time scale of variation of the functions  $\alpha_{nk}(t)$ is long with respect to $\lambda_{nk}^{-1}$ so that Eq. (23) gives: $\beta_{nk} \simeq \lambda_{nk} \alpha_{nk}$. Considering Eq. (30), we find:

$$% \omega_{\rm m}(t) = - \int_{V} {\vec \tau}_{1} \cdot \left\... ..._{nk}} \nabla \left[ \zeta_{nk} P_{n}^{(1,1)} \right]\right\}.$$ (64)

Equation (64) can be used to find a lower limit to the average of $\vert{\vec \tau}_{1}\vert$ over the convection zone volume. If we indicate by $M_{\rm w}$ the maximum of the function:

$$% \left\vert \sum_{n}^{N_{\rm w}} \sum_{k}^{K_{\rm w}} \frac{... ...k}} \nabla \left[ \zeta_{nk} P_{n}^{(1,1)} \right] \right\vert$$ (65)

over the volume of the convection zone, we find:

$$% \langle \vert {\vec \tau}_{1}(t)\vert \rangle \equiv \frac{... ...{V} \geq \frac{\vert\omega_{\rm m}(t)\vert}{M_{\rm w} V}\cdot$$ (66)

3 Application to the Sun

   
3.1 Interior model, eigenvalues and eigenfunctions

A model of the solar interior can be used to specify the functions $\rho(r)$ and $\eta_{\rm t}(r)$ that appear in our equations. While the density stratification can be determined with an accuracy better than 0.5%, the turbulent dynamical viscosity is uncertain by at least one order of magnitude and it is estimated from the mixing-length theory according to the formula:

$$% \eta_{t} = \frac{1}{3} \alpha_{\rm ML} \rho u_{\rm c} H_{\rm p},$$ (67)

where $\alpha_{\rm ML}$ is the ratio of the mixing-length to the pressure scale height $H_{\rm p}$ and $u_{\rm c}$ is the convective velocity given by:

$$% u_{\rm c} = \left( \frac{\alpha_{\rm ML} L_{\odot} }{40 \pi r^{2} \rho}\right)^{\frac{1}{3}},$$ (68)

where $L_{\odot}$ is the luminosity of the Sun. In our computations we adopt $\alpha_{\rm ML} = 1.5$ and assume the solar model S for the interior quantities (Christensen-Dalsgaard et al. 1996). In that model, the base of the convection zone is at $r = 0.713~R_{\odot}$. We consider also the effect of an overshoot layer extending between $r= 0.673~R_{\odot}$ and the base of the convection zone within which the turbulent dynamical viscosity is assumed to increase linearly from zero up to the value at the base of the convection zone. The density and the turbulent viscosity are plotted in Fig. 1 where their values have been normalized to the values at the base of the convection zone (i.e., $\rho_{0} = 0.1875$ g cm^-3 and $\eta _{0}=2.56$ $\times$ 10¹² g cm^-1 s^-1, respectively).

	Figure 1: The ratio of the density to the density at the base of the convection zone ( $\rho / \rho _{0}$ - solid line) and the ratio of the turbulent dynamical viscosity to the turbulent viscosity at the base of the convection zone ( $\eta _{\rm t} /\eta _{0}$ - dotted) versus the fractionary radius  $r/R_{\odot }$ in our solar interior model.
Open with DEXTER

The basic equations of our model (i.e., 12, 13 and 14) can be made nondimensional by adopting as the unit of length the solar radius $R_{\odot }$, as the unit of density $\rho_{0}$, i.e., the density at the base of the solar convection zone, and as a unit of time $t_{0} = \rho_{0} R_{\odot}^{2}/\eta_{0}$, where $\eta_{0}$ is the turbulent viscosity at the base of the convection zone. Indeed, the value of the diffusion coefficients estimated from the mixing-length theory leads to a too short period for the solar cycle in mean-field dynamo models. Covas et al. (2004) adopted a turbulent magnetic diffusivity $\nu_{\rm t} = 3$ $\times$ 10¹¹ cm² s^-1 to get a sunspot cycle of $\sim$11 yr. This implies $\eta_{0} \sim \rho \nu_{\rm t} = 5.62$ $\times$ 10¹⁰ g cm^-1 s^-1 in our model.

The radial eigenfunctions $\zeta_{nk}$ and the Jacobian polynomials P_n^(1,1) have been computed from the respective Sturm-Liouville problem equations by means of the Fortran 77 subroutine sleign2.f (Bailey et al. 2001). For the radial eigenfunctions  $\zeta_{nk}$, the Sturm-Lioville problem has been solved with Neumann boundary conditions at both ends, set at $0.675~R_{\odot}$ and $0.99~R_{\odot}$ to avoid divergence at the surface. For the Jacobian polynomials, limit point boundary conditions have been adopted at $\mu = \pm 1$.

Note that, to be rigorous, it would be better to use the helioseismic estimate of $\frac{\partial \Omega}{\partial r}$ at r=0.99 where we fixed our outer boundary for the computation of the radial eigenfunctions instead of the stress-free boundary condition that is valid only at the surface. However, the differences are confined to the outermost layer of the solar convection zone, the moment of inertia of which is so small that there are no pratical consequences.

The eigenvalues $\lambda_{nk}$ and the eigenfunctions  $\zeta_{nk}$ computed by sleign2.f have been compared with those computed by means of the code introduced by Lanza (2006b). The relative differences in the eigenvalues and in the eigenfunctions are lower than 1.5% for $n \leq 14$, $k \leq 10$. However, some problems of convergence of the numerical algorithm used by sleign2.f have been found for $k\geq 20$, particularly for $n \geq 30$, so we decided to limit its application up to k = 19.

The Jacobian polynomials and the eigenvalues computed by sleign2.f are very good up to n =30, as it has been found by comparison with their analytic expressions up to n=10 and their asymptotic expressions for $n\geq 12$. We conclude that for $n \geq 30$ and $k\geq 20$ it is better to use the asymptotic formulae (22) and (18) instead of the numerically computed $\zeta_{nk}$ and P_n^(1,1).

The eigenvalue $\lambda_{nk}$ gives the inverse of the characteristic timescale of angular momentum transfer of the mode corresponding to $\zeta_{nk}$ under the action of the turbulent viscosity. The longest timescale corresponds to the lowest eigenvalue, i.e., $\lambda_{20}^{-1} = 0.078$ in nondimensional units. It corresponds to a timescale of 0.86 yr with the turbulent viscosity given by the mixing-length theory and to 39.5 yr with $\eta _{0}=5.62$ $\times$ 10¹⁰ g cm^-1 s^-1.

   
3.2 Localization functions for the source of the torsional oscillations

The available data on the torsional oscillations are displayed with a typical radial resolution of 0.05 $R_{\odot }$ and a latitudinal resolution of $15^{\circ}$ in Howe et al. (2006,2005). The rotational inversion kernels of Schou et al. (1998) show a higher radial resolution close to the surface, but, given the small amplitude of the torsional oscillations, the choice of a uniform resolution of 0.05 $R_{\odot }$ seems to be better.

We have found that the best results of the localization of the source term with the method outlined in Sect. 2.4 are obtained with localization functions that depend on r or $\mu$ only and that have a smooth derivative. As a typical function to probe the radial localization, we adopt:

$$% f(r)= \left\{ \begin{array}{ll} 0 & \mbox{for $r \leq r_{1... ...$ ,} \\ 2 & \mbox{for $r \geq r_{2}$ .} \end{array} \right.$$ (69)

This function has zero derivative, except in the interval  ]r₁, r₂[ where its derivative varies smoothly reaching a maximum in the mid of the interval. In Fig. 2 we plot the modulus of the derivative $\left\vert \frac{{\rm d}f}{{\rm d}r} \right\vert$, as given by the sum of the series of the radial eigenfunctions truncated at k =19, for three different intervals centered at 0.725, 0.85 and 0.965 $R_{\odot }$, all with an amplitude of r₂ - r₁ = 0.05 $R_{\odot }$. The derivative is well approximated by the truncated series, with small sidelobes the amplitude of which increases toward the surface of the Sun because the eigenfunctions scale as $(\rho \eta_{\rm t})^{-\frac{1}{4}}$, according to the asymptotic expression (22).

Figure 2: The modulus of the radial localization kernel versus the relative radius for the function (69), as obtained by truncating the series of the radial eigenfunctions at $K_{\rm m }= 19$, for three different intervals centered at r=0.725 ( upper panel), 0.85 ( middle panel) and 0.965 $R_{\odot }$ ( lower panel), all with an amplitude of r2 - r1 = 0.05 $R_{\odot }$. The modulus of the derivative has been normalized to its maximum value.

Open with DEXTER

The performance of the localization method introduced in Sect. 2.4 has been tested with simulated data in the absence of noise. The results of two such tests are plotted in Fig. 3 where we consider the case of a purely radial perturbation ${\vec \tau}_{1} = \tau_{1}\hat{\vec r}$ localized within an interval of amplitude 0.05 $R_{\odot }$. Its time dependence is assumed to be purely cosinusoidal and the corresponding coefficients  $\beta_{nk}$ are computed by means of Eq. (30) up to the orders n=38 and k=19. From the $\beta_{nk}$, the coefficients  $\alpha_{nk}$ are computed by means of Eq. (26) and the simulated angular velocity perturbation follows from Eq. (16).

When the interval in which ${\vec \tau}_{1}$ is localized coincides with one of the inversion intervals  [r₁, r₂], the lower limit for $\vert{\vec \tau}_{1}\vert$ turns out to be $\sim$80%-90% of the value assumed in the simulation (Fig. 3, upper panel). The agreement increases up to 90%-95% if we consider a case in which $\tau_{1 r}$ has a constant sign and put the value of the derivative $\frac{{\rm d}f}{{\rm d}r}$ in the denominator of Eq. (46) instead of its modulus. This happens because the modulus of the derivative increases the effects of the sidelobes by increasing the value of the denominator in Eq. (46). When the interval in which the assumed ${\vec \tau}_{1}$ is localized does not coincide with an interval of the inversion grid, the inversion method still performs well distributing the contributions among neighbour intervals nearly in the correct proportion (Fig. 3, lower panel). The case of simulations including a noise component is straightforward to treat thanks to the linear character of our inversion method. The inverted value turns out to be the sum of the inverted noiseless value and of the contribution coming from the inversion of the noise. Their amplitude ratio is equal to the ratio of the amplitudes of the input noise to the input signal (for constant relative signal errors), as discussed in Sect. 2.4.

Figure 3: Upper panel: test case of the application of the inversion method from Sect. 2.4. An input profile with ${\vec \tau}_{1} = 1.0 \hat{\vec r}$ in nondimensional units between 0.80 and 0.85 $R_{\odot }$ (solid line) is used to simulate a noiseless profile of angular velocity perturbation. The reconstructed lower-limit profile according to Eq. (46) is plotted as a dotted line. Lower panel: same as in the upper panel, but with an input profile localized between 0.775 and 0.825 $R_{\odot }$.

Open with DEXTER

The localization in latitude can be sampled by means of a localization function of the kind:

$$% f(\mu)= \left\{ \begin{array}{ll} 0 & \mbox{for $ -1 < \mu... ...0 & \mbox{for $\mu_{2} \leq \mu < 1 $ .} \end{array} \right.$$ (70)

Such a function is symmetric with respect to the equator so that its derivative $\frac{{\rm d}f}{{\rm d}\mu}$ is antisymmetric, as it is the source term  ${\vec \tau}_{1}$ leading to a symmetric perturbation of the angular velocity. The localization function has a derivative always equal to zero except in two intervals  $]\!-\!\mu_{2}, -\mu_{1}[$ and $]\mu_{1}, \mu_{2}[$, symmetric with respect to the solar equator. Its representation by means of the Jacobian polynomials with degree up to $N_{\rm } = 30$ is given in Fig. 4. Note that the maximum of $\left\vert \frac{{\rm d}f}{{\rm d}\mu} \right\vert$is reached in two sharp peaks at $\mu = \pm 1$. However, they are so narrow that their contribution to the integral in Eq. (46) is modest in comparison to those of the broader peaks in the intervals  $[-\mu_{2}, -\mu_{1}]$ and $[\mu_{1}, \mu_{2}]$. Several tests have been performed, as in the case of the radial localization, to assess the performance of the proposed method. The results are similar to those obtained for the radial case and are not discussed here.

Figure 4: The modulus of the latitudinal localization kernel versus $\mu$ for the function (70), obtained by truncating the series of the Jacobian polynomials at $N_{\rm m }= 30$, for three different intervals of $\mu$ in the Northern hemisphere, i.e., [0, 0.26] ( upper panel), [0.5, 0.7] ( middle panel) and [0.87, 0.99] ( lower panel). Note the symmetry of the modulus of the kernel with respect to the equator. The value of $\left\vert \frac{{\rm d}f}{{\rm d}\mu} \right\vert$is normalized at its maximum at $\mu = \pm 1$.

Open with DEXTER

We conclude that our choice of $N_{\rm m }= 30$ and $K_{\rm m }= 19$ is perfectly adequate to invert the available data on the solar torsional oscillations to derive information on the location of the perturbation term within the convection zone.

   
3.3 Results on the sources of the torsional oscillations

The data plotted in Fig. 4 of Howe et al. (2006) can be used for an illustrative application of the inversion methods introduced in Sect. 2.4. They are given with a sampling of $15^{\circ}$ between the equator and $60^{\circ}$ of latitude, i.e., in the range in which the rotational inversion techniques perform better (Schou et al. 1998). The features distinguishable on the plots indicate an actual radial resolution of $\sim$0.05 $R_{\odot }$, in agreement with the sampling adopted in Fig. 3 of Howe et al. (2005). An average statistical error of about 30% can be assumed for the amplitude, whereas the phase errors become very large below 0.80 $R_{\odot }$, especially at low latitudes, because of the uncertainty in the reconstruction of the signal in the deep layers. Note that the error intervals reported in Fig. 4 of Howe et al. (2006) indicate only how the particular method solution (here an OLA inversion of SoHO/MDI data) would vary with a different realization of the input data affected by a randon Gaussian noise. Unfortunately, they do not give the statistical ranges in which the true values of the amplitude and phase are likely to lie. This is not a major limitation in the context of the present study because we aim at illustrating the capabilities of the proposed approach rather than derive definitive conclusions.

To perform our inversion, we interpolate linearly the values of the amplitude and phase over the grid used to compute the radial eigenfunctions and the Jacobian polynomials. We assume that the amplitude is zero at the poles and increases linearly toward $60^{\circ}$ of latitude whereas the phase is constant poleward of $60^{\circ}$ of latitude.

The divergence of the angular momentum flux perturbation  ${\vec \tau}_{1}$ can be obtained from Eq. (40). We define its amplitude and phase as:

$$% A_{\tau} \equiv \sqrt{[\nabla \cdot {\vec \tau}_{1}^{(c)}]^{2} + [\nabla \cdot {\vec \tau}_{1}^{(s)}]^{2}},$$ (71)

$$% \Phi_{\tau} \equiv \tan^{-1} \left( \frac{\nabla \cdot {\ve... ...u}_{1}^{(c)}}{\nabla \cdot {\vec \tau}_{1}^{(s)}} \right)\cdot$$ (72)

They are plotted in Figs. 5 and 6, respectively, in the case of $\eta _{0}=2.56$ $\times$ 10¹² g cm^-1 s^-1. A suitable smoothing has been applied to have a resolution of $\sim$0.05 $R_{\odot }$ in the radial coordinate and $\approx$ $15^{\circ}$ in latitude. In order to show the effect of a smaller value of the turbulent viscosity, we plot in Figs. 7 and 8 the isocontours of $A_{\tau }$ and $\Phi_{\tau}$ for $\eta_{\rm t} = 5.62$ $\times$ 10¹⁰ g cm^-1 s^-1. The amplitude of the perturbation term is higher close to the equator because the data in Fig. 4 of Howe et al. (2006) mainly sample the low-latitude branch of the torsional oscillations. The relative maxima of $A_{\tau }$ are reached close to the base of the convection zone and at a radius of $\sim$0.85 $R_{\odot }$ and their locations show only a minor dependence on the value of $\eta_{0}$. Conversely, the value of $\eta_{0}$ significantly affects  $\Phi_{\tau}$, as it follows from Eqs. (39). When $\eta_{0}$ is large, the timescales for angular momentum exchange, as given by the inverse of the lowest eigenvalues, are significantly shorter than the eleven-year cycle and the perturbations are almost in phase with the angular velocity variations. When $\eta_{0}$ is sufficiently low, the timescales for angular momentum exchange become comparable or longer than the eleven-year cycle so the torsional oscillations lag behind the perturbations. Considering Fig. 6, we see that the phase is in agreement with that in Fig. 4 of Howe et al. (2006), except for $r \ga 0.93$ $R_{\odot }$. The disagreement in those layers is due to the small values of $\nabla \cdot {\vec \tau}_{1}$ close to the surface that makes the corresponding phases uncertain (cf. Fig. 5). The phase lag is apparent by comparing Fig. 8 with Fig. 6, especially close to the equator for $0.72 \leq r/R_{\odot} \leq 0.92$. It is less evident near the base of the convection zone and in the upper layers due to the smaller values of $\nabla \cdot {\vec \tau}_{1}$ in those regions.

	Figure 5: The isocontours of the function $A_{\tau } (r, \mu )$ as defined by Eq. (71) in the case of $\eta _{0}=2.56$ $\times$ 1012 g cm-1 s-1. The scale on the left indicates the ranges corresponding to the different colors and is in units of 105 g cm-1 s-2. The relative statistical uncertainty of $A_{\tau }$ is $\sim$30%, as it follows from the relative uncertainty of the data.
Open with DEXTER

	Figure 6: The isocontours of $\Phi _{\tau } (r, \mu )$ as defined by Eq. (72) in the case of $\eta _{0}=2.56$ $\times$ 1012 g cm-1 s-1. The phase ranges from $-\frac{\pi}{2}$ to $\frac{\pi}{2}$. The scale on the left indicates the phase ranges corresponding to the different colors and is in radians.
Open with DEXTER

	Figure 7: As for Fig. 5 in the case of $\eta _{0}=5.62$ $\times$ 1010 g cm-1 s-1. The scale on the left is in units of 103 g cm-1 s-2. The relative statistical uncertainty of $A_{\tau }$ is about 30%.
Open with DEXTER

	Figure 8: As for Fig. 6 in the case of $\eta _{0}=5.62$ $\times$ 1010 g cm-1 s-1.
Open with DEXTER

It is interesting to note that the dependence of the amplitude and phase of the torsional oscillations on the turbulent viscosity can lead to its estimate in the framework of mean-field models (e.g., Rüdiger et al. 1986; Rüdiger 1989). Specifically, Rempel (2007) finds that a mean turbulent kinematic viscosity about one order of magnitude smaller than the mixing-length estimate is needed to reproduce the polar branch of the torsional oscillations.

Lower limits for the modulus of the angular momentum flux vector  $\vert{\vec \tau}_{1}\vert$ can be derived from Eq. (46). From the lower limits on $\vert{\vec \tau}_{1}^{(c)}\vert$ and $\vert{\vec \tau}_{1}^{(s)}\vert$, a lower limit on $\vert{\vec \tau}_{1}\vert= \sqrt{[{\vec \tau}_{1}^{(c)}]^{2} + [{\vec \tau}_{1}^{(s)}]^{2}}$ can be derived and it is plotted in Figs. 9 and 10 for the radial and latitudinal localizations described in Sect. 3.2, respectively. Given the uncertainty in the penetration depth of the torsional oscillations, in addition to the results for a penetration down to the base of the convection zone, we also plot those for penetration depths of 0.80 and 0.90 $R_{\odot }$, respectively. They are obtained simply by assuming that the oscillation amplitude has the value in Fig. 4 of Howe et al. (2006) above the penetration depth and drops to zero immediately below it.

When the oscillations penetrate down to the base of the convection zone, the maximum of the perturbation term is reached between 0.80 and 0.85 $R_{\odot }$. It is mainly localized into two latitude zones, i.e., within $\pm$ $15^{\circ}$ from the equator and between $30^{\circ}$ and $60^{\circ}$. Note that our data refer mainly to the low-latitude branch of the oscillations, so the possible source at latitude > $60^{\circ}$ cannot be detected. When the turbulent dynamical viscosity is reduced with respect to the mixing-length value, the amplitude of the perturbation term drops, but its radial and latitudinal localizations are not greatly affected, except for a shift of the nearly equatorial band towards higher latitudes.

According to the hypothesis that the Lorentz force is the only source of the torsional oscillations, the amplitude of the Maxwell stress $B_{\rm r}B_{\phi}$ can be estimated from the lower limit of $\vert{\vec \tau}_{\rm 1r}\vert$ and is $\sim$2.7 $\times$ 10⁵ G² in the case of $\eta _{0}=2.56$ $\times$ 10¹² g cm^-1 s^-1. Considering that the poloidal field $B_{\rm r}$ is of the order of 1-10 G, this leads to very high toroidal fields in the bulk of the convection zone, which would be highly unstable because of magnetic buoyancy. On the other hand, with $\eta _{0}=5.62$ $\times$ 10¹⁰ g cm^-1 s^-1, the Maxwell stress is $B_{\rm r}B_{\phi} \simeq 8.1$ $\times$ 10³ G² which leads to a toroidal field intensity of the order of 10³ G. It may be stably stored for timescales comparable to the solar cycle thanks to the effects of the downward turbulent pumping in the convection zone (Brandenburg 2005, and references therein).

Note that the maximum radial stress $\vert B_{\rm r} B_{\phi} \vert$ is reduced by a factor of $\sim$30 when $\eta_{\rm t}$ is decreased from 2.56 $\times$ 10¹² to 5.62 $\times$ 10¹⁰ g cm^-1 s^-1, whereas the maximum of $\vert B_{\theta} B_{\phi}\vert$ is reduced only by a factor of $\sim$4 (Figs. 9 and 10). This mainly reflects the predominance of the radial gradient over the latitudinal gradient of the angular velocity perturbation in the deeper layers of the solar convection zone.

The average phase lags $\Pi _{\rm r}$ and $\Pi _{\theta }$ between $B_{\rm r}$ and $B_{\theta }$ and the azimuthal field $B_{\phi }$, as derived by the method in Sect. 2.6, are plotted in Figs. 11 and 12, respectively. It is interesting to note that $B_{\rm r} B_{\phi} > 0$ below $\sim$0.85 $R_{\odot }$ when $\eta_{\rm t} = 5.62$ $\times$ 10¹⁰ g cm^-1 s^-1, in agreement with the finding of most dynamo models in which the azimuthal field is produced by the stretching of the radial field in the low-latitude region, where $\frac{\partial \Omega}{\partial r} > 0$ for r < 0.95 $R_{\odot }$(e.g., Schlichenmaier & Stix 1995; Rüdiger & Hollerbach 2004). Above $\sim$0.85 $R_{\odot }$, the phase relationship between the radial and the azimuthal fields leads to $B_{\rm r} B_{\phi} < 0$ with a phase lag of $\sim$$\pi$, in agreement with the early finding that the photospheric zone equatorward of the activity belt is rotating faster than that at higher latitudes (Rüdiger 1989). Note that $\frac{\partial \Omega}{\partial r}$ becomes negative for r > 0.95 $R_{\odot }$ at low latitudes, leading to a reversal of the phase relationship between the two field components in the outer layers of the solar convection zone. Our limited spatial resolution and the uncertainty of the measurements of the torsional oscillations inside the Sun might explain why we find the transition from a mostly positive to a negative $B_{\rm r}B_{\phi}$ at $\sim$0.85 $R_{\odot }$.

Figure 9: Upper panel: the lower limit of the amplitude of the perturbation term  $\vert{\vec \tau}_{\rm 1r}\vert$ averaged over spherical shells of thickness 0.05 $R_{\odot }$ versus the relative radius; different linestyles and colors refer to the depth at which the torsional oscillations are assumed to vanish: black solid line - oscillations extending down to the base of the convection zone; green dotted line - oscillations extending down to 0.80 $R_{\odot }$; red dashed line - oscillations extending down to 0.90 $R_{\odot }$. Results are obtained assuming a turbulent dynamical viscosity at the base of the convection zone $\eta _{0}=2.56$ $\times$ 1012 g cm-1 s-1. Lower panel: as in the upper panel, but for $\eta _{0}=5.62$ $\times$ 1010 g cm-1 s-1. The relative statistical uncertainties are in all cases about 30%, as follows from the uncertainties of the data.

Open with DEXTER

Figure 10: Upper panel: the lower limit of the amplitude of the perturbation term  $\vert{\vec \tau}_{1\theta}\vert$ averaged over different latitude zones versus their average value of $\mu$. Different linestyles and colors are used to indicate the results for different penetration depths of the torsional oscillations, as in Fig. 9. Results are obtained assuming a turbulent dynamical viscosity at the base of the convection zone $\eta _{0}=2.56$ $\times$ 1012 g cm-1 s-1. Lower panel: as in the upper panel, but for $\eta _{0}=5.62$ $\times$ 1010 g cm-1 s-1. The relative statistical uncertainties are in all cases about 30%, as follows from the uncertainties of the data.

Open with DEXTER

The phase lag between $B_{\theta }$ and $B_{\phi }$ depends remarkably on the colatitude for $\eta_{\rm t} = 2.56$ $\times$ 10¹² g cm^-1 s^-1, whereas it is almost constant for $\eta_{\rm t} = 5.62$ $\times$ 10¹⁰ g cm^-1 s^-1, leading in the latter case mostly to $B_{\theta} B_{\phi} < 0$. This, together with $\frac{\partial \Omega}{\partial \theta} > 0$, suggests that the toroidal field is mainly produced by the stretching of the radial field.

On the other hand, if the angular momentum transport leading to the torsional oscillations is produced only by a perturbation of the meridional flow, as in the thermal wind model, then a very small perturbation follows from the lower limit of $\vert{\vec \tau}_{1}\vert$. For $\eta _{0}=2.56$ $\times$ 10¹² g cm^-1 s^-1, we find a minimum amplitude of the meridional flow component oscillating with the eleven-year cycle of $\approx$3 cm s^-1 at a depth of 0.85 $R_{\odot }$ and a latitude of $45^{\circ}$. Note, however, that if we estimate $\vert{\vec \tau}_{1}\vert$ from its divergence and the typical lengthscale of its variations in Fig. 5, we find a value about one order of magnitude larger, that is in agreement with the estimate by Rempel (2007). A similar argument applies also to the estimate of the Maxwell stresses made above.

Figure 11: Upper panel: the phase lag $\Pi _{\rm r}$ between $B_{\rm r}$ and $B_{\phi }$ as derived by Eqs. (54) averaged over spherical shells of thickness 0.05 $R_{\odot }$ versus the relative radius. Different linestyles are used to indicate the results for different penetration depths of the torsional oscillations, as in Fig. 9. Results are obtained assuming a turbulent dynamical viscosity at the base of the convection zone $\eta _{0}=2.56$ $\times$ 1012 g cm-1 s-1. Lower panel: as in the upper panel, but for $\eta _{0}=5.62$ $\times$ 1010 g cm-1 s-1.

Open with DEXTER

Figure 12: Upper panel: the phase lag $\Pi _{\theta }$ between $B_{\theta }$ and $B_{\phi }$, as derived by equations analogous to (54), averaged over different latitude zones versus their average value of $\mu$. Different linestyles are used to indicate the results for different penetration depths of the torsional oscillations, as in Fig. 9. Results are obtained assuming a turbulent dynamical viscosity at the base of the convection zone $\eta _{0}=2.56$ $\times$ 1012 g cm-1 s-1. Lower panel: as in the upper panel, but for $\eta _{0}=5.62$ $\times$ 1010 g cm-1 s-1.

Open with DEXTER

The average kinetic energy variation associated with the torsional oscillations, as computed from Eq. (49), is $\langle {\cal T}\rangle _{\rm c} = 4.88$ $\times$ 10²⁸ erg, whereas the average dissipated power is 6.7 $\times$ 10²⁶ erg s^-1, when $\eta _{0}=2.56$ $\times$ 10¹² g cm^-1 s^-1, and only 1.5 $\times$ 10²⁵ erg s^-1, when $\eta_{0} = 5.6$ $\times$ 10¹⁰ g cm^-1 s^-1.

When the torsional oscillations are assumed to be confined to shallower and shallower layers, the location of the perturbation term is shifted closer and closer to the surface and it becomes more and more uniformly distributed in latitude. Its amplitude shows a remarkable decrease because of the smaller moment of inertia of the surface layers.

4 Application to solar-like stars

The results of Sect. 2.7 can be applied to the variation of the surface differential rotation observed in, e.g., LQ Hya between 1996.99 and 2000.96 (see Table 2 of Donati et al. 2003). Assuming an internal structure analogous to the solar one and that the differential rotation variations observed at the surface are representative of those at a depth of 0.99 R, we can estimate a lower limit for $\vert \nabla \cdot {\vec \tau}_{1}\vert$ in the overshoot layer from Eq. (59). This is used to estimate the Maxwell stresses by assuming that: $r_{\rm o} \frac{B_{\rm p}B_{\phi}}{\tilde{\mu}} \simeq \delta r \vert \nabla \cdot {\vec \tau}_{1}\vert$, where $r_{\rm o} = 0.67$ R is the radius of the overshoot layer and $\delta r=0.04$ R its thickness. As such, the minimum magnetic field strength is: $B_{\rm min} = \sqrt{B_{\rm p}B_{\phi} } \sim 2500$ G. However, if we assume a poloidal field strength of $\sim$100 G, as indicated by the Zeeman Doppler imaging, we find an azimuthal field of $B_{\phi} \sim 6.2$ $\times$ 10⁴ G, which is in the range of the values estimated by Lanza (2006a) in the framework of the Taylor-Proudman hypothesis. Note that the result is independent of $\eta_{0}$ in the limit $t_{\rm DR} \ll \lambda_{nk}^{-1}$. Although this application is purely illustrative, it suggests that our method can be applied to derive estimates of the internal magnetic torques (as well as other sources of angular momentum transport) when asteroseismic results are available, in combination with surface rotation measurements, to further constrain rotation variations in solar-like stars.

5 Discussion and conclusions

We introduced a general solution of the angular momentum transport equation that takes into account the density stratification and the radial dependence of the turbulent viscosity  $\eta_{\rm t}$. Its main limitation is due to the uncertainty of the turbulent viscosity in stellar convection zones. It is interesting to note that the method of the separation of variables to solve Eq. (3) can be applied also when $\eta_{\rm t}$ is the product of a function of the radius by one of the latitude. However, when $\eta_{\rm t}$ depends on the latitude, the angular eigenfunctions are no longer Jacobian polynomials.

Our formalism can be applied to compute the response of a turbulent convection zone to prescribed time-dependent Lorentz force and meridional circulation. From the mathematical point of view, it is a generalization of that of Rüdiger et al. (1986) and can be easily compared with it by considering that: $P_{n}^{(1,1)} (\mu) = \frac{2}{n+2} \frac{{\rm d}P_{n+1}}{{\rm d}\mu} = - \frac{2 \sin \theta}{n+2} P_{n}^{1}(\theta)$, where P_n is the Legendre polynomial of degree n and $P_{n}^{1} = -\frac{{\rm d}P_{n}}{{\rm d}\theta}$. Moreover, our method can be used to estimate the torques leading to the angular momentum redistribution within the solar (or stellar) convection zone, thus generalizing the approach suggested by Komm et al. (2003). The main limitation, in addition to the uncertainty of the turbulent viscosity, comes from the low resolution and the limited accuracy of the present data on solar torsional oscillations, particularly in the deeper layers of the convection zone. In fact, these layers are the most important because torsional oscillations with an amplitude of $\sim$0.5%-1% of the solar angular velocity extending down to the base of the convection zone lead to Maxwell stresses with an intensity of at least $\approx$8 $\times$ 10³ G² around 0.85 $R_{\odot }$ or a perturbation of the order of several percents of the meridional flow speed at the same depth (e.g., Rempel 2007). If the torsional oscillations are due solely to the Maxwell stresses associated with the mean field of the solar dynamo, we can also estimate the phase relationships between $B_{\rm r}$, $B_{\theta }$ and $B_{\phi }$. Our preliminary results indicate that $B_{\rm r} B_{\phi} > 0$ in the layers below $\approx$0.85 $R_{\odot }$ and $B_{\rm r} B_{\phi} < 0$ in the outer layers which, together with helioseismic measurements of the internal angular velocity, suggest that the toroidal field is mainly produced by the stretching of the poloidal field by the radial shear.

Future helioseismic measurements may improve our knowledge of the torsional oscillations, essentially by extending the time series of the data or by means of space-borne instruments, such as those foreseen for the Solar Dynamic Observatory (e.g., Howe et al. 2006). On the other hand, asteroseismic measurements may open the possibility of investigating similar phenomena in solar-like stars, particularly in those young, rapidly rotating objects showing variations of the angular velocity one or two orders of magnitude larger than the Sun.

Acknowledgements

The author wishes to thank an anonymous Referee for valuable comments and Professor G. Rüdiger for interesting discussion. Solar physics and active star research at INAF-Catania Astrophysical Observatory and the Department of Physics and Astronomy of Catania University is funded by MIUR (Ministero dell'Università e della Ricerca), and by Regione Siciliana, whose financial support is gratefully acknowledged.
This research has made use of the ADS-CDS databases, operated at the CDS, Strasbourg, France.

Appendix A: Proof of convergence of the Green function series

The convergence of the Green function series in Eqs. (34) can be studied by considering the asymptotic formulae for $\zeta_{nk}$ and the Jacobian polynomials, i.e., for $n \gg 1$ and $k \gg 1$. In the asymptotic limit, those series can be written as:

 

$$% \sum_{n} \sum_{k} \frac{F_{n}}{\lambda_{nk}} \zeta_{nk}(r) \zet... ...}^{(1,1)} (\mu^{\prime}) = {\sum_{n} h_{n} ~ \sqrt{F_{n}} P_{n}^{(1,1)} (\mu)},$$ (A.1)

where:

$$% h_{n} \equiv \sum_{k} \frac{\sqrt{F_{n}}}{\lambda_{nk}} \zeta_{nk}(r) \zeta_{nk}(r^{\prime}) P_{n}^{(1,1)} (\mu^{\prime}).$$ (A.2)

From the asymptotic formulae it follows that for a point in the domain $[r_{\rm b}, R_{\odot}]$ $\times$ $[r_{\rm b}, R_{\odot}]$ $\times$ ]-1, 1[, the quantity $\left\vert\! \sqrt{F_{n}} \zeta_{nk}(r) \zeta_{nk}(r^{\prime}) P_{n}^{(1,1)} (\mu^{\prime}) \right\vert$ is limited with an upper bound independent of n and k. Therefore, the series (A.2) that defines the coefficient h_n is uniformly convergent in that domain if the series $\sum_{k} \frac{1}{\lambda_{nk}}$ converges.

This can be proven by considering the inequalities in (21) and the formula:

$$% \sum_{k=1}^{\infty} \frac{z}{a k^{2} + y} = \frac{1}{2} \le... ...ft[ g\left( \frac{y}{a}\right) \sqrt{\frac{y}{a}} - 1 \right],$$ (A.3)

where z, $a \not= 0$ and y are real numbers and the function g(x) is defined as:

$$% g(x) \equiv \pi \left[ \frac{\exp(\pi x) + \exp(-\pi x)}{\exp(\pi x) - \exp(-\pi x)} \right]\cdot$$ (A.4)

Equation (A.3) follows from the equality (1.217.1) of Gradshteyn & Ryzhik (1994). In this way, we find:

 

$$\frac{m}{2n(n+3)Q}\left\{ g\left[ \frac{n(n+3) Q l^{2}}{\pi^{2} P... ... q l^{2}}{\pi^{2} p}\right] \sqrt{\frac{n(n+3) q l^{2}}{\pi^{2}p}} -1 \right\}.$$ (A.5)

This result indicates that $h_{n} \sim \frac{1}{\sqrt{n(n+3)}}$ for $n \gg 1$. Since the eigenfunctions $\sqrt{F_{n}} P_{n}^{(1,1)}(\mu)$ form an orthonormal set, the convergence of the series in (A.1) follows by the Riesz-Fisher Theorem given that the series $\sum_{n} h_{n}^{2} \sim \sum_{n} \frac{1}{n(n+3)}$ converges. As such, the uniform convergence of the series in Eqs. (34) in the domain $[r_{\rm b}, R_{\odot}] \times [r_{\rm b}, R_{\odot}] \times ]-1, 1[ \times ]-1, 1[$ is proven.

References

• Bailey, P. B., Everitt, W. N., & Zettl, A. 2001, ACM Trans. Math. Software, 21, 143 (In the text)
• Basu, S., & Antia, H. M. 2001, ApJ, 559, L67 [NASA ADS] [CrossRef]
• Basu, S., & Antia, H. M. 2003, ApJ, 585, 553 [NASA ADS] [CrossRef] (In the text)
• Brandenburg, A. 2005, ApJ, 625, 539 [NASA ADS] [CrossRef] (In the text)
• Christensen-Dalsgaard, J., Däppen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [NASA ADS] [CrossRef] (In the text)
• Corbard, T., & Thompson, M. J. 2005, Sol. Phys., 205, 211 [NASA ADS] [CrossRef] (In the text)
• Covas, E., Moss, D., & Tavakol, R. 2004, A&A, 416, 775 [NASA ADS] [CrossRef] [EDP Sciences]
• Covas, E., Moss, D., & Tavakol, R. 2005, A&A, 429, 657 [NASA ADS] [CrossRef] [EDP Sciences]
• Donati, J.-F., Collier Cameron, A., & Petit, P. 2003, MNRAS, 345, 1187 [NASA ADS]
• Gradshteyn, I. S., & Ryzhik, I. M. 1994, Table of Integrals, Series, Products, 5th edn., ed. A. Jeffrey (San Diego: Academic Press) (In the text)
• Howard, R., & LaBonte, B. 1980, ApJ, 239, L33 [NASA ADS] [CrossRef] (In the text)
• Howe, R., Christensen-Dalsgaard, J., Hill, F., et al. 2000, Science, 287, 2456 [NASA ADS] [CrossRef] (In the text)
• Howe, R., Christensen-Dalsgaard, J., Hill, F., et al. 2005, ApJ, 634, 1405 [NASA ADS] [CrossRef]
• Howe, R., Rempel, M., Christensen-Dalsgaard, J., et al. 2006, ApJ, 649, 1155 [NASA ADS] [CrossRef]
• Jeffers, S. V., Donati, J.-F., & Collier Cameron, A. 2007, MNRAS, 375, 567 [NASA ADS]
• Kitchatinov, L. L., Pipin, V. V., Makarov, V. I., & Tlatov, A. G. 1999, Sol. Phys., 189, 227 [NASA ADS] [CrossRef] (In the text)
• Kichatinov, L. L., Pipin, V. V., & Rüdiger, G. 1994, Astron. Nachr., 315, 157 [NASA ADS] (In the text)
• Komm, R., Howe, R. Durney, B. R., & Hill, F. 2003, ApJ, 586, 650 [NASA ADS] [CrossRef] (In the text)
• Küker, M., Rüdiger, G., & Pipin, V. V. 1996, A&A, 312, 615 [NASA ADS] (In the text)
• Lanza, A. F. 2006a, MNRAS, 373, 819 [NASA ADS] (In the text)
• Lanza, A. F. 2006b, MNRAS, 369, 1773 [NASA ADS] (In the text)
• Lochard, J., Samadi, R., & Goupil, M. J. 2005, A&A, 438, 939 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics (New York: McGraw-Hill) (In the text)
• Parker, E. N. 1993, ApJ, 408, 707 [NASA ADS] [CrossRef] (In the text)
• Rempel, M. 2006, ApJ, 647, 662 [NASA ADS] [CrossRef]
• Rempel, M. 2007, ApJ, 655, 651 [NASA ADS] [CrossRef]
• Rüdiger, G. 1989, Differential rotation and stellar convection: Sun and solar-type stars (Gordon & Breach Sci. Publ.)
• Rüdiger, G., & Hollerbach, R. 2004, The magnetic Universe (Weinheim: Wiley-VCH)
• Rüdiger, G., & Kichatinov, L. L. 1990, A&A, 236, 503 [NASA ADS] (In the text)
• Rüdiger, G., Tuominen, I., Krause, F., & Virtanen, H. 1986, A&A, 166, 306 [NASA ADS]
• Schlichenmaier, R., & Stix, M. 1995, A&A, 302, 264 [NASA ADS] (In the text)
• Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [NASA ADS] [CrossRef]
• Schüssler, M. 1981, A&A, 94, L17 [NASA ADS] (In the text)
• Smirnov, V. I. 1964a, A Course of Higher Mathematics, Vol. III, Part II (Oxford: Pergamon Press) (In the text)
• Smirnov V. I., 1964b, A Course of Higher Mathematics, Vol. IV (Oxford: Pergamon Press) (In the text)
• Spruit, H. C. 2003, Sol. Phys., 213, 1 [NASA ADS] [CrossRef] (In the text)
• Toomre, J., Christensen-Dalsgaard, J., Howe, R., et al. 2000, Sol. Phys., 192, 437 [NASA ADS] [CrossRef]
• Vorontsov, S. V., Christensen-Dalsgaard, J., Schou, J., Strakhov, V. N., & Thompson, M. J. 2002, Science, 296, 101 [NASA ADS] [CrossRef] (In the text)
• Yoshimura, H. 1981, ApJ, 247, 1102 [NASA ADS] [CrossRef] (In the text)