Issue 
A&A
Volume 547, November 2012



Article Number  A75  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201219548  
Published online  01 November 2012 
Pulsations of rapidly rotating stars
I. The ACOR numerical code
^{1} LESIA, UMR 8109, Université Pierre et Marie Curie, Université Denis Diderot, Observatoire de Paris, 92195 Meudon Cedex, France
email: rhitamaria.ouazzani@obspm.fr
^{2} Institut d’Astrophysique et de Géophysique de l’Université de Liège, Allée du 6 Août 17, 4000 Liège, Belgium
Received: 4 May 2012
Accepted: 21 September 2012
Context. Very high precision seismic space missions such as CoRoT and Kepler provide the means of testing the modeling of transport processes in stellar interiors. For some stars, such as solarlike and red giant stars, a rotational splitting is measured. However, to fully exploit these splittings and constrain the rotation profile, one needs to be able to calculate them accurately. For some other stars, such as δ Scuti and Be stars, for instance, the observed pulsation spectra are modified by rotation to such an extent that a perturbative treatment of the effects of rotation is no longer valid.
Aims. We present here a new twodimensional nonperturbative code called ACOR (adiabatic code of oscillation including rotation) that allows us to compute adiabatic nonradial pulsations of rotating stars without making any assumptions on the sphericity of the star, the fluid properties (i.e., baroclinicity) or the rotation profile.
Methods. The 2D nonperturbative calculations fully take into account the centrifugal distortion of the star and include the full influence of the Coriolis acceleration. The numerical method is based on a spectral approach for the angular part of the modes and a fourthorder finite differences approach for the radial part.
Results. We test and evaluate the accuracy of the calculations by comparing them with those coming from the TOP (twodimensional oscillation program) for the same polytropic models. We illustrate the effects of rapid rotation on stellar pulsations through the phenomenon of avoided crossings.
Conclusions. As shown by the comparison with the TOP for simple models, the code is stable, and gives accurate results up to nearcritical rotation rates.
Key words: asteroseismology / methods: numerical / stars: oscillations / stars: rotation
© ESO, 2012
1. Introduction
Rotation plays a key role in the evolution of stars across the HertzsprungRussell diagram, which shows the distribution of stars in the luminosity versus effective temperature plane. For instance, the centrifugal distortion breaks the thermal equilibrium and provokes largescale currents called meridional circulation (Eddington 1925). These currents as well as shear and baroclinic instabilities modify the angular momentum distribution (Zahn 1992; Mathis et al. 2004), and thereby the rotation profile inside the stars. Meanwhile, these processes transport chemical elements and change the evolution of the stars. That is the reason why determining rotational profiles inside stars is crucial for modeling stellar structure and evolution.
Asteroseismology is currently the only tool that allows such a determination. But rotation also changes stellar pulsations. The centrifugal force distorts the resonant cavity of the pulsations, and the Coriolis force modifies the modes’ dynamic. Usually, for slow rotators, the effects of rotation are taken into account as a perturbation of pulsations (see for instance Ledoux 1951, for firstorder effects, Gough & Thompson 1990; and Dziembowski & Goode 1992, for secondorder effects; and Soufi et al. 1998, for thirdorder effects). But these methods, although elegant and simple to use, have shown their limits (Reese et al. 2006; Suárez et al. 2010). On the one hand these perturbative methods approximate the effects of the centrifugal force on acoustic pulsations in stars that show very high surface velocities, such as δ Scuti and Be stars. On the other hand, they approximate the impact of the Coriolis force on gravity modes in stars whose surface rotates slowly, but in which the pulsation periods are of the same order as their rotation period (P_{rot} ∼ P_{osc}), such as SPB stars or γ Doradus. Finally, perturbative methods may also fail in modeling pulsations of cooler stars, such as subgiant or red giant stars, which have very low surface velocities but rapidly rotating cores, and in which all pulsation modes are of a mixed nature, i.e. gravity close to the core and acoustic in the envelope (see for example Beck et al. 2012).
This concerns many stars in the CoRoT and Kepler observation fields. Hence, if one aims at correctly extract the rotation profile from seismic observations, they need to correctly apprehend the effects of rotation on pulsations.
This work aims at presenting a model that accurately takes into account the nonperturbative effects of rotation on oscillation spectra. A new two dimensional nonperturbative code is presented. The 2D nonperturbative calculations fully take into account the centrifugal distortion of the star and include the full influence of the Coriolis acceleration. This 2D nonperturbative code is useful for studying pulsational spectra of highly distorted evolved models of stars, as well as stars presenting highly differential rotation profiles. Section 2 introduces the basic pulsation equations in spheroidal geometry using a coordinate system adapted to the star’s shape. Section 3 is dedicated to the numerical method, which is based on a spectral approach for the angular part of the modes, and on a finite difference method for the radial part. In Sect. 4, we test these calculations and evaluate their accuracy. Finally, in Sect. 5, the results are validated by comparing them with those of Reese et al. (2006), and an example of application is given in Sect. 6.
2. Basic equations in spheroidal geometry
When dealing with the subject of computing the pulsations of rotating stars, one has to face two main problems. Firstly, rotation, through the centrifugal force, distorts the resonant cavity of the pulsations. If a solenoidal rotation profile is assumed (around a northsouth axis), the azimuthal symmetry is conserved, whereas the spherical one is broken: the star takes on a spheroidal geometry. This centrifugal distortion, if it is strong enough, has to be treated by a twodimensional approach. Secondly, when rotation is accounted for in the pulsation equations, the Coriolis acceleration enters the momentum equation and modifies the dynamics of the modes. If weak enough, this Coriolis effect can be approximated by perturbative methods. But for moderate to high rotational velocities, as well as for high order g modes when P_{rot} ∼ P_{osc}, the perturbative treatment is no longer relevant, and a nonperturbative approach has to be adopted.
2.1. Spheroidal geometry
Because of the distorted shape of a rotating star, we chose a new coordinate system that defines the star’s surface at a constant pseudoradial coordinate. To do so, we adopted a multidomain approach, which consists in dividing the physical space into domains whose boundaries correspond to the model’s discontinuities (e.g. Canuto et al. 1988): one domain V_{1}, which encloses the star, and one external domain, V_{2}, bounded by the stellar surface and a sphere of twice the equatorial radius. Following Bonazzola et al. (1998), we introduced a coordinate system that extends from spherical symmetry at the center to a spheroidal shape at the stellar surface, and back to a spherical geometry at the external boundary of V_{2}. In this system, the radial coordinate, ζ, is no longer r, the distance to the center. However, for a fixed colatitude θ, r and ζ are related thanks to the following continuous and bijective function: In domain V_{1}: (1)where ζ ranges from 0 to 1, ϵ = 1 − R_{pol} / R_{eq} is the flatness, and r(ζ = 1,θ) = R_{s}(θ) the stellar surface. In domain V_{2}: (2)where ζ ranges from 1 to 2, ζ = 2 corresponding to a spherical surface that encompasses the star (see Fig. 1). The use of such a coordinate system helps significantly with establishing the boundary conditions.
Then one has to define a new vector basis. We chose the following orthogonal spheroidal basis: (3)where

is tangential to isoζ surfaces in the meridional plane;

is directly orthogonal to in the meridional plane;

,
and where r_{θ} = ∂_{θ}r. In a forthcoming paper, this method will be applied to a realistic stellar model using three domains: one that encloses the convective core, the second the radiative envelope, and the third an external domain, which allows us to avoid discontinuities at the top of the convective zone.
Fig. 1 Coordinate system used in ACOR. V_{1} extends from the center to the stellar surface, and V_{2} encompasses the star. 
2.2. Basic equations
Here, we consider the stellar interior to be an inviscid selfgravitating fluid, where the magnetic field is neglected. Therefore its physics is governed by the conservation of mass, momentum and energy, the energy transfer equation, and Poisson’s equation for the gravitational potential (see for example Kippenhahn & Weigert 1994).
Using the Eulerian formalism, we compute the oscillation modes as the adiabatic response of the structure to small perturbations – i.e., of the density, pressure, gravitational potential and velocity field. As in Unno et al. (1989), we perturb the hydrodynamics equations around the equilibrium state. Considering that the velocity field in the equilibrium state is only due to rotation, the linearized equations in the inertial frame are given by where primed quantities denote Eulerian perturbations, whereas quantities with the subscript “0” correspond to the stationary state. The symbol e_{i} stands for the spherical basis vectors. Note that the energy conservation equation was replaced by the adiabatic relation. Given that the equilibrium state is stationary and axisymmetric, the time and azimuthal dependences are of the form e^{i(ωt+mϕ)}, where ω is the pulsation frequency, and m the azimuthal order of the pulsation mode.
We put the system of equations in nondimensional form using the following transformations: where Ω_{k} stands for the Keplerian critical velocity, i.e., the rotation velocity at which the centrifugal force compensates gravity at the equator.
From now on, we work in terms of nondimensional variables and drop the tilded notation.
Rather than projecting the motion equation onto the basis vectors, we chose to decompose it into one radial, one poloidal, and one toroidal field. This decomposition allows us to obtain separate partial differential equations for the radial and angular coordinates, and helps to reduce the number of unknowns, as will be shown below.
Moreover, we apply the change of variable π^{′} = p^{′} / ρ to avoid singularity problems at the surface of polytropic models, and we introduce an auxiliary variable dΦ^{′}, as well as the equation relating it to Φ^{′}: dΦ^{′} = ∂Φ / ∂ζ = ∂_{ζ}Φ^{′}, thereby reducing Poisson’s equation from a secondorder differential equation to two firstorder ones.
2.3. Spectral expansion
Thanks to an appropriate choice of the coordinate system, the behavior of the eigenfunctions on isoζ surfaces is smooth. Therefore, the angular behavior of pulsation modes is well described in terms of an expansion on the basis of spherical harmonics. This spectral method is known to be very wellsuited in fluid dynamics. For instance, for a function of class , the numerical error decreases as e^{−aM}, where M is the number of spherical harmonics in the expansion and a a constant (Canuto et al. 1988).
Therefore, we perform a spectral expansion of the unknowns in terms of the spherical harmonics (Rieutord 1987). Any vector field can be decomposed into a radial, a poloidal, and a toroidal part. Therefore, the components of the pulsation velocity field are expressed as follows in the vector basis (a_{ζ},a_{θ},e_{ϕ}) (8)All other scalar unknowns, namely Φ^{′}, dΦ^{′}, and ρ^{′}, are expanded in the same way as the scaled pressure perturbation: (9)
2.4. Symmetries and mode classification
Because of the symmetries of the equilibrium model with respect to the rotation axis and to the equator, one obtains a separate eigenvalue problem for each value of the azimuthal order, m, and each parity, par, with respect to the equator. Thus, for a given value of m, there are two independent sets of ordinary differential equations (ODE) coupling the spectral coefficients, one with ℓ of the same parity as m, and the other with opposite parities. This means that for a given m, when including M terms in the spectral expansion, ∀j ∈ [1,M] ,with par = 0 if ℓ is of the same parity as m (even mode), and par = 1 otherwise (odd mode). We then obtain a system of ODE of the variable ζ for the coefficients of the spherical harmonic expansion , and .
2.5. Projections
After expanding the unknowns on the basis of spherical harmonics, the second step consists of projecting the equation system onto the spherical harmonics basis, which is truncated to M terms, as is the spectral expansion. In this subsection, the unknowns are generically designated by X_{i}(ζ,θ,ϕ). Each partial differential equation of the form (12)is replaced by a system of M differential equations in ζ, obtained by projecting these equations on a basis of M spherical harmonics: (13)in which ℓ_{1} and ℓ_{2} take on the values defined in Eqs. (10) or (11).
The equilibrium quantities, which are expressed as functions of ζ and θ, are implicitly contained in the operator E.
The impact of the basis dimension (M) on the precision of the computations will be discussed in detail in Sect. 4.1.
2.6. Boundary conditions
To complete the eigenvalue problem, it is necessary to specify a number of boundary conditions. The system of equations contains four sets of firstorder ODEs. Acordingly, we impose four boundary conditions. As the system is solved simultaneously for all layers, two boundary conditions are imposed at the center, one at the stellar surface and one on the external spherical surface of V_{2} (see Fig. 1).
Taking boundary conditions at the center is a delicate problem because of the coordinate system. It consists of imposing the regularity of the solutions at the center. To do so, we take the limits of the equations as ζ goes to zero. The different scalar unknowns behave as follows and the components of the velocity field obey This results in two algebraic relations between the unknowns. A detailed explanation of these central boundary conditions is presented in Appendix B.
At the stellar surface, a stressfree condition is imposed: δp = 0. The stellar surface is assumed to be on an isobar, at ζ = 1, thus the boundary condition corresponds to (14)where the subscript “0” stands for equilibrium quantities.
The external condition for the gravitational potential consists in imposing that it vanishes at infinity. This can be achieved by matching the gravitational potential to a vacuum potential at ζ = 2, i.e., on the spherical external surface V_{2}. The advantage in using a spherical boundary is that the Laplace equation is separable for solutions of the form . This gives very simple conditions for a continuous match at the ζ = 2 spherical boundary. Specifically, for ζ ≥ 2, this equation is decomposed over the spherical harmonic basis, and each component solved separately. For each ℓ, this yields two independent solutions: , which diverges at infinity and is therefore discarded, and , which vanishes. Hence, the corresponding boundary condition is (15)where r_{ζ} = ∂_{ζ}r.
3. Numerical method
The spectral form of Eqs. (A.1), (A.6), (A.9), (A.12), (A.13), (A.15), and (A.16) constitute a firstorder ordinary differential system of 7 × M equations, in terms of the coordinate ζ. Of these, 4 × M are ordinary differential equations for the spectral coefficients u_{ℓ}, , and d (i.e., Eqs. (A.1), (A.13), (A.15), and (A.16)), and the remaining 3 × M equations are algebraic equations for the spectral coefficients v_{ℓ}, w_{ℓp} and (i.e., Eqs. (A.6), (A.9), and (A.12)).
We chose to solve this system by an Newtonlike method, which consists of taking a guess at the pulsation frequency σ_{0} and looking for the pulsation mode with the closest frequency to this guess. Solving the system yields a deviation, δσ, from the initial guess, σ = σ_{0} + δσ is taken as a new guess for the next iteration, and this process is iterated till convergence (three steps are generally enough).
We therefore isolate the terms proportional to δσ, and obtain the following system of equations, which can be put under the form of a matrix: where y_{1} and z_{1} are column vectors containing the spectral coefficients of the unknowns (18)and where ℓ =  m  + 2(j − 1) + par, ∀j ∈ [1,M] .
A_{11}, A_{12}, A_{21} and A_{22} correspond to the following equations (19)whereas B_{11}, B_{12}, B_{21} and B_{12} correspond to the algebraic equations (20)Thanks to the three last equations, the nondifferentiated unknowns can be expressed in terms of the differentiated ones with the help of Eq. (17). To do so, the matrix (B_{21} + δσ B_{22}), which is a real square matrix of rank 3M, has to be inverted. It is straightforward to show that if δσ is small enough, (21)Thus, the matrix system can be written in the form:
3.1. Radial discretization
In the radial direction, structural quantities, and as a consequence eigenmodes, undergo sharp variations (for instance at the top of a convective core, or at the star’s surface). Therefore, in the radial direction, spectral methods are inappropriate when dealing with realistic stellar models. We chose to discretize the differential equations using a finite differences approach. The continuous domain of integration is replaced by a discrete set of N_{r} radial points. The number of points in the radial grid is an important factor that affects the numerical precision, as discussed in detail in Sect. 4.1. Another characteristic of the discretization scheme, which comes into play in the numerical precision, is the estimation of the derivatives. With classical finite differences methods, the more precise you get, the less stable are the computations.
We adopted a fourthorder difference scheme, which relies on the following identity (24)where the primes denote derivatives with respect to ζ, h is given by h = ζ_{i + 1} − ζ_{i}, and the subscripts “i” and “i + 1” denote the layer indexes. The great advantage of this scheme is that it is accurate up to h^{4}, while retaining high numerical stability, as it only involves two consecutive grid points. This finite differences scheme has already been implemented by Scuflaire et al. (2008) in the Liège Oscillation Code, which has been proven to be very stable and accurate for the calculations of all types of pulsation modes in all kinds of stars.
From now on, we drop the subscripts “1” from y_{1}. Thanks to Eq. (22), it is possible to express the derivatives of y. The identity (24) can then be valid as long as the matrix coefficients in A and A_{δσ} are continuous and have continuous derivatives. We then obtain the following matrix equation at each layer i(25)where Here, I_{d} is the square 4M × 4M identity matrix. The system of equations over the whole stellar interior can then be written in the canonical form: (26)where the vector has been introduced: (27) and being block diagonal matrices.
We impose boundary conditions that can be expressed as algebraic relations, as explained in detail in Appendix B.
We note that Eq. (26) is a generalization of the classical eigenvalue problem Ay = λy.
3.2. Inverse iteration algorithm
To calculate the eigenvectors, , and eigenvalues, δσ, we generalized the classical inverse iteration algorithm to the more general problem formulated in Eq. (26), as developed by Dupret (2001) in the stellar pulsations context.
The estimate at step k + 1 of the eigenvector, , is obtained from the estimate at step k and the initial guess, δσ_{0}, of the eigenvalue by the formula (28)If we assume that is inversible and that is diagonalizable, it is easy to prove that this algorithm has the same geometrical convergence rate as the classical inverse iteration algorithm. The inverse matrix is not explicitly calculated but we solve the system: (29)To do so, we perform a LU factorization of the system with partial pivoting. Then, we iterate solving the two triangular systems at each step. Note that if the initial guess for the eigenvalue is good, the algorithm converges quickly toward the solution even with a poor eigenvector as an initial estimate.
The corresponding eigenvalue is then determined by the generalization of the Rayleigh ratio (30)where and are the Hermitian conjugates of and , respectively. It can be easily shown that δσ, given by Eq. (30), minimizes: (31)
4. Tests and accuracy
As mentioned in the introduction, to be integrated into a stellar modeling chain for massive computations, the program has been developed with a constant concern for simplicity and rapidity of use. In this section, we assess the role of numerical parameters in the convergence process toward an oscillation mode and establish the computational performances with respect to these parameters.
What takes up the most computing resources in ACOR are the integrations over θ of the equation coefficients (Eq. (13)). These coefficients are evaluated at each radial layer (for N_{r} layers) by projecting the equations onto the spherical harmonics basis (of dimension M). Therefore, by assessing the role of the two parameters M and N_{r}, it is possible to define the optimal values for a good compromise between computation time, memory resources and accuracy of the results.
Note that there is no automatic method that allows us to identify the modes. In this work, we followed the progression of the modes, as we gradually increased the rotation rate from zero to a high value. We used the kinetic energy distribution in the meridional plane, such as in Fig. 7, to correctly select the mode at each step.
All tests presented in this paper were made assuming uniformly rotating polytropes of polytropic index N = 3 (polytropic exponent γ = 4 / 3) with an adiabatic index of Γ_{1} = 5 / 3.
4.1. Convergence tests
In the ideal case where equilibrium quantities would be , the numerical errors would decrease as e^{−aM}. Rotation induces nonspherical profiles for the equilibrium quantities and causes the eigenfunctions to depart from a single spherical harmonic. Therefore, the higher the rotation rate, the stronger the deviations from sphericity, and the larger the spherical harmonic basis has to be, as illustrated indeed in Fig. 2. The convergence calculations illustrated in Fig. 2 show that convergence is reached for 7 terms for Ω = 18.5%Ω_{k}, 16 terms for Ω = 37.9%Ω_{k}, and 25 terms for Ω = 58.9%Ω_{k} in the spectral expansion. This figure also shows that the convergence reaches machine precision.
Fig. 2 Convergence as a function of the number of spherical harmonics, taken as the relative frequency difference between computations using consecutive spectral resolutions, M − 1 and M, for an n = 3 mode dominated by an ℓ = 1, m = 0 component, and for three different rotational velocities: 18.5% Ω_{k} in blue, 38% Ω_{k} in orange, and 59% Ω_{k} in red. 
Fig. 3 Convergence as a function of the radial resolution, taken as the relative frequency differences between computations using consecutive radial resolutions, r and (r + 1) points, at three different rotational velocities: 18.5% Ω_{k} in green, 37.9% Ω_{k} in light blue, and 58.9% Ω_{k} in dark blue. In both panels, the mode is dominated by the ℓ = 2, m = 0 component. The upper panel corresponds to an n = 3 mode and uses irregular grids, whereas the lower panel shows an n = 9 mode calculated with even distributions of points. 
Numerical resources – i.e., time in minutes and memory in MB or GB – used by ACOR with a spectral resolution M and a radial one N_{r}.
Concerning the convergence with respect of the radial resolution, we present in Fig. 3 the worst (bottom) and the best case (top). The higher the radial order n, the more numerous the nodes in the eigenfunction and the higher the radial resolution has to be. These plots also show that using a nonregular radial grid, whose number of points increases in an outward direction, allows us to increase the accuracy of computations. This is because the computed modes are acoustic modes, with consequently a short wavelength in the outer layers.
4.2. Numerical resources
Of all operations performed by the code, calculating the projection integrals (Eq. (13)) is the most demanding in terms of computational time, whereas inverting the two matrices and requires the most memory. In Table 1 we indicate the memory and time resources needed by the computations for the parameters M and N_{r}. Note that and are bloc matrices, their size corresponds to the number of nonzero elements they contain.
5. Comparison with Reese et al. (2006)
After the numerical tests presented in Sect. 4, the aim of this section is to validate ACOR’s results by comparing them with those from the TOP code. The TOP code has been developed by Reese et al. (2006) for two dimensional polytropes as a first step. The approach is based on a twodimensional spectral method that uses Tchebichev polynomials for the radial dependence. We present here the comparison between TOP and ACOR results for identical polytropic models.
Roughly, a polytrope is a simplified model for which one assumes an ad hoc relation between density (ρ_{0}) and pressure (p_{0}): (32)where K is the polytropic constant and γ is called the polytropic exponent, which is taken to be 4 / 3 here. The detailed method used to compute the rotating polytropic models is given in Rieutord et al. (2005). This method is an iterative scheme based on a spectral decomposition using Tchebichev polynomials for the radial dependence, and spherical harmonics for the horizontal one. We subsequently interpolate these models onto a radial grid that is appropriate for finite differences.
Concerning the parameters used for the calculations within the two codes, the angular spectral resolutions were fixed to reach convergence. It therefore depends on the rotation velocity and varies from 10 terms in the harmonics expansion for low rotation rates to 25 for the highest ones. For TOP, the spectral resolution in terms of the Tchebichev polynomials varies from 50 to 80 terms. For ACOR, the radial resolution is fixed at N_{r} = 2000 grid points.
5.1. Eigenfunctions comparison
The major effect of centrifugal distortion is the loss of spherical symmetry, which results in the coupling of spherical harmonics of different degrees to describe the horizontal behavior of a single mode. In Appendix C we provide the contributions of dominant spherical harmonics in the spectral expansion of two modes: an odd mode dominated by an ℓ = 1, m = 0 component (see Fig. C.1) and an even one dominated by ℓ = 2, m = 0 (see Fig. C.2). The solid lines correspond to the calculations made with ACOR and the dotted lines to those made with TOP. These plots clearly show that the evolution of the angular behavior of the modes with respect to rotation obtained by the two codes is very similar. This allows us to validate the analytical calculations as well as the inverse iteration algorithm, which converges onto eigenfunctions, while eigenfrequencies are computed a posteriori through the minimization of Eq. (31).
5.2. Eigenfrequencies comparison
Concerning the comparisons of eigenfrequencies, Fig. 4 shows the frequency differences between the two codes for odd and even eigenmodes in the low (top) and high (bottom) frequency regimes. Generally, the discrepancies between results from the two codes are of the order of 0.0001% − 0.08% (3 × 10^{3} Ω_{k} at most) for p modes, and 0.5% for g modes (1 × 10^{2} Ω_{k} at most), which seems very satisfying considering that the two codes rely on different and independent computations, and in particular on a different treatment for the central boundary conditions. Once more, this makes numerical programing mistakes unlikely in our calculations.
Fig. 4 Discrepancies between eigenfrequencies computed by ACOR and TOP as a function of the rotation angular velocity for five axisymmetric modes: top: for two n = 3 modes, corresponding to ℓ = 1 (light blue) and ℓ = 2 (orange) in the nonrotating case and one n = 1 g mode with ℓ = 2; bottom: n = 9 modes, corresponding to ℓ = 1 (dark blue) and ℓ = 2 (red). 
Fig. 5 Behavior of the eigenfrequencies computed by ACOR and TOP with respect of the rotational angular velocity for two multiplets; top: centroid modes with dominant component ℓ = 1, m = 0, and ℓ = 2, m = 0; bottom: sectorial and tesseral modes, dominated by components ℓ = 1, m = ± 1 and ℓ = 2, m = ± 1 and ± 2. 
One of the observational characteristics of rotation in seismology is the rotational generalized splitting, i.e., the frequency difference between two modes with the same radial order and harmonic degree, but with opposite azimuthal orders. Figure 5 shows the evolution of the structure of two multiplets with respect to the rotation velocity. The plots in the bottom panel show that the two codes find the same structure for the multiplets, regardless of the rotation rate (from 0 to 60% Ω_{k}) and the symmetry class of modes. The impact of rotation on centroid modes (top panel) is also the same with the two methods.
The agreement between the two oscillation programs developed separately, not only on the eigenfunctions and on frequencies, but also on the structure of the spectra, allows us to validate the approach adopted by ACOR.
Fig. 6 Avoided crossing, illustrated by plotting the frequency with respect to the rotation rate for two n = 2 modes (referred to as mode #1 in red and mode #2 in blue), which, in the nonrotating case, are ℓ = 1 and ℓ = 5. The harmonic degree given in the graph is the degree of the dominant component for each mode. 
Fig. 7 Spatial distribution of the pressure perturbation in a meridional plane for two low frequency modes (radial order n = 2, azimuthal order m = 0) involved in an avoided crossing occurring around Ω = 40%Ω_{k}, as modeled by ACOR. Left: at low rotation velocities, mode #1 is dominated by an ℓ = 1 component and changes nature to become dominated by ℓ = 5 as the velocity increases. Right: at low rotation velocities, mode #2 is dominated by an ℓ = 5 component and changes nature to become dominated by ℓ = 1 as the velocity increases. 
6. Illustration: avoided crossing
In quantum mechanics, an avoided crossing (or level repulsion) occurs, for instance, in a two level system (  1 ⟩ and  2 ⟩ ), placed in a magnetic field that acts differently on the two levels (see for example CohenTannoudji et al. 1973). When the two states are coupled, the levels repulse each other, since the system’s energy cannot be degenerated.
Accordingly, in a rotating star, pulsation mode frequencies tend to cross because the modes are not affected the same way by rotation (particularly by the centrifugal force, see Lignières et al. 2006). As two modes of the same parity cannot have the same frequency, an avoided crossing occurs during which the two modes exchange their characteristics.
This is illustrated in Fig. 6 by the evolution of the frequencies of two coupled modes with respect to the rotation rate. The two modes are of the same symmetry with m = 0, and ℓ = 1 or 5. As explained above, their frequencies cannot be degenerated, therefore the crossing is avoided, and as shown in Fig. 7, they progressively exchange angular characteristics. With the help of the distribution of the pressure perturbation in the meridional plane, we show that the mode with geometry dominated by ℓ = 1 at moderate rotation rates (mode #1), ends up with a dominant ℓ = 5 component at high rotation rates.
Therefore, this work (see also Reese et al. 2009) shows that modes in rapidly rotating stars can no longer be identified by one angular degree ℓ only. Indeed, when rotation increases, the different components in the spectral expansion are more and more coupled by the nonspherical terms of the system of equations. As a consequence, each mode is composed of a mixture of spherical harmonics of the same symmetry, and it is not even possible to follow a mode considering its dominant component, as it can change during an avoided crossing.
7. Conclusion
A new oscillation code that computes nonradial adiabatic pulsations in rotating stars has been developed. The accuracy of the calculations was achieved with a hybrid method based on a spectral expansion on the spherical harmonics basis and a fourthorder finite differences scheme. The code was tested and validated in the present study for polytropic models, but no assumptions were made in the implementation concerning the structure of the model used as input to ACOR.
Although we limited ourselves to barotropic models in this paper, it must be emphasized that our code is fully able to handle nonbarotropic rapidly rotating models, as will be presented in a forthcoming paper. This is an important point for studying the pulsations of realistic stellar models such as those including rotationally induced transport processes according to Zahn (1992). Indeed, these processes lead to shellular rotation profiles in radiative zones, which are as a consequence nonbarotropic. Moreover, the radial treatment based on a finite differences method that is accurate up to the fourthorder is particularly wellsuited for stellar models that present sharp variations of the structural quantities.
Acknowledgments
R.M.O. is indebted to the “Fédération WallonieBruxelles – Fonds Spéciaux pour la Recherche/Crédit de démarrage – Université de Liège” for financial support. D.R.R. is financially supported through a postdoctoral fellowship from the “Subside fédéral pour la recherche 2011”, University of Liège. The authors would like to thank J. Ballot for providing g modes frequencies. R.M.O. also thanks MarieJo Goupil for financial support and significant scientific contribution.
References
 Beck, P. G.,Montalban, J.,Kallinger, T., et al. 2012, Nature, 481, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Bonazzola, S.,Gourgoulhon, E., &Marck, J.A. 1998, Phys. Rev. D, 58, 104020 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, C., Hussaini, M. Y., Quarteroni, A., & Zang, T. A. 1988, Spectral methods in Fluid Dynamics, eds. R. Glowinski, M. Holt, P. Hut, H. B. Keller, J. Killeen, S. A. Orsag, & V. V. Rusanov (SpringerVerlag), 1 [Google Scholar]
 CohenTannoudji, C., Dui, B., & Laloe, F. 1973, Collection Enseignement des Sciences (Paris: Herman) [Google Scholar]
 Dupret, M. A. 2001, A&A, 366, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dziembowski, W. A., &Goode, P. R. 1992, ApJ, 394, 670 [NASA ADS] [CrossRef] [Google Scholar]
 Eddington, A. S. 1925, The Observatory, 48, 73 [NASA ADS] [Google Scholar]
 Gough, D. O., &Thompson, M. J. 1990, MNRAS, 242, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., &Weigert, A. 1994, Stellar Structure and Evolution, XVI, 468, also Astronomy and Astrophysics Library (Berlin, Heidelberg, New York: SpringerVerlag), 192 [Google Scholar]
 Ledoux, P. 1951, ApJ, 114, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Lignières, F.,Rieutord, M., &Reese, D. 2006, A&A, 455, 607 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S.,Palacios, A., &Zahn, J.P. 2004, A&A, 425, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D.,Lignières, F., &Rieutord, M. 2006, A&A, 455, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. R.,Thompson, M. J.,MacGregor, K. B., et al. 2009, A&A, 506, 183 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 1987, Geophys. Astrophys. Fluid Dyn., 39, 163 [Google Scholar]
 Rieutord, M., Corbard, T., Pichon, B., Dintrans, B., & Lignières, F. 2005, SF2A2005: Semaine de l’Astrophysique Francaise, 759 [Google Scholar]
 Scuflaire, R.,Montalbán, J.,Théado, S., et al. 2008, Ap&SS, 316, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Soufi, F.,Goupil, M. J., &Dziembowski, W. A. 1998, A&A, 334, 911 [NASA ADS] [Google Scholar]
 Suárez, J. C.,Goupil, M. J.,Reese, D. R., et al. 2010, ApJ, 721, 537 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars, 2nd edn. (Tokyo: University of Tokyo Press) [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
Appendix A: Equation system
Rather than projecting the motion equation onto the basis vectors, we chose to decompose it into one radial, one poloidal, and one toroidal field. We obtain one equation without derivatives relative to the angular coordinates by projecting the motion equation onto . Then by taking the divergence, and the pseudoradial componant of the curl of the motion equation, we obtain two equations without radial derivatives.
A.1. Radial motion
(A.1)with, (A.2)where the colatitude θ is replaced by μ = cosθ, so r_{θ} = − sinθ r_{μ}. As stated in Sect. 2.2, we applied the change of variable π^{′} = p^{′} / ρ to avoid singularity problems at the surface of polytropic models.
A.2. Poloidal motion
Here we calculate the divergence of the horizontal component of the motion equation. To do so, one has to apply the operator r∇_{ ⊥ }·, where for any vector : (A.3)Once applied to the motion equation, this gives (A.4)where g is given by , and the operator L^{2} is the angular momentum operator: (A.5)whose eigenfunctions are the spherical harmonics: . Thus, the equation becomes (A.6)with, (A.7)
Central behavior of the first terms in the spectral expansion in different cases.
A.3. Toroidal motion
The curl of the motion equation is (Unno et al. 1989) (A.8)where we introduce the vorticity vector: ω = ∇ × v. To find the toroidal component of the motion, we take the pseudoradial component of the above equation (i.e., along ): (A.9)where
t is given in Eq. (A.7) and a, b, c, d, e, f and h are defined by
A.4. Adiabatic relation
We define quantities analogous to the Schwarzschild discriminant in the two dimensional case: Therefore, the adiabatic relation given in Eq. (6) can be written in the following form: (A.12)where we have introduced , and .
A.5. Conservation of mass
After linearization of Eq. (4), one obtains With the help of Eq. (6), the equation can be rewritten which finally gives (A.13)with
A.6. Poisson equation
Computing the Laplacien of Φ^{′} in the new coordinate system, one obtains a dimensionless version of the Poisson equation: (A.14)As mentioned in Sect. 2.2, we introduce a new variable and add an auxiliary equation to have a firstorder differential system: with, (A.17)Equations (A.1), (A.6), (A.9), (A.12), (A.13), (A.15), and (A.16) are the seven equations that make up the firstorder differential equation system, which yields 2D pulsations of any rotating model.
In the second domain V_{2}, i.e., in the vacuum, only Eqs. (A.15) and (A.16) remain with ρ^{′} = 0 for the Poisson equation.
Appendix B: Central boundary conditions
As explained in Sect. 2.6, we chose to impose two boundary conditions at the center. The goal of this section is to find two algebraic equations to replace the differential ones at the center of the star.
Table B.1 shows the central behavior of the first terms in the spectral expansion of the unknowns. As can be seen in the table, the parity of the w_{ℓ} coefficients is different from the parity of other spectral coefficients.
The global parity of the mode is defined by the parity of u_{ℓ}, v_{ℓ}, , and . For an odd axisymmetric mode, w_{0} has no meaning since the spherical harmonic Y_{ℓ = 0, m = 0} is spherically symmetric and has no toroidal component. A similar argument applies to the v_{0} component in even axisymmetric modes.
To emphasize the dominant terms in the equations near the center, we propose to express the spectral components in the following form: (B.1)where the tilded quantities have a constant behavior at the center. For the derivative of the gravitational potential, we obtain .
Introducing this variable change into the system of equations and taking the limit when ζ goes to zero, we obtain
Radial motion (Eq. (A.1))
Where ℳ_{ρ} = M_{ρ} / ζ and ℳ_{π} = M_{π} / ζ.
Poloidal motion (Eq. (A.6)) (B.3)Toroidal motion (Eq. (A.9)) (B.4)Adiabatic relation (Eq. (A.12)) (B.5)Where .
Continuity equation (Eq. (A.13)) (B.6)Where .
Poisson’s equation (Eq. (A.16)) (B.7)where is given by (B.8)and where ϵ is the flatness, given in Sect. 2.1.
We note, first of all, that the equations of radial motion, of continuity, and the Poisson equation are singular at the center. To enforce the regularity of those equations, the boundary conditions that are imposed are given by the fact that singular terms go to zero when ζ → 0. Therefore the radial motion and the Poisson equations lead to
This case is specific because, near the center, . Thus, the radial motion equation is not singular, whereas the continuity equation is (B.11)Therefore, the algebraic equation to impose near the center when  m  = 0, par = 0 and ℓ = 0 is (B.12)The Poisson equation is also singular in this specific case, and Eq. (B.10) is still relevant.
It seems then that we obtained the algebraic boundary conditions necessary to impose regular behavior of the unknowns near the center. However, the condition given in Eq. (B.9) is redundant with a linear combination of Eqs. (B.3) and (B.4). Hence, the system is underdetermined at the center. It is not possible to enforce the correct behavior on the unknowns near the center with the tilded variables. We therefore have to directly impose central conditions on the nontilded spectral components. Here again, different cases have to be explored:
The general case is composed of the following symmetry cases:

 m  = 0 and par = 1,

 m  ≠ 0 and par = 0,
as well as all the components after the first one for

 m  = 0 and par = 0, i.e., ℓ = 2,4,6,...

 m  ≠ 0 and par = 1, i.e., ℓ =  m  + 3,  m  + 5,  m  + 7,...
The regularity of the velocity field leads to The regularity of the motion equation leads to (B.15)The regularity of Poisson’s equation leads to (B.16)
In this specific case, w_{ℓ − 1} has no meaning, we thus impose where the condition over the radial velocity component is obtained from the continuity equation.
The continuity equation, the toroidal and radial motion equations, and the Poisson equation give, respectively, The components w_{ℓ} are neglected at the center in the remaining cases. This is an approximation, but after verifications, this method is the most efficient and numerically suitable to enforce regular behavior of the solutions at the center.
Appendix C: Comparison of eigenfunctions
Fig. C.1 Radial parts of the eigenfunction for a low frequency mode dominated by an ℓ = 1, m = 0 component, for an N = 3 polytropic model uniformly rotating at 18.5%, 46%, 63.5%, and 83.7% of the Keplerian breakup veloci ty. 
Fig. C.2 Radial parts of the eigenfunction for a low frequency mode dominated by an ℓ = 2, m = 0 component, for an N = 3 polytropic model uniformly rotating at 3.7%, 18.5%, 46%, and 64.5% of the Keplerian breakup velocity. 
All Tables
Numerical resources – i.e., time in minutes and memory in MB or GB – used by ACOR with a spectral resolution M and a radial one N_{r}.
Central behavior of the first terms in the spectral expansion in different cases.
All Figures
Fig. 1 Coordinate system used in ACOR. V_{1} extends from the center to the stellar surface, and V_{2} encompasses the star. 

In the text 
Fig. 2 Convergence as a function of the number of spherical harmonics, taken as the relative frequency difference between computations using consecutive spectral resolutions, M − 1 and M, for an n = 3 mode dominated by an ℓ = 1, m = 0 component, and for three different rotational velocities: 18.5% Ω_{k} in blue, 38% Ω_{k} in orange, and 59% Ω_{k} in red. 

In the text 
Fig. 3 Convergence as a function of the radial resolution, taken as the relative frequency differences between computations using consecutive radial resolutions, r and (r + 1) points, at three different rotational velocities: 18.5% Ω_{k} in green, 37.9% Ω_{k} in light blue, and 58.9% Ω_{k} in dark blue. In both panels, the mode is dominated by the ℓ = 2, m = 0 component. The upper panel corresponds to an n = 3 mode and uses irregular grids, whereas the lower panel shows an n = 9 mode calculated with even distributions of points. 

In the text 
Fig. 4 Discrepancies between eigenfrequencies computed by ACOR and TOP as a function of the rotation angular velocity for five axisymmetric modes: top: for two n = 3 modes, corresponding to ℓ = 1 (light blue) and ℓ = 2 (orange) in the nonrotating case and one n = 1 g mode with ℓ = 2; bottom: n = 9 modes, corresponding to ℓ = 1 (dark blue) and ℓ = 2 (red). 

In the text 
Fig. 5 Behavior of the eigenfrequencies computed by ACOR and TOP with respect of the rotational angular velocity for two multiplets; top: centroid modes with dominant component ℓ = 1, m = 0, and ℓ = 2, m = 0; bottom: sectorial and tesseral modes, dominated by components ℓ = 1, m = ± 1 and ℓ = 2, m = ± 1 and ± 2. 

In the text 
Fig. 6 Avoided crossing, illustrated by plotting the frequency with respect to the rotation rate for two n = 2 modes (referred to as mode #1 in red and mode #2 in blue), which, in the nonrotating case, are ℓ = 1 and ℓ = 5. The harmonic degree given in the graph is the degree of the dominant component for each mode. 

In the text 
Fig. 7 Spatial distribution of the pressure perturbation in a meridional plane for two low frequency modes (radial order n = 2, azimuthal order m = 0) involved in an avoided crossing occurring around Ω = 40%Ω_{k}, as modeled by ACOR. Left: at low rotation velocities, mode #1 is dominated by an ℓ = 1 component and changes nature to become dominated by ℓ = 5 as the velocity increases. Right: at low rotation velocities, mode #2 is dominated by an ℓ = 5 component and changes nature to become dominated by ℓ = 1 as the velocity increases. 

In the text 
Fig. C.1 Radial parts of the eigenfunction for a low frequency mode dominated by an ℓ = 1, m = 0 component, for an N = 3 polytropic model uniformly rotating at 18.5%, 46%, 63.5%, and 83.7% of the Keplerian breakup veloci ty. 

In the text 
Fig. C.2 Radial parts of the eigenfunction for a low frequency mode dominated by an ℓ = 2, m = 0 component, for an N = 3 polytropic model uniformly rotating at 3.7%, 18.5%, 46%, and 64.5% of the Keplerian breakup velocity. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.