A&A 402, 683-692 (2003)
DOI: 10.1051/0004-6361:20030239

Asteroseismic determination of stellar angular momentum

F. P. Pijpers1

Theoretical Astrophysics Center (TAC), Department of Physics and Astronomy, Aarhus University, Ny Munkegade, 8000 Århus C, Denmark

Received 29 November 2002 / Accepted 14 February 2003

The total angular momentum of stars is an important diagnostic for testing theories of formation of stars. Angular momentum can also play a decisive role in the evolution of stars, in particular towards the higher masses. Asteroseismology can play an important role in providing measurement of this parameter. The spectrum of stellar oscillations is affected by rotation in a way similar to the splitting of spectral lines through the magnetic Zeeman effect. It is shown here that this can be used to determine accurately the total angular momentum of stars, taking full differential rotation into account, and without any dependence of inclination angle of rotation axis with respect to the line of sight. All aspects of this technique are addressed, including the issue of mode identification, determination of the moment of inertia, and the influence of the reference model on the seismic inference. It is a realistic prospect to do asteroseismogyrometry with planned space satellite missions such as MOST, MONS, COROT, and Eddington.

Key words: stars: oscillations - stars: evolution - stars: rotation

1 Introduction

The aim of asteroseismology is to use the properties of the oscillations of stars to place constraints on their internal structure. The frequencies of stellar oscillations can be used to determine parameters of stars that are either inaccessible to classical observation or can not be determined with an accuracy that is sufficient for many purposes. This paper concentrates on one property of stars only accessible through seismic observations, which is the internal rotation rate. An important constraint for theories of mechanisms mediating angular momentum loss during star formation would be the measurement of the angular momentum of stars for a range of masses (cf. Bouvier et al. 1997; see also Bodenheimer 1995 and references therein). Also, for the post-main sequence evolution of in particular high-mass stars the loss of angular momentum during evolution is similarly important (cf. Maeder & Meynet 2002). Spectroscopic observations have access to only the projected surface rotation velocity.

Previous asteroseismic work on rotation has concentrated on the possibility of measuring differential rotation in $\delta$ Scuti stars (Goupil et al. 1996) and in white dwarf stars (Kawaler et al. 1999). $\delta$ Scuti stars are multimode pulsators. In principle two separate driving mechanisms could operate in these stars (cf. Samadi et al. 2002); the opacity driving mechanism that also operates in large amplitude pulsators such as Cepheids and RR Lyrae variables, and the turbulent driving by convection that is known from the Sun. Because of this Goupil et al. (1996) assumed that modes with l up to 4 or 6 might be detectable in these stars, where l is the number of node lines of the pulsation mode across the stellar surface. In stars more similar to the Sun, where only the turbulent driving mechanism is expected to operate, the oscillation amplitudes are likely to be smaller because of stronger damping effects. This means that the cancellation in the signal, due to averaging over the stellar surface, for the higher l values may limit the range in l for which a signal is measurable to $l\leq4$ or even $l\leq 2$. In the best studied oscillating white dwarf so far, PG 1159 (Kawaler et al. 1999), only l=1 modes could be used for the determination of an internal rotation rate, with a marginal detection of radially differential rotation. It is therefore of interest to investigate what measures of internal rotation can be determined given more pessimistic assumptions than those of Goupil et al. (1996).

The angular momentum of stars is one important example of a parameter that is of astrophysical interest and that can be determined with reasonable accuracy even from relatively little data. This paper presents a method for the measurement of the total angular momentum of stars from asteroseismology: asteroseismogyrometry.

Section 2 of this paper gives a brief overview of the relationship between oscillation frequencies and rotation, shows how the angular momentum can be inferred, and presents tests of this method using solar data. Section 3 presents a method for determining the moment of inertia, and also discusses issues concerning mode identification and the choice of reference model.

2 The relationship between oscillation frequencies and rotation

2.1 Asymptotic relationship

The oscillation modes of stars are fully characterised by the radial node number n and the degree l and order m of spherical harmonics which specify the pattern of node lines over surfaces concentric with the star. For stars which are nearly spherically symmetric, the oscillation frequencies of the modes are distributed in multiplets labelled by n and l, with the individual frequencies within each multiplet corresponding to different values of m. Each central frequency of an oscillation multiplet of a star represents a weighted average of internal properties of the star, usually represented by two thermodynamic variables such as sound speed $c_{\rm s}$ and density $\rho$. The weighting function of the average depends uniquely on the eigenfunction of the mode of oscillation which is fully characterised by n, l and m. However, in a simple approximate relationship which improves in precision asymptotically for large values of n, the frequency of p-mode (sound-wave like) oscillations can be expressed as:
$\displaystyle \nu_{nlm} = \left(n + {l\over 2} + {1\over 4} + \alpha\right)\Del...
+ m{\displaystyle\langle\Omega\rangle_{nl}\over\displaystyle 2\pi},$     (1)

where $\Delta\nu$ is known as the large separation and the small separation is defined as $\nu_{nl}-\nu_{n-1,~l+2}\approx (4l+6)D_0 $. The large and small separation can be used as measures of mass and evolutionary age of the star (cf. Christensen-Dalsgaard 1993b). The rotational term $\langle\Omega\rangle_{nl}/2\pi$is responsible for the splitting up of multiplets: modes with different m for the same n and l. For the Sun this term is in the range of $400{-}450~{\rm nHz}$. The factors $\alpha$ and $\delta$ are introduced to account for near-surface effects. These cause a frequency shift that depends on the frequency $\nu$and can be filtered out before further analysis of the frequencies (cf. Pérez Hernández & Christensen-Dalsgaard 1994; Basu et al. 1996).

2.2 Perturbation analysis

In a spherically symmetric, non-rotating star, the frequency of an eigenmode would be independent of m and thus there would be multiplets of (2l+1) modes with identical frequencies, each multiplet corresponding to an (n,l) pair. Rotation lifts this (2l+1)-fold degeneracy. The difference in frequency between modes in the same multiplet is called the (rotational) frequency splitting. Since the oscillations propagate through large parts of the interior, the frequency splittings are determined by the rotation rate inside the star and can be used to probe this internal rotation by solving an inverse problem. In particular, they enable one to perform 2-D inversions for the rotation rate as a function of radius and latitude (cf. Thompson et al. 1996). The oscillation displacement eigenfunction is

 \begin{displaymath}\left( \xi_{nl}(r), L^{-1}\eta_{nl}{\partial\over\partial\the...\partial\phi}\right)
P_l^m(\cos\theta){\rm e}^{{\rm i}m\phi}
\end{displaymath} (2)

with respect to spherical polar coordinates $(r,\theta,\phi)$: here $\xi_{nl}$ and $\eta_{nl}$ are calculable functions, given a stellar model;  $L = {\sqrt{l(l+1)}}$; and Plm is the associated Legendre polynomial of degree l and order m. For the Sun and for stars with $\langle\Omega\rangle/2\pi < 2 D_0$ linear perturbation theory provides accurate oscillation frequencies. According to linear perturbation theory the difference between frequencies of modes with the same values of n and l, but opposite m, is given in terms of the eigenfunctions (2) of the non-rotating star by

\begin{displaymath}2\pi\left(\nu_{nlm} - \nu_{nl~-m}\right) = 2m {\rm d}_{nlm},
\end{displaymath} (3)


 \begin{displaymath}d_{nlm} = \int\limits_{-1}^1 \int\limits_0^1 K_{nlm}(r,\theta)
\Omega(r,\theta) {\rm d}r~{\rm d}\cos\theta.
\end{displaymath} (4)

Here and in the following, the radial variable r is made dimensionless by dividing it by the surface radius. The kernels Knlm are given by:

 \begin{displaymath}K_{nlm} (r, \theta) =\left( F_1^{nl}(r) G_1^{lm}(\theta)\;+\;
F_2^{nl}(r) G_2^{lm}(\theta)\right)/ I_{nl},
\end{displaymath} (5)

                                Fnl1 $\textstyle \equiv$ $\displaystyle \rho(r) r^2 \left[ \xi_{nl}^2(r) - 2L^{-1} \xi_{nl}(r)
\eta_{nl} (r) + \eta_{nl}^2(r) \right]$  
Fnl2 $\textstyle \equiv$ $\displaystyle \rho(r) r^2 \left[ \eta_{nl}^2(r)L^{-2} \right] ,$ (6)

$\rho$ being the density and $\xi_{nl}$, $\eta_{nl}$ the components of the displacement eigenfunction in the non-rotating star, and

 \begin{displaymath}I_{nl} = \int\limits_0^1 {\rm d}r~ \rho r^2 \left(\xi_{nl}^2 +
\eta_{nl}^2 \right);
\end{displaymath} (7)

                          Glm1 $\textstyle \equiv$ $\displaystyle { (l - \vert m \vert)! \over (l + \vert m \vert)!} (l+1/2)
\left[ P_l^m(u) \right]^2$  
Glm2 $\textstyle \equiv$ $\displaystyle {1\over 2} (1-u^2) {{\rm d}^2 G^{lm}_1\over{\rm d}u^2}$ (8)


\begin{displaymath}u \equiv \cos\theta .
\end{displaymath} (9)

A different but equivalent way of describing the effect of rotation on the oscillation frequencies is through so-called a-coefficients. The a-coefficients are constructed by fitting polynomials in m to ridges of power in an $(m, \nu)$ diagram for each n and l. Details of their use can be found in e.g. Schou et al. (1994). The relationship between a-coefficients and mode frequencies is:

\begin{displaymath}\nu_{nlm} = \sum_{s=0}^{s_{\rm max}} a_{nl s} {\cal P}_{s}^{(l)} (m).
\end{displaymath} (10)

Here the polynomials ${\cal P}_{s}^{(l)}$ have degree s and are often normalised so that ${\cal P}_{s}^{(l)} (l) = l$. For the purposes of measuring the rotation rate only the odd a-coefficients are of interest since the even order a-coefficients are sensitive only to effects that are antisymmetric around the equator. As noted by Schou et al. (1994), who follow the projection proposed by Ritzwoller & Lavely (1991), these properties together with the orthogonality condition

\begin{displaymath}\sum_{m=-l}^l {\cal P}_{s}^{(l)} (m) {\cal P}_{s'}^{(l)} (m) =0 \quad
{\rm for} s\neq s'
\end{displaymath} (11)

specify the polynomials ${\cal P}_{s}^{(l)}$ completely. As noted in Appendix A of Schou et al. (1994) the polynomials are usually constructed in an iterative process. However explicit expressions exist as shown by Pijpers (1997), and the equivalent expression of Eq. (4) is:

 \begin{displaymath}2\pi a_{nl~ 2s+1} = \int\limits_{0}^{1} {\rm d}r \int\limits_{-1}^1
{\rm d}u~ K_{nl s}(r, u) \Omega (r, u),
\end{displaymath} (12)

where the "mode kernels'' are now given by:

 \begin{displaymath}K_{nls} (r, u) =\left( F_1^{nl}(r) G_1^{ls}(u)\;+\;
F_2^{nl}(r) G_2^{ls}(u)\right)/ I_{nl} ,
\end{displaymath} (13)

in which the F functions are given by Eq. (6), the factor Inl is given by Eq. (7), and the G functions are defined by:
                            G1ls $\textstyle \equiv$ $\displaystyle -{1\over 2}{(4s+3) v^{(l)}_{2s+1}\over
(2s+2)(2s+1)} \left(1-u^2\right)^{1\over 2} P_{2s+1}^1 (u)$  
G2ls $\textstyle \equiv$ $\displaystyle {1\over 4}v^{(l)}_{2s+1} (4s+3) \left(1-u^2\right)^{1\over 2}
P_{2s+1}^1 (u).$ (14)

The factors v(l)2s+1 are given by:
                                      $\displaystyle v^{(l)}_{2s+1} = (-1)^{s} {1\over l}{\displaystyle (2l+1)! (2s+2)!
(l+s+1)! \over\displaystyle s!(s+1)! (l-s-1)! (2l+2s+2)!}$  
    $\displaystyle {\rm for}~ s \leq l-1 .$ (15)

For s=0 this reduces to v(l)1 = 1. With these equations a 2-D inverse problem is defined which can be solved using methods identical to those used to invert for the solar rotation rate from individual frequency splittings (cf. Pijpers & Thompson 1992; Pijpers & Thompson 1994; Pijpers 1996).

2.3 The angular momentum

It was shown by Pijpers (1998) that it is not necessary to first obtain the resolved rotation rate $\Omega(r,\theta)$ in order to measure, among other quantities, the total angular momentum of stars. Instead it can be determined directly from the data. In this paper it is shown that it is possible to do the same for other stars, even though significantly less data are available. For convenience the appropriate equations from Pijpers (1998) are repeated here. The total angular momentum H of a star is related to its internal rotation rate by:

\begin{displaymath}H\equiv \int\limits_0^1 {\rm d}r \int\limits_{-1}^1 {\rm d}u~ {\cal I}\Omega(r,u)
\end{displaymath} (16)

with the moment of inertia kernel ${\cal I}$:

\begin{displaymath}{\cal I}\equiv 2\pi R^5\rho r^4 \left(1-u^2\right),
\end{displaymath} (17)

where R is the stellar radius and $\rho$ is the density inside the star.

The Subtractive Optimally Localised Averages (SOLA) algorithm (Pijpers & Thompson 1992; Pijpers & Thompson 1994) provides measures of the rotation rate by constructing appropriate linear combinations of the kernels, and the data, from expression (5) or from expression (13). For resolving the rotation rate as a function of rand $\theta$ one would construct a linear combination of kernels that would resemble as closely as possible for instance a two-dimensional Gaussian. For measuring the total angular momentum one instead constructs a linear combination of kernels, that resembles as closely as possible the kernel ${\cal I}$.

Recalling the expression for the associated Legendre polynomial (cf. Gradshteyn & Ryzhik 1994):

\begin{displaymath}P_1^1(u) = -\left(1-u^2\right)^{1/2}
\end{displaymath} (18)

it can be seen that the expressions (14) for s=0 reduce to:

G_1^{l0} &\equiv \displaystyle\frac{3}{4} \...
...equiv -\displaystyle\frac{3}{4} \left(1-u^2\right),
\end{array}\end{displaymath} (19)

which means that they have the same latitudinal (u) dependence as ${\cal I}$. Given that the functions Gls correspond to an orthogonal polynomial expansion one deduces that only the lowest order a-coefficients a1 are required for any l in order to match exactly the u-dependence of ${\cal I}$. The a1 are also the most robust and straightforward to obtain from stellar data. They correspond to the linear least squares slope of $\nu_{nl~
m}-\nu_{nl~ 0}$ as a function of m.

It is convenient at this point to normalise the density by the mean density of the star which defines $\varrho$

\begin{displaymath}\varrho\equiv {4\pi R^3\over 3 M}\rho;
\end{displaymath} (20)

where M is the total stellar mass. It is also useful to rewrite ${\cal I}$, using explicitly normalised radial and latitudinal parts for the integration kernel:

\begin{displaymath}{\cal I}\equiv {2\over 5} M R^2 I_{\rm M} {\cal T}(r) {3\over 4}\left(1-u^2\right),
\end{displaymath} (21)

where $I_{\rm M}$ is defined as

 \begin{displaymath}I_{\rm M}\equiv \int\limits_0^1{\rm d}r~ 5 \varrho r^4,
\end{displaymath} (22)

and ${\cal T}$ defined as

 \begin{displaymath}{\cal T}(r)\equiv \varrho r^4 \left/\int\limits_0^1{\rm d}r~ \varrho
\end{displaymath} (23)

The normalisation of these kernels refers to both  ${\cal T}(r)$ and  ${3\over 4} (1-u^2)$ giving unity when integrated over their respective domains. For a homogeneous gas sphere $I_{\rm M} = 1$, for real stars $I_{\rm M} < 1$. Because the two G functions for s=0 in (19) are identical except for a sign difference, the expression (12) can be reformulated as follows

 \begin{displaymath}2\pi {a_{nl~ 1}\over\beta_{nl}} = \int\limits_0^1{\rm d}
r {\...
...limits_{-1}^1{\rm d}u {3\over 4} \left(1-u^2\right)\Omega(r,u)
\end{displaymath} (24)

in which the ${\cal K}_{nl}$ are new, normalised kernels, that are functions of r only:

\begin{displaymath}{\cal K}_{nl}\equiv {\varrho r^2 \left[\xi_{nl}^2 - {2\over
...i_{nl}\eta_{nl}+\left(1-{1\over L^2}\right)\eta_{nl}^2\right]}
\end{displaymath} (25)

and the factors $\beta_{nl}$ are constants of order unity, defined by:

\begin{displaymath}\beta_{nl}\equiv {\int\limits_0^1{\rm d}r \varrho r^2 \left[\...
...rm d}r \varrho r^2 \left[\xi_{nl}^2 + \eta_{nl}^2\right]}\cdot
\end{displaymath} (26)

The problem of finding the SOLA estimate of the total angular momentum H is reduced to finding linear coefficients cnl, such that the integrated squared difference:

 \begin{displaymath}\int\limits_0^1\!{\rm d}r \left[\sum\limits_{nl}c_{nl}{\cal K...
\end{displaymath} (27)

is minimised. The term $\mu E_{nl~n'l'}$ is added in (27) because there are measurement errors on the a-coefficients, for which  Enl n'l' is the variance-covariance matrix. If the parameter $\mu$ were to be set to 0 (no error weighting) the error on the resulting value for H propagated from the measurement errors on the a-coefficients would be unacceptably large. By adjusting $\mu$ a balance can be found between on the one hand the mismatch to the kernel ${\cal T}$, with the associated uncertainty on H, and on the other hand the estimated statistical error in H propagated from measurement errors. Since the ${\cal K}_{nl}$ and ${\cal T}$ are normalised, one can explicitly add the constraint that $\sum_{nl}c_{nl} =
1$ so that the cnl correspond to averaging weights. In this way an estimate of H can be obtained by taking the corresponding weighted average of the data:

\begin{displaymath}H \approx {4\pi\over 5} M R^2 I_{\rm M} \sum\limits_{nl}c_{nl}
{a_{nl~ 1}\over\beta_{nl}}\cdot
\end{displaymath} (28)

In order to estimate the error for H for stars several factors must be accounted for. For the Sun the errors on M, R, and $I_{\rm M}$ are all negligible since the first two are directly measured to high accuracy and the solar models are sufficiently good that the uncertainty on the internal density, and therefore $I_{\rm M}$, is small. For single stars the mass and radius can be constrained, either from stellar oscillations as well or from other measurements, but the associated uncertainty is no longer negligible. This implies that the internal density distribution and therefore $I_{\rm M}$ is also no longer a negligible source of uncertainty. For stars the measurement errors for the a-coefficients are likely to be larger than for the Sun. Furthermore there is always a mismatch between the best matching averaging kernel $\sum c_{nl}{\cal K}_{nl}$ and the kernel ${\cal T}$, which introduces a source of error as well, which is smaller for the Sun than it is for other stars. The range of l values for which data can be obtained for stars is restricted to $l\leq4$, whereas solar data are available for l up to 250 or even higher with present-day instrumentation. The mismatch between the kernel ${\cal T}$ and the best matching averaging kernel $\sum c_{nl}{\cal K}_{nl}$ increases if the number of available modes decreases, and so does the associated uncertainty in the estimate of H. The magnitude of the error introduced in this way can be estimated as outlined in Appendix A of Pijpers & Thompson (1994).

2.3.1 Special case: Solid body rotation
If the rotation rate $\Omega$ is independent of radius, such as for solid body rotation with a rotation rate $\Omega_{\rm S}$, then Eq. (24) reduces to

\begin{displaymath}{a_{nl~1}\over \beta_{nl}} = {\Omega_{\rm S}\over 2\pi},
\end{displaymath} (29)

because of the normalisation of the kernels. In this case there is no need to attempt to construct a linear combination of kernels that matches ${\cal T}$ and instead one can use weights resulting in the minimal variance for H. This corresponds to the limiting case $\mu\rightarrow \infty$ in (27). If the error estimates for the individual anl 1 are given by  $\sigma_{nl}$, and the errors are uncorrelated, then the linear coefficients should be set to:

\begin{displaymath}c_{nl} = {\beta_{nl}^2 /\sigma_{nl}^2 \over\sum\limits_{nl} \beta_{nl}^2
\end{displaymath} (30)

In this case the contribution to the error estimate for H from the a-coefficients $\sigma_{\rm a}(H)$ is:

\begin{displaymath}\sigma_{\rm a}(H) = \left[\sum\limits_{nl} \beta_{nl}^2/\sigma_{nl}^2\right]^{-1/2}.
\end{displaymath} (31)

It is usually not known a-priori whether a star rotates as a solid body. However if there is no trend or correlation of the anl 1 with $\nu_{nl}$ or l, and the spread of values among the anl 1is consistent with a distribution coming entirely from measurement errors, then this procedure will lead to a value for H that is consistent with the data, and which has a tighter error bound than what would be obtained from minimization of (27) with finite $\mu$. One should note that the averaging kernel $\sum_{nl} c_{nl} {\cal
K}_{nl}$ corresponding to this choice of cnl will usually increase strongly towards the surface. If near-surface layers of the star are in fact rotating at a different rate than the interior, an estimate of H obtained in this way evidently will be biased.

2.4 Tests using solar data

\end{figure} Figure 1: SOLA coefficients cnl as a function of mode turning point. Modes for l=1, 2, 3 and 4 are identified by crosses, asterisks, circles, and squares respectively. Coefficients for the higher values of l are indicated by points.

Application of the method to helioseismic data from the Global Oscillations Network Group (GONG) network and also data from the Solar and Heliospheric Observatory (SoHO) satellite yields a solar total angular momentum of: $H = (190.0 \pm 1.5)\times 10^{39}~{\rm kg~ m^2~ s^{-1}}$(Pijpers 1998). In Fig. 1 are plotted the SOLA coefficients cnl for the latter dataset as a function of the turning point radius $r_{\rm t}$ of the mode in the Sun defined by:

\begin{displaymath}{c_{\rm s}(r_{\rm t})\over r_{\rm t}} = {2\pi\nu_{nl}\over L},
\end{displaymath} (32)

where $c_{\rm s}$ is the internal solar sound speed. Although this set included modes with l up to 250 the figure shows that the weights cnl for l>4 rapidly decrease in value as the turning point radius of the mode increases. One would therefore expect that for a mode set restricted to $l\leq4$ or even $l\leq 2$ it should still be possible to obtain a reasonably accurate estimate of H.

The same solar data set after excluding all modes for which l>4contains 39 modes. The value of H obtained is $H = (199. \pm
2.4)\times 10^{39}~{\rm kg~ m^2~ s^{-1}}$, where the error is just the statistical error propagated from the data. There is an extra contribution, not included in this error estimate, which comes from the mismatch between the averaging kernel ${\cal K} = \sum
c_{nl} {\cal K}_{nl}$ and the target kernel ${\cal T}$ which is of course rather larger here than it is for the full solar set.

If instead all modes for which l>2 are excluded, only 14 modes remain and $H = (191. \pm 5.)\times 10^{39}~{\rm kg~ m^2~ s^{-1}}$. The mismatch between the averaging kernel and the target kernel is substantial, and the close correspondence between this H and the value obtained for the full set is in part through fortuitous cancellation: the averaging kernel oscillates around zero with a large amplitude and this happens to cancel out in the integration.

Evidently it is possible to obtain an accurate estimate of the solar angular momentum even from small sets of data covering a limited range in l. It is realistic to expect to observe in stars with the upcoming satellite projects as many multiplets or more over the same limited range of l, with a comparable precision of the frequencies and a1 coefficients. Thus stellar angular momenta can be measured with a precision comparable to the one for the Sun shown here.

2.5 The algorithm

The minimization of Eq. (27) for the coefficients cnlreduces to solving a set of linear equations. The general case can be found in the papers presenting the method (Pijpers & Thompson 1992, 1994). In the specific case being considered here the linear equations to solve, using vector and matrix notation, are:

\begin{displaymath}{\vec v} = {\cal C}\cdot {\vec c},
\end{displaymath} (33)

where ${\vec c}$ is the vector composed of all the coefficients cnl, which is to be solved for. The matrix ${\cal C}$ is composed of two components; ${\cal Y}$ and E. Replacing for convenience of notation the combined nl indices with a single index i, and similarly n'l' with j, the individual matrix elements are:

{\cal C}_{ij} &= {\cal Y}_{ij} + \mu E_{ij}...
...&= \int\limits_0^1{\rm d}r {\cal K}_i{\cal K}_j \\
\end{array}\end{displaymath} (34)

where E is the error (co)-variance matrix. The elements of the vector ${\vec v}$ are:

\begin{displaymath}{\vec v}_i = \int\limits_0^1{\rm d}r {\cal K}_i {\cal T}.
\end{displaymath} (35)

Adding the constraint that the sum of the coefficients be unity is implemented by augmenting ${\cal C}$ with a row and a column filled with 1's except for the diagonal element which is 0. Further the vector ${\vec v}$ is augmented with a single element 1.

3 Seismic determination of the moment of inertia

There can be substantial variation in the value of $I_{\rm M}$ between stars of different masses, and also for the same star during its life time. Because the variation in $I_{\rm M}$ is so large it is worthwhile to attempt to measure it seismically, in order to reduce the overall uncertainty for the measurement of the angular momentum.

For a small set of stellar models computed using the stellar evolution code of Christensen-Dalsgaard (cf. Christensen-Dalsgaard 1993a) the scaled angular momentum is calculated and shown in Fig. 3 and the appropriate evolutionary tracks in Fig. 2. Every model has been evolved from Zero Age Main Sequence (ZAMS) until a time that is roughly its main sequence life time. The same He and heavy element abundance is used for all models, as well as the same properties for the convection which is parameterised using a standard mixing length recipe.

\end{figure} Figure 2: Evolutionary tracks in a theoretical HR diagram $(T_{\rm eff}, \log L/L_\odot)$ for stellar models of  $0.85~M_\odot, 1.0~M_\odot, 1.5~M_\odot, 2.0~M_\odot, 3.0~M_\odot, 4.0~M_\odot$. The tracks are plotted from ZAMS up to the ages indicated, roughly corresponding to the total main sequence life time.

\end{figure} Figure 3: The scaled moment of inertia $I_{\rm M}$ normalised to the solar value  $I_{\rm M,\odot }$ versus the mean stellar density in solar units, for the stellar models of Fig. 2. For all models the mean density and $I_{\rm M}$ decrease as a function of time while on the main sequence.

Beyond the main sequence, as these stars evolve towards and up the giant branch, their mean density will tend to continue to decrease. However, the moment of inertia at some point increases quite sharply. Stellar models with masses below $\sim$ $1.2~M_{\odot}$have a convective envelope, models with masses higher than that value have little or no convective envelope but instead a convective core that increases in fractional mass with the total mass. In Fig. 2 the evolution for the lower mass stars first is directed towards higher  $T_{\rm eff}$ before turning around whereas the higher mass stars evolve towards lower  $T_{\rm eff}$ throughout this evolutionary stage. For each evolutionary track in Fig. 2 the corresponding track is shown in Fig. 3: the lower mass tracks appear to correspond to a single relation between mean density and $I_{\rm M}$ whereas the higher mass tracks lie roughly parallel to each other. For the purposes of this paper the important point to notice is that the value of $I_{\rm M}$ between stars of different masses, and also for the same star during its life time, changes by as much as 50%.

3.1 The inverse problem

Since the value of $I_{\rm M}$ of an observed star can be substantially different from the value  $I_{\rm M~ {\rm ref}}$ of a model for that star it is worthwhile to quantify that difference using seismic data. It is possible to do this given a reference model and the observed central frequencies  $\nu_{nl 0}$ of multiplets. The starting point is a stellar model that matches within the measurement errors the effective temperature, luminosity, and metallicity of the star; the reference model. The structure of this model is used to compute eigenfrequencies of oscillations  $\nu_{nl~{\rm ref}}$. From this point it is implied that these frequencies are for m=0 and the subscript 0 is omitted. For a model that is sufficiently close in structure to the star, the difference between the observed frequencies $\nu_{nl}$ and the reference model frequencies  $\nu_{nl~{\rm ref}}$ can be related to the structural differences as follows (cf. Dziembowski et al. 1990):
$\displaystyle {\nu_{nl}-\nu_{nl~{\rm ref}}\over\nu_{nl~{\rm ref}}} =
\rho}(r)\right] +Q_{nl}^{-1}{\cal G}(\nu_{nl}) + \epsilon_{nl}.$     (36)

Here $c_{\rm s}^2$ and $\rho$ are the square of the sound speed and the density, and  $\delta_{\rm r}$ refers to the difference in these quantities between the star and the reference model, which are functions of r. The sets of functions $K^{nl}_{c_{\rm s}^2,\rho}$ and $K^{nl}_{\rho,c_{\rm s}^2}$ can be calculated from the reference model (cf. Gough & Thompson 1991). As pointed out before, other combinations of two variables can be chosen. For the purposes of this paper it is convenient to use the density $\rho$ as one of the two, the sound speed is chosen as the other variable merely for definiteness. The term ${\cal G}(\nu_{nl})$ accounts for the near-surface effects on the frequencies: the precise boundary conditions for the oscillations are determined by the pressure and density structure of stellar atmospheres which are not well known. However it can be shown that these effects can be accounted for through appropriate subtraction of a low-order polynomial in frequency ${\cal
G}$ (cf. Pérez Hernández & Christensen-Dalsgaard 1994). The $\epsilon$ are the measurement errors on the data.

In the same spirit of linearisation around a reference model the aim is now to measure the difference  $\delta I_{\rm M}$ between the moment of inertia of the star and the reference model. Using Eq. (22) an expression for the relative difference in moment of inertia is:

\begin{displaymath}{\displaystyle\delta I_{\rm M}\over\displaystyle I_{\rm M}} =...
...l T}(r){\displaystyle\delta\varrho\over\displaystyle
\end{displaymath} (37)

where the kernel ${\cal T}(r)$ is the same as before, given by Eq. (23). Using the SOLA method the task is now to find coefficients fnl such that the linear combination of density kernels matches  ${\cal T}(r)$, while at the same time minimising the effect of the measurement errors, the effect of the surface terms, and also ensuring that the linear combination of the sound speed kernels is as close to 0 as possible everywhere. This is achieved by minimising for the coefficients fnl the quantity:
$\displaystyle \int\limits_0^1 {\rm d}r \Big[\sum f_{nl}K^{nl}_{\rho,c_{\rm s}^2...} K^{nl}_{c_{\rm s}^2,\rho}(r)\Big]^2 + \mu_2 \sum f_{nl}f_{n'l'}E_{nl~ n'l'}$     (38)

with the additional constraint of unimodularity:

 \begin{displaymath}\sum f_{nl}\int\limits_0^1 K^{nl}_{\rho,c_{\rm s}^2}(r) =1.
\end{displaymath} (39)

The $\mu_1, \mu_2$ in Eq. (38) are weights that must be chosen to balance the various contributions to the function to be minimised. This method has been used for inversions for resolving solar structure (cf. Dziembowski et al. 1994; Basu et al. 1996; Rabello-Soares et al. 1999; Di Mauro et al. 2002). It should be noted also that the minimization is usually done under the additional constraint that the total mass for the star and the reference model be identical:

 \begin{displaymath}\int\limits_0^1{\rm d}r \varrho r^2
{\displaystyle\delta_{\rm r}\rho\over\displaystyle \rho} = 0.
\end{displaymath} (40)

It can be seen that Eq. (40) is of the same form as Eq. (36) but instead of a frequency difference, there is a single datum 0, and the kernels are:
$\displaystyle K^{\rm mass}_{\rho,c_{\rm s}^2}(r)$ = $\displaystyle \varrho r^2$  
$\displaystyle K^{\rm mass}_{c_{\rm s}^2,\rho}(r)$ = 0 . (41)

In this case the measurement error for the datum 0 is the relative error on the stellar mass. For the Sun this relative error is $\sim$10-6, for single stars this error is a few percent or more. The surface term ${\cal
G}$ is 0 for this constraint.

There are various ways of carrying out the surface term filtering. One can add constraints in the same way as is done for Eqs. (39) and (40) using for instance a set of orthogonal polynomials in frequency  $\Phi_\lambda$ of low order. Thus one adds constraints of the following form:

\begin{displaymath}\sum f_{nl}Q_{nl}^{-1} \Phi_\lambda(\nu_{nl}) = 0\hskip 1cm \lambda=1,..,\Lambda.
\end{displaymath} (42)

Equivalently one can apply an appropriate linear transform to both the data and the kernels (cf. Basu et al. 1996) in Eq. (36), so that a new linear inverse relation is generated very similar to Eq. (36) except that the surface contribution disappears.

It is known from helioseismology that, given enough data, it is possible to accommodate all of these constraints, even if instead of ${\cal T}$ a Gaussian function with a small width is chosen (cf. Basu et al. 1997). Generally the uncertainty in the localised estimate of the density difference or sound speed difference increases as the width of the Gaussian target kernel is reduced. Since the function  ${\cal T}(r)$ is quite broad, the smaller expected signal-to-noise and number of nl-multiplets obtainable for stars compared to the Sun, should not prevent a seismic determination of  $\delta I_{\rm M}$. However, for each additional surface term to fit, the number of degrees of freedom is reduced. Since there is less data available for other stars it is likely to be more difficult than it is for the Sun to accommodate all constraints including extensive surface term filtering. Because of this it is worth noting that the particular combination of $\rho$ and $c_{\rm s}$ as variables for which to invert may not be optimal. If the equation of state is assumed to be known, then it is possible to reformulate Eq. (36) in terms of differences in $\rho$ and in Helium abundance Y. The kernels for Yhave large amplitudes only near the surface and the influence of Ycan therefore be considered as a surface effect to be filtered out. In that case the term analogous to the term with $\mu_1$ can be removed from Eq. (38). Since the number of constraints on the inversion is reduced by removing this term, one will obtain a better match to the kernel ${\cal T}(r)$ and/or a smaller propagated statistical error. Detailed testing is deferred at this stage.

Once the coefficients fnl are obtained, by solving a set of linear equations similar to those in Sect. 2.5, the angular momentum is then estimated as follows:

\begin{displaymath}H = {\displaystyle 4\pi\over\displaystyle 5} M_{\rm ref}R_{\r...
{\displaystyle a_{nl~ 1}\over\displaystyle\beta_{nl}}\right),
\end{displaymath} (43)

where the subscript ${\rm ref}$ refers to the reference model. The issue of choosing the most appropriate reference model is closely linked with the issue of mode identification which is discussed in the next section.

3.2 The mode identification problem

In order to be able to apply the inverse techniques discussed in the previous sections, it is necessary to know with which n and l to label any given measured frequency or a1 coefficient. Otherwise one does not know which coefficient cnl to combine with that datum. This is generally referred to as the mode identification problem. In the following sections techniques are proposed to determine l and n from observations. In principle such methods must be subjected to the same type of numerical tests that are presented in Sect. 2.4. Extensive testing is in progress in the more general context of asteroseismic inversions for structure, where the same problems occur, but this falls beyond the scope of the present paper. The results of such tests will be reported in a subsequent paper.

3.2.1 Determination of l
The first problem that arises in attempts to apply helioseismic techniques to other stars is due to the inability to resolve the stellar surface when obtaining asteroseismic time series. For the Sun appropriate spatial filters can be designed to isolate the signal for any specified l roughly up to the number of resolution elements over the solar disk. For an unresolved stellar disk the number of resolution elements is 1 and therefore another method needs to be employed to provide the required filter. Otherwise the correspondence between frequencies and their appropriate labels of n, l and m remains uncertain, which precludes effective use of inverse methods.

One aid in assigning the appropriate l value to a multiplet for a rotating star is the structure of the multiplet itself. Since mcannot be higher than l, the number of peaks in the multiplet determines the minimum possible value of l. The oscillatory signal corresponding to some m values within a multiplet can be suppressed, or too weak to measure given the noise, due to for instance projection effects because of the angle of inclination of the rotation axis with the line of sight. This means that for a given oscillation spectrum multiplet the value of l is higher than indicated by the number of peaks present in the multiplet.

A second aid in assigning l is the use of colour information in the oscillations. The number of nodelines over the surface of a star increases with l. Because of averaging over regions which have alternating excess or deficit local brightness (or velocity) the signal that remains for unresolved observations of the stellar disk reduces quickly with increasing l (cf. Christensen-Dalsgaard & Gough 1982). In this averaging procedure the limb darkening of the stellar surface has to be taken into account, which is wavelength dependent. By comparing the amplitude of a signal in (at least) two different wavelength bands it is possible to identify l. The limb darkening effectively functions as a (non-optimal) spatial filter for l (cf. Heynderickx et al. 1994; Bedding et al. 1996). The factor by which the amplitude of intensity variations has to be multiplied in order to obtain the signal for unresolved observations in a spectral band X of an oscillation with degree l and for m=0 is:

 \begin{displaymath}S_X(l) = 2\sqrt{2l+1}\int{\rm d}\lambda
\int\limits_0^1{\rm d}u P_l(u) I_X(\lambda, u) u,
\end{displaymath} (44)

where Pl is the Legendre polynomial of degree l, and $I_X(\lambda,
u)$ is the centre-to-limb variation of intensity for band X: the limb darkening profile. For convenience of notation the polar coordinate axis for the rotation has been chosen to be along the line of sight in which case the only non-zero contribution is for m=0. This orientation would preclude measurement of the rotational splittings. However, for other orientations of the axis the response for the different m can be obtained by multiplying the SX(l) with an appropriate factor  $\Gamma_{ml}$ (cf. Christensen-Dalsgaard & Gough 1982). An equivalent analysis then applies after summing over m.

The spectral response of the filter for the band X is assumed to be accounted for in  $I_X(\lambda,
u)$. For observations in velocity or equivalent width of spectral lines, similar relations hold (cf. Bedding et al. 1996). If one has time series in several bands, or in a combination of intensity, velocity, and equivalent width, construction of a filter to select a given degree l0 can be considered as an inverse problem (cf. Christensen-Dalsgaard 1984; Toutain & Kosovichev 2000). Ideally one would combine the observed relative amplitudes $S_X = \sum_l S_X(l)$ with weights gX(l0) so that:

 \begin{displaymath}\sum\limits_X g_X(l_0) \left(\sum\limits_{l}S_X(l)\right) =
N_{l_0} \delta_{l~ l_0},
\end{displaymath} (45)

where Nl0 is a normalisation factor, and $\delta_{l~ l_0}$ is the Kronecker delta symbol. In practice only a finite amount of data is available, with a finite signal-to-noise ratio, so this ideal cannot be attained. Solving this problem using SOLA could be achieved by minimising for the gX(l0) the following:
                                         $\displaystyle \sum\limits_X \Big[g_X(l_0)\sum\limits_l 2\sqrt{2l+1}\int\limits_0^1{\rm d}u
\int\limits{\rm d}\lambda~ P_l(u) u I_X(\lambda, u)$  
    $\displaystyle - {\cal T}_l\left(\vert l-l_0\vert\right)\Big]^2 +
\mu\sum\limits_{X~ X'} g_X(l_0)g_{X'}(l_0)E_{X~ X'},$ (46)

where $\mu$ is a weighting term analogous to those used in Eq. (38), and EX  X' is the (co-)variance matrix of the errors in the observables. In implementations of the SOLA method for other inverse problems it is common to use a Gaussian form for the filter  ${\cal T}_l$ around l0 in which the width is a free parameter to be set as small as feasible given the measurement errors and the amount of data available. In the case at hand a Gaussian might be inappropriate. To illustrate this consider a set of weights gX(l0) such that:

 \begin{displaymath}\sum\limits_X g_X(l_0) \int\limits{\rm d}\lambda~ u I_X(\lambda, u) =
\end{displaymath} (47)

If it were possible to do this and if the entire stellar surface were to have been observed, the orthogonality of Legendre polynomials would ensure that Eq. (45) would be satisfied exactly. This can be seen by multiplying in Eq. (44) left- and right-hand sides with the gX(l0), summing over X, and integrating over u from -1 to 1. For real observations covering only half the stellar surface the integration is from 0 to 1, as written in Eq. (44), and Eq. (45) is replaced by (cf. Gradshteyn & Ryzhik 1994):

\sum\limits_X &g&_X(l_0) \left(\sum\limi...
... = {\rm even}(l,l_0),~ l_{\rm o} ={\rm odd}(l,l_0)
\end{array}\end{displaymath} (48)

where the function ${\rm even}(l,l_0)$ selects the even element of the pair l, l0 and ${\rm odd}(l,l_0)$ selects the odd element. Although it is possible that Eq. (47) does not correspond to the best attainable filter, it is not unlikely that it comes close. In particular it has the desirable feature of producing 0 cross-talk between modes differing in l by an even number, because the linear combination of Eq. (48) is 0 for $\vert l - l_0\vert = {\rm even}$. From Eq. (1) it can be seen that modes differing in l by an odd number are separated in frequency by an odd multiple of $\sim \Delta\nu/2$, and therefore well separated. Modes (n1, l1) and (n2, l2) differing in l by an even number can nearly coincide in frequency if (n1 - n2)=(l2-l1)/2. It is therefore exactly modes with $\vert l - l_0\vert = {\rm even}$ for which one wishes to achieve a minimal cross-talk with the mode identification filter. Therefore minimization of Eq. (46) can be replaced with minimization for the gX(l0) of:
$\displaystyle \int\limits_0^1{\rm d}u~ \left[\sum\limits_X g_X(l_0)
... u) - P_{l_0}(u)\right]^2
+\mu\sum\limits_{X~ X'} g_X(l_0)g_{X'}(l_0)E_{X~ X'}.$     (49)

In practice this minimization procedure usually will not achieve the exact equality of Eq. (47) and thus not achieve an exactly 0 cross-talk between modes for which $\vert l - l_0\vert = {\rm even}$. Better results can be achieved for increasing numbers of bands (or other independent observables) X, and for increasing signal-to-noise ratio in each band. The advantage compared to minimization of Eq. (46) is that an infinite summation over l, which would have to be truncated in practice, is avoided altogether. In either case it is essential to have accurate knowledge of limb darkening profiles as a function of wavelength for target stars of asteroseismic campaigns.

3.2.2 Determination of n
Having accounted for l to the extent that this is possible given the available data, what remains is to assign the appropriate radial node number n to any given frequency. Often recourse is taken to Eq. (1). If a regular pattern of oscillation frequencies is identified, a large separation can be determined. The large separation together with a (tentative) identification of l can provide n, but some assumption for the surface term $\alpha$ is required. A way to partially circumvent this problem is to use frequency separations as data instead of the frequencies. Define the frequency separation $\Delta_n^{(l)}\nu$ as:

\begin{displaymath}\Delta_n^{(l)}\nu\equiv \nu_{n+1~ l}-\nu_{n~ l}.
\end{displaymath} (50)

An integral relation equivalent to Eq. (36) can be obtained by subtracting Eq. (36) for n from the same equation for n+1. This yields:

 \begin{displaymath}{\Delta_n^{(l)}\nu - \Delta_n^{(l)~ {\rm ref}}\nu \over
...l}^{-1}{\cal G}^{(\Delta )}(\nu_{(\Delta )l}) + \epsilon_{nl},
\end{displaymath} (51)

where the appropriate kernels are defined as:

K^{(\Delta)l}_{c_{\rm s}^2,\rho} &= \left[\...
\Delta_n^{(l)~{\rm ref}}\nu\right..\\
\end{array}\end{displaymath} (52)

By doing this one achieves elimination of the lowest order surface term compared to Eq. (36). The other surface terms ${\cal
G}^{(\Delta )}$ are still functions of frequency only. If necessary, further filtering can be carried out using the same procedure as outlined in Sect. 3.1. However, it is expected, based on experience in helioseismology, that the remaining surface terms are negligible up to some frequency  $\nu_{\rm max}$. At low frequencies the asymptotic relation Eq. (1) loses accuracy, but there is an intermediate frequency range where the $\Delta_n^{(l)~{\rm ref}}\nu$ are constant and equal to the large separation $\Delta\nu$ from asymptotic theory. $\Delta\nu$ is determined by the structure of the star through the crossing time for sound waves $\tau_0$:

\begin{displaymath}\Delta\nu^{-1} = \tau_0 = 2R\int\limits_0^1{\rm d}r {1\over c_{\rm s}(r)}\cdot
\end{displaymath} (53)

If in the intermediate frequency range the individual relative differences of separations between model and star are not consistent with 0, this implies that the mean sound crossing time of the star is different from the reference model. In this case it would be beneficial to use a different reference model for which the sound crossing time is identical to that of the star. As soon as the appropriate value of n for one frequency or frequency separation is established, all other modes are identified as well by virtue of either the asymptotic relation (1) or more detailed calculations of eigenfrequencies. If Eq. (1) were satisfied exactly for all n given any l, ambiguity would remain in the determination of n, since only the sum $n+\alpha$ could be determined. It is exactly because at low frequencies there are departures from the asymptotic relation, that it becomes possible to assign a most likely value of n for each frequency, given a reference model. A useful strategy is to minimise the quantity  $\chi_{n'\Delta}$:

 \begin{displaymath}\chi_{n'\Delta} \equiv\sum\limits_{nl} {1\over\sigma_{n+n'~
...}\nu -
\Delta_n^{(l)~ {\rm ref}}\nu -\Delta\Delta\nu\right]^2
\end{displaymath} (54)

for the two variables n' and $\Delta\Delta\nu$. Here n' relates the index of the measured frequencies to the n provided by oscillation models, which is the mode identification sought. $\Delta\Delta\nu$ is the difference in the asymptotic large separation between star and model. As stated above this difference should ideally be 0 and thus a value significantly different from 0 indicates that it is desirable to use a different reference model, and also what sound crossing time that reference model should have. The  $\sigma_{n+n'~ l}$ are the measurement errors on the data. The summation should be done for modes with oscillation frequencies $\nu<\nu_{\rm max}$ in order to avoid biases due to remaining surface term effects.

At this point it should be noted that there may be a range of models that have the same sound crossing time as the star, and thus $\Delta\Delta\nu =0$ in the minimum of Eq. (54), but for which n' is not necessarily the same. The minimum value of  $\chi_{n'\Delta}$ for all of these different models should then indicate the most appropriate reference model and the associated mode identification. The implication is that mode identification performed in this way is model dependent. In practice one would perform inversions using reference models spanning the allowed range, with n' also being allowed to vary within an allowed range. Linear inversion methods, of which SOLA is one, are fast. Therefore no significant computational effort is required in estimating the influence of these uncertainties on the uncertainty in the moment of inertia.

3.3 Constructing the reference model

As stated above, the mass and radius of the Sun are known to high precision. This implies that it is possible to construct a reference model for the Sun, using standard physics and a stellar evolution code, that is very close in internal structure to that determined from seismology. One can therefore be confident that the linearisation procedure is valid, and tests have confirmed this. For other stars the mass and radius are not nearly as well determined, even for those with the best spectroscopic determinations of their classical parameters: effective temperature, luminosity, heavy element abundance (cf. Pijpers 2003). This may mean that the linearisation procedure, on which Eq. (36) is based, is no longer valid for all reference models that are consistent with these classical parameters. In this case inversions are not guaranteed to provide accurate estimates for the true differences between star and model.

It may be necessary to repeat the same inversion using several different reference models, all matching within the uncertainties the classical parameters, and also the sound crossing time of the star as discussed in Sect. 3.2.2. That model for which the differences between observed and measured frequencies are least is the model most likely to be an adequate reference model. By doing this one can account for effects from propagation of measurement errors on the oscillation frequencies, as well as minimising systematic effects due to linearisation around a reference model that may be too dissimilar from the star.

4 Conclusions

In this paper an end-to-end method is outlined for the asteroseismic determination of stellar angular momentum for stars with low or moderate rotation rates: asteroseismogyrometry. It is shown how to determine the appropriate weighted average of the rotation rate from observations of rotationally split oscillation frequency multiplets, and how to determine the moment of inertia from the central frequencies. In both cases an inverse problem needs to be solved which is implemented with the SOLA framework (cf. Pijpers & Thompson 1994).

A method is outlined as well for identifying oscillation modes, i.e. assigning appropriate values of n and l to measured frequencies. In part this is achieved through solving an inverse problem as well. It is shown that an essential tool in this process are accurate wavelength dependent limb-darkening profiles for the target stars of asteroseismic campaigns. In order to obtain these it is necessary to construct appropriate 3D stellar atmosphere models and/or measure stellar limb darkening using interferometric techniques. In order to prepare for data to be gathered with upcoming satellite missions it is important to commence such work now.


The author thanks Jørgen Christensen-Dalsgaard and Teresa Teixeira for useful comments and discussions when writing the paper. The author also thanks the Theoretical Astrophysics Center, a collaborative centre between Copenhagen University and Aarhus University funded by the Danish Research Foundation, for support of this work. GONG is managed by NSO, a division of NOAO that is operated by the AURA under co-operative agreement with NSF. The GONG data were acquired by instruments operated by the BBSO, HAO, Learmonth, Udaipur, IAC and CTIO. The MDI project operating the SOI/MDI experiment on board the SoHO spacecraft is supported by NASA contract NAG5-3077 at Stanford University.


Copyright ESO 2003