A&A 455, 607-620 (2006)
DOI: 10.1051/0004-6361:20065015

Acoustic oscillations of rapidly rotating polytropic stars

I. Effects of the centrifugal distortion

F. Lignières - M. Rieutord - D. Reese

Laboratoire d'Astrophysique de Toulouse et Tarbes, UMR 5572, Observatoire Midi-Pyrénées, Université Paul Sabatier - Toulouse 3, 14 avenue É. Belin, 31400 Toulouse, France

Received 13 February 2006 / Accepted 31 March 2006

Abstract
Aims. A new non-perturbative method to compute accurate oscillation modes in rapidly rotating stars is presented.
Methods. The effect of the centrifugal force is fully taken into account while the Coriolis force is neglected. This assumption is valid when the time scale of the oscillation is much shorter than the inverse of the rotation rate and is expected to be suitable for high radial order p-modes of $\delta$ Scuti stars. Axisymmetric p-modes have been computed in uniformly rotating polytropic models of stars.
Results. In the frequency and rotation range considered, we found that as rotation increases (i) the asymptotic structure of the non-rotating frequency spectrum is first destroyed then replaced by a new form of organization (ii) the mode amplitude tends to concentrate near the equator (iii) differences to perturbative methods become significant as soon as the rotation rate exceeds about fifteen percent of the Keplerian limit. The implications for the seismology of rapidly rotating stars are discussed.

Key words: stars: oscillations - stars: rotation

1 Introduction

Since helioseismology revolutionized our knowledge of the solar interior, great advances in stellar structure and evolution theory are expected from asteroseismology. Major efforts including space missions are under way to detect pulsation frequencies with unprecedented accuracy across the HR diagram (Catala et al. 1995; Walker et al. 2003). To draw information from the observed frequencies, seismology relies on the theoretical computation of eigenmodes for a given model of a star. Yet, except for slowly rotating stars, the effect of rotation on the gravito-acoustic modes is not fully taken into account in the present theoretical calculations (e.g. Rieutord 2001).

Rotational effects have been mostly studied through perturbative methods. In this framework, both  $\Omega/\omega$, the ratio of the rotation rate $\Omega$ to the mode frequency $\omega$, and  $\Omega /\! \sqrt{GM/R^3}$, the square root of the ratio of the centrifugal force to the gravity at equator are assumed to be small and of the same order. Solutions valid up to the first, second, and even third order in  $\Omega/\omega$have been obtained by Ledoux (1951), Saio (1981) and Soufi et al. (1998). The first order analysis proved fully adequate to match the observed acoustic frequency of the slowly rotating sun (Dziembowski & Goode 1992). At the other extreme, the perturbative methods are not expected to be correct for stars approaching the Keplerian limit $\Omega_{\rm K} = \sqrt{GM/R_{\rm e}^3}$, where $R_{\rm e}$ is the equatorial radius. Achernar is a spectacular example of such star since interferometric observations showed a very important distortion of its surface, the equatorial radius $R_{\rm e}$ being at least one and a half times larger than the polar radius $R_{\rm p}$ (Domiciano de Souza et al. 2003). In the context of Roche models, such a flattening occurs at the Keplerian limit  $\Omega_{\rm K}$. For intermediate rotation rates, second or third order perturbative methods might be used, but the main problem is that the limit of validity of the perturbative methods is unknown. Departures from the perturbative results would impact the values of individual frequency but also other properties that are commonly used to analyze the spectrum of observed frequencies. This concerns in particular the rotational splitting, the asymptotic large and small frequency separations or the mode visibility.

New methods able to compute accurate eigenmodes in rotating stars are therefore needed to allow progress in the seismology of rapidly rotating stars. Such methods would also assess the limit of validity of perturbative analysis. The main difficulty comes from the fact that, except in the special cases of spherically symmetric media and uniform density ellipsoids, the eigenvalue problem of gravito-acoustic resonances in arbitrary axially symmetric cavities is not separable in the radial and meridional variables. For self-gravitating and rotating stars, a two-dimensional eigenvalue problem has to be solved.

Clement (1998,1981) made the first attempts to solve this eigenvalue problem for gravito-acoustic modes, investigating various numerical schemes. However, the accuracy of his calculations is generally difficult to estimate. Moreover, the different numerical schemes could not converge for low frequency g-modes when  $\Omega/\omega$ exceeds about 0.5. Since then, eigenmodes in this frequency range have been successfully calculated by Dintrans & Rieutord (2000) using spectral methods. These authors however did not consider the effect of the centrifugal acceleration in their model. The search for unstable modes in neutron stars also triggered the development of numerical schemes able to solve the two-dimensional eigenvalue problem. But only surface gravity modes (f-modes) and some inertial modes (r-modes) have been determined in this context (Yoshida et al. 2005). Espinosa et al. (2004) reported calculations of adiabatic acoustic modes in MacLaurin spheroids of uniform density neglecting the Coriolis force, the potential perturbation and the Brunt-Väisälä frequency.

In this paper, we present a new method to compute accurate eigenmodes in rotating stars. For the first application of the method, we only consider the effect of the centrifugal force through its impact on the equilibrium state of the star; we thus neglect the Coriolis force. This assumption is valid when the time scale of the oscillation is much shorter than the inverse of the rotation rate and is expected to be suitable for high radial order p-modes of $\delta$ Scuti stars. The problem is further simplified by using uniformly rotating polytropes as equilibrium models and assuming adiabatic perturbations as well as the Cowling approximation. Low frequency axisymmetric p-modes have been computed for rotation rates varying from $\Omega =0$ up to $\Omega / \Omega _{\rm K} = 0.59$, this ratio corresponding to a typical $\delta$ Scuti star ( $M=1.8~M_{\odot}, R=2~R_{\odot}$) with an equatorial velocity of 240  ${\rm km~s}^{-1}$. The centrifugal force modifies the effective gravity in two ways: it makes it smaller and aspherical. Decreasing the effective gravity should affect sound waves by reducing the sound speed inside the star and by increasing the star's volume, thus potentially the volume of the resonant cavity. The physical consequences of the non-spherical geometry are unknown.

In the following, the formalism and the numerical method are presented along with accuracy tests. Then, the parameter range of the calculations is given together with the method used to label the eigenmodes. The structure of the frequency spectrum, some properties of the eigenfunctions and the differences with perturbative methods are further analyzed as a function of the rotation rate. These results are discussed in the last section.

2 Formalism

Accurate numerical solutions of 2D eigenvalue problems require a careful choice of the numerical method and the mathematical formalism. In this section we explain the choices that have been made for the variables, the coordinate system, the numerical discretization, and the method to solve the resulting algebraic eigenvalue problem. All play a role in the accuracy of the eigenfrequency determinations that will be presented at the end of this section.

2.1 Equilibrium model

The equilibrium model is a self-gravitating uniformly rotating polytrope. It is therefore governed by a polytropic relation, the hydrostatic equilibrium in a rotating frame, and Poisson's equation for the gravitational potential:

 

$$P_0 = K \rho_0^{1+1/N}$$ (1)

$$0 = - \vec{\nabla} P_0 - \rho_0 \vec{\nabla} \left( \psi_0 -\Omega^2 s^2/2 \right)$$ (2)

$$\Delta \psi_0 = 4\pi G\rho_0$$ (3)

where P₀ is the pressure, $\rho_0$ the density, K the polytropic constant, N the polytropic index, $\psi_0$ the gravitational potential, s the distance to the rotation axis and G the gravitational constant.

The polytropic relation and uniform rotation ensure that the fluid is barotropic. A pseudo-enthalpy can then be introduced $h_0=\int {\rm d}P_0/\rho_0 = (1+N)P_0/\rho_0$ and the integration of the hydrostatic equation reads:

$$% h_0 = h_{\rm c} -(\psi_0 - \psi_{\rm c}) + \frac{1}{2}\Omega^2 s^2$$ (4)

where the subscript "c'' denotes the value in the center of the polytropic model. Equation (4) is then inserted into Poisson's equation to yield:

$$% \Delta \psi_{\rm o} = 4 \pi G \rho_{\rm c} \left( 1 - \fra... ...h_{\rm c}} + \frac{\Omega^2 s^2}{2 h_{\rm c}} \right)^{N}\cdot$$ (5)

Equation (5) is solved numerically, using an iterative scheme. Since the shape of the star is not spherical, a system of coordinates  $(\zeta,\theta,\phi)$ based on Bonazzola et al. (1998) is used, such that $\zeta=1$corresponds to the surface of the spheroid (more details on the coordinate system are given in Sect. 2.3). This enables the use of spectral methods both for the radial coordinate $\zeta$ and the angular ones. The numerical method is further detailed in Rieutord et al. (2005).

2.2 Perturbation equations and boundary conditions

Neglecting the Coriolis force, the linear equations governing the evolution of small amplitude adiabatic perturbations read:

$$% {\partial}_t \rho + \vec{\nabla} \cdot (\rho_{0} \vec{v})=0,$$ (6)

$$% \rho_{0} {\partial}_t \vec{v}= - \vec{\nabla} P + \rho \vec{g}_0 - \rho_{0} \vec{\nabla} \psi,$$ (7)

$$% {\partial}_t P + \vec{v}\cdot \vec{\nabla} P_0 = c_0^{2} \;... ...{\partial}_t \rho + \vec{v}\cdot \vec{\nabla} \rho_0 \right),$$ (8)

$$% \Delta \psi = 4 \pi G \rho$$ (9)

where $\vec{v}$, $\rho$, P and $\psi$, are the perturbations of velocity, density, pressure and gravitational potential. The sound speed is $c_0 = \sqrt{{\Gamma}_{1,0} P_0 / \rho_{0}}$, ${\Gamma}_{1,0}$ being the first adiabatic exponent of the gas, and the effective gravity $\vec{g}_0 = - \vec{\nabla} \left( \psi_0 -\Omega^2 s^2/2 \right)$ has been introduced. In the framework of polytropic models of stars, the polytropic relation (1) is assumed to give a reasonably good approximation of the relation between the pressure and the density of the equilibrium state. As the first adiabatic exponent  ${\Gamma}_{1,0}$ is obtained from the equation of state of the gas, ${\Gamma}_{1,0}$ is in general not equal to 1+1/N.

We simplified Eqs. (6)-(9) following two constraints: first, the governing equations should be written for general coordinate systems because we shall use a surface-fitting non-orthogonal coordinate system. Second, they should take the form $\cal{M} \vec{X} = \lambda \cal{Q} \vec{X}$ where $\lambda$ is the eigenvalue, $\vec{X}$ is the eigenfunction, and $\cal{M}$ and $\cal{Q}$ are linear differential operators. Indeed, the method that we shall use to solve the algebraic eigenvalue problem obtained after discretization works for problem reading $[M] X = \lambda [Q] X$, where X is the discretized eigenvector and, [M] and [Q] are matrices. Taking the time derivative of Eqs. (7) and (9) and using Eqs. (6) and (8) to eliminate the pressure and density perturbations, we obtain two equations for the velocity and the gravitational potential perturbation:

$$% {\partial}^{2}_{tt} \vec{v}= \vec{\nabla} \;\left(c_0^2 \ch... ...v}\cdot \vec{g}_0 - {\partial}_t \psi \right) - \chi \vec{A}_0$$ (10)

$$% \Delta {\partial}_t \psi = - 4 \pi G \;\left(\vec{v}\cdot \vec{\nabla} {\rho}_0 + {\rho}_0 \chi \right)$$ (11)

where $\chi = \vec{\nabla} \cdot \vec{v}$ is the divergence of the velocity. The vector $\vec{A}_0$ characterizes the stratification of the equilibrium model:

$$% \vec{A}_0 = c^2_0 \;\left(\frac{1}{{\Gamma}_{1,0}}\frac{ \v... ... \right)= \frac{c_0^2 N_0^2}{\Vert\vec{g}_0\Vert} \vec{n}_0,$$ (12)

where N₀ is the Brunt-Väisälä frequency and $\vec{n}_0$ is the unit vector in the direction opposite to the effective gravity defined by:

$$% \vec{g}_0 = - \Vert\vec{g}_0\Vert \vec{n}_0.$$ (13)

Note that the barotropicity of the fluid has been used to obtain (10).

In addition to the gravitational potential perturbation, the right hand sides of Eqs. (10) and (11) only involve the divergence of the velocity and the scalar product of the velocity with vectors parallel to gravity. Then, the scalar product of Eq. (10) with gravity,

$$% {\partial}^2_{tt} \vec{v}\cdot \vec{g}_0 = \vec{g}_0 \cdot ... ..._0 - {\partial}_t \psi \right)- \chi \vec{g}_0 \cdot \vec{A}_0$$ (14)

and the divergence of Eq. (10),

$$% {\partial}^2_{tt} \chi = \Delta (c_0^2 \chi + \vec{v}\cdot ... ...{g}_0 - \partial_t \psi) - \vec{\nabla} \cdot (\chi \vec{A}_0)$$ (15)

yield, together with Eq. (11), a complete set of differential equations for the variables $\vec{v}\cdot \vec{g}_0$, $\chi$and $\psi$. Pekeris (1938) who studied the oscillations of spherically symmetric polytropes considered similar variables but, instead of Eq. (10), used a combination of Eqs. (14) and (10) which does not involve second order derivative with respect to the radial coordinate. For general system of coordinate as well, the order of the differential system can be lowered replacing Eq. (10) by the following one:

 

$$% {\partial}^2_{tt} \left[\chi - g^{11} {\partial}_1 \!\! \;\left... ...{\partial}_1 \!\! \;\left(\frac{\chi \vec{g}_0 \cdot \vec{A}_0 }{g_0^1} \right)$$ (16)

where the linear operator ${\cal L}$, defined by,

$$% {\cal L}(F) = \Delta F - g^{11} {\partial}_1 \!\! \;\left(\frac{\vec{g}_0 \cdot \vec{\nabla} F}{g_0^1} \right)$$ (17)

does not contain second order derivatives with respect to the first coordinate x¹. In this expression, g₀¹ is the first contravariant component of the gravity in the natural basis $(\vec{E}_1, \vec{E}_2, \vec{E}_3)$ defined by $\vec{E}_i = \partial \vec{OM} / \partial x^i$, and g¹¹ is the first contravariant component of the metric tensor.

The equations are non-dimensionalized using the equatorial radius, $R_{\rm e}$, as length unit, the density at the center of the polytrope, $\rho_{\rm c}$, as density unit and $T_0 = \;\left(4 \pi G \rho_{\rm c} \right)^{-1/2}$ as time unit. As we look for harmonic solutions in time, the variable are written $F = \hat{F} \exp~(i \omega t)$. Dropping the hat and denoting dimensionless quantities as previous dimensional ones, the governing equations read:

   

$$\lambda W = \vec{g}_0 \cdot \vec{\nabla} \;\left(c_0^2 \chi + W + \Psi \right)- c_0^2 N_0^2 \chi$$ (18)

$$\lambda \left[\chi - g^{11} {\partial}_1 \!\! \;\left(\frac{W}{g_... ...A}_0) + g^{11} {\partial}_1 \!\! \;\left(\frac{c_0^2 N_0^2 \chi}{g_0^1} \right)$$ (19)

$$0 = \Delta \Psi - d_0 W - {\rho}_0 \chi$$ (20)

where $\lambda = - \omega^2$, $W = \vec{v}\cdot \vec{g}_0$, $\Psi = - i \omega \psi$ and d₀ denotes

$$% d_0 = \frac{\Vert \vec{\nabla} \rho_0\Vert}{\Vert\vec{g}_0\Vert}\cdot$$ (21)

Another form of these equations may be obtained replacing W by a new variable $U = c_0^2 \chi + W + \Psi$. The set of equations then reads:

$$% \vec{g}_0 \cdot \vec{\nabla} U - c_0^2 N_0^2 \chi = \lambda (U - \Psi - c_0^2 \chi)$$ (22)

 

$${\cal L}(U) - \vec{\nabla} \cdot (\chi \vec{A}_0) + g^{11} {\part... ...chi + g^{11} {\partial}_1 \!\! \;\left(\frac{c_0^2 \chi}{g_0^1} \right) \right]$$ (23)

$$% - d_0 U + (d_0 c^2_0 - {\rho}_0) \chi + d_0 \Psi + \Delta \Psi = 0.$$ (24)

As in Pekeris (1938), the boundary conditions are that the gravitational potential vanishes at infinity and that U and $\chi$ be regular everywhere.

2.3 Coordinates choice

The choice of the coordinate system has been guided by two considerations. First, for the accuracy of the numerical method, it seems preferable to apply the boundary conditions on a surface of coordinate. This imposes that the stellar surface is described by an equation $\zeta = {\rm const.}$, where $\zeta$ is one of the coordinates. Second, when using spherical harmonic expansions, the regularity conditions at the center have a simple form for spherical coordinates only. Therefore, the coordinate system should become spherical near the center. If  $(r, \theta, \phi)$ denotes the usual spherical coordinates and $r= S(\theta)$ describes the surface, families of coordinates  $(\zeta,\theta',\phi')$ verifying both conditions have been proposed by Bonazzola et al. (1998):

$$% \left\{\begin{array}{l} r = r(\zeta, \theta') \\ \theta = \theta' \\ \phi = \phi', \end{array}\right.$$ (25)

where

$$% r(\zeta, \theta) = \alpha \zeta + A(\zeta) \left[S(\theta) - \alpha \right].$$ (26)

The polynomial $A(\zeta) = (5 \zeta^3 - 3 \zeta^5)/2$ ensures that $r \propto \zeta$ near the center and that $\zeta=1$ describes the surface $r= S(\theta)$. The scalar $\alpha$ is chosen so that the transformation $(r, \theta, \phi) \mapsto (\zeta, \theta, \phi)$ is not singular and the numerical convergence is fast. We find that $\alpha = 1 - \epsilon$ is a convenient choice, $\epsilon = 1 - R_{\rm p}/R_{\rm e}$ being the flatness of the star surface. In the following, we shall refer to $\zeta$ as the pseudo-radial coordinate.

To express the governing equations in this non-orthogonal coordinate system, we use the covariant and contravariant components of the corresponding metric tensor. The non-vanishing components read:

$$% \begin{array}{ll} g_{11} = r_{\zeta}^{2} & g_{12} = g_{21} ... ...{22} = 1/r^{2} & g^{33} = 1/(r^{2} \sin^2\!\theta), \end{array}$$ (27)

and the square-root of the absolute value of the metric tensor determinant is:

$$% \sqrt{\mid g \mid} = r^{2} r_{\zeta} \sin~ \theta.$$ (28)

In Appendix A, the linear operators involved in Eqs. (22)-(24) are expressed in terms of the partial derivatives of  $r(\zeta,\theta)$. Note that, for vectorial operators, we used the natural basis defined above.

2.4 The numerical method

The method follows and generalizes the one presented in Rieutord & Valdettaro (1997). The numerical discretization is done with spectral methods, spherical harmonics for the angular coordinates $\theta$ and $\phi$and Chebyshev polynomials for the pseudo-radial coordinate $\zeta$. The variables U, $\Psi$ and $\chi$ are expanded into spherical harmonics:

$$% U(\zeta,\theta,\phi) = \sum_{\ell=0}^{L}\sum_{m=-\ell}^{+\ell} U^\ell_m (\zeta) Y^m_\ell (\theta,\phi),$$ (29)

and equivalent expressions for $\Psi$ and $\chi$, where $\ell$ and m are respectively the degree and the azimuthal number of the spherical harmonic  $Y^m_\ell (\theta,\phi)$. Then, the governing equations are projected onto spherical harmonics to obtain a system of ordinary differential equations (ODE) of the variable $\zeta$ for the coefficients of the spherical harmonic expansion $U^{\ell}_m(\zeta), \chi^{\ell}_m(\zeta), \Psi^{\ell}_m(\zeta)$. This system is then discretized on the collocation points of a Gauss-Lobatto grid associated with Chebyshev polynomials. It results in an algebraic eigenvalue problem $[M] X = \lambda [Q] X$, where X is the eigenvector of  $L \times N_r$ components and [M] and [Q] are square matrices of $L \times N_r$ lines, L and N_r being respectively the truncation orders on the spherical harmonics and Chebyshev basis. The algebraic eigenvalue problem is solved using either a QZ algorithm or an Arnoldi-Chebyshev algorithm. The QZ algorithm provides the whole spectrum of eigenvalues while the iterative calculation based on the Arnoldi-Chebyshev algorithm computes a few eigenvalues around a given value of the frequency.

Because of the symmetries of the equilibrium model with respect to the rotation axis and the equator, one obtains a separated eigenvalue problem for each absolute value of the azimuthal number  $\mid \! m \! \mid$and each parity with respect to the equator. Thus, for a given $m \geq 0$, we have two independent sets of ODE coupling the coefficients of the spherical harmonic expansion having respectively even and odd degree numbers, that is:

$$% \begin{array}{lll} \hat{U}^{+} = U^{m+2k}_m (\zeta) & \hat{... ...} (\zeta) & \hat{\Psi}^{-}=\Psi^{m+2k+1}_m (\zeta), \end{array}$$ (30)

where $0 \leq k < +\infty$. The solutions of these two ODE systems are either symmetric or antisymmetric with respect to the equator.

The two sets of ODE can be written in the form:

$$% \;\left({\cal A} \frac{{\rm d}^2 }{{\rm d}\zeta^2} + {\cal... ... D} \frac{{\rm d} }{{\rm d}\zeta} + {\cal E} \right)\vec{\Xi},$$ (31)

where $\vec{\Xi}$ denotes

$$% \vec{\Xi}^{+} = \left\vert \begin{array}{l} \hat{U}^{+} \\ \ha... ...y}{l} \hat{U}^{-} \\ \hat{\chi}^{-} \\ \hat{\Psi}^{-} \\ \end{array} \right.$$ (32)

and where the matrices are defined by blocks as follows:

$$% {\cal A} = \;\left(\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0... ...box{B}_{22} & 0 \\ 0 & 0 & \mbox{B}_{33} \end{array} \right)$$ (33)

$$% {\cal C} = \;\left(\begin{array}{ccc} \mbox{C}_{11} & \mbo... ...ox{C}_{31} & \mbox{C}_{32} & \mbox{C}_{33} \end{array} \right)$$ (34)

$$% {\cal D} = \;\left(\begin{array}{ccc} 0 & 0 & 0 \\ \mbox... ...x{E}_{22} & -\mbox{E}_{21} \\ 0 & 0 & 0 \end{array} \right).$$ (35)

Equivalently, one can write Eq. (31) as:

$$\begin{array}{l}% \;\left(\mbox{B}_{11} \frac{{\rm d} }{{\rm... ...rm d}\zeta} + \mbox{C}_{33} \right)\hat{\Psi} = 0. \end{array}$$ (36)

Each sub-matrices can be expressed in terms of the two following functionals:

$$% {\mbox{I}}^m_{ \ell \ell '}({\rm G})= 2\pi \int_0^{\pi}~{\r... ...(\theta) \hat{Y}^m_{\ell '}(\theta) \sin \theta {\rm d} \theta$$ (37)

$$% {\mbox{J}}^m_{ \ell \ell '}({\rm G}) = 2\pi \int_0^{\pi}~{\... ..._{\ell '}}{\partial \theta}(\theta) \sin \theta {\rm d} \theta$$ (38)

where $\hat{Y}^m_\ell(\theta) = Y^m_\ell (\theta,\phi) {\rm e}^{-{\rm i}m\phi}$ is a normalized Legendre polynomial.

In Appendix B, all the coefficient of the sub-matrices are made explicit in terms of the function  $r(\zeta,\theta)$ and its first and second order derivatives as well as in terms of the enthalpy of the equilibrium model, its first and second order derivatives.

In the following, we consider the Cowling approximation thus neglecting the gravitational potential perturbation. The ODE system (31) is simplified accordingly and in particular reduces to the first order.

3 Tests and accuracy

The formalism and the numerical method presented in the previous section have been tested and the accuracy of the frequency determinations has been estimated.

3.1 Tests

A first test of the method has been performed in the case of axisymmetric ellipsoids of uniform density. We choose this configuration because the eigenvalue problem is fully separable using the oblate spheroidal coordinates  $(\xi, \eta, \phi)$defined as $(x = a \cosh \xi \sin \eta \sin \phi, \; y = a \cosh \xi \sin \eta \cos \phi, \; z = a \sinh \xi \cos \eta)$, where $0\leq \xi <+ \infty$, $0\leq \eta \leq \pi$ et $0\leq \phi \leq 2 \pi$. The eigenfrequencies obtained with this method were compared with the eigenfrequencies computed with our general method, $S(\theta)$ describing an ellipse. We found the same frequencies with a high degree of accuracy for arbitrary values of the ellipsoid flatness between 0 and 0.5. Moreover, as the flatness goes to zero, the frequencies were found to converge towards the values given by a first order perturbative analysis in terms of flatness. More details about this test are given in Lignières et al. (2001) and Lignières & Rieutord (2004).

The frequencies of axisymmetric p-modes in a self-gravitating uniformly rotating N=3 polytrope that will be presented in the following sections have been also tested. As shown in the previous section, the method involves lengthy analytical calculations of the coefficients of the ODE system (31). Terms involved in the non-rotating case have been tested by comparing our result with the p-modes frequencies in a non-rotating self-gravitating N=3 polytrope published in Christensen-Dalsgaard & Mullan (1994). The relative error is smaller than 10^-7 for the $\ell =0$ to 3, n=1 to 10 modes. In the rotating case, we compared our results with the ones obtained by solving the same problem but using a different form of the starting equations. This alternative system of equations aims at including the Coriolis force; thus the variables and the resulting ODE systems to be solved are different. We verified that when the terms involving the Coriolis force are omitted from the equations, eigenfrequencies presented in the next section are recovered with a very high precision (Reese et al. 2006a). For instance, the maximum relative error on the eigenfrequency for all the frequencies computed at $\Omega / \Omega _{\rm K} = 0.59$ has been found of the order of 10^-6, for a given set of numerical parameters. This last test gives us strong confidence in the method and its implementation.

	Figure 1: Evolution of the frequency relative error, $E(L) = \vert\omega (L) - \omega (L^{{\rm max}})\vert / \omega (L^{{\rm max}})$, as the spatial resolution in latitude is increased. Two modes labeled ( $\ell=1,\;n=1$) and ( $\ell=6,\;n=8$) are considered at the rotation rate $\Omega /\Omega _{\rm K} = 0.46$ with $L^{{\rm max}} = 80$ and Nr = 61.
Open with DEXTER

3.2 Accuracy

As pointed out by Clement (1998,1981), accurate numerical solutions of the 2D eigenvalue problem are not easy to obtain. It is therefore important to estimate the accuracy of our method. In the following we first investigate the influence of the spatial resolution on the eigenfrequencies and then consider other sources of errors.

A relative spectral error E is defined as the absolute value of the relative difference between the frequency computed at a given resolution and the frequency obtained at the maximum resolution considered. Let us first consider the effects of the angular resolution. Figure 1 displays $E(L) = \vert\omega (L) - \omega (L^{{\rm max}})\vert / \omega (L^{{\rm max}})$ as a function of L, the truncation order of the spherical harmonic expansion, for two axisymmetric modes labeled ( $\ell=1,\;n=1$) and ( $\ell=6,\;n=8$) whose spatial structures are dominated by large and small length scales, respectively (the labeling of the mode will be described in the next section). The maximum angular resolution is $L^{{\rm max}} = 80$, the resolution in the pseudo-radial coordinate is fixed to N_r = 61 and the rotation rate is $\Omega /\Omega _{\rm K} = 0.46$. In the same way, Fig. 2 illustrates the effects of the pseudo-radial resolution by showing $E(N_r) = \vert\omega (N_r) - \omega (N_r^{{\rm max}})\vert / \omega (N_r^{{\rm max}})$ as a function of N_r, the truncation order of the Chebyshev polynomial expansion, for the same modes and rotation rate. The maximum radial resolution is $N_r^{{\rm max}} = 61$ and the latitudinal resolution is fixed to L=62. In both figures, we observe that the error first decreases and then reaches a plateau which means that a better approximation of the eigenfrequency cannot be obtained by increasing the spatial resolution. The plateau are significantly higher for the ( $\ell=6,\;n=8$) mode than for the ( $\ell=1,\;n=1$) mode. We verified that this difference is due to the presence of smaller radial length scales (rather than to smaller angular length scales).

Even when the spatial resolution is sufficient, two other sources of numerical errors can indeed limit the accuracy of eigenfrequency determination. First, the component of the matrix L and M being computed numerically, they are determined with a certain error. Second, even when this error is reduced to round-off errors, the accuracy of the algebraic eigenvalue solver, the Arnoldi-Chebyshev algorithm, remains limited.

	Figure 2: Evolution of the frequency relative error, $E(N_r) = \vert\omega (N_r) - \omega (N_r^{{\rm max}})\vert / \omega (N_r^{{\rm max}})$as the resolution in the radial coordinate is increased. Two modes labeled ( $\ell=1,\;n=1$) and ( $\ell=6,\;n=8$) are considered at a rotation rate $\Omega /\Omega _{\rm K} = 0.46$ with $N_r^{{\rm max}} = 61$ and L = 62.
Open with DEXTER

Errors on the matrix component that arise from quadratures (see Eqs. (37) and (38)) can approach round-off errors using weighted Gauss-Lobatto quadratures. The other source of error in the matrix components comes from the computation of equilibrium quantities. Indeed, the accuracy of the enthalpy, its first and second derivatives and the surface shape, is at best limited by the effect of round-off errors on the convergence of the algorithm used to compute the polytropic stellar models. The effect of these errors on the eigenfrequencies have been investigated and appears to be smaller than the effect of the Arnoldi-Chebyshev algorithm itself which is now described.

As any solver in linear algebra, the Arnoldi-Chebyshev algorithm amplifies the round-off error that affect the matrix components. Thus, the error on the eigenvalue and the associated eigenvector is usually much larger than the round-off error of double precision arithmetic. The accuracy of the Arnoldi-Chebyshev algorithm has been studied in details by Valdettaro et al. (2006) in the context of inertial modes in a spherical shell where the matrix component are known analytically. Theoretically, it can be estimated by computing the spectral portrait of the eigenvalue problem $[M] X = \lambda [Q] X$, which shows how small variations of [M] and [Q]affects the determination of each eigenfrequencies. In fact, as the iterative Arnoldi-Chebyshev algorithm requires an initial guess of the eigenfrequency, a practical alternative to measure the accuracy of a frequency determination is to compute frequencies for slightly different values of the initial guess. This has been done for a large number (100) of initial guess values randomly distributed around the eigenfrequency of the ( $\ell=1,\;n=1$) and ( $\ell=6,\;n=8$) modes. The histogram in Fig. 3 shows the resulting frequencies distribution around a most probable mean eigenfrequency. The width of the histogram determined by the standard deviation of the distribution provides a measure of the algorithm accuracy. The standard deviation $\sigma$ is equal to 5.6 $\times$ 10^-6 for the ( $\ell=6,\;n=8$) mode and to 6.2 $\times$ 10^-10for the ( $\ell=1,\;n=1$) mode. The error thus grows with the radial order of the mode, this trend being general in our results (as in Valdettaro et al. 2006). Moreover, the width of the histogram does not depend on the amplitude of the initial guess perturbation provided it is sufficiently small.

We also observe that, except for the dependence of the ( $\ell=1,\;n=1$) frequency on the angular resolution, the levels of the plateau shown in Figs. 1 and 2 are of the same order as the error of the algorithm. It means that, in these cases, the changes in the matrix component and size associated with the modification of the resolution have a similar effect on the frequency as varying the initial guess of the algorithm. However, the convergence of the ( $\ell=1,\;n=1$) frequency at a 10^-14 level, much lower than the 6.2 $\times$ 10^-10 accuracy of the algorithm, shows that it is not always true and that the spectral error can underestimate the true error.

Figure 3: Histogram of 100 frequencies obtained for 100 different values of the initial guess of the Arnoldi-Chebyshev algorithm randomly chosen in a small interval around $\omega _m = 12.0547$. The standard deviation of this frequency distribution $\sigma = 5.6$ $\times$ 10-6 measures the algorithm relative error on the frequency for this particular mode labeled ( $\ell=6,\;n=8$). The spatial resolution is Nr = 61 and L = 62 and the rotation rate is $\Omega /\Omega _{\rm K} = 0.46$.

Open with DEXTER

Although it is too demanding to compute a global accuracy by repeating the statistical study on the initial guess for all eigenfrequencies, the relative accuracy on all tested frequencies is always better than 2 $\times$ 10^-5 using double precision arithmetic.

Note that another potential source of error will be discussed below when describing avoided crossings between modes.

Figure 4: Evolution of all the computed p-modes frequencies from $\Omega =0$ to $\Omega / \Omega _{\rm K} = 0.59$. The frequencies have been adimensionalized by  $(GM/R_{\rm p}^3)^{1/2}$ because we expect that the polar radius $R_{\rm p}$ does not change much as the rotation of the star increases. Non-rotating $\ell = (0,...,7)$, $n= (1,...,n_{\rm max})$ p-modes have been followed by progressively increasing the rotation. This mode tracking requires special care when an avoided crossing occurs between two modes of the same equatorial parity. The figure on the left shows an overview of the frequency evolution while the two right figures display zooms to illustrate avoided crossings between the $\ell =1, \; n=6$ and $\ell =5, \; n=5$ modes and the $\ell =0, \; n=4$ and $\ell =4, \; n=3$ modes, respectively. Although the two "interacting'' modes have a mixed character near the closest frequency approach, their original properties are recovered after the crossing which enables to unambiguously follow and label the modes. This is illustrated in Fig. 5 by considering the spectra of Legendre expansion components of the $\ell =0, \; n=4$ and $\ell =4, \; n=3$ modes at the rotation rates marked by an arrow. Note that in the above figures crossings do occur between equatorially symmetric and anti-symmetric modes. In the global view, there are two examples of discontinuous frequency changes due to avoided crossing with modes which frequency is not represented on the figure. Actually, the $\ell = 8$, n= (1,2,3) modes have been displayed in this view to avoid more discontinuous changes.

Open with DEXTER

4 Results

The parameter range of the calculations is first presented. Then, we describe the method used to label the eigenmodes, the structure of the frequency spectrum, some properties of the eigenfunctions and the differences with perturbative methods.

4.1 Parameter range

Self-gravitating uniformly rotating polytropes of index N=3 and specific heat ratio ${\Gamma}_{1,0} = 5/3$ have been computed for rotation rates varying from $\Omega =0$ up to $\Omega / \Omega _{\rm K} = 0.59$. In this range, the flatness of the star's surface $\epsilon = 1 - R_{\rm p}/R_{\rm e}$ increases from 0 to 0.15.

Low frequency axisymmetric p-modes have been computed for each polytropic model. We started with the non-rotating model and computed the $\ell = (0,...,7), \; n= (1,...,n_{\rm max})$ axisymmetric p-modes, the largest radial order depending on the degree $\ell$: $n_{\rm max}=10$ for $\ell=(0,1)$, $n_{\rm max}=9$ for $\ell=(2,3,4)$and $n_{\rm max}=8$ for $\ell=(5,6,7)$. All these 71 modes were then calculated at higher rotation rates by progressively increasing the rotation of the polytropic model. In the next section, we explain how we could track and label them from zero rotation to $\Omega / \Omega _{\rm K} = 0.59$.

Figure 5: Evolution of the Legendre components $C(\ell ) = {\rm max}_{n_r} \vert U(\ell ,n_r )\vert$ of the $\ell =0, \; n=4$ and $\ell =4, \; n=3$ modes during the avoided crossing shown in Fig. 4. Prior to ( $\Omega / \Omega _{\rm K} = 0.38$) and after ( $\Omega / \Omega _{\rm K} = 0.55$) the avoided crossing, the spectrum of Legendre components peaks at a given degree while, near the closest approach ( $\Omega /\Omega _{\rm K} = 0.46$), the double peaks of the spectra show the mixed character of the eigenmodes.

Open with DEXTER

4.2 Mode labeling

In the absence of rotation, modes are identified and classified by the three "quantum'' numbers $n, \ell, m$ characterizing their radial, latitudinal and azimuthal structure respectively. Because of separability, independent 1D eigenvalue problems are solved for each pair ($\ell$, m) and it is then an easy task to order the computed frequencies, the order n additionally indicating the number of radial nodes of the mode. By contrast, in the presence of rotation, independent 2D eigenvalue problems are solved for a given $\mid \! m \! \mid$ and a given equatorial parity. The computed modes are then obtained without a priori information about their latitudinal and radial structures. Therefore, an important issue is whether it is possible to define a meaningful classification of these modes. In this work, we investigate the possibility of associating unambiguously each mode with a non-rotating mode thus identifying it by the three quantum numbers  $n, \ell, m$of the non-rotating mode. Similarly, Clement (1986) followed some equatorially symmetric acoustic modes to high rotation rates but in a limited frequency range and using low spatial resolution calculations.

In practice, instead of backtracking modes towards zero rotation, we started at zero rotation with a mode we are interested in and tried to follow it by progressively increasing the rotation. We managed to track all the $\ell = (0,...,7)$, $n= (1,...,n_{\rm max})$ axisymmetric p-modes from $\Omega =0$ to $\Omega / \Omega _{\rm K} = 0.59$, a global view of the eigenfrequencies evolution being displayed in Fig. 4 (left panel). As explained below, the main difficulty comes from avoided crossings between modes of the same equatorial parity.

Zooms in the $\omega - \Omega$ plane displayed in Fig. 4 (right panels) provide two examples of avoided crossings respectively between odd ( $\ell =1, \; n=6$ and $\ell =5, \; n=5$) and even ( $\ell =0, \; n=4$ and $\ell =4, \; n=3$) modes. Modes tends to cross because their frequencies are not affected in the same way by the centrifugal force but, as two eigenstates of the same parity cannot be degenerate, an avoided crossing takes place during which the two eigenfunctions exchange their characteristics. This exchange of property is illustrated in Fig. 5 in the case of the ( $\ell =0, \; n=4$), ( $\ell =4, \; n=3$) crossing. A mean Legendre spectrum is displayed before, near the closest frequency separation and after the avoided crossing. The mean Legendre spectrum of a field U is defined as $C (\ell) = {\rm max}_{n_r} \vert U(\ell,n_r )\vert/ {\rm max} \vert U(\ell,n_r )\vert$, where $U(\ell,n_r)$ are the components of the Legendre/Chebyshev expansion, n_r being the degree of the Chebyshev polynomial. The quantity $C(\ell)$ thus represents the largest Chebyshev component for a given value of $\ell$ normalized by the maximum over all spectral components. The mean Legendre spectra peak at one characteristic degree before and after the avoided crossing, thus showing that the modes recover their original properties after the crossing and therefore can be unambiguously recognized. Up to the fastest rotation considered, the $\ell =0{-}7, \; n=1{-}10, \; m=0$, p-modes undergo a limited number of avoided crossing and could be followed unambiguously.

It remains that near the crossing the labeling is somewhat ambiguous. First, it is difficult to define a criterion to assign a label. Here, we mostly use the degree at which the mean Legendre spectrum reaches a maximum. But it occurred that the two interacting modes peak at the same degree in which case we determined the location of the smallest frequency separation. Second, as shown by Fig. 5, a more fundamental problem is that a single label cannot reflect the mixed nature of the eigenfunction.

	Figure 6: Frequency spectrum of $\ell =0{-}7, \; n=1-n_{\rm max}, m=0$ modes at $\Omega =0$ ( top panel) and $\Omega / \Omega _{\rm K} = 0.59$ ( bottom panel). The degree number $\ell$ associated with the frequency is shown by the height of the vertical bar.
Open with DEXTER

Another issue related to avoided crossings concerns their influence on the accuracy of the eigenfunction computation. Indeed, if a large-scale, well resolved eigenfunction undergoes an avoided crossing with a small-scale unresolved mode, the accuracy of the eigenfunction determination will be affected. The effect on the frequency accuracy should be small as the frequency gap induced by the avoided crossing of two modes of well separated length scales is small. But, at the closest approach, the eigenfunctions will be much affected. At zero rotation, the highest degree mode present in our frequency range is $\ell = 51 , n=1$. Thus, if one of the low degree modes that we computed undergoes an avoided crossing with a mode of such a high degree, the high degree mode should be resolved to ensure an accurate determination of the eigenfunction of the low degree mode.

4.3 The structure of the frequency spectrum

The effect of the centrifugal force on the acoustic frequency spectrum of axisymmetric modes is investigated. The mean modifications of the spectrum are discussed then we investigate how regularities in the frequency spacings evolve with rotation. Finally, differences between equatorially symmetric and anti-symmetric modes are outlined.

4.3.1 Global spectrum evolution

Figure 6 compares the frequency spectrum of the $\ell =0{-}7, n=1-n_{{\rm max}}, m=0$ modes at $\Omega =0$(upper panel) and at $\Omega / \Omega _{\rm K} = 0.59$ (lower panel), the height of the vertical bars corresponding to the degree $\ell$ of the mode. It appears that the centrifugal force induces a mean contraction of the frequency spectrum. This is expected as the decrease of the sound speed and the increase of the star volume induced by the centrifugal force both tend to lessen the frequency of acoustic modes.

To illustrate this effect, consider a spherically symmetric decrease of the effective gravity. In a homologous series of spherical models of increased volume V, the decreasing rate of the frequencies  $\Delta \omega/\omega$ is $-(1/2) (\Delta V/V)$, as the normalized frequencies  $\omega/(GM/R^3)^{1/2}$ remain constant. For non-homologous spherically symmetric changes, $\Delta \omega/\omega$ is asymptotically equal to $- \Delta (\ln \int_0^{R_{\rm p}} {\rm d}r/c)$for high order modes verifying the following asymptotic formula valid for low degree and high order p-modes (Tassoul 1980):

$$% \omega = \frac{\pi}{\int_0^R \frac{{\rm d}r}{c}} (n + (\ell + 1/2)/2 + \alpha)$$ (39)

where 1/ $\int_0^R \frac{{\rm d}r}{c}$ is the sound travel time along a stellar radius and $\alpha$ is a constant. When, as in these two previous cases, $\Delta \omega/\omega$ does not depend on the frequency, the concentration of the frequency spectrum is homothetic.

Figure 7: Regularities in the frequency spacings of axisymmetric (m=0) modes. The large frequency separation between modes of consecutive order $\Delta _n = \omega _{n \ell } - \omega _{n-1 \ell }$, the frequency separation between $\ell +2$ and $\ell$ modes, $\Delta _{2,\ell }= \omega _{n,\ell +2} - \omega _{n,\ell }$, and the small frequency separation $\delta = \Delta _n - \Delta _{2,\ell }$ are displayed as a function of the radial order nfor four different rotation rates, a)  $\Omega =0$, b)  $\Omega / \Omega _{\rm K} =0.32$, c)  $\Omega /\Omega _{\rm K} = 0.46$ and d)  $\Omega / \Omega _{\rm K} = 0.59$. We plotted the opposite of the small frequency separation $-\delta$ for clarity. Continuous lines have been drawn between frequencies of the same degree $\ell$.

Open with DEXTER

This is clearly not the case here since the frequencies cross each other (see Fig. 4). But there is still an average contraction rate which is of the order of  $-(1/2) (\Delta V/V)$, where now V is the volume of the centrifugally distorted star. In addition, the contraction rates of individual frequencies appears to be comprised between the logarithmic derivative of the sound travel times computed respectively along the polar and equatorial radii:

$$% \partial_{\Omega} \left(\ln \int_0^{R_{\rm p}} \frac{{\rm d... ...a} \left(\ln \int_0^{R_{\rm e}} \frac{{\rm d}r}{c}\right)\cdot$$ (40)

Another interesting property is that, at small rotation rates, say $\Omega/\Omega_{K} \leq 0.05$, the contraction rate $\partial_{\Omega} (\ln \omega)$ tends to be independent of $\ell$ and n for the large degree modes $\ell \geq 3$. This suggests that an asymptotic regime exists for modes with horizontal wavelengths smaller than the dominant length scales of the centrifugal distortion. In this regime, the contraction rate has a constant value that is not equal to  $-(1/2) (\Delta V/V)$. We already found such behaviour in the case of homogeneous ellipsoids (Lignières et al. 2001) where a perturbative analysis shows that the contraction rate of axisymmetric modes is constant for high $\ell$ and n and that it can be related to the increase of the ellipse perimeter.

Nevertheless, for the low degree modes $\ell \leq 2$ below $\Omega/\Omega_{K} \approx 0.05$, and for all modes at higher rotation rates, $\partial_{\Omega} (\ln \omega)$ depends on $\ell$ and n. This differential effect modifies the structure of the frequency spectrum as the rotation increases.

4.3.2 Regular frequency spacings

In a non-rotating star, the frequency spectrum presents some regular frequency spacings which can be accounted for by an asymptotic theory in the high frequency limit $\omega \rightarrow \infty$. The asymptotic formula (39), valid for low degree and high order modes, shows that the large frequency separation between modes of consecutive order n, $\Delta_n = \omega_{n+1, \ell} - \omega_{n, \ell}$, does not depend on $\ell$ and n and is equal to $\pi/$ $\int_0^R \frac{{\rm d}r}{c}$. A more detailed asymptotic analysis also shows how the so-called small frequency separation $\delta = \omega_{n+1, \ell} - \omega_{n, \ell+2}$ vanishes as a function of the frequency. Although our calculations are restricted to the low frequency part of the acoustic spectrum, we observe a clear tendency towards these asymptotic behaviors in the non-rotating case. We can therefore investigate whether these properties are modified by rotation.

Figure 7 presents the large frequency separation  $\Delta_{n}$and the frequency separation between consecutive modes of the same order and parity:

$$% \Delta_{2,\ell}= \omega_{n,\ell+2} - \omega_{n,\ell},$$ (41)

as a function of the radial order n for four different rotation rates, (a)  $\Omega =0$, (b)  $\Omega / \Omega _{\rm K} =0.32$, (c)  $\Omega /\Omega _{\rm K} = 0.46$ and (d)  $\Omega / \Omega _{\rm K} = 0.59$. As in the previous figures, the frequencies are adimensionalized by  $(GM/R_{\rm p}^3)^{1/2}$. Continuous lines have been drawn between frequencies of the same degree $\ell$. We first observe that the large frequency separation tends to be independent of n and $\ell$ at all rotation rates. In accordance with the mean contraction of the frequency spectrum mentioned above, the large frequency separation decreases with rotation. It is always between $\pi/$ $\int_0^{R_{\rm p}} {\rm d}r/c$ and $\pi/$ $\int_0^{R_{\rm e}} {\rm d}r/c$.

The dispersion of the large frequency separations around their mean value also has an interesting evolution with rotation. In the non-rotating case, the dispersion reflects a regular departure from the asymptotic limit. It is larger for high degrees and monotonically decreases with frequency (see Fig. 7a). In the rotating cases, the dispersion is not as regular. The largest departures, some of which are most clearly visible in Fig. 7c, can be attributed to an ongoing avoided crossing. The residual dispersion is irregular and decreases with rotation. At  $\Omega / \Omega _{\rm K} = 0.59$, if we exclude all n < 4 values from our sample, the mean large frequency separation  $\langle \Delta_n \rangle$ is equal to  $1.095 (GM/R_{\rm p}^3)^{1/2}$ and its standard deviation is  $0.017 \langle \Delta_n \rangle$.

We now consider the small frequency separation $\delta = \Delta _n - \Delta _{2,\ell }$. As expected, in the absence of rotation the small frequency separation tends to vanish as n increases. But, at $\Omega / \Omega _{\rm K} =0.32$, the small frequency separation no longer decreases with n for some values of $\ell$and for the higher rotation rates it becomes nearly constant. At the same time, the $\Delta_{2,\ell}$ separation becomes more and more uniform as rotation increases. As shown in Fig. 7b, $\Delta_{2,\ell}$ becomes approximatively constant with n first for low degree modes while it still increases with frequency for high degree modes. In addition, equatorially antisymmetric modes reach this new regime at a lower rotation rate than the equatorially symmetric modes of similar degree. This is illustrated in Fig. 7d by the $\Delta_{2,\ell=4}$ curve which still remains above the mean $\Delta_{2,\ell}$ value while the $\Delta_{2,\ell=5}$separation collapses with the other curves. At $\Omega / \Omega _{\rm K} = 0.59$, if we exclude all n < 4 values from our sample, the mean frequency separation  $\Delta_{2,\ell}$ is equal to  $0.387 (GM/R_{\rm p}^3)^{1/2}$ and its standard deviation is $0.033 \langle \Delta_n \rangle$.

As a consequence of the near uniformity of $\Delta_{n}$ and  $\Delta_{2,\ell}$, the frequencies of low degree and high order modes can be approximated by the following expressions:

$$% \tilde{\omega}_{n \ell} = \left\{\begin{array}{l} n \delta_... ...ll} + \alpha^- \;\mbox{if} \; \ell = 2p + 1 \end{array}\right.$$ (42)

where $\delta_n=$ $\langle \Delta_n \rangle$, $\delta_{\ell}=$ $\langle \Delta_{2,\ell} \rangle$, $\alpha^+$ and $\alpha^-$ only depend on the equilibrium model. Using a reference frequency to determine the $\alpha$ constants (the $\ell=0, \; n=8$ frequency for $\alpha^+$ and the $\ell=1, \; n=8$ frequency for $\alpha^-$), we computed the root mean square error $\sqrt{1/N \sum ({\tilde{\omega}} - \omega)^2}$ and the maximum error made in using the approximate expressions (42). For a frequency subset containing the n>4 and $\ell<5$ modes, the rms error is $0.017 \delta_n$ while the maximum error amounts to  $0.05 \delta_n$. Both errors are a very small fraction of the large separation which shows that Eq. (42) yields useful approximations of the frequency spectrum.

4.3.3 Equatorially symmetric versus anti-symmetric frequency spectra

We have seen that the regular frequency spacings $\Delta_{n}$ and  $\Delta_{2,\ell}$ have similar values for symmetric and anti-symmetric modes with respect to the equator. The evolution of the equatorially symmetric and anti-symmetric frequency spectra are nevertheless quite different. Indeed, considering two modes of similar frequency but of opposite equatorial parity, the frequency of the symmetric mode generally decreases faster with rotation than the frequency of the antisymmetric modes. The consequence is that the frequency separation between modes of consecutive degree (and thus of opposite parity) $\Delta_l = \omega_{n,\ell+1} - \omega_{n,\ell}$tends to increase when $\ell$ is even and to decrease when $\ell$ is odd. The frequency separation $\Delta_l$ can even become negative which implies that, contrary to the non-rotating case, frequencies of a given order n do not increase monotonically with the degree $\ell$. This striking modification of the usual frequency ordering is apparent in Fig. 6 where the $(\ell =2 , n)$ frequencies are smaller than the $(\ell =1 , n)$ frequencies for all the order n that we calculated, that is n= (1,...,10). In the same way, the $(\ell =4 , n)$ frequencies are smaller than the $(\ell =3 , n)$ frequencies if $n \geq 3$, and again the $(\ell =6 , n)$ frequencies are smaller than the $(\ell =5 , n)$ frequencies if $n \geq 5$.

Figure 8: Isocontours of the $\ell =4, n=4$ mode amplitude in a meridional plane as a function of the rotation a)  $\Omega =0$, b)  $\Omega / \Omega _{\rm K} =0.32$, c)  $\Omega /\Omega _{\rm K} = 0.46$ and d)  $\Omega / \Omega _{\rm K} = 0.59$. The amplitude is normalized to the maximum of its absolute value. Continuous lines correspond to positive amplitudes, dashed lines to the zero amplitude and dotted lines to negative amplitudes. At zero rotation, the angular distribution is given by the $\hat{Y}^0_4 (\theta)$ Legendre polynomial while the radial distribution is characterized by the surface concentration of the amplitude and the presence of n=4 nodes in the inner part. For larger rotation rates, the largest amplitudes concentrates toward the equator. We also note that the number of radial nodes decreases along the polar radius while it increases along the equatorial radius.

Open with DEXTER

4.4 Equatorial concentration

In this section, we focus on the most notable effect of the centrifugal force on the eigenmodes, namely the equatorial concentration and consider its consequences on the mode visibility. Note that Clement (1981) also reported an equatorial concentration of the equatorially symmetric modes that he calculated.

Figure 8 shows this effect on the ( $\ell =4, n=4$) mode. Contours of the amplitude of the Lagrangian pressure perturbation are plotted in a meridional plane for increased rotation rates, (a)  $\Omega =0$, (b)  $\Omega / \Omega _{\rm K} =0.32$, (c)  $\Omega /\Omega _{\rm K} = 0.46$ and (d)  $\Omega / \Omega _{\rm K} = 0.59$. We observe that the number of nodes increases along the equatorial radius and decreases along the polar one. Along the surface, the number of nodes remains equal to $\ell$ before $\Omega / \Omega _{\rm K} = 0.59$ where additional nodes appear. The equatorial concentration is clearly seen in the outermost layers.

In Fig. 9, the equatorial concentration is shown for other modes including the lowest and highest degree modes of our sample as well as symmetric and anti-symmetric modes. The latitudinal variation of the mode amplitude is displayed at the surface for the following modes, (a)  $\ell =0, n=1$, (b)  $\ell =1, n=1$, (c)  $\ell =6, n=8$, (d)  $\ell =7, n=8$. In each case, the equatorial concentration grows with rotation. At the largest rotation rate, symmetric modes are maximal at the equator while anti-symmetric modes peak at small latitudes since they must vanish at the equator. The contrast between these maxima and the polar amplitude is strong.

Figure 9: The mode amplitude at the surface of the polytropic model as a function of the rotation rate a)  $\ell =0, n=1$, b)  $\ell =1, n=1$, c)  $\ell =6, n=8$, d)  $\ell =7, n=8$. The amplitude is normalized to the maximum of its absolute value. While the angular distribution is given by the corresponding Legendre polynomial  $\hat{Y}^0_\ell (\theta)$ in the absence of rotation, the oscillation amplitude progressively concentrates towards the equator ( $\theta = \pi /2$) as rotation increases.

Open with DEXTER

The equatorial concentration reveals a modification of the resonant cavity of the acoustic waves. In particular, the reduction of the volume of the resonant cavity should tend to increase the frequency. The equatorial concentration seems also to be associated with the near-uniformity of the frequency separation  $\Delta_{2,\ell}$. At low rotation rates, the concentration is not completed and  $\Delta_{2,\ell}$ is clearly not constant. At the largest rotation rate, all modes are concentrated near the equator and  $\Delta_{2,\ell}$ is nearly uniform.

Figure 10: The disk-averaging factor D(i) is shown as a function of the inclination angle i for various axisymmetric modes at two different rotation rates $\Omega =0$ a) and $\Omega / \Omega _{\rm K} = 0.59$. The degree of the modes varies from $\ell =0$ to $\ell =7$. In the rotating case, the surface distribution also depends on the order of the mode. Two values n=1 b) and n=8 c) have been considered at $\Omega / \Omega _{\rm K} = 0.59$.

Open with DEXTER

Besides its effect on the frequency spectrum, the equatorial concentration of eigenmodes should induce a profound modification of the mode visibility as compared to the non-rotating case. The photometric mode visibility is determined by the integration over the visible part of the star's perturbed surface of the radiation intensity perturbations associated with a particular pulsation mode. Rigorous calculations of photometric visibilities are beyond the scope of the present paper as they require non-adiabatic calculations of the oscillation modes and stellar atmosphere models (e.g. Daszynska-Daszkiewicz et al. 2002). But we can still determine the effects of averaging the perturbations over the visible surface which have a direct impact on the visibility. The disk-averaging factor is defined as:

$$% D(i)= \frac{1}{\pi R_{\rm e}^2 \delta T_0} \int\!\!\!\!\int_{S_v} \delta T (\theta,\phi) {\rm d} {\bf S} \cdot {\bf e}_i$$ (43)

where i is the inclination angle between the line-of-sight and the rotation axis, ${\bf e}_i$ is a unit vector in the observer's direction and $\delta T$ is the spatial part of the Lagrangian temperature perturbation at the stellar surface, $\delta T$ being proportional to the velocity divergence $\chi$ in the approximation of adiabatic perturbations. The mode amplitude is normalized by  $\delta T_0$ the root mean square of the perturbation over the whole stellar surface

$$% \delta T_0 = \;\left(\int\!\!\!\!\int_{S} \delta T^2 (\theta,\phi) {\rm d} S \right)^{1/2}$$ (44)

and the visible surface S_v has been normalized by  $\pi R_{\rm e}^2$, the visible surface of a star seen pole-on. With these normalizations the disk-averaging factor of a radial mode seen pole-on is unity.

In the absence of rotation, the surface distribution of modes is determined by a unique spherical harmonic and the disk-averaging factor takes a simple analytical form (Dziembowski 1977). For even degree and for $\ell =1$, the disk-averaging factor varies with the inclination angle as the Legendre polynomial  $\hat{Y}^m_\ell(i)$ while it vanishes altogether for odd degree $\ell \ge 3$. For rotating stars, the method of the calculation is detailed in Appendix C. Note that for modes that are equatorially anti-symmetric and axisymmetric, the disk-averaging factor also has a simple dependency on the inclination angle as it is proportional to $\cos(i)$.

Figure 10 shows the disk-averaging factor of various axisymmetric modes as a function of the inclination angle. The non-rotating case is displayed in Fig. 10a where $\ell =0$ to $\ell =7$ modes are considered. We recall that, at $\Omega =0$, modes of different radial orders but same $\ell$ and m have the same surface distribution. This is not true in the rotating case and, at $\Omega / \Omega _{\rm K} = 0.59$, Figs. 10b and c present the disk-averaging factor for modes of the same degree numbers but for two different radial orders n=1 and n=8, respectively. Note also that the disk-averaging factor was allowed to take a negative value for clarity of the figure although it is its absolute value that is relevant for the mode's visibility. Figure 10 shows that rotation strongly modifies the dependency of the disk-averaging factor on the inclination angle as well as on the degree number.

Figure 11: The variation of the disk-averaging factor as a function of the degree $\ell$ is shown at two fixed values of the inclination angle (i = 0 and $i= \pi /2$). The three curves in each figure correspond to the non-rotating case and, at $\Omega / \Omega _{\rm K} = 0.59$, to two different orders, n=1 and n=8. The absolute value of the disk-averaging factor has been rescaled by its maximum value among the different degree considered to outline its dependency on the degree number. In sharp contrast to the non-rotating case, the disk-averaging factor at $\Omega / \Omega _{\rm K} = 0.59$ does not show a strong decrease with $\ell$.

Open with DEXTER

Figure 10c shows that for all n=8 equatorially symmetric modes the absolute value of the disk-averaging factor tends to increase with the inclination angle. This is due to the equatorial concentration of these modes (see for example the surface distribution of the $(\ell=6, n=8)$ mode at $\Omega / \Omega _{\rm K} = 0.59$ shown in Fig. 9c). This tendency is less pronounced for the n=1 symmetric modes shown in Fig. 10b (except for the $\ell=2$ mode) although these modes are also equatorially concentrated. This is due to a cancellation effect between positive and negative perturbations concentrated near the equator as illustrated by the surface distribution of the $(\ell=0, n=1)$ mode in Fig. 9a. The non-rotating case strongly differs since the absolute value of the disk-averaging factor for even degree $\ell > 0$ modes does not vary monotonically with the inclination. Indeed, they have $\ell/2$ nodes between 0 and $\pi/2$. For odd $\ell$ modes, the disk-averaging factor is also modified by rotation since it no longer vanishes for $\ell \ge 3$. This occurs because the projected elementary surfaces  ${\rm d} {\bf S} \cdot {\bf e}_i$ are no longer symmetric with respect to the observer's direction and because the projection of the eigenmode surface distribution onto the Legendre polynomial  ${\hat{Y}^0_1}$ is not zero for $\ell \ge 3$ modes.

In non-rotating stars, the cancellation effect between positive and negative perturbations results in a rapid decrease of the disk-averaging factor as the degree $\ell$ of the mode increases. Consequently, modes above a certain degree $\ell \geq 3{-}4$ are not expected to be detectable with photometry and are therefore not included when trying to identify the observed frequencies. As shown in Fig. 11, this property must be reconsidered for rapidly rotating stars. The absolute disk-averaging factor normalized by its maximum value over the degree considered $0 \leq \ell \leq 7$, $\vert D(i)\vert/{\rm max}_{\ell}\vert D(i)\vert$, is plotted as a function of $\ell$ for two fixed values of the inclination angle, i = 0 in Fig. 11a and $i= \pi /2$in Fig. 11b. The three curves correspond to $\Omega =0$ and to $\Omega / \Omega _{\rm K} = 0.59$ for the n=1 and n=8 modes, respectively. In contrast to the non-rotating case, the disk-averaging factor has no tendency to decrease above $\ell=2$. Again, this can be explained by the equatorial concentration as modes of different degree have a similar surface distribution.

4.5 Comparison with perturbative methods

According to the perturbative analysis, centrifugal effects appear at second order in $\Omega$ (Saio 1981). To determine the second-order perturbative coefficient from our complete calculations, we performed a series of calculations for small rotation rates ( $\Omega= 0,1.8$ $\times$ 10^-3, 1.8 $\times$ 10^-2, 4.6 $\times$ 10^-2, 0.09 $\times$ 10^-2, ... in units of  $\Omega_{\rm K}$). From them, we determined the second-order perturbative coefficient, denoted $\omega_1$, as the limit of the ratio $(\omega(\Omega) - \omega_0)/ \Omega^2$, where $\omega_0$ denotes a non-rotating eigenfrequency. Thus the approximate frequencies valid up to the second order in $\Omega$ read $\omega_{{\rm pert}} = \omega_0 + \omega_1 \Omega^2$, where the frequencies are in units of  $(GM/R_{\rm p}^3)^{1/2}$ and the rotation is in units of  $\Omega_{\rm K}$. To assess the range of validity of the second order perturbative approach, we compared these approximate frequencies to the "exact'' frequencies. In Fig. 12, the relative differences between the two calculations, $(\omega - \omega_{{\rm pert}})/\omega$, is plotted as a function of the rotation rate for the $\ell =0{-}2, \; n=1{-}10$ modes. The departures computed for the other modes, $\ell = (3,...,7)$, $n= (1,...,n_{\rm max})$, are smaller than the extremal differences shown in Fig. 12 and are not displayed for clarity. The relative differences are generally larger for low degree modes and, for small rotation rates, are a monotonic function of the radial order n (an increasing function for the $\ell = 0{-}2$ modes shown in Fig. 12). As mentioned before the low degree modes seem to be sensitive to the precise form of the distortion that occurs at similar lengthscales. As rotation increases, it appears that higher than second order effects are important to describe the effect of the centrifugal distortion on these modes. The second order approximation is much better for large $\ell$ modes which are sensitive to global distortion properties.

	Figure 12: The relative difference between exact frequencies and their second order perturbative approximation (second order in terms of the small parameter  $\Omega / \Omega _{\rm K}$), namely $\delta \omega / \omega$ where $\delta \omega = \omega - \omega _{{\rm pert}}$is displayed as a function of the rotation rate for the $\ell =0{-}2, \; n=1{-}10, \; m=0$ modes.
Open with DEXTER

As compared to the observational uncertainties on the frequency determinations, the error made in using second order perturbative methods becomes rapidly significant as rotation increases. For a ratio $\Omega/\Omega_{\rm K} = 0.24$, corresponding to a typical $\delta$ Scuti star with an equatorial velocity of 100  ${\rm km~s}^{-1}$, the maximum absolute difference is $11~\mu {\rm Hz}$ which is much larger than typical observational uncertainties.

The fact that in our frequency sample the absolute difference increases with frequency suggests that departure from the perturbative approach could be detectable for moderately rotating stars pulsating in high order modes. In our data limited to $n \leq 10$, the largest relative difference is 2.87 $\times$ 10^-3for a  $\Omega/\Omega_{\rm K} =0.139$ rotation corresponding to a solar-type star with an equatorial velocity of 60  ${\rm km~s}^{-1}$. If we assume that the same relative difference holds for high order p-modes in the range of $2000~\mu {\rm Hz}$generated in convective envelopes of these stars, the absolute difference is $5.7~\mu {\rm Hz}$ for a typical $2000~\mu {\rm Hz}$ frequency. This would be also easily detectable given the observational uncertainties (Bouchy et al. 2005). However, a firm conclusion should await a direct comparison between the perturbative approach and the complete calculation for p-modes generated in the convective envelope of rotating solar-type stars.

5 Discussion and conclusion

A new non-perturbative method to compute accurate oscillation modes in rotating stars was presented. The accuracy of the computed frequencies has been obtained by testing the effect of the different parameters of the numerical method. Then, the effects of the centrifugal force on low frequency axisymmetric acoustic oscillation modes were investigated in uniformly rotating polytropic models of stars. Seventy-one low degree $\ell \leq 7$ and low order $n \leq 10$ modes were first determined at zero rotation and then tracked at higher rotation rates up to $\Omega / \Omega _{\rm K} = 0.59$.

In the frequency and rotation ranges considered in this paper, the zero rotation quantum numbers $\ell$ and n were used to label the modes. This labeling turned out to be meaningful since we found regular frequency spacings between modes of the same degree and consecutive orders, $\Delta_{n}$, and, within the subsets of modes of the same equatorial parity, between modes of the same order and consecutive degree, $\Delta_{2,\ell}$. We noted however that near avoided crossings, when the eigenfunction is a mix of the two "interacting'' modes, a unique label cannot reflect the actual eigenfunction structure. Although successful in the frequency and rotation ranges considered, it remains to be proved that this labeling can be performed in practice at higher rotation and at higher frequencies. Indeed, the main difficulty of the labeling procedure arises from the avoided crossing between modes of the same equatorial parity and such crossings will be more frequent as the eigenfrequency density increases with the frequency. The coupling between modes is also stronger at higher rotation rates. It might then be necessary to investigate tools other than the mean Legendre power spectrum to characterize the modes.

The study of the frequency spectrum showed a quite unexpected result, namely that, at the highest rotation rates, a new form of organization sets in after the zero rotation asymptotic structure of the spectrum has been destroyed. In the absence of rotation, the asymptotic theory is directly related to the spherical symmetry of the stars and ultimately to the integrability of the underlying ray dynamics. In the presence of rotation, the eigenvalue problem is not fully separable and the underlying acoustic ray dynamics is most probably not integrable. The regular spacings observed at high rotation rates were not expected. They might be the sign of a near-integrable ray dynamic rather than a chaotic system. These aspects will be investigated in a ray dynamic study of rotating polytropic models of stars.

Most importantly for asteroseismology, the existence of regular spacings in the spectrum can potentially provide tools for the mode identification in rapidly rotating pulsating stars. A complete acoustic frequency spectrum including $m \ne 0$ modes and the effects of the Coriolis acceleration should however be computed and analyzed to assess the practical relevance of these regular spacings.

Apparently, there is a relation between the new spectrum structure and the equatorial concentration of the mode amplitudes. A consequence would be that this spectrum structure does not apply to the whole spectrum. Indeed, sufficiently high degree modes should still be of the whispering gallery type (e.g. Rieutord 2001). Then, being so different from equatorially concentrated modes, they are not expected to follow the same regular spacings.

Another interesting issue is the difference between the modes of different equatorial symmetry. We have seen that although the structure of the symmetric and anti-symmetric spectra is similar, the frequency spectrum of the symmetric modes as a whole seems to evolve independently from the anti-symmetric spectrum. The equatorial symmetry also influences the "strength'' of the avoided crossings measured by the frequency separation at the closest frequency approach. As illustrated in Fig. 4 (right panels), avoided crossings between symmetric modes are always stronger than avoided crossings between anti-symmetric modes since they remain further apart.

Modes undergoing an avoided crossing are particular because they have close frequencies and similar eigenfunctions. As a consequence, both can be excited to observable levels by some excitation mechanism. They are therefore good candidates to explain the occurrence of close frequencies in observed spectra (Breger & Pamyatnykh 2006) as well as the associated amplitude variations induced by beating between the two close frequencies.

The most striking effect of the centrifugal force on the eigenfunction is the equatorial concentration of the mode amplitude. Again, the study of the ray dynamics should help specify the conditions in which the sound waves stay focused in the equatorial region. As compared to the non-rotating case, the equatorial concentration strongly modifies the integrated light visibility and in particular its variation with respect to the mode degree and the inclination angle. Our results showing a global increase of the disk-integration factor as the star is seen equator-on are compatible with observations of $\delta$ Scuti pulsations which also suggest an increase of the pulsation amplitudes with i (Suárez et al. 2002). Another finding of practical interest is that, for rapidly rotating stars, the cancellation effect of the disk averaging no longer sharply decreases with the degree of the mode and also varies with the order of the mode. Realistic calculations of the mode visibility including non-adiabatic calculations of the oscillation modes, stellar atmosphere models as well as the gravity and limb darkening effect will however be needed to draw firm observational conclusions.

The omission of the Coriolis force did not allow a complete treatment of the rotational effects. However, the effect of the Coriolis force vanishes for sufficiently large frequency (as the time scale of the Coriolis acceleration $1/\Omega$ becomes much larger than the pulsation period) while the modification of the equilibrium model by the centrifugal force affects all frequencies. Therefore, the results presented here should be useful for the high frequency part of the acoustic spectrum in rotating stars. In a companion paper (Reese et al. 2006b), we extend the present results by taking into account the Coriolis acceleration which, among other things, allows us to specify the domain of validity of perturbative calculations.

Acknowledgements

We thank L. Valdettaro for his contribution to the numerical part of this work and B. Georgeot for fruitful discussions. We also thank the referee for constructive comments. Numerical simulations have been performed with the computing resources of Institut du Développement et des Ressources en Informatique Scientifique (IDRIS, Orsay, France) and of CALcul en MIdi-Pyrénées (CALMIP, Toulouse, France) which are gratefully acknowledged.

References

• Bonazzola, S., Gourgoulhon, E., & Marck, J. A. 1998, Phys. Rev. D, 58, 104020 [NASA ADS] [CrossRef] (In the text)
• Bouchy, F., Bazot, M., Santos, N. C., Vauclair, S., & Sosnowska, D. 2005, A&A, 440, 609 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Breger, M., & Pamyatnykh, A. A. 2006, MNRAS, in press (In the text)
• Catala, C., Mangeney, A., Gautier, D., et al. 1995, GONG 1994. Helio- and Astro-Seismology from the Earth and Space, ed. R. K. Ulrich, E. J. Rhodes, & W. Dappen, ASP Conf. Ser., 76, 426
• Christensen-Dalsgaard, J., & Mullan, D. J. 1994, MNRAS, 270, 921 [NASA ADS] (In the text)
• Clement, M. J. 1981, ApJ, 249, 746 [NASA ADS] [CrossRef]
• Clement, M. J. 1986, ApJ, 301, 185 [NASA ADS] [CrossRef] (In the text)
• Clement, M. J. 1998, ApJS, 116, 57 [NASA ADS] [CrossRef]
• Daszynska-Daszkiewicz, J., Dziembowski, W. A., Pamyatnykh, A. A., & Goupil, M.-J. 2002, A&A, 392, 151 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Dintrans, B., & Rieutord, M. 2000, A&A, 354, 86 (In the text)
• Domiciano de Souza, A., Kervella, P., Jankov, S., et al. 2003, A&A, 407, L47 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Dziembowski, W. 1977, Acta Astronomica, 27, 203 [NASA ADS] (In the text)
• Dziembowski, W. A., & Goode, P. R. 1992, ApJ, 394, 670 [CrossRef] (In the text)
• Edmonds, A. R. 1960, Angular Momentum in Quantum Mechanics, second edn. (Princeton Univ. Press) (In the text)
• Espinosa, F., Pérez Hernández, F., & Roca Cortès, T. 2004, Proceedings of the SOHO 14 GONG 2004 Workshop, ed. D. Danesy, ESA SP-559, 424 (In the text)
• Ledoux, P. 1951, ApJ, 114, 373 [NASA ADS] [CrossRef] (In the text)
• Lignières, F., & Rieutord, M. 2004, IAU Symp., 215, 414 [NASA ADS] (In the text)
• Lignières, F., Rieutord, M., & Valdettaro, L. 2001, in SF2A 2001: Semaine de l'Astrophysique Française (EDP Sciences), ed. F. Combes et al., 127 (In the text)
• Pekeris, C. L. 1938, ApJ, 88, 189 [NASA ADS] [CrossRef] (In the text)
• Reese, D., Lignières, F., & Rieutord, M. 2006a, Commun. Asteroseismol., 147, 65 [NASA ADS] (In the text)
• Reese, D., Lignières, F., & Rieutord, M. 2006b, A&A, 455, 621 [EDP Sciences] [CrossRef] (In the text)
• Rieutord, M. 2001, in Stellar Astrophysical Fluid Dynamics, ed. M. J. Thompson, & J. Christensen-Dalsgaard (Cambridge Univ. Press), 99 (In the text)
• Rieutord, M., & Valdettaro, L. 1997, J. Fluid Mech., 341, 77 [NASA ADS] [CrossRef] (In the text)
• Rieutord, M., Dintrans, B., Lignières, F., Corbard, T., & Pichon, B. 2005, in SF2A 2005: Semaine de l'Astrophysique Française (EDP Sciences), ed. F. Combes et al., 759 (In the text)
• Saio, H. 1981, ApJ, 244, 299 [NASA ADS] [CrossRef] (In the text)
• Soufi, F., Goupil, M. J., & Dziembowski, W. A. 1998, A&A, 334, 911 [NASA ADS] (In the text)
• Suárez, J. C., Michel, E., Pérez Hernández, F., et al. 2002, A&A, 390, 523 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Tassoul, M. 1980, ApJS, 43, 469 [NASA ADS] [CrossRef] (In the text)
• Valdettaro, L., Rieutord, M., Braconnier, T., & Frayssé, V. 2006, [arXiv:physics/0604219] (In the text)
• Walker, G., Matthews, J., Kuschnig, R., et al. 2003, PASP, 115, 1023 [NASA ADS] [CrossRef]
• Yoshida, S., Yoshida, S., & Eriguchi, Y. 2005, MNRAS, 356, 217 [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

Appendix A: Linear operators in spheroidal coordinates

Let us now express the linear operators involved in Eqs. (22)-(24), using the spheroidal coordinates given by Eq. (25). We need the general expression of the divergence

$$% \vec{\nabla} \cdot \vec{V} = \frac{1}{\sqrt{\mid g \mid}} {\partial}_i \;\left(\sqrt{\mid g \mid} V^i \right)$$ (A.1)

and the Laplacian:

$$% \Delta \Phi = \frac{1}{\sqrt{\mid g \mid}} {\partial}_i \left( \sqrt{\mid g \mid} g^{in} {\partial}_n \Phi \right)$$ (A.2)

where $\vec{V} = V^1 \vec{E}_1 + V^2 \vec{E}_2 + V^3 \vec{E}_3 = V_1 \vec{E}^1 + V_2 \vec{E}^2 + V_3 \vec{E}^3$ is written in the natural basis $\vec{E}_i = \partial \vec{OM} / \partial x^i$ or the conjugated basis $\vec{E}^i$ verifying $\vec{E}_i \cdot \vec{E}^j = \delta_{ij}$, gⁱⁿ are the components of the metric tensor and $\mid$g$\mid$ is the absolute value of metric tensor determinant.

From these expressions, we derived the form of the following operators:

$$% \vec{g}_0 \cdot \vec{\nabla} \equiv (g^{1}_0 \vec{E}_1 + g^... ...^i \equiv g_0^1 {\partial}_{\zeta} + g_0^2 {\partial}_{\theta}$$ (A.3)

$$% r^2 \Delta \equiv h_1 {\partial}_{\zeta \zeta}^2 - 2 h_2 {\... ...theta \zeta}^2 + h_4 {\partial}_{\zeta} + \Delta_{\theta \phi}$$ (A.4)

$$% r^2 {\cal L} \equiv - \;\left(h_1 \frac{g_0^2}{g_0^1} + 2 h... ...0^1} \right) \right]{\partial}_{\theta} + \Delta_{\theta \phi}$$ (A.5)

$$% r^2 \vec{\nabla} \cdot (\; \bullet \; \vec{A_0}) \equiv (r^... ...rtial}_{\theta} + r^2 \vec{\nabla} \cdot \vec{A_0} [\mbox{id}]$$ (A.6)

where $\Delta_{\theta \phi}$ represent the horizontal part of the Laplacian in spherical coordinates:

$$% \Delta_{\theta \phi} \equiv {\partial}_{\theta \theta}^2 + ... ...\theta} + \frac{1}{\sin^2\!\theta} {\partial}_{\phi \phi}^{2},$$ (A.7)

and

$$% h_1 = \frac{r^2 + r_{\theta}^2}{r^2_{\zeta}}$$ (A.8)

$$% h_2 = \frac{r_{\theta}}{r_{\zeta}}$$ (A.9)

$$% h_3 = \frac{r}{r_{\zeta}}$$ (A.10)

$$% h_4 = \frac{1}{r_{\zeta}} \left[{\partial}_{\zeta} \!\! \;\... ...n \theta} {\partial}_{\theta}(r_{\theta} \sin \theta) \right].$$ (A.11)

We recall that:

$$% r = (1-\epsilon) \zeta + A(\zeta) \;\left(S(\theta) - 1+\epsilon \right)$$ (A.12)

where $S(\theta)$ describes the stellar surface.

Appendix B: Coupling matrix

The components of the sub-matrices which define the ODE system (36) are specified below using the functionals  ${\mbox{I}}^m_{ \ell \ell '}$ and  ${\mbox{J}}^m_{ \ell \ell '}$ defined in Eqs. (37) and (38):

$$% \mbox{A}_{33} \qquad {\mbox{I}}^m_{ \ell \ell '}(h_1)$$ (B.1)

$$% \mbox{B}_{11} \qquad {\mbox{I}}^m_{ \ell \ell '}(r^2 g_0^1)$$ (B.2)

$$% \mbox{B}_{21} \qquad {\mbox{I}}^m_{ \ell \ell '}(h_4) - {\m... ...{ \ell \ell '} \;\left(2 h_2 + h_1 \frac{g_0^2}{g_0^1} \right)$$ (B.3)

$$% \mbox{B}_{22} \qquad {\mbox{I}}^m_{ \ell \ell '} \;\left(h_1 \frac{c_0^2 N^2_0}{g_0^1} - r^2 A_0^1 \right)$$ (B.4)

$$% \mbox{B}_{33} \qquad {\mbox{I}}^m_{ \ell \ell '}(h_4) - {\mbox{J}}^m_{ \ell \ell '}(2 h_2)$$ (B.5)

$$% \mbox{C}_{11} \qquad {\mbox{J}}^m_{ \ell \ell '}(r^2 g_0^2)$$ (B.6)

$$% \mbox{C}_{12} \qquad -{\mbox{I}}^m_{ \ell \ell '}(r^2 c_0^2 N^2_0)$$ (B.7)

$$% \mbox{C}_{21} \qquad - \ell(\ell+1)\delta_{ \ell \ell '} - ... ...tial}_{\zeta} \!\! \;\left(\frac{g_0^2}{g_0^1} \right) \right]$$ (B.8)

$$% \mbox{C}_{22} \qquad {\mbox{I}}^m_{ \ell \ell '} \left[h_1 ... ...cdot \vec{A}_0 \right]- {\mbox{J}}^m_{ \ell \ell '}(r^2 A_0^2)$$ (B.9)

$$% \mbox{C}_{31} \qquad - {\mbox{I}}^m_{ \ell \ell '}(r^2 d_0)$$ (B.10)

$$% \mbox{C}_{32} \qquad {\mbox{I}}^m_{ \ell \ell '} \left[r^2 (d_0 c_0^2 - \rho_0) \right]$$ (B.11)

$$% \mbox{C}_{33} \qquad - \ell(\ell+1)\delta_{ \ell \ell '} + {\mbox{I}}^m_{ \ell \ell '}(r^2 d_0)$$ (B.12)

$$% \mbox{D}_{21} \qquad - {\mbox{I}}^m_{ \ell \ell '} \;\left(\frac{h_1}{g_0^1} \right)$$ (B.13)

$$% \mbox{D}_{22} \qquad {\mbox{I}}^m_{ \ell \ell '} \;\left(\frac{h_1 c_0^2}{g_0^1} \right)$$ (B.14)

$$% \mbox{E}_{11} \qquad {\mbox{I}}^m_{ \ell \ell '}(r^2)$$ (B.15)

$$% \mbox{E}_{12} \qquad -{\mbox{I}}^m_{ \ell \ell '}(r^2 c_0^2)$$ (B.16)

$$% \mbox{E}_{21} \qquad -{\mbox{I}}^m_{ \ell \ell '} \left[h_1 {\partial}_{\zeta} \!\! \;\left(\frac{1}{g_0^1} \right) \right]$$ (B.17)

$$% \mbox{E}_{22} \qquad {\mbox{I}}^m_{ \ell \ell '} \left[r^2 ... ...tial}_{\zeta} \!\! \;\left(\frac{c_0^2}{g_0^1} \right) \right]$$ (B.18)

where $\ell = m+ 2k$, $\ell' = m+ 2k'$ when applied to the $\vec{\Xi}_{m}^{+}$ vector and $\ell = m+ 2k+1$, $\ell' = m+ 2k'+1$for  $\vec{\Xi}_{m}^{-}$.

For a polytropic model of index N, the quantities describing the equilibrium can be expressed in terms of the dimensionless enthalpy H as follows:

$$% \begin{array}{ll} \vec{g}_0 = \vec{\nabla} H & \vec{A}_0 = ... ...\rho_0 = \Lambda^N H^N & d_0 = N \Lambda^N H^{N-1}, \end{array}$$ (B.19)

where $\Lambda$ is such that

$$% \Lambda = \frac{4 \pi G \rho_{\rm c} R_{\rm e}^2}{h_{\rm c}}$$ (B.20)

where $h_{\rm c}$ and $\rho_{\rm c}$ are the dimensional enthalpy and density at the center of the polytropic model.

The components of the ODE, given by Eqs. (B.1) to (B.18), can then be expressed in terms of the enthalpy and its derivatives, $H_{\zeta}, H_{\theta}, H_{\zeta \theta}, H_{\zeta \zeta} H_{\theta \theta}$. This has been done in order to minimize the numerical error in the calculation of these components. The most useful expressions are:

$$% g_0^1 = \frac{h_1 H_{\zeta} - h_2 H_{\theta}}{r^2}$$ (B.21)

$$% g_0^2 = \frac{- h_2 H_{\zeta} + H_{\theta}}{r^2}$$ (B.22)

$$% {\Vert \vec{\nabla} H\Vert}^2 = \frac{h_1 H_{\zeta}^2 -2 h_2 H_{\zeta} H_{\theta} + H_{\theta}^2}{r^2}$$ (B.23)

$$% h_1 \frac{c_0^2 N_0^2}{g_0^1} - r^2 A_0^1 = \;\left(1 - \fr... ...ma}_{1,0}}{N+1} \right) \frac{H_{\theta}^2}{r_{\zeta}^2 g_0^1}$$ (B.24)

 

$$% \partial_{\zeta} \!\! \;\left(\frac{g_0^2}{g_0^1} \right) \frac{1}{ \;\left(r^2 g_0^1 \right)^2} \left[h_3^2 \;\left(H_{\ze... ...eta \zeta} H_{\theta} \right)+ \;\left(h_2 \partial_{\zeta} h_1 \right. \right.$$ =

$$\left. \left. -~ h_1 \partial_{\zeta} h_2 \right)H_{\zeta}^2 - \p... ...al_{\zeta} h_1 H_{\zeta} H_{\theta} + \partial_{\zeta} h_2 H_{\theta}^2 \right]$$ (B.25)

$$% h_1 {\partial}_{\zeta} \!\! \;\left(\frac{c_0^2 N^2_0}{g_0^... ...\;\left(1 - \frac{N{\Gamma}_{1,0}}{N+1} \right)r^2 {\cal L}(H)$$ (B.26)

 

$$% r^2 {\cal L}(H) -h_1 \left[\frac{g_0^2}{g_0^1} H_{\theta \zeta} + \partial_{\zeta... ...! \;\left(\frac{g_0^2}{g_0^1} \right)H_{\theta} \right] -2 h_2 H_{\theta \zeta}$$ =

$$+~ h_4 H_{\zeta} + \Delta_{\theta \phi} H$$ (B.27)

$$% d_0 c_0^2 - \rho_0 = - \;\left(1 - \frac{N{\Gamma}_{1,0}}{N+1} \right)\Lambda^N H^N$$ (B.28)

 

$$% \partial_{\zeta} \!\! \;\left(\frac{1}{g_0^1} \right) - \frac{r^2}{ \;\left(r^2 g_0^1 \right)^2} \left[h_1 H_{\zeta \zeta} -h_2 H_{\zeta \theta} + \;\left(\partial_{\zeta} h_1 \right. \right.$$ =

$$\left. \left. - ~2 h_1 / h_3 \right)H_{\zeta} + \;\left(2 h_2 / h_3 - \partial_{\zeta} h_2 \right)H_{\theta} \right]$$ (B.29)

 

$$% \partial_{\zeta} \!\! \;\left(\frac{c_0^2}{g_0^1} \right)= \fra... ...l_{\zeta} \!\! \;\left(\frac{1}{g_0^1} \right)+ \frac{H_{\zeta}}{g_0^1} \right]$$ (B.30)

where

$$% \partial_{\zeta} h_1 = 2 \;\left(\frac{r_{\theta \zeta}}{r_... ...a}} h_2 - \frac{r_{ \zeta \zeta}}{r_{\zeta}} h_1 + h_3 \right)$$ (B.31)

$$% \partial_{\zeta} h_2 = \frac{r_{\theta \zeta}}{r_{\zeta}} - \frac{r_{ \zeta \zeta}}{r_{\zeta}} h_2$$ (B.32)

Appendix C: Calculation of the disk-integration factor

According to the definition of the disk-integration factor, Eq. (43), we are led to calculate integrals of the following form:

 

$$\int\!\!\!\int_{S_v} F(\theta,\phi) {\rm d} {\bf S} \cdot {\bf e}_i$$                     I = (C.1)

$$\int\!\!\!\int_{S_v} G(\theta,\phi,i) {\rm d}\mu {\rm d}\phi$$ = (C.2)

where $\mu = \cos~(\theta)$ and $F(\theta,\phi)$ is the surface distribution of the eigenfunction obtained in the coordinate system (25) in which the polar axis is the rotation axis. The integral is most simply calculated in the coordinate system in which the polar axis is aligned with the direction of the observer. This coordinate system results from a rotation of angle i around the y axis of the original coordinate system, the new angular variables being denoted $\theta'$ and $\phi'$. To express G in these coordinates, we use the formula relating the spherical harmonics in both systems:

$$% Y^m_\ell (\theta,\phi) = \sum_{m' =-\ell}^{+\ell} d_{m m'}^{\ell}(i) Y^{m'}_\ell (\theta ',\phi ')$$ (C.3)

where $d_{m m'}^{\ell}(i)$ do not generally have a simple form (Edmonds 1960). Then, using the spherical harmonic expansion of G, we obtain:

 

$$\sum_{\ell=0}^{+\infty}\sum_{m=-\ell}^{+\ell} G^\ell_m(i) Y^m_\ell (\theta,\phi)$$                          G = (C.4)

$$\sum_{\ell=0}^{+\infty}\sum_{m=-\ell}^{+\ell} \sum_{m'=-\ell}^{+\ell} G^\ell_m(i) d_{m m'}^{\ell}(i) Y^{m'}_\ell (\theta',\phi').$$ = (C.5)

Then, integrating over the longitude $\phi'$, from 0 to $2 \pi$, the terms involving $Y^{m'}_\ell (\theta',\phi')$ vanish if $m' \neq 0$. It follows that

$$% I = 2 \pi \sum_{\ell=0}^{+\infty}\sum_{m=-\ell}^{+\ell} J_{\ell} G^\ell_m(i) \hat{Y}^m_\ell(i)$$ (C.6)

where we used the following relations,

$$% d_{m 0}^{\ell}(i) = \sqrt{\frac{4 \pi}{2 \ell + 1}} \hat{Y}^m_\ell(i)$$ (C.7)

$$% {\rm d}\mu {\rm d}\phi = {\rm d}\mu' {\rm d}\phi' \;\; \mbox{where $\mu' = \cos \theta'$ }$$ (C.8)

and defined $J_{\ell}$ as,

 

$$% J_{\ell} \sqrt{\frac{4 \pi}{2\ell+1}} \int_0^{1} \hat{Y}^0_\ell(\mu') {\rm d}\mu'$$ = (C.9)

$$\left\{ \begin{array}{lll} 0 & \mbox{if $\ell$\space is even and ... ... and $\ell \neq 1$ }, \\ \frac{1}{2} & \mbox{if $\ell=1$ }.\end{array} \right.$$ = (C.10)

Because of axial axial symmetry, the function to integrate reads

$$% F(\theta,\phi) = W(\theta) {\rm e}^{{\rm i}m\phi}.$$ (C.11)

Then, from the expression of the vector  ${\rm d}\vec{S}$ at the star's surface:

$$% {\rm d}\vec{S} = \partial_{\theta} \vec{OM} \times \partial... ...ta {\rm d}\phi = \sqrt{g} \vec{E}^1 {\rm d}\theta {\rm d}\phi$$ (C.12)

we deduce that

$$% G = r A(\theta,\phi,i) W(\theta) {\rm e}^{{\rm i}m\phi} \\$$ (C.13)

where

 

$$% A(\theta,\phi,i) r r_{\zeta} \vec{E^1} \cdot \vec{e}_i = \sqrt{r^2 + r^2_{\theta}} \vec{e^s} \cdot \vec{e}_i$$ = (C.14)

$$r \;\left(\sin\theta\cos\phi\sin i + \cos\theta \cos i \right)$$ =

$$+ r_{\theta} \;\left(\sin\theta \cos i - \cos \theta \cos\phi \sin i \right)$$ (C.15)

$$\cos i \frac{{\rm d}}{{\rm d} \theta} \;\left(r \sin\theta \right)-\sin i \cos \phi \frac{{\rm d}}{{\rm d} \theta} \;\left(r \cos \theta \right)$$ = (C.16)

where $\vec{e^s}$ denotes the unit vector perpendicular to the surface, r and  $r_{\theta}$ are calculated at the star surface $\zeta=1$. Thus the dependency of G on i, $\phi$ and $\theta$ can be specified as follows:

 

$$A(\theta) \cos i e^{{\rm i}m\phi} - B(\theta) \sin i \cos \phi e^{{\rm i}m\phi}$$ (C.17)

where

 

$$A = r \frac{{\rm d}}{{\rm d} \theta} \;\left(r \sin\theta \right)W(\theta)$$ (C.18)

$$B = r \frac{{\rm d}}{{\rm d} \theta} \;\left(r \cos\theta \right)W(\theta).$$ (C.19)

It follows that

$$% G^\ell_k = 0 \;\; \mbox{if $k \neq m-1,m,m+1$ }$$ (C.20)

so that the integral now reads:

 

$$I/2\pi = I_{m-1}+ I_m + I_{m+1} \;\; \mbox{where}$$ (C.21)

$$I_m = \cos i \hat{A}_m(i)$$ (C.22)

$$I_{m-1} = - \frac{\sin i}{2} \hat{B}_{m-1}(i)$$ (C.23)

$$I_{m+1} = - \frac{\sin i}{2} \hat{B}_{m+1}(i)$$ (C.24)

where $\hat{A}_m$ denotes:

 

$$\hat{A}_m(i) = \sum_{\ell=\vert m\vert}^{+\infty} J_\ell A^\ell_m \hat{Y}^m_\ell(i)$$ (C.25)

$$A^\ell_m = 2\pi \int_0^{\pi} A(\theta) \hat{Y}^m_\ell(\theta) \sin \theta {\rm d} \theta$$ (C.26)

the $\hat{B}_m$ terms being defined accordingly.

Note that for modes which are equatorially anti-symmetric and axisymmetric (m=0), $\hat{A}_0(i)=J_0 A^0_0 \hat{Y}^0_0(i)$ and $\hat{A}_1(i) = \hat{A}_{-1}(i) = 0$, thus the integral I reduces to:

 

$$% I = 4 \pi \sqrt{\pi} A^0_0 \cos(i).$$ (C.27)