© 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 (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 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 P_rot ∼ P_osc, 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 V₁, which encloses the star, and one external domain, V₂, 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₂. 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₁: $\begin{array}{ccc}\mathit{r}\mathrm{\left(}\mathit{\zeta ,\theta }\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\hspace{0.17em}\mathit{\zeta }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\frac{\mathrm{5}{\mathit{\zeta }}^{\mathrm{3}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{3}{\mathit{\zeta }}^{\mathrm{5}}}{\mathrm{2}}\left({\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{ϵ}\right)\mathit{,}& & \end{array}$(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₂: $\mathit{r}\mathrm{\left(}\mathit{\zeta ,\theta }\mathrm{\right)}\mathrm{=}\mathrm{2}\hspace{0.17em}\mathit{ϵ}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\hspace{0.17em}\mathit{\zeta }\mathrm{+}\mathrm{(}\mathrm{2}{\mathit{\zeta }}^{\mathrm{3}}\mathrm{-}\mathrm{9}{\mathit{\zeta }}^{\mathrm{2}}\mathrm{+}\mathrm{12}\mathit{\zeta }\mathrm{-}{\mathrm{4}}^{\mathrm{)}}\left({\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{ϵ}\right)\mathit{,}$(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: $\begin{array}{ccc}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\zeta }}\end{array}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{{\mathit{r}}_{\mathit{\theta }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}}\hspace{0.17em}\left(\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{r}}\end{array}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\hspace{0.17em}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{\theta }}\end{array}\right)\mathit{,}\\ \begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\theta }}\end{array}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{{\mathit{r}}_{\mathit{\theta }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}}\hspace{0.17em}\left(\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{\theta }}\end{array}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\hspace{0.17em}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{r}}\end{array}\right)\mathit{,}\\ \begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\varphi }}\end{array}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{\varphi }}\end{array}\mathit{,}\end{array}$(3)where

• $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\theta }}\end{array}$ is tangential to iso-ζ surfaces in the meridional plane;

• $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\zeta }}\end{array}$ is directly orthogonal to $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\theta }}\end{array}$ in the meridional plane;

• $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\varphi }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{\varphi }}\end{array}$,

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. V1 extends from the center to the stellar surface, and V2 encompasses the star.

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 $\begin{array}{ccc}& & {\mathit{\partial }}_{\mathit{t}}{\mathit{\rho }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{\nabla }\mathrm{·}\left({\mathit{\rho }}_{\mathrm{0}}{{v}}^{\mathbf{\prime }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{\rho }}^{\mathrm{\prime }}{v}\mathrm{0}\right)\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}\mathit{,}\\ & & \mathrm{[}\mathrm{(}{\mathit{\partial }}_{\mathit{t}}\mathrm{+}\mathrm{\Omega }{\mathit{\partial }}_{\mathit{\varphi }}\mathrm{)}{{\mathit{v}}_{\mathit{i}}^{\mathrm{\prime }}}^{\mathrm{]}}{e}\mathit{i}\mathrm{+}\mathrm{2}\mathbf{\Omega }\mathrm{\times }{{v}}^{\mathbf{\prime }}\mathrm{+}\left({{v}}^{\mathbf{\prime }}\mathrm{·}\mathrm{\nabla \Omega }\right)\mathit{r}\sin \mathit{\theta }{e}\mathit{\varphi }\mathrm{=}\\ & & \mathrm{-}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{\nabla }{\mathit{p}}^{\mathrm{\prime }}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{\rho }}^{\mathrm{\prime }}}{{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{\nabla }{\mathit{p}}_{\mathrm{0}}\mathit{,}\\ & & \left({\mathit{\partial }}_{\mathit{t}}\mathrm{+}{v}\mathrm{0}\mathrm{·}\mathrm{\nabla }\right)\left(\frac{{\mathit{\rho }}^{\mathrm{\prime }}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{-}\frac{{\mathit{p}}^{\mathrm{\prime }}}{{\mathrm{\Gamma }}_{\mathrm{1}}{\mathit{p}}_{\mathrm{0}}}\right)\mathrm{+}{{v}}^{\mathbf{\prime }}\mathrm{·}\left(\frac{\mathrm{\nabla }{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{-}\frac{\mathrm{\nabla }{\mathit{p}}_{\mathrm{0}}}{{\mathrm{\Gamma }}_{\mathrm{1}}{\mathit{p}}_{\mathrm{0}}}\right)\mathrm{=}\mathrm{0}\mathit{,}\\ & & {\mathrm{\nabla }}^{\mathrm{2}}{\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}^{\mathrm{\prime }}\mathit{,}\end{array}$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 non-dimensional form using the following transformations: $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ \mathit{\rho }\end{array}& \mathrm{=}& \left(\frac{\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{eq}}^{\mathrm{3}}}{\mathit{M}}\right)\mathit{\rho }\begin{array}{}\mathrm{˜}\\ \mathit{p}\end{array}\mathrm{=}\left(\frac{\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{eq}}^{\mathrm{4}}}{\mathit{G}{\mathit{M}}^{\mathrm{2}}}\right)\mathit{p}\begin{array}{}\mathrm{˜}\\ \mathit{\phi }\end{array}\mathrm{=}\left(\frac{{\mathit{R}}_{\mathrm{eq}}}{\mathit{GM}}\right)\mathit{\phi }\\ \mathit{\sigma }& \mathrm{=}& \frac{\mathit{\omega }}{{\mathrm{\Omega }}_{\mathrm{k}}}\begin{array}{}\mathrm{˜}\\ \mathrm{\Omega }\end{array}\mathrm{=}\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_{\mathrm{k}}}\mathit{,}\end{array}$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 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 ${\mathit{Y}}_{\mathit{\ell }}^{\mathit{m}}$ (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_ϕ) $\begin{array}{ccc}{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}& \mathrm{=}& \mathrm{i}\sum _{\mathit{\ell }\mathrm{\ge }\mathrm{\left|}\mathit{m}\mathrm{\right|}}^{\mathrm{+}\mathrm{\infty }}{\mathit{u}}_{\mathit{\ell }}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}{\mathit{Y}}_{\mathit{\ell ,m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\mathit{,}\\ {\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}& \mathrm{=}& \mathrm{i}\sum _{\mathit{\ell }\mathrm{\ge }\mathrm{\left|}\mathit{m}\mathrm{\right|}}^{\mathrm{+}\mathrm{\infty }}\left({\mathit{v}}_{\mathit{\ell }}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\hspace{0.17em}{\mathit{\partial }}_{\mathit{\theta }}{\mathit{Y}}_{\mathit{\ell ,m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\mathrm{+}{\mathit{w}}_{\mathit{\ell }}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\frac{\mathit{m}}{\sin \mathit{\theta }}{\mathit{Y}}_{\mathit{\ell ,m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\right)\mathit{,}\\ {\mathit{v}}_{\mathit{\varphi }}^{\mathrm{\prime }}& \mathrm{=}& \end{array}$(8)All other scalar unknowns, namely Φ^′, dΦ^′, and ρ^′, are expanded in the same way as the scaled pressure perturbation: $\begin{array}{ccc}{\mathit{\pi }}^{\mathrm{\prime }}\mathrm{=}\sum _{\mathit{\ell }\mathrm{\ge }\mathrm{\left|}\mathit{m}\mathrm{\right|}}^{\mathrm{+}\mathrm{\infty }}{\mathit{\pi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}{\mathit{Y}}_{\mathit{\ell ,m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\mathit{.}& & \end{array}$(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] ,$\begin{array}{ccc}& & \mathit{\ell }\mathrm{=}\mathrm{\left|}\mathit{m}\mathrm{\right|}\mathrm{+}\mathrm{2}\mathrm{(}\mathit{j}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{+}\hspace{0.17em}\mathit{par,}\hspace{0.17em}\mathrm{for}\hspace{0.17em}{\mathit{\pi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}{\mathit{\phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}\mathrm{d}{\mathit{\phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}{\mathit{u}}_{\mathit{\ell }}\mathit{,}\hspace{0.17em}{\mathit{v}}_{\mathit{\ell }}\mathit{,}{\mathit{\rho }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\\ & & {\mathit{\ell }}_{\mathit{p}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{\left|}\mathit{m}\mathrm{\right|}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{2}\mathrm{(}\mathit{j}\mathrm{-}\mathrm{1}\mathrm{)}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{par}\mathrm{for}{\mathit{w}}_{{\mathit{\ell }}_{\mathit{p}}}\mathit{,}\end{array}$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 ${\mathit{u}}_{\mathit{\ell }}\mathit{,}\hspace{0.17em}{\mathit{v}}_{\mathit{\ell }}\mathit{,}\hspace{0.17em}{\mathit{w}}_{{\mathit{\ell }}_{\mathit{p}}}\mathit{,}\hspace{0.17em}{\mathit{\pi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}{\mathit{\rho }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}{\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}$, and $\mathit{d}{\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}$.

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 $\begin{array}{ccc}\mathrm{E}\left({\mathit{X}}_{\mathit{i}}\mathrm{\left(}\mathit{\zeta ,\theta ,\varphi }\mathrm{\right)}\right)\mathrm{=}\mathrm{E}\left(\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathit{X}}_{\mathit{i,}{\mathit{\ell }}_{\mathrm{2}}}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\hspace{0.17em}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{2}}\mathit{,}\mathrm{m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\right)\mathrm{=}\mathrm{0}& & \end{array}$(12)is replaced by a system of M differential equations in ζ, obtained by projecting these equations on a basis of M spherical harmonics: $\begin{array}{ccc}\mathrm{\int }\frac{\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\varphi }}{\mathrm{4}\mathit{\pi }}\hspace{0.17em}\mathrm{E}\left(\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathit{X}}_{\mathit{i,}{\mathit{\ell }}_{\mathrm{2}}}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}\hspace{0.17em}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{2}}\mathit{,}\mathrm{m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\right)\hspace{0.17em}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{1}}\mathit{,}\mathrm{m}}^{\mathrm{\ast }}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(13)in which ℓ₁ and ℓ₂ 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 V₂ (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 $\begin{array}{ccc}{\mathit{\pi }}_{\mathit{\ell }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{𝒪}\mathrm{(}{\mathit{\zeta }}^{\mathit{\ell }}\mathrm{)}{\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{=}\mathrm{𝒪}\mathrm{(}{\mathit{\zeta }}^{\mathit{\ell }}\mathrm{)}{\mathit{\rho }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{=}\mathrm{𝒪}\mathrm{(}{\mathit{\zeta }}^{\mathit{\ell }}\mathrm{)}& & \end{array}$and the components of the velocity field obey $\begin{array}{ccc}{\mathit{u}}_{\mathit{\ell }}& \mathrm{=}& \mathrm{𝒪}\mathrm{(}{\mathit{\zeta }}^{\mathit{\ell }\mathrm{-}\mathrm{1}}\mathrm{)} \mathrm{and}{\mathit{u}}_{\mathrm{0}}\mathrm{=}\mathrm{𝒪}\left(\mathit{\zeta }\right)\mathit{,}\\ {\mathit{v}}_{\mathit{\ell }}& \mathrm{=}& \mathrm{𝒪}\mathrm{(}{\mathit{\zeta }}^{\mathit{\ell }\mathrm{-}\mathrm{1}}\mathrm{)} \mathrm{and}{\mathit{v}}_{\mathrm{0}}\mathrm{=}\mathrm{0}\mathit{,}\\ {\mathit{w}}_{\mathit{\ell }}& \mathrm{=}& \mathrm{𝒪}\mathrm{(}{\mathit{\zeta }}^{\mathit{\ell }}\mathrm{)}\mathit{.}\end{array}$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 $\begin{array}{ccc}\mathrm{i}\mathrm{(}\mathit{m}\mathrm{\Omega }\mathrm{+}\mathit{\sigma }\mathrm{)}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\frac{{\mathit{\partial }}_{\mathit{\zeta }}{\mathit{p}}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}\hspace{0.17em}{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}\mathit{,}& & \end{array}$(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₂. The advantage in using a spherical boundary is that the Laplace equation is separable for solutions of the form ${\mathit{X}}_{\mathit{i}}\mathrm{\left(}\mathit{r,\theta ,\varphi }\mathrm{\right)}\mathrm{=}{\mathit{X}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{Y}}_{\mathit{\ell }}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}$. 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: ${\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{=}\mathit{A}\hspace{0.17em}{\mathit{r}}^{\mathit{\ell }}$, which diverges at infinity and is therefore discarded, and ${\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{=}\mathit{B}\hspace{0.17em}\mathit{/}\hspace{0.17em}{\mathit{r}}^{\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}$, which vanishes. Hence, the corresponding boundary condition is $\begin{array}{ccc}\mathrm{d}{\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}\frac{{\mathrm{r}}_{\mathit{\zeta }}}{\mathrm{r}}\hspace{0.17em}{\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}& & \end{array}$(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_ℓ, ${\mathit{\pi }}_{\mathit{\ell }}^{\mathrm{\prime }}$, ${\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}$ and d${\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}$ (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 ${\mathit{\rho }}_{\mathit{\ell }}^{\mathrm{\prime }}$ (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 σ₀ and looking for the pulsation mode with the closest frequency to this guess. Solving the system yields a deviation, δσ, from the initial guess, σ = σ₀ + δσ 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: $\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{y}}_{\mathrm{1}}}{\mathrm{d}\mathit{\zeta }}\mathrm{=}\mathrm{(}{\mathit{A}}_{\mathrm{11}}\mathrm{+}\mathit{\delta \sigma }\hspace{0.17em}{\mathit{A}}_{\mathrm{12}}\mathrm{)}{\mathit{y}}_{\mathrm{1}}\mathrm{+}\mathrm{(}{\mathit{A}}_{\mathrm{21}}\mathrm{+}\mathit{\delta \sigma }{\mathit{A}}_{\mathrm{22}}\mathrm{)}\hspace{0.17em}{\mathit{z}}_{\mathrm{1}}\mathit{,}\\ & & \end{array}$where y₁ and z₁ are column vectors containing the spectral coefficients of the unknowns ${\mathit{y}}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\left[\begin{array}{c}\\ {\mathit{\pi }}_{\mathit{\ell }}^{\mathrm{\prime }}\\ {\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\\ \mathrm{d}{\mathrm{\Phi }}_{\mathit{\ell }}^{\mathrm{\prime }}\\ {\mathit{u}}_{\mathit{\ell }}\end{array}\right]\mathrm{and}{\mathit{z}}_{\mathrm{1}}\mathrm{=}\left[\begin{array}{c}\\ {\mathit{v}}_{\mathit{\ell }}\\ {\mathit{w}}_{{\mathit{\ell }}_{\mathit{p}}}\\ {\mathit{\rho }}_{\mathit{\ell }}^{\mathrm{\prime }}\end{array}\right]\mathit{,}$(18)and where ℓ = | m | + 2(j − 1) + par,   ∀j ∈ [1,M] .

A₁₁, A₁₂, A₂₁ and A₂₂ correspond to the following equations $\left[\begin{array}{c}\\ \hspace{0.17em}\mathrm{pseudo}\mathrm{-}\mathrm{radialmotion}\mathrm{(}\mathrm{Eq}\mathrm{.} \mathrm{(}\mathrm{A}\mathrm{.1}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\\ \hspace{0.17em}\mathrm{definitionofd}{\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{(}\mathrm{Eq}\mathrm{.} \mathrm{(}\mathrm{A}\mathrm{.15}\mathrm{)}\hspace{0.17em}\\ \hspace{0.17em}\mathrm{Poisson}\mathrm{(}\mathrm{Eq}\mathrm{.} \mathrm{(}\mathrm{A}\mathrm{.16}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\\ \hspace{0.17em}\mathrm{continuity}\mathrm{(}\mathrm{Eq}\mathrm{.} \mathrm{(}\mathrm{A}\mathrm{.13}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\end{array}\right]\mathit{,}$(19)whereas B₁₁, B₁₂, B₂₁ and B₁₂ correspond to the algebraic equations $\left[\begin{array}{c}\\ \hspace{0.17em}\mathrm{poloidalmotion}\mathrm{\left(}\mathrm{Eq}\mathrm{.}\mathrm{\right(}\mathrm{A}\mathrm{.6}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\\ \hspace{0.17em}\mathrm{toroidalmotion}\mathrm{\left(}\mathrm{Eq}\mathrm{.}\mathrm{\right(}\mathrm{A}\mathrm{.9}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\\ \hspace{0.17em}\mathrm{adiabaticrelation}\mathrm{\left(}\mathrm{Eq}\mathrm{.}\mathrm{\right(}\mathrm{A}\mathrm{.12}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\end{array}\right]\mathit{.}$(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 (B₂₁ + δσ   B₂₂), which is a real square matrix of rank 3M, has to be inverted. It is straightforward to show that if δσ is small enough, ${\mathrm{(}{\mathit{B}}_{\mathrm{21}}\mathrm{+}\mathit{\delta \sigma }\hspace{0.17em}{\mathit{B}}_{\mathrm{22}}\mathrm{)}}^{-1}\mathrm{=}{\mathit{B}}_{\mathrm{21}}^{-1}\mathrm{-}\mathit{\delta \sigma }\hspace{0.17em}{\mathit{B}}_{\mathrm{21}}^{-1}\hspace{0.17em}{\mathit{B}}_{\mathrm{22}}\hspace{0.17em}{\mathit{B}}_{\mathrm{21}}^{-1}\mathrm{+}\mathrm{o}\mathrm{\left(}\mathit{\delta }{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\right)}\mathit{.}$(21)Thus, the matrix system can be written in the form: $\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{y}}_{\mathrm{1}}}{\mathrm{d}\mathit{\zeta }}\mathrm{=}\mathrm{(}\mathit{A}\mathrm{+}\mathit{\delta \sigma }\hspace{0.17em}{\mathit{A}}_{\mathit{\delta \sigma }}\mathrm{)}\hspace{0.17em}{\mathit{y}}_{\mathrm{1}}\\ & & {\mathit{z}}_{\mathrm{1}}\mathrm{=}\mathrm{(}\mathit{B}\mathrm{+}\mathit{\delta \sigma }\hspace{0.17em}{\mathit{B}}_{\mathit{\delta \sigma }}\mathrm{)}\hspace{0.17em}{\mathit{y}}_{\mathrm{1}}\mathit{.}\end{array}$

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 fourth-order difference scheme, which relies on the following identity ${\mathit{y}}_{\mathit{i}}\mathrm{+}\frac{\mathit{h}}{\mathrm{2}}{\mathit{y}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{12}}{\mathit{y}}_{\mathit{i}}^{\mathrm{\prime \prime }}\mathrm{=}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{-}\frac{\mathit{h}}{\mathrm{2}}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{12}}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{+}\mathit{o}\mathrm{\left(}{\mathit{h}}^{\mathrm{5}}\mathrm{\right)}\mathit{,}$(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⁴, 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₁. 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${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{+}}{\mathit{y}}_{\mathit{i}}\mathrm{-}{\mathit{\alpha }}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{-}}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{=}\mathit{\delta \sigma }\hspace{0.17em}\mathrm{[}{\mathit{\beta }}_{\mathit{i}}^{\mathrm{+}}{\mathit{y}}_{\mathit{i}}\mathrm{+}{\mathit{\beta }}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{-}}{\mathit{y}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{]}\mathrm{+}\mathit{o}\mathrm{\left(}{\mathit{h}}^{\mathrm{5}}\mathrm{\right)}\mathit{,}$(25)where $\begin{array}{ccc}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{+}}& \mathrm{=}& {\mathit{I}}_{\mathit{d}}\mathrm{+}\frac{\mathit{h}}{\mathrm{2}}\hspace{0.17em}{\mathit{A}}_{\mathit{i}}\mathrm{+}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{12}}\hspace{0.17em}\mathrm{(}{\mathit{A}}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{{\mathit{A}}_{\mathit{i}}^{\mathrm{\prime }}}^{\mathrm{)}}\mathit{,}\\ {\mathit{\beta }}_{\mathit{i}}^{\mathrm{+}}& \mathrm{=}& \frac{\mathit{h}}{\mathrm{2}}\hspace{0.17em}{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}}\mathrm{+}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{12}}\hspace{0.17em}\mathrm{(}{\mathit{A}}_{\mathit{i}}\hspace{0.17em}{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}}\mathrm{+}{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}}\hspace{0.17em}{\mathit{A}}_{\mathit{i}}\mathrm{+}{{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}}^{\mathrm{\prime }}}^{\mathrm{)}}\mathit{,}\\ {\mathit{\alpha }}_{\mathit{i}}^{\mathrm{-}}& \mathrm{=}& {\mathit{I}}_{\mathit{d}}\mathrm{-}\frac{\mathit{h}}{\mathrm{2}}\hspace{0.17em}{\mathit{A}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{12}}\hspace{0.17em}\mathrm{(}{\mathit{A}}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{A}}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{\prime }}}^{\mathrm{)}}\mathit{,}\\ {\mathit{\beta }}_{\mathit{i}}^{\mathrm{-}}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathit{h}}{\mathrm{2}}\hspace{0.17em}{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}\\ & & \mathrm{+}\hspace{0.17em}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{12}}\hspace{0.17em}\mathrm{(}{\mathit{A}}_{\mathit{i}\mathrm{+}\mathrm{1}}\hspace{0.17em}{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{+}{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}\hspace{0.17em}{\mathit{A}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{+}{{\mathit{A}}_{\mathit{\delta \sigma ,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{\prime }}}^{\mathrm{)}}\mathit{.}\end{array}$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: $\mathrm{𝒜}\hspace{0.17em}\mathrm{𝒴}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\delta \sigma }\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}\hspace{0.17em}\mathrm{𝒴}\mathit{,}$(26)where the vector has been introduced: $\mathrm{𝒴}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\left[\begin{array}{c}\\ {\mathit{y}}_{\mathrm{1}}\\ {\mathit{y}}_{\mathrm{2}}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{y}}_{\mathit{N}}\end{array}\right]\mathrm{(}\mathrm{layers}\mathrm{1}\mathit{,}\hspace{0.17em}\mathit{...}\mathit{,}{\mathit{N}}_{\mathit{r}}\mathrm{)}$(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, δσ₀, of the eigenvalue by the formula ${\mathrm{𝒴}}_{\mathit{k}\mathrm{+}\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{(}\mathrm{𝒜}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\delta }{\mathit{\sigma }}_{\mathrm{0}}\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}{\mathrm{)}}^{-1}\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}\hspace{0.17em}{\mathrm{𝒴}}_{\mathit{k}}\mathit{.}$(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: $\mathrm{(}\mathrm{𝒜}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\delta }{\mathit{\sigma }}_{\mathrm{0}}\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}\mathrm{)}\hspace{0.17em}{\mathrm{𝒴}}_{\mathit{k}\mathrm{+}\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}\hspace{0.17em}{\mathrm{𝒴}}_{\mathit{k}}\mathit{.}$(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 $\begin{array}{ccc}\mathit{\delta \sigma }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathrm{𝒴}}^{\mathrm{\ast }}\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}^{\mathrm{\ast }}\hspace{0.17em}\mathrm{𝒜}\hspace{0.17em}\mathrm{𝒴}}{{\mathrm{𝒴}}^{\mathrm{\ast }}\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}^{\mathrm{\ast }}\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}\hspace{0.17em}\mathrm{𝒴}}\hspace{0.17em}\mathit{,}& & \end{array}$(30)where and ${\mathrm{𝒜}}_{\mathit{\delta \sigma }}^{\mathrm{\ast }}$ are the Hermitian conjugates of and , respectively. It can be easily shown that δσ, given by Eq. (30), minimizes: $\begin{array}{ccc}{\mathit{S}}^{\mathrm{2}}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}{\mathrm{𝒴}}^{\mathrm{\ast }}\hspace{0.17em}\mathrm{(}{\mathrm{𝒜}}^{\mathrm{\ast }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\delta \sigma }\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}^{\mathrm{\ast }}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{𝒜}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\delta \sigma }\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}\mathrm{)}\hspace{0.17em}\mathrm{𝒴}\\ & \mathrm{=}& \hspace{0.17em}\mathrm{\left|}\mathrm{\right|}\mathrm{(}\mathrm{𝒜}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\delta \sigma }\hspace{0.17em}{\mathrm{𝒜}}_{\mathit{\delta \sigma }}\mathrm{)}\hspace{0.17em}\mathrm{𝒴}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}\mathit{.}\end{array}$(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 Γ₁  = 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 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.

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.

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 N_r. 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 (ρ₀) and pressure (p₀): $\begin{array}{ccc}{\mathit{p}}_{\mathrm{0}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{K}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}^{\mathit{\gamma }}\mathit{,}& & \end{array}$(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 non-rotating 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 non-rotating 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 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 $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{r}}\end{array}$. 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

$\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\pi }}^{\mathrm{\prime }}}{\mathit{\partial \zeta }}\hspace{0.17em}& \mathrm{=}& \mathrm{-}\hspace{0.17em}\mathrm{i}\hspace{0.17em}\left(\mathit{m}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\sigma }\right)\hspace{0.17em}{\mathit{M}}_{\mathit{\zeta }}\hspace{0.17em}{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{i}\hspace{0.17em}\left(\mathit{m}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\sigma }\right)\hspace{0.17em}\sin \mathit{\theta }\hspace{0.17em}{\mathit{M}}_{\mathit{\mu }}\hspace{0.17em}\stackrel{}{{{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}}_{\mathrm{˜}}}\\ & & \mathrm{-}\hspace{0.17em}\sin \mathit{\theta }\hspace{0.17em}{\mathit{M}}_{\mathit{\phi }}\hspace{0.17em}{\mathit{v}}_{\mathit{\phi }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathit{\partial }{\mathrm{\Phi }}^{\mathrm{\prime }}}{\mathit{\partial \zeta }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{M}}_{{\mathit{\rho }}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{\rho }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{M}}_{{\mathit{\pi }}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{\prime }}\mathit{,}\end{array}$(A.1)with, $\begin{array}{ccc}{\mathit{M}}_{\mathit{\zeta }}& \mathrm{=}& {\mathit{r}}_{\mathit{\zeta }}\mathit{g}\mathrm{\left(}\mathit{\zeta ,\mu }\mathrm{\right)}\mathit{,}\\ {\mathit{M}}_{\mathit{\mu }}& \mathrm{=}& \mathrm{-}{\mathit{r}}_{\mathit{\zeta }}\mathit{g}\mathrm{\left(}\mathit{\zeta ,\mu }\mathrm{\right)}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\mathit{,}\\ {\mathit{M}}_{\mathit{\phi }}& \mathrm{=}& \mathrm{-}\mathrm{2}\mathrm{\Omega }{\mathit{r}}_{\mathit{\zeta }}\mathit{,}\\ {\mathit{M}}_{{\mathit{\rho }}^{\mathrm{\prime }}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{2}}}{\frac{\mathit{\partial p}}{\mathit{\partial \zeta }}}_{\mathit{\theta ,\varphi }}\mathit{,}\\ {\mathit{M}}_{{\mathit{\pi }}^{\mathrm{\prime }}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}{\frac{\mathit{\partial \rho }}{\mathit{\partial \zeta }}}_{\mathit{\theta ,\varphi }}\mathit{,}\end{array}$(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 ${X}\mathrm{=}{\mathit{X}}_{\mathit{\zeta }}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\zeta }}\end{array}\mathrm{+}{\mathit{X}}_{\mathit{\theta }}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\theta }}\end{array}\mathrm{+}{\mathit{X}}_{\mathit{\varphi }}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{e}}_{\mathit{\varphi }}\end{array}$: $\begin{array}{ccc}\mathbf{\nabla }\mathrm{\perp }\mathrm{·}{X}\mathrm{=}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{r}\sin \mathit{\theta }}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial \theta }}{\left(\sin \mathit{\theta }{\mathit{X}}_{\mathit{\theta }}\right)}_{\mathit{\zeta }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{r}\sin \mathit{\theta }}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial \varphi }}\mathrm{(}{\mathit{X}}_{\mathit{\varphi }}\mathrm{)}\mathit{.}& & \end{array}$(A.3)Once applied to the motion equation, this gives $\begin{array}{ccc}& & \frac{\mathrm{i}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}{\left(\sin \mathit{\theta }\mathrm{(}\mathit{m}\mathrm{\Omega }\mathrm{+}\mathit{\sigma }\mathrm{)}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\right)}_{\mathit{\zeta }}\mathrm{-}\frac{\mathrm{2}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}(\mathrm{\Omega }\sin \mathit{\theta }(\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\sin \mathit{\theta }\\ & & { \mathrm{+}\cos {\mathit{\theta }}^{\mathrm{)}}{{\mathit{v}}_{\mathit{\phi }}^{\mathrm{\prime }}}^{\mathrm{)}}}_{\mathit{\zeta }}\mathrm{-}\frac{\mathit{m}\hspace{0.17em}\mathrm{(}\mathit{m}\mathrm{\Omega }\mathrm{+}\mathit{\sigma }\mathrm{)}}{\sin \mathit{\theta }}{\mathit{v}}_{\mathit{\phi }}^{\mathrm{\prime }}\\ & & \mathrm{+}\mathrm{2}\mathrm{i}\mathit{m}\mathrm{\Omega }\mathit{g}\left(\cot \mathit{\theta }\mathrm{+}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\right){\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{+}\mathrm{i}\mathit{mg}\\ & & \mathrm{\times }\left[\mathrm{2}\mathrm{\Omega }\left(\mathrm{1}\mathrm{-}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\cot \mathit{\theta }\right)\mathrm{+}\frac{\mathit{r}}{{\mathit{r}}_{\mathit{\zeta }}}\mathrm{(}\mathrm{1}\mathrm{-}\frac{{\mathit{r}}_{\mathit{\theta }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{)}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \zeta }}\mathrm{-}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \theta }}\right]{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}\\ & & \mathrm{+}\mathrm{i}\mathit{mg}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \theta }}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{=}\frac{\mathrm{1}}{\mathit{r}}{\mathit{L}}^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}{\mathit{L}}^{\mathrm{2}}{\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{r}}_{\mathit{\theta }}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }{\mathit{\pi }}^{\mathrm{\prime }}}{\mathit{\partial \theta }}\mathrm{+}\frac{{\mathit{r}}_{\mathit{\theta }}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }{\mathrm{\Phi }}^{\mathrm{\prime }}}{\mathit{\partial \theta }}\\ & & \mathrm{-}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}{\left(\frac{\sin \mathit{\theta }}{\mathit{r\rho }}\frac{\mathit{\partial \rho }}{\mathit{\partial \theta }}{\mathit{\pi }}^{\mathrm{\prime }}\right)}_{\mathit{\zeta }}\\ & & \mathrm{+}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}{\left(\frac{\sin \mathit{\theta }}{\mathit{r}{\mathit{\rho }}^{\mathrm{2}}}\frac{\mathit{\partial p}}{\mathit{\partial \theta }}{\mathit{\rho }}^{\mathrm{\prime }}\right)}_{\mathit{\zeta }}\mathit{,}\end{array}$(A.4)where g is given by $\mathit{g}\mathrm{\left(}\mathit{\zeta ,\mu }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}{\mathit{r}}_{\mathit{\mu }}^{\mathrm{2}}\mathit{/}{\mathit{r}}^{\mathrm{2}}\mathrm{)}$, and the operator L² is the angular momentum operator: $\begin{array}{ccc}{\mathit{L}}^{\mathrm{2}}\mathrm{\Psi }\hspace{0.17em}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\hspace{0.17em}{\sin \mathit{\theta }\hspace{0.17em}\frac{\mathit{\partial }\mathrm{\Psi }}{\mathit{\partial \theta }})}_{\mathit{\zeta }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{1}}{{\sin }^{\mathrm{2}}\mathit{\theta }}\hspace{0.17em}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Psi }}{\mathit{\partial }{\mathit{\varphi }}^{\mathrm{2}}}\\ & \mathrm{=}& \hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}{\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Psi }}{\mathit{\partial }{\mathit{\mu }}^{\mathrm{2}}})}_{\mathit{\zeta }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{2}\mathit{\mu }{\frac{\mathit{\partial }\mathrm{\Psi }}{\mathit{\partial \mu }})}_{\mathit{\zeta }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Psi }}{\mathit{\partial }{\mathit{\varphi }}^{\mathrm{2}}}\mathit{,}\end{array}$(A.5)whose eigenfunctions are the spherical harmonics: ${\mathit{L}}^{\mathrm{2}}\hspace{0.17em}{\mathit{Y}}_{\mathit{\ell }}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\ell }\hspace{0.17em}\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{Y}}_{\mathit{\ell }}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}$. Thus, the equation becomes $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \mathrm{(}\mathit{m}\mathrm{\Omega }\mathrm{+}\mathit{\sigma }\mathrm{)}\left[\frac{\mathrm{i}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\left(\sin \mathit{\theta }\hspace{0.17em}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\right)\mathrm{-}\frac{\mathit{m}}{\sin \mathit{\theta }}\hspace{0.17em}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{\prime }}\right]\\ & & \mathrm{-}{\mathit{D}}_{\mathit{\varphi \mu }}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\mathrm{\left(}\sin \mathit{\theta }{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{\right)}\\ & & \mathrm{-}\mathrm{i}{\mathit{D}}_{\mathit{\zeta }}{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}\mathrm{-}\mathrm{i}\sin \mathit{\theta }{\mathit{D}}_{\mathit{\mu }}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\sin \mathit{\theta }{\mathit{D}}_{\mathit{\varphi }}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{r}}\hspace{0.17em}{\mathit{L}}^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{r}}\hspace{0.17em}{\mathit{L}}^{\mathrm{2}}{\mathrm{\Phi }}^{\mathrm{\prime }}\\ & & \mathrm{+}\hspace{0.17em}{\mathit{D}}_{\mathrm{\Phi }\mathit{\mu }}\frac{\mathit{\partial }{\mathrm{\Phi }}^{\mathrm{\prime }}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{D}}_{\mathit{\rho }}{\mathit{\rho }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{D}}_{\mathit{\rho \mu }}\frac{\mathit{\partial }{\mathit{\rho }}^{\mathrm{\prime }}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{D}}_{\mathit{\pi }}{\mathit{\pi }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{D}}_{\mathit{\pi \mu }}\frac{\mathit{\partial }{\mathit{\pi }}^{\mathrm{\prime }}}{\mathit{\partial \mu }}\mathit{,}\end{array}$(A.6)with, $\begin{array}{ccc}& & {\mathit{D}}_{\mathit{\varphi \mu }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}\mathit{t}\\ & \mathrm{where}\mathit{t}\mathrm{\left(}\mathit{\zeta ,\mu }\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{2}\mathrm{(}\mathit{\mu }\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}{\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}}^{\mathrm{)}}\mathit{,}& \\ & & {\mathit{D}}_{\mathit{\zeta }}\mathrm{=}\mathrm{-}\hspace{0.17em}\mathit{mg}\hspace{0.17em}\left[\mathrm{2}\mathrm{\Omega }\hspace{0.17em}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\hspace{0.17em}\cot \mathit{\theta }\right)\hspace{0.17em}\mathrm{+}\hspace{0.17em}\frac{\mathit{r}}{{\mathit{r}}_{\mathit{\zeta }}}\mathrm{(}\mathrm{1}\mathrm{-}\frac{{\mathit{r}}_{\mathit{\theta }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{)}\hspace{0.17em}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \zeta }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \theta }}\right]\mathit{,}\\ & & {\mathit{D}}_{\mathit{\mu }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{m}\hspace{0.17em}\left[\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \mu }}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{g}\mathrm{)}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{\Omega }\hspace{0.17em}\mathit{g}\hspace{0.17em}\mathit{t}}{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\right]\mathit{,}\\ & & {\mathit{D}}_{\mathit{\varphi }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\frac{\mathit{\partial }\mathrm{\left(}\mathrm{\Omega }\hspace{0.17em}\mathit{t}\mathrm{\right)}}{\mathit{\partial \mu }}\mathit{,}\\ & & {\mathit{D}}_{\mathit{\rho }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}{\mathit{\rho }}^{\mathrm{2}}}\frac{\mathit{\partial p}}{\mathit{\partial \mu }}\right)\mathit{,}\\ & & {\mathit{D}}_{\mathit{\rho \mu }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{\mathit{r}{\mathit{\rho }}^{\mathrm{2}}}\frac{\mathit{\partial p}}{\mathit{\partial \mu }}\mathit{,}\\ & & {\mathit{D}}_{\mathit{\pi }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r\rho }}\frac{\mathit{\partial \rho }}{\mathit{\partial \mu }}\right)\mathit{,}\\ & & {\mathit{D}}_{\mathit{\pi \mu }}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\left(\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathit{\partial \rho }}{\mathit{\partial \mu }}\mathrm{-}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\right)\mathit{,}\\ & & {\mathit{D}}_{\mathrm{\Phi }\mathit{\mu }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}{\mathit{r}}_{\mathit{\mu }}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{·}\end{array}$(A.7)

A.3. Toroidal motion

The curl of the motion equation is (Unno et al. 1989) $\begin{array}{ccc}& & \left[\left(\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{+}\mathrm{\Omega }\frac{\mathit{\partial }}{\mathit{\partial \phi }}\right){\mathit{\omega }}_{\mathit{i}}^{\mathrm{\prime }}\right]{e}\mathit{i}\mathrm{+}\left({\mathit{\omega }}^{\mathrm{\prime }}\mathrm{·}\mathrm{\nabla }\mathbf{\Omega }\right)\mathit{r}\sin \mathit{\theta }{e}\mathit{\phi }\mathrm{+}\left({{v}}^{\mathbf{\prime }}\mathrm{·}\mathrm{\nabla }\right)\mathrm{2}\mathbf{\Omega }\\ & & \mathrm{-}\left(\mathrm{2}\mathbf{\Omega }\mathrm{·}\mathrm{\nabla }\right){{v}}^{\mathbf{\prime }}\mathrm{+}\left(\mathrm{\nabla }\mathrm{·}{{v}}^{\mathbf{\prime }}\right)\mathrm{2}\mathbf{\Omega }\mathrm{+}\mathrm{\nabla }\mathrm{\times }\left(\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{\nabla }\mathit{p}\right)\mathrm{=}\mathbf{0}\mathit{,}\end{array}$(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 $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ {\mathit{a}}_{\mathit{\zeta }}\end{array}$): $\begin{array}{ccc}\mathrm{0}& \mathrm{=}& \mathrm{i}\mathrm{(}\mathit{m}\mathrm{\Omega }\mathrm{+}\mathit{\sigma }\mathrm{)}\left[\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\mathrm{\left(}\sin \mathit{\theta }{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{-}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}}{\mathit{\partial \varphi }}\right]\\ & & \mathrm{+}\mathrm{i}\mathrm{(}\mathit{m}\mathrm{\Omega }\mathrm{+}\mathit{\sigma }\mathrm{)}\frac{{\mathit{r}}_{\mathit{\theta }}}{\mathit{r}}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{+}\mathrm{i}\sin \mathit{\theta }{\mathit{R}}_{\mathit{\varphi }}\hspace{0.17em}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{+}{\mathit{R}}_{\mathit{\zeta }}\hspace{0.17em}{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}\mathrm{+}\sin \mathit{\theta }{\mathit{R}}_{\mathit{\mu }}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\\ & & \mathrm{+}{\mathit{R}}_{\mathit{\zeta \mu }}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}}{\mathit{\partial \mu }}\mathrm{+}{\mathit{R}}_{\mathit{\mu \mu }}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\mathrm{\left(}\sin \mathit{\theta }{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{+}\mathrm{i}{\mathit{R}}_{\mathit{\pi }}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{\prime }}\mathrm{+}\mathrm{i}{\mathit{R}}_{\mathit{\rho }}\hspace{0.17em}{\mathit{\rho }}^{\mathrm{\prime }}\end{array}$(A.9)where $\begin{array}{ccc}& & {\mathit{R}}_{\mathit{\zeta }}\mathrm{=}\mathit{a}\hspace{0.17em}\mathit{t}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\hspace{0.17em}\mathit{c}\mathrm{+}\mathrm{\Omega }\hspace{0.17em}\mathit{t}\hspace{0.17em}\mathrm{(}\mathit{d}\mathrm{+}\mathit{f}\mathrm{)}\\ & & {\mathit{R}}_{\mathit{\mu }}\mathrm{=}\mathrm{-}\mathit{t}\hspace{0.17em}\mathit{b}\mathrm{+}\mathit{c}\mathrm{-}\mathrm{2}\mathrm{\Omega }\mathrm{+}\mathrm{\Omega }\mathrm{(}\mathit{e}\mathrm{+}\mathit{h}\mathrm{)}\mathit{t}\\ & & {\mathit{R}}_{\mathit{\zeta \mu }}\mathrm{=}\mathrm{-}\mathrm{\Omega }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\left(\mathrm{2}\mathrm{+}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\hspace{0.17em}\mathit{g}\hspace{0.17em}\mathit{t}\right)\\ & & {\mathit{R}}_{\mathit{\mu \mu }}\mathrm{=}\mathrm{\Omega }\hspace{0.17em}\mathit{t}\hspace{0.17em}\mathit{g}\\ & & {\mathit{R}}_{\mathit{\varphi }}\mathrm{=}\frac{\mathit{m}\mathrm{\Omega }\mathit{t}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}\\ & & {\mathit{R}}_{\mathit{\pi }}\mathrm{=}\mathrm{-}\frac{\mathit{m}}{\mathit{r\rho }}\frac{\mathit{\partial \rho }}{\mathit{\partial \mu }}\\ & & {\mathit{R}}_{\mathit{\rho }}\mathrm{=}\frac{\mathit{m}}{\mathit{r}{\mathit{\rho }}^{\mathrm{2}}}\frac{\mathit{\partial p}}{\mathit{\partial \mu }}\end{array}$

t is given in Eq. (A.7) and a, b, c, d, e, f and h are defined by $\begin{array}{ccc}\mathit{a}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathit{r}}{{\mathit{r}}_{\mathit{\zeta }}}\hspace{0.17em}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \zeta }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \mu }}\\ \mathit{b}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial \mu }}\\ \mathit{c}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}\left(\frac{\mathit{t}}{{\mathit{r}}_{\mathit{\zeta }}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{r}}_{\mathit{\mu }}}{\mathit{\partial \zeta }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\frac{\mathrm{2}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{\mathit{r}}\frac{\mathit{\partial }{\mathit{r}}_{\mathit{\mu }}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{2}\hspace{0.17em}\mathit{\mu }\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\right)\\ \mathit{d}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\mathrm{2}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\hspace{0.17em}\left(\mathrm{-}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\frac{\mathrm{1}}{{\mathit{r}}_{\mathit{\zeta }}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{r}}_{\mathit{\mu }}}{\mathit{\partial \zeta }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\mu }\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\right)\\ \mathit{e}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{1}}{{\mathit{r}}_{\mathit{\zeta }}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{r}}_{\mathit{\mu }}}{\mathit{\partial \zeta }}\\ \mathit{f}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2}\mathit{g}\mathrm{)}}{\mathit{r}}\hspace{0.17em}\left(\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}{\mathit{r}}_{\mathit{\mu }}\mathit{t}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{r}}_{\mathit{\mu }}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{2}\hspace{0.17em}\mathit{\mu }\hspace{0.17em}{\mathit{r}}_{\mathit{\mu }}\right)\\ \mathit{h}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\mathrm{2}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\hspace{0.17em}\mathit{g}\hspace{0.17em}\left[\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}\left(\frac{\mathrm{1}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{r}}_{\mathit{\mu }}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\mu }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\right)\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\mu }\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{\mu }}}{\mathit{r}}\right]\mathrm{·}\end{array}$

A.4. Adiabatic relation

We define quantities analogous to the Schwarzschild discriminant in the two dimensional case: $\begin{array}{ccc}{\mathit{A}}_{\mathit{\zeta }}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{1}}{{\mathrm{\Gamma }}_{\mathrm{1}}}{\left(\frac{\mathit{\partial }\ln \mathit{p}}{\mathit{\partial \zeta }}\right)}_{\mathit{\mu }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}{\left(\frac{\mathit{\partial }\ln \mathit{\rho }}{\mathit{\partial \zeta }}\right)}_{\mathit{\mu }}\\ {\mathit{A}}_{\mathit{\mu }}\hspace{0.17em}& \mathrm{=}& \end{array}$Therefore, the adiabatic relation given in Eq. (6) can be written in the following form: $\begin{array}{ccc}\mathrm{i}\mathrm{(}\mathit{m}\mathrm{\Omega }\mathrm{+}\mathit{\sigma }\mathrm{)}\hspace{0.17em}\left(\frac{\mathrm{1}}{\mathit{\rho }}\hspace{0.17em}{\mathit{\rho }}^{\mathrm{\prime }}\mathrm{-}\frac{\mathit{\rho }}{{\mathrm{\Gamma }}_{\mathrm{1}}\hspace{0.17em}\mathit{p}}\hspace{0.17em}{\mathit{\pi }}^{\mathrm{\prime }}\right)\mathrm{-}\begin{array}{}\mathrm{˜}\\ {\mathit{A}}_{\mathit{\zeta }}\end{array}{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}\mathrm{+}\sin \mathit{\theta }\hspace{0.17em}\begin{array}{}\mathrm{˜}\\ {\mathit{A}}_{\mathit{\mu }}\end{array}\hspace{0.17em}\hspace{0.17em}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(A.12)where we have introduced $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{A}}_{\mathit{\mu }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{g}\hspace{0.17em}{\mathit{A}}_{\mathit{\mu }}\hspace{0.17em}\mathit{/}\mathit{r}$, and $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{A}}_{\mathit{\zeta }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathit{A}}_{\mathit{\zeta }}\mathit{/}{\mathit{r}}_{\mathit{\zeta }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}{\mathit{r}}_{\mathit{\mu }}\hspace{0.17em}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{A}}_{\mathit{\mu }}\end{array}\mathit{/}\mathit{r}$.