Free Access
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/0004-6361/201219548
Published online 01 November 2012

© ESO, 2012

1. Introduction

Rotation plays a key role in the evolution of stars across the Hertzsprung-Russell 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 large-scale 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 first-order effects, Gough & Thompson 1990; and Dziembowski & Goode 1992, for second-order effects; and Soufi et al. 1998, for third-order 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 (Prot ∼ Posc), 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 non-perturbative effects of rotation on oscillation spectra. A new two dimensional non-perturbative code is presented. The 2D non-perturbative calculations fully take into account the centrifugal distortion of the star and include the full influence of the Coriolis acceleration. This 2D non-perturbative 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 north-south 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 two-dimensional 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 Prot ∼ Posc, the perturbative treatment is no longer relevant, and a non-perturbative 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 pseudo-radial 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 V1, which encloses the star, and one external domain, V2, 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 V2. 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 V1: (1)where ζ ranges from 0 to 1, ϵ = 1 − Rpol / Req is the flatness, and r(ζ = 1) = Rs(θ) the stellar surface. In domain V2: (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.

thumbnail Fig. 1

Coordinate system used in ACOR. V1 extends from the center to the stellar surface, and V2 encompasses the star.

Open with DEXTER

2.2. Basic equations

Here, we consider the stellar interior to be an inviscid self-gravitating 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 ei 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 ei(ωt+), where ω is the pulsation frequency, and m the azimuthal order of the pulsation mode.

We put the system of equations in non-dimensional 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 non-dimensional 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 second-order differential equation to two first-order 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 well-suited in fluid dynamics. For instance, for a function of class , the numerical error decreases as eaM, 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 Xi(ζ,θ,ϕ). 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 first-order 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 V2 (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 stress-free 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 V2. 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 first-order 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, wp and (i.e., Eqs. (A.6), (A.9), and (A.12)).

We chose to solve this system by an Newton-like 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 y1 and z1 are column vectors containing the spectral coefficients of the unknowns (18)and where  = | m | + 2(j − 1) + par,   ∀j ∈ [1,M] .

A11, A12, A21 and A22 correspond to the following equations (19)whereas B11, B12, B21 and B12 correspond to the algebraic equations (20)Thanks to the three last equations, the non-differentiated unknowns can be expressed in terms of the differentiated ones with the help of Eq. (17). To do so, the matrix (B21 + δσ   B22), 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 Nr 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 fourth-order 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 h4, 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 y1. 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, Id 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 Nr layers) by projecting the equations onto the spherical harmonics basis (of dimension M). Therefore, by assessing the role of the two parameters M and Nr, 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 eaM. Rotation induces non-spherical 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.

thumbnail 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.

Open with DEXTER

thumbnail 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.

Open with DEXTER

Table 1

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 Nr.

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 non-regular 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 Nr. Note that and are bloc matrices, their size corresponds to the number of non-zero 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 two-dimensional 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 (p0): (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 Nr = 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.

thumbnail 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 non-rotating case and one n = 1 g mode with  = 2; bottom: n = 9 modes, corresponding to  = 1 (dark blue) and  = 2 (red).

Open with DEXTER

thumbnail 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.

Open with DEXTER

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.

thumbnail 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 non-rotating case, are  = 1 and  = 5. The harmonic degree given in the graph is the degree of the dominant component for each mode.

Open with DEXTER

thumbnail 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.

Open with DEXTER

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 Cohen-Tannoudji 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 non-spherical 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 non-radial 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 fourth-order 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 non-barotropic 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 non-barotropic. Moreover, the radial treatment based on a finite differences method that is accurate up to the fourth-order is particularly well-suited for stellar models that present sharp variations of the structural quantities.

Acknowledgments

R.-M.O. is indebted to the “Fédération Wallonie-Bruxelles – 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 Marie-Jo Goupil for financial support and significant scientific contribution.

References

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 pseudo-radial 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 L2 is the angular momentum operator: (A.5)whose eigenfunctions are the spherical harmonics: . Thus, the equation becomes (A.6)with, (A.7)

Table B.1

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 pseudo-radial 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 first-order 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 first-order differential equation system, which yields 2D pulsations of any rotating model.

In the second domain V2, 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, w0 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 v0 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))

(B.2)

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 under-determined 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 non-tilded 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

thumbnail 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 break-up veloci ty.

Open with DEXTER

thumbnail 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 break-up velocity.

Open with DEXTER

All Tables

Table 1

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 Nr.

Table B.1

Central behavior of the first terms in the spectral expansion in different cases.

All Figures

thumbnail Fig. 1

Coordinate system used in ACOR. V1 extends from the center to the stellar surface, and V2 encompasses the star.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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 non-rotating case and one n = 1 g mode with  = 2; bottom: n = 9 modes, corresponding to  = 1 (dark blue) and  = 2 (red).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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 non-rotating case, are  = 1 and  = 5. The harmonic degree given in the graph is the degree of the dominant component for each mode.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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 break-up veloci ty.

Open with DEXTER
In the text
thumbnail 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 break-up velocity.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.