© D. R. Reese et al. 2021

1. Introduction

Much effort has gone into producing realistic models of rapidly rotating stars. This includes the pioneering works by Roxburgh et al. (1965), Ostriker & Mark (1968), and Jackson (1970) and continues on in the present with various 1D codes (e.g. Palacios et al. 2003; Eggenberger et al. 2008; Marques et al. 2013) as well as 2D codes such as the one from the Evolution STEllaire en Rotation (ESTER) project (Rieutord & Espinosa Lara 2009; Espinosa Lara & Rieutord 2013; Rieutord et al. 2016). An extensive monograph on the effects of rotation on stellar structure and evolution has also recently been published (Maeder 2009). In parallel, much work has gone into calculating pulsation spectra in such models in order to interpret observations from recent space missions such as CoRoT (Baglin et al. 2009; Auvergne et al. 2009), Kepler (Borucki et al. 2009), and TESS (Ricker et al. 2015). Some of the most recent works include Lovekin & Deupree (2008), Lovekin et al. (2009), Lignières & Georgeot (2008, 2009), Ballot et al. (2010), Reese et al. (2009, 2013), and Ouazzani et al. (2015, 2017). Of these works, only Ouazzani et al. (2015) addresses pulsations in baroclinic stellar models, that is, models in which surfaces of constant pressure, temperature, or density do not coincide. This is a major ingredient of realistic models, as rotating stars are expected to be baroclinic (e.g. Zahn 1992). The work by Ouazzani et al. (2015) used stellar models from Roxburgh (2006) in which the rotation profile is imposed beforehand rather than being calculated in a self-consistent way using energy conservation. In contrast, the ESTER code deduces the rotation profile in a self-consistent way when constructing stellar models. Hence, it is important to study pulsation modes in such models.

One of the first signatures of rotation on stellar pulsations is rotational splittings, the frequency differences between consecutive modes with the same radial order and harmonic degree but different azimuthal orders. At slow rotation rates, rotational splittings can be used to invert 1D or 2D rotation profiles using a first-order perturbative approach (e.g. Deheuvels et al. 2014; Schou et al. 1998; Thompson et al. 2003). At high rotation rates, higher-order effects come into play and must be addressed before meaningful information on the rotation profile can be deduced (e.g. Soufi et al. 1998; Suárez et al. 2009). In this context, a particularly interesting quantity to investigate is the generalised rotational splitting, namely the frequency difference between prograde modes and their retrograde counterparts. In particular, Ouazzani & Goupil (2012) showed that it is possible to distinguish between third-order effects of rotation and latitudinal differential rotation in such splittings. At higher rotation rates, Reese et al. (2009) showed that such splittings are a weighted integral of the rotation profile, provided the degree of differential rotation is not too large. This would potentially provide the basis for carrying out rotation inversions in such stars. This work, however, was restricted to cylindrical rotation profiles and furthermore neglected the influence of the Coriolis force in the integrals. This raises the open questions of whether such weighted integrals can be generalised to general 2D rotation profiles, and if so, how accurate they are.

Another important consideration concerns frequency separations. Indeed, a number of recent studies have shown that the pulsation frequencies of low-degree acoustic modes of rapidly rotating stars follow an asymptotic formula. Such a formula was first explored on an empirical basis (Lignières et al. 2006; Reese et al. 2008, 2009) before being justified using ray dynamics (Lignières & Georgeot 2008, 2009; Pasek et al. 2011, 2012). Reese et al. (2017) studied theoretical pulsation spectra with realistic mode visibilities in rapidly rotating 1.8 and 2 M_⊙ stellar models based on the self-consistent field (SCF) method (Jackson et al. 2005; MacGregor et al. 2007). They showed that it may be possible, depending on the configuration, to detect the rotating counterpart to the large frequency separation, or half its value, as well as frequency spacings corresponding to multiples of the rotation rate. More recently, Mirouh et al. (2019) set up a machine learning algorithm to automatically identify to which class a given mode belongs. They went on to characterise the large frequency separation in a large set of models at different rotation rates and with different core compositions (thus mimicking the effects of stellar evolution), and showed a tight scaling relation between it and the stellar mean-density. From an observational point of view, recurrent frequency spacings have been detected in a number of δ Scuti stars (Mantegazza et al. 2012; Suárez et al. 2014; García Hernández et al. 2009, 2013; Paparó et al. 2016), including a very recent study involving interferometry, spectroscopy, and space photometry (Bouchaud et al. 2020), and interpreted as the large frequency separation or half its value. García Hernández et al. (2015) studied a number of δ Scuti pulsators in binary systems, for which independent estimates of the mass and radius are available, and have shown that this separation scales with the mean density, as expected based on the calculations in Reese et al. (2008). Ensemble asteroseismology has recently been applied to CoRoT δ Scuti stars by Michel et al. (2017) who also found regular patterns related to the large separation, although we note that Bowman & Kurtz (2018) applied a similar strategy to Kepler δ Scuti stars without the same degree of success. Finally, in the very recent work by Bedding et al. (2020), the pulsation spectra of 57 δ Scuti stars observed by TESS and three by Kepler were matched to axisymmetric ℓ = 0 and ℓ = 1 modes from non-rotating models via echelle diagrams. Such modes were shown to be relatively invariant as a function of rotation rate up to ∼0.5Ω_K using pulsation calculations in SCF models, apart from a scale factor related to the mean density, thus justifying the use of non-rotating models.

However, it is unclear to what extent the asymptotic formula would hold in the presence of discontinuities within the stellar model. Based on results previously obtained in non-rotating models with sharp gradients (e.g. Monteiro et al. 1994), one can expect the asymptotic formula to still apply albeit with a supplementary oscillatory component. However, it is not clear how strong this component is, how it behaves in the presence of rapid rotation, and whether it can hinder the interpretation of observed oscillation spectra in rapidly rotating stars, as discussed in Breger et al. (2012). The recent works by Bouabid et al. (2013) and Ouazzani et al. (2017) have shown, using the traditional approximation and full 2D pulsation calculations, respectively, how a sharp gradient around the core of rotating γ Dor stars affects g-modes. In particular, they show the presence of a periodic component in the period spacings of such modes, analogous to what was found in the non-rotating case (Miglio et al. 2008), in agreement with observations from the Kepler mission (e.g. Van Reeth et al. 2015). Likewise, similar observations in SPB stars have also revealed an oscillatory behaviour in the period spacing of their g-modes (e.g. Pápics et al. 2017). A similar study is needed for acoustic modes in rapidly rotating stars.

In order to address the above questions, we investigate low-degree acoustic modes in rapidly rotating stellar models from the ESTER code. One of the advantages of the ESTER code is its multi-domain spectral approach, ideal for introducing discontinuities while retaining a high numerical accuracy. The pulsation modes are calculated using a multi-domain spectral version of the Two-dimensional Oscillation Program (TOP, Reese et al. 2006, 2009). The article is organised as follows: the following section describes stellar models based on the ESTER code. This is then followed by a description of the pulsation calculations as well as the variational principle, with a particular emphasis on the effects of discontinuities. Section 4 deals with generalised rotational splittings. Section 5 then goes on to describe the effects of discontinuities, both on the pulsation frequencies and on the eigenmodes. This is then followed by the conclusion.

2. Stellar models based on the ESTER code

The aim of the ESTER project is to produce and evolve self-consistent stellar models of rapidly rotating stars. Consequently, a fully 2D approach is used in order to solve the relevant fluid equations while taking into account energy conservation when modelling the stationary structure of the star. This leads to centrifugal deformation of the stellar structure, as well as more subtle effects, namely differential rotation and meridional circulation, resulting from baroclinicity. Consequently, the rotation profile depends on both the radial coordinate and colatitude, and the isobars, isochores, and isotherms are distinct.

In terms of microphysics, it is possible to apply various equations of state (EOS) in ESTER. These include: the ideal gas law with or without radiation pressure, the OPAL EOS (Rogers & Nayfonov 2002), and FreeEOS¹ (Irwin 2012). In what follows, we applied the ideal gas law (without radiation pressure) in the discontinuous models (see below) and one of the continuous models, in order to avoid introducing numerical errors coming from a tabulated EOS, and the OPAL EOS in the other continuous model for the sake of realism. In terms of opacities, there are two options currently implemented: Kramer’s opacities and OPAL opacities (Iglesias & Rogers 1996). We used Kramer’s opacities in conjunction with the ideal gas law to rely entirely on analytical expressions thus reducing numerical errors, and OPAL opacities with the OPAL EOS for the sake of realism and consistency. Models with Kramer’s opacity have significantly larger radii and hence lower mean densities.

Currently, the ESTER code has some limitations. Firstly, it is unable to simulate convective envelopes. Indeed, applying a strong entropy diffusion as is done in the convective core is too approximate for the envelope. Various numerical difficulties have so far prevented the code from converging to a convective solution in such regions. Accordingly, ESTER is currently not suitable for stars with masses below ∼1.6 M_⊙. Secondly, the ESTER code is unable to simulate time evolution using a full chain of nuclear reactions. However, it is possible to alter the core composition in order to mimic the effects of stellar evolution or to include a rudimentary implementation of hydrogen combustion.

From a numerical point of view, the star is divided into multiple domains in the radial direction. There are two main reasons for doing this. First, this allows us to overcome the limitations inherent to using a spectral approach with its imposed collocation grid. In particular, it enables us to have a high resolution near the surface where it is needed. The second reason is that one can place a discontinuity between two domains without losing spectral accuracy. This goes hand in hand with the use of a dedicated coordinate system, (ζ, θ, ϕ), where ζ is a surface-fitting radial coordinate that is constant across the stellar surface and across the surfaces which delimit the boundary between consecutive domains (see Rieutord et al. 2016, for more details).

In this study, we use 2 M_⊙ stellar models at 70% of the Keplerian break-up rotation rate. We note that this value is not too far from the rotation rates of Rasalhague (α Oph) for which Ω ∼ 0.64Ω_K (see e.g. Deupree 2011; Mirouh et al. 2017, and references therein) and Altair for which Ω = 0.74Ω_K (Bouchaud et al. 2020), two well-studied δ Scuti stars with photometric observations from the space missions MOST and WIRE respectively. These models use a spectral approach based on Chebyshev polynomials in the radial direction, and spherical harmonics in the horizontal directions. The radial direction is subdivided into eight domains, the resolution in each domain being 30, 55, 45, 40, 40, 50, 70, and 70, that is, a total of 400 radial points. In the horizontal directions, 22 or 32 points on half of a Gauss-Legendre collocation grid are used depending on the model, thus corresponding to 22 or 32 spherical harmonics with even ℓ values. Density discontinuities are achieved by modifying the hydrogen content abruptly. Table 1 gives the characteristics of the five models (Mreal, M, M6, M7, and M7b) involved in this study. In all of the models, Z = 0.02 everywhere. Figure 1 shows where the discontinuity is located in model M6, and Fig. 2 gives the density and sound velocity profiles in the M, M6, and M7 models.

Fig. 1. Meridional cross-section of model M6 showing where the discontinuity is located. The other models have a discontinuity closer to the surface.

Fig. 2. Density (upper panel) and sound velocity (lower panel) profiles along the equator in models M, M6, and M7.

Models Mreal and M are our most realistic models, and serve as a reference, since they do not feature ad hoc discontinuities. Models M6 and M7 include a drop in density near the surface for two different radii, while M7b includes an increase in density near the surface. In realistic models, such as Mreal, such discontinuities are not expected. Instead, more subtle phenomena, such as dips in the Γ₁ profile due to the hydrogen and helium ionisation zones, occur near the stellar surface and can lead to a glitch pattern in the frequencies. However, it is still useful to test models with discontinuities as they exaggerate the phenomena we wish to study, namely acoustic glitches, and should thus make it easier to detect its signature in the pulsation spectrum. Furthermore, one can easily modify the different parameters related to the discontinuity such as depth and intensity in order to study its impact on the frequencies. Finally, the lack of radiation pressure in these models leads to flat Γ₁ profiles meaning that the only glitch signatures expected are those arising from the discontinuities, thus simplifying the subsequent analysis. Nonetheless, the realistic model also allows us to test acoustic glitches related to the Γ₁ profile.

We do note that stars can have a discontinuity around the core due to the depletion of hydrogen by nuclear reactions. However, acoustic modes are sensitive to the near-surface layers of stars and are thus not the most suitable for studying such discontinuities. This is particularly true of island modes as the ray trajectory orbits around which they are concentrated remain far away from the convective core for Ω ≳ 0.2Ω_K (at least for the models in this study). In order to probe such discontinuities, it is more useful to look at gravity-mode glitches (e.g. Miglio et al. 2008; Ouazzani et al. 2017), which is beyond the scope of the present article.

In model M7b, the denser layer is on top. At first sight, this may seem unrealistic, but density inversions can occur in the near-surface layers of stars such as the one shown in Fig. 3 for a 2 M_⊙ non-rotating main sequence model from grid B of Marques et al. (2008). Such density inversions typically occur as a result of a sharp temperature drop in low density regions of the star (Marques, private communication). Furthermore, including a model with a denser layer on top allows us to test the Snell-Descartes law in different situations.

Fig. 3. Near-surface density profile in a 2 M⊙ non-rotating main sequence model from grid B of Marques et al. (2008).

3. Pulsation calculations

The pulsation modes are calculated using the Two-dimensional Oscillation Program (TOP, Reese et al. 2006, 2009). This program fully takes into account the centrifugal deformation and has been set up to apply a multi-domain spectral approach, in accordance with the models from ESTER. The next subsections describe the set of pulsation equations, the interface conditions that apply between different domains, the boundary conditions, and the numerical approach.

3.1. Pulsation equations

The following set of equations are used to calculate pulsation modes. They are, respectively, the continuity equation, Euler’s equation, the adiabatic relation, and Poisson’s equation:

$$\begin{array}{cc}& 0=\frac{\delta \rho }{{\rho }_0}+\mathbf{\nabla }·\mathit{\xi },\hfill \end{array}$$(1)

$$\begin{array}{cc}& 0={(\omega +m\mathrm{\Omega })}^2\mathit{\xi }-2i(\omega +m\mathrm{\Omega })\mathbf{\Omega }\times \mathit{\xi }-\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathit{\xi })\hfill \\ \hfill & -\mathit{\xi }·\mathbf{\nabla }\left(s{\mathrm{\Omega }}^2{\mathit{e}}_s\right)-\frac{P_0}{{\rho }_0}\mathbf{\nabla }\left(\frac{\delta p}{P_0}\right)+\frac{\mathbf{\nabla }{P}_0}{{\rho }_0}(\frac{\delta \rho }{{\rho }_0}-\frac{\delta p}{P_0})-\mathbf{\nabla }\mathrm{\Psi }\hfill \\ \hfill & +\mathbf{\nabla }\left(\frac{\mathit{\xi }·\mathbf{\nabla }{P}_0}{{\rho }_0}\right)+\left[\frac{(\mathit{\xi }·\mathbf{\nabla }{P}_0)\mathbf{\nabla }{\rho }_0-(\mathit{\xi }·\mathbf{\nabla }{\rho }_0)\mathbf{\nabla }{P}_0}{{\rho }_0^2}\right],\hfill \end{array}$$(2)

$$\begin{array}{cc}& 0=\frac{\delta p}{P_0}-{\mathrm{\Gamma }}_1\frac{\delta \rho }{{\rho }_0},\hfill \end{array}$$(3)

$$\begin{array}{cc}& 0=\mathrm{\Delta }\mathrm{\Psi }-\mathrm{\Lambda }({\rho }_0\frac{\delta \rho }{{\rho }_0}-\mathit{\xi }·\mathbf{\nabla }{\rho }_0),\hfill \end{array}$$(4)

where quantities with the subscript “0” are equilibrium quantities, those with “δ” in front Lagrangian perturbations, ρ the density, P the pressure, Ψ the Eulerian gravitational potential perturbation, ξ the Lagrangian displacement, Ω the rotation profile (which depends on ζ, the surface-fitting radial coordinate, and θ, the co-latitude), Λ = 4πG, G the gravitational constant, s the distance from the rotation axis, and e_s the associated unit vector. The term in square brackets (last line of Eq. (2)) does not cancel, since the stellar model is not barotropic. The Lagrangian density perturbation is eliminated in favour of the Lagrangian pressure perturbation using Eq. (3). In the above set of equations, we have neglected the meridional circulation, given that it is expected to have a negligible effect on the pulsation modes.

3.2. Non-dimensionalisation

The following reference length, pressure, and density scales are used:

$$\begin{array}{c}\hfill R_{\mathrm{ref}}=R_{\mathrm{eq}},P_{\mathrm{ref}}=\frac{G{M}^2}{R_{\mathrm{eq}}^4},{\rho }_{\mathrm{ref}}=\frac{M}{R_{\mathrm{eq}}^3},\end{array}$$(5)

where R_eq is the equatorial radius and M the mass. As a result of this choice of reference scales, the frequencies are non-dimensionalised by the inverse of the dynamic time scale:

$$\begin{array}{c}\hfill {\omega }_{\mathrm{ref}}={\mathrm{\Omega }}_K=\sqrt{\frac{\mathit{GM}}{R_{\mathrm{eq}}^3}}.\end{array}$$(6)

Using this non-dimensionalisation leads to the same set of pulsation equations as previously (Eqs. (1)–(4)) except that Λ is now equal to 4π.

3.3. Interface conditions

Given that ESTER models are calculated over multiple domains, interface conditions are needed to describe the relation between various quantities on either side of the different boundaries. Furthermore, care is needed when expressing these conditions given that some of the models contain discontinuities. The first condition is simply that the fluid domain is continuous. In other words, the deformation of the boundary must be the same on either side. This yields the following first order expression:

$$\begin{array}{c}\hfill {\mathit{\xi }}_{-}·\mathit{n}={\mathit{\xi }}_{+}·\mathit{n},\end{array}$$(7)

where n is the normal to the unperturbed surface, and the subscripts “−” and “+” denote quantities below and above the boundary. This condition does allow the fluid to “slip” along the boundary. A more detailed derivation is given in Appendix A.1.

A second condition is that the pressure remains continuous across the perturbed boundary. This condition is simply expressed as follows (see Appendix A.2):

$$\begin{array}{c}\hfill \delta p_{-}=\delta p_{+}.\end{array}$$(8)

The third condition is the continuity of Ψ and its gradient across the perturbed boundary. This is enforced by the following conditions (see Appendix A.3):

$$\begin{array}{cc}& {\mathrm{\Psi }}_{-}={\mathrm{\Psi }}_{+},\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\partial }_{\zeta }{\mathrm{\Psi }}_{-}+\frac{\mathrm{\Lambda }{\rho }_{-}{\zeta }^2{r}_{\zeta }}{r^2+r_{\theta }^2}{\xi }^{\zeta }={\partial }_{\zeta }{\mathrm{\Psi }}_{+}+\frac{\mathrm{\Lambda }{\rho }_{+}{\zeta }^2{r}_{\zeta }}{r^2+r_{\theta }^2}{\xi }^{\zeta }.\hfill \end{array}$$(10)

3.4. Boundary conditions

As usual, various boundary conditions are needed to complete the system. In the centre, the solutions need to be regular. At the surface, we apply the simple mechanical boundary condition δp = 0. We note that in Reese et al. (2013), a more complex condition was imposed in order to have a non-zero value for δT/T₀ at the surface, useful for mode visibility calculations. However, with such a condition, the pulsation equations do not derive from a variational principle. Here, since we are seeking to obtain accurate frequencies, we prefer the simpler boundary condition (δp = 0), so that we can then apply the variational principle as a supplementary check on the accuracy. Finally, the gravitational potential must match a vacuum potential at infinity. This is achieved by extending the gravitational potential, thanks to Eqs. (9) and (10), into an external domain which encompasses the star and has a spherical outer boundary. The outer boundary condition is then (see Reese et al. 2006):

$$\begin{array}{c}\hfill \frac{1}{r_{\zeta }}\frac{\mathrm{d}{\mathrm{\Psi }}_m^{\ell }}{\mathrm{d}\zeta }+\frac{\ell +1}{r_{\mathrm{ext}}}{\mathrm{\Psi }}_m^{\ell }=0,\end{array}$$(11)

where r_ζ = ∂_ζr = 1 − ε, r_ext = 2, and where we have used a harmonic decomposition of Ψ, ℓ being the spherical harmonic degree.

3.5. Numerical approach

The above system of equations, as well as boundary and interface conditions, are discretised using the spherical harmonic basis for the angular coordinates (θ, ϕ), and using Chebyshev polynomials in the radial direction. This leads to a generalised matrix eigenvalue problem of the form Ax = λBx. This problem is modified using a shift-invert approach to target frequencies around a given shift, σ, before being solved through the Arnoldi-Chebyshev approach (e.g. Braconnier 1993; Chatelin 2012).

The multi-domain spectral approach used in the radial direction leads to matrices A and B which are block tri-diagonal. The matrix A − σB (which intervenes in the shift-invert approach) can be efficiently factorised using successive factorisations of the diagonal blocks (including a corrective term from the non-diagonal blocks).

3.6. Accuracy of the pulsation calculations

3.6.1. Various numerical resolutions

In order to check the accuracy of the frequencies, it is useful to recalculate the pulsation modes using different radial resolutions or numbers of spherical harmonics. Accordingly, we recalculated 28–30 axisymmetric (m = 0) modes in three of the models, using various resolutions. Table 2 gives the maximum relative differences on the pulsation frequencies. The first column corresponds to a 50% increase in the number of spherical harmonics in the pulsation calculations, that is, the pulsation modes are calculated with N_θ = 60 rather N_θ = 40 spherical harmonics. The second column corresponds to a ∼50% increase of the radial resolution in the pulsation calculations (after having interpolated the model). Specifically, the resolutions in the eight domains are 45, 85, 70, 60, 60, 75, 105, and 105, that is, a total of N_r = 605 points. Finally, the third column corresponds to ∼50% increase of the radial resolution both in the model (that is, the model is calculated with ESTER using an increased radial resolution rather than being interpolated) and pulsation calculations.

Two trends can be seen in Table 2. First, modifying the resolution in both the model and the pulsation calculations has a greater impact than only modifying the resolution of the pulsation calculations. This is expected as the higher resolution will be taken into account in the ESTER convergence process when calculating the model in the former case. We note that no value is provided in the last column of Table 2 for model M6 since ESTER was unable to converge in that situation. Secondly, modifying the number of spherical harmonics in the pulsation calculations has a greater impact than modifying the radial resolution. This probably simply illustrates the need for a sufficient harmonic resolution to resolve the intricate island mode geometry, particularly in model M6. Overall, these differences remain small (especially bearing in mind these are the maximal differences), except possibly for the differences related to the harmonic resolution in model M6.

3.6.2. Variational principle

Another way of checking the accuracy of the pulsation calculations consists in comparing the numerical frequencies with those obtained using a variational formula. Such a formula is an integral relation between the frequencies and their associated eigenfunctions. According to the variational principle, the error on the variational frequency scales as the square of the error on the eigenfunctions (e.g. Christensen-Dalsgaard 1982). A general formulation of the variational principle in differentially rotating bodies has previously been obtained by Lynden-Bell & Ostriker (1967). However, the formulation of some of the terms, notably the use of Green’s theorem for the gravitational potential, is not the most suitable for numerical implementation. Previous, numerically-friendly expressions similar to those in Unno et al. (1989), have been obtained in Reese et al. (2006) and Reese et al. (2009), but these expressions were only obtained for uniform or cylindrical rotation profiles, assumed that the star is barotropic, and did not include the effects of discontinuities. In Appendix B, we give a full derivation for baroclinic models with 2D rotation profiles and discontinuities. The final expression is:

$$\begin{array}{cc}& 0={\displaystyle \sum _i\underset{V_i}{\int }\{{(\omega +m\mathrm{\Omega })}^2{\rho }_0\mathit{\xi }·{\mathit{\eta }}^{\ast }-2i(\omega +m\mathrm{\Omega }){\rho }_0\mathbf{\Omega }·(\mathit{\xi }\times {\mathit{\eta }}^{\ast })}\hfill \\ \hfill & -{\rho }_0(\mathbf{\Omega }·\mathit{\xi })(\mathbf{\Omega }·{\mathit{\eta }}^{\ast })+{\rho }_0{\mathrm{\Omega }}^2\mathit{\xi }·{\mathit{\eta }}^{\ast }-{\mathit{\eta }}^{\ast }·[\mathit{\xi }·\mathbf{\nabla }\left(\mathbf{\nabla }{P}_0\right)]\hfill \\ \hfill & -{\rho }_0{\mathit{\eta }}^{\ast }·[\mathit{\xi }·\mathbf{\nabla }\left(\mathbf{\nabla }{\mathrm{\Psi }}_0\right)]-\frac{{\pi }^{\ast }P}{{\mathrm{\Gamma }}_1{P}_0}+\frac{(\mathit{\xi }·\mathbf{\nabla }{P}_0)({\mathit{\eta }}^{\ast }·\mathbf{\nabla }{P}_0)}{{\mathrm{\Gamma }}_1{P}_0}\}\mathrm{d}V\hfill \\ \hfill & +{\displaystyle \sum _i{\int }_{S_i}\mathit{\xi }·(\mathbf{\nabla }{P}_0^{-}-\mathbf{\nabla }{P}_0^{+}){\mathit{\eta }}^{\ast }·\mathbf{dS}+{\int }_{V_{\infty }}\frac{\mathbf{\nabla }\mathrm{\Psi }·\mathbf{\nabla }{\mathrm{\Phi }}^{\ast }}{\mathrm{\Lambda }}\mathrm{d}V,}\hfill \end{array}$$(12)

where Λ is 4πG or 4π in the dimensionless case, V_i are the different domains over which the stellar model is continuous, S_i are the surfaces of the discontinuities (including the stellar surface), the subscripts “+” and “−” represent quantities right above and below the discontinuities, respectively ($\mathbf{\nabla }{P}_0^{+}=\mathbf{0}$ at the stellar surface), and V_∞ is infinite space (including the star).

Figure 4 shows the relative differences between the numerical and variational frequencies for our models. In each case a set of 168 modes with quantum numbers $\stackrel{\sim }{n}=19$ to 30, $\stackrel{\sim }{\ell }=0$ to 1, m = −3 to 3, was used. We recall that $\stackrel{\sim }{n}$ is the number of nodes along an island mode’s orbit whereas $\stackrel{\sim }{\ell }$ the number of nodes parallel to it. These are related to the usual quantum numbers, (n, ℓ, m), of pulsation modes in the non-rotating case via the relations $\stackrel{\sim }{n}=2n+\epsilon$ and $\stackrel{\sim }{\ell }=\frac{\ell -|m|-\epsilon }{2}$, where $\epsilon \equiv \ell +m\mathrm{mod}2\equiv \stackrel{\sim }{n}\mathrm{mod}2$ corresponds to the mode’s parity, that is, symmetry with respect to the equatorial plane (Reese 2008). As can be seen in the figure, relative differences range from 10⁻¹¹ to 10⁻⁴, apart from an outlier in model M6². This compares quite favourably with the typical accuracy obtained with space missions. For instance, Kepler observations spanned up to four years during the main mission thus leading to a Rayleigh resolution of 0.008 μHz. The frequency at maximum amplitude of δ Scutis can reach approximately 700 μHz (e.g. Bowman & Kurtz 2018) thus leading to a relative precision as low as 10⁻⁵ in the best cases. This is higher than the errors on most of the variational frequencies, except for model M6.

Fig. 4. Relative differences between numerical and variational frequencies.

The very high accuracy which is reached in a number of cases is due to the use of spectral methods. Such an accuracy was not reached straight away but rather by repeating the calculation using the numerical frequency as the shift, σ, in the second calculation and refining the solution through supplementary iterations. Factors that decrease the accuracy (even in the second calculations) are the presence of a discontinuity in the stellar model, especially if it is sharp, and the occurrence of avoided crossings³ which lead to island modes which are “polluted” by contributions from neighbouring modes.

4. Rotation profile

In order to understand the approximate effects of differential rotation on pulsation frequencies, we consider a prograde acoustic mode and its retrograde counterpart. The azimuthal orders of these modes will be denoted −m and m, respectively⁴. Furthermore, the subscript “+” will designate the prograde mode and “–” the retrograde mode. The variational principle can be expressed in the following approximate form for these two modes:

$$\begin{array}{c}\hfill 0\simeq ({\omega }_{\pm }\mp |m|{\mathrm{\Omega }}_{\pm }^{\mathrm{eff}}{)}^2+2({\omega }_{\pm }\mp |m|{\mathrm{\Omega }}_{\pm }^{\mathrm{eff}}){\mathcal{C}}_{\pm }+{\mathrm{rest}}_{\pm },\end{array}$$(13)

where we have used the following definitions/approximations:

$$\begin{array}{cc}& {\mathrm{\Omega }}_{\pm }^{\mathrm{eff}}=\frac{{\int }_V\mathrm{\Omega }{\rho }_0{\Vert {\mathit{\xi }}_{\pm }\Vert }^2\mathrm{d}V}{{\int }_V{\rho }_0{\Vert {\mathit{\xi }}_{\pm }\Vert }^2\mathrm{d}V},\hfill \end{array}$$(14)

$$\begin{array}{cc}& {({\mathrm{\Omega }}_{\pm }^2)}^{\mathrm{eff}}=\frac{{\int }_V{\mathrm{\Omega }}^2{\rho }_0{\Vert {\mathit{\xi }}_{\pm }\Vert }^2\mathrm{d}V}{{\int }_V{\rho }_0{\Vert {\mathit{\xi }}_{\pm }\Vert }^2\mathrm{d}V}\simeq {\left({\mathrm{\Omega }}_{\pm }^{\mathrm{eff}}\right)}^2,\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\mathcal{C}}_{\pm }=\frac{i{\int }_V{\rho }_0\mathbf{\Omega }·({\mathit{\xi }}_{\pm }^{\ast }\times {\mathit{\xi }}_{\pm })\mathrm{d}V}{{\int }_V{\rho }_0{\parallel {\mathit{\xi }}_{\pm }\parallel }^2\mathrm{d}V},\hfill \end{array}$$(16)

$$\begin{array}{cc}& {\left({\mathcal{C}}_{\pm }{\mathrm{\Omega }}_{\pm }\right)}^{\mathrm{eff}}=\frac{i{\int }_V{\rho }_0\mathrm{\Omega }\mathbf{\Omega }·({\mathit{\xi }}_{\pm }^{\ast }\times {\mathit{\xi }}_{\pm })\mathrm{d}V}{{\int }_V{\rho }_0{\parallel {\mathit{\xi }}_{\pm }\parallel }^2\mathrm{d}V}\simeq {\mathcal{C}}_{\pm }{\mathrm{\Omega }}_{\pm }^{\mathrm{eff}}.\hfill \end{array}$$(17)

If the two modes are of sufficiently high frequency so that the Coriolis force only has a small impact, and if the rotation profile is not too differential, then the two modes will be close to symmetric. This means that by taking the difference between Eq. (13) applied to the prograde mode, and the same equation applied to the retrograde mode, the terms “rest₊” and “rest₋” nearly cancel. Neglecting the difference between these two terms leads to the following equation:

$$\begin{array}{cc}& ({\omega }_{+}-|m|{\mathrm{\Omega }}_{+}^{\mathrm{eff}}{)}^2+2({\omega }_{+}-|m|{\mathrm{\Omega }}_{+}^{\mathrm{eff}}){\mathcal{C}}_{+}\hfill \\ \hfill & \simeq ({\omega }_{-}+|m|{\mathrm{\Omega }}_{-}^{\mathrm{eff}}{)}^2+2({\omega }_{-}+|m|{\mathrm{\Omega }}_{-}^{\mathrm{eff}}){\mathcal{C}}_{-}.\hfill \end{array}$$(18)

This can be re-expressed as:

$$\begin{array}{cc}& ({\omega }_{+}-|m|{\mathrm{\Omega }}_{+}^{\mathrm{eff}}{)}^2[1+\frac{2{\mathcal{C}}_{+}}{({\omega }_{+}-|m|{\mathrm{\Omega }}_{+}^{\mathrm{eff}})}]\hfill \\ \hfill & \simeq ({\omega }_{-}+|m|{\mathrm{\Omega }}_{-}^{\mathrm{eff}}{)}^2[1+\frac{2{\mathcal{C}}_{-}}{({\omega }_{-}+|m|{\mathrm{\Omega }}_{-}^{\mathrm{eff}})}].\hfill \end{array}$$(19)

Taking the square-root of both sides and assuming ${\mathcal{C}}_{\pm }\ll ({\omega }_{\pm }\mp |m|{\mathrm{\Omega }}_{\pm }^{\mathrm{eff}})$ leads to:

$$\begin{array}{cc}& ({\omega }_{+}-|m|{\mathrm{\Omega }}_{+}^{\mathrm{eff}})[1+\frac{{\mathcal{C}}_{+}}{({\omega }_{+}-|m|{\mathrm{\Omega }}_{+}^{\mathrm{eff}})}]\hfill \\ \hfill & \simeq ({\omega }_{-}+|m|{\mathrm{\Omega }}_{-}^{\mathrm{eff}})[1+\frac{{\mathcal{C}}_{-}}{({\omega }_{-}+|m|{\mathrm{\Omega }}_{-}^{\mathrm{eff}})}].\hfill \end{array}$$(20)

This equation can finally be rearranged to yield:

$$\begin{array}{c}\hfill \frac{{\omega }_{+}-{\omega }_{-}}{2|m|}\simeq \frac{{\mathrm{\Omega }}_{+}^{\mathrm{eff}}+{\mathrm{\Omega }}_{-}^{\mathrm{eff}}}{2}+\frac{-{\mathcal{C}}_{+}+{\mathcal{C}}_{-}}{2|m|}.\end{array}$$(21)

This equation is particularly interesting because it provides a linear relation between the generalised rotational splitting, which only depends on the frequency of the modes, and the rotation profile. The weighting function that intervenes in the integral is known as the rotation kernel and only depends on the eigenfunctions. If the Coriolis force is neglected, this equation reduces to the linearised version of Eq. (32) from Reese et al. (2009).

Figure 5 shows what a typical rotation kernel will look like for an island mode. As can be seen, the rotation kernel closely follows the geometry of the island mode much like in Reese et al. (2009). Accordingly, these modes are especially sensitive to the rotation rate in this region, in particular near the surface at mid-latitudes.

Fig. 5. Rotation kernel for an island mode in model Mreal. The island mode is an m = 1 mode which is symmetric with respect to the equator and has a frequency of ν = 760.6 μHz.

In Figs. 6 and 7, we compare the generalised splittings with the right-hand sides of Eq. (21) for models M and Mreal, respectively. The latter is for a much more extensive set of modes. As can be seen, a good agreement is obtained in most cases, but there are some notable exceptions. Such exceptions typically occur for avoided crossings. Indeed, the geometry of the modes changes rapidly as a function of the rotation rate during avoided crossings thereby causing prograde modes and their retrograde counterparts to be at different parts of their avoided crossings and to have different geometric structures. Figure 8 provides an example of such modes. As a result, the terms “rest₊” and “rest₋” do not cancel each other out. This interpretation is confirmed in Table 3 which provides a detailed comparison between modes in this situation (Solutions 3 and 4) and those which are not undergoing an avoided crossing (Solutions 1 and 2). By including the difference between the terms “rest₊” and “rest₋”, it is possible to correct Eq. (21) and improve the agreement by a factor of 20 for Solutions 3 and 4. The supplementary rows in this Table also show that the approximations given in Eqs. (15) and (17) are well justified. Hence, apart from the cases involving avoided crossings, the agreement between the generalised splittings and the weighted integrals of the rotation profile (that is, the right-hand side of Eq. (21)) is excellent thus potentially providing the basis for probing the rotation profile via inversions.

Fig. 6. Comparison between generalised rotational splittings and the corresponding weighted integrals of the rotation profile for model M (see Eq. (21)).

Fig. 7. Same as Fig. 6 but for model Mreal and a more extensive set of modes.

Fig. 8. Meridional cross-sections of a prograde mode with m = −6 and its retrograde counterpart. The frequencies of these modes are, respectively, 610.2 and 358.8 μHz. These modes are involved in avoided crossings with other modes (not shown). As a result of being at different stages of the avoided crossing, their geometry is different and applying Eq. (21) yields less accurate results.

5. Acoustic glitches

We now turn our attention to pulsations in the discontinuous models and focus on acoustic glitches. We recall that glitches are regions in the star with a strong gradient or near discontinuity, which can lead to an oscillatory behaviour in the pulsation spectrum (e.g. Monteiro et al. 1994).

5.1. Frequencies

Figure 9 shows the pulsation frequencies obtained for the various models for modes with $\stackrel{\sim }{n}=19$ to 30, $\stackrel{\sim }{\ell }=0$ to 1, and m = −3 to 3. As can be seen, these frequencies follow fairly closely the asymptotic formula given in Reese et al. (2009). However, a closer look reveals irregularities in the pulsation spectra of the discontinuous models. This is brought out more clearly with the frequency separations ${\mathrm{\Delta }}_{\stackrel{\sim }{n}}={\omega }_{\stackrel{\sim }{n}+1,\stackrel{\sim }{\ell },m}-{\omega }_{\stackrel{\sim }{n},\stackrel{\sim }{\ell },m}$. In Fig. 10, we plot averaged large separations, $⟨{\mathrm{\Delta }}_{\stackrel{\sim }{n}}⟩=⟨{\omega }_{\stackrel{\sim }{n}+1,\stackrel{\sim }{\ell }}⟩-⟨{\omega }_{\stackrel{\sim }{n},\stackrel{\sim }{\ell }}⟩$, where $⟨{\omega }_{\stackrel{\sim }{n},\stackrel{\sim }{\ell }}⟩$ is the pulsation frequency averaged over the azimuthal orders m = −3 to 3. This is done in order to reduce the effects of avoided crossings which tend to be more numerous in the discontinuous models and tend to mask the frequency variations caused by the glitch. Even then, the averaged large separations in the discontinuous models are more irregular than in the continuous model. This raises the question whether these variations can be explained by glitch theory.

Fig. 9. Corotating pulsation frequencies in models M, M6,M7,M7b, and Mreal. These frequencies are fairly well described by the empirical formula from Reese et al. (2009).

Fig. 10. Averaged pseudo large frequency separations, $\langle {\mathrm{\Delta }}_{\stackrel{\sim }{n}}\rangle$, for the different models.

5.2. Glitch analysis and ray dynamics

In order to investigate the behaviour of the frequencies in a more detailed way, we carried out a simplified ray dynamics analysis. We used the following dispersion relation, valid for axisymmetric modes in the high frequency limit:

$$\begin{array}{c}\hfill {\omega }^2=c_0^2{k}^2,\end{array}$$(22)

where k is the norm of the wave-vector. A simple reflection was used at the stellar surface, rather than a more realistic but complex approach involving the cut-off frequency (e.g. Lignières & Georgeot 2009). Furthermore, we applied the Snell-Descartes refraction law at the discontinuity:

$$\begin{array}{c}\hfill \frac{sin{\vartheta }_{+}}{c_{+}}=\frac{sin{\vartheta }_{-}}{c_{-}},\end{array}$$(23)

where ϑ is the angle between the surface normal and the wave-vector, c the local sound velocity, and the subscripts “+” and “−” the upper and lower domains at the discontinuity. We neglect the partial wave reflection at the discontinuity, since we are only searching for the island mode periodic orbit. A more complete description of the ray dynamics is provided in Appendix C. Figure 11 shows the periodic orbit for island modes superimposed on an island mode in model M6. As can be seen, the orbit reproduces very well the location of the mode.

Fig. 11. Island mode with periodic orbit superimposed, in model M6. The mode is axisymmetric with a frequency of 309.2 μHz. The dotted line shows the location of the discontinuity. The acoustic travel times τT and τ1 along the ray trajectory are illustrated.

Figure 12 then shows the sound velocity, density, and perturbed pressure ($\delta p/\sqrt{P_0}$) profiles calculated along the periodic orbit, both as a function of distance along the profile and acoustic travel time. As expected, a sharp transition in wavelength occurs at the discontinuity. Furthermore, when plotted as a function of acoustic travel time, the wave takes on a nearly sinusoidal behaviour as indicated by the comparison with the simple sine curve, apart from a phase shift at the discontinuity and a variable amplitude.

Fig. 12. Sound velocity (first row), density (second row), and perturbed pressure ($\delta p/\sqrt{P_0}$, third row) profiles as a function of distance (first and second columns) and acoustic travel time (third column) along island periodic orbit. The second column is a zoom of the first column around the first discontinuity. The vertical light blue solid lines indicate the discontinuity and the vertical light blue dotted lines correspond to the equator. The dotted curve in the lower right panel is a simple sine curve with the same periodicity as the mode.

These observations provide the basis for a simple toy model which is described in Appendix D. According to this model, the frequencies are given to first order by:

$$\begin{array}{c}\hfill \omega =\frac{1}{2{\tau }_T}[n\pi +ϵsin\left(n\pi \frac{{\tau }_1}{{\tau }_T}\right)],\end{array}$$(24)

where ${\tau }_T={\int }_{\mathrm{surf}.}^{\mathrm{eq}.}\frac{\mathrm{d}r}{c}$ is the acoustic travel time from the surface to the equator along the ray path, and ${\tau }_1={\int }_{\mathrm{surf}.}^{\mathrm{disc}.}\frac{\mathrm{d}r}{c}$, the acoustic travel time from the surface to the discontinuity, as illustrated in Fig. 11. The quantity ϵ is given by the relation:

$$\begin{array}{c}\hfill ϵ=\frac{k_{-}}{k_{+}}-1=\frac{c_{+}}{c_{-}}-1\end{array}$$(25)

and is treated as a small parameter. As shown in Appendix D, even for ϵ = 0.39 (for model M7), Eq. (24) gives an accurate estimate of the glitch period and a rough idea of its amplitude. However, we do not expect the toy model to give an accurate idea of the phase of the glitch pattern on the oscillation frequencies as it would require fully treating surface effects.

Table 4 provides the acoustic travel times τ_T and τ₁ for the different models in our study. Although model Mreal is continuous, we included the τ₁ value for the He II ionisation zone. Indeed, the Γ₁ profile undergoes a dip in that region, as illustrated in Fig. 13. Based on these values, Fig. 14 compares the predictions from the toy model with the $(\stackrel{\sim }{l},m)=(0,0)$ frequencies minus a second or third-order polynomial fit in order to isolate the glitch pattern. Indeed, using the large separations rather than the frequencies would tend to amplify the impact of avoided crossings thus making it harder to see the glitch pattern. We note that a second rather than third-order polynomial fit was used for model M7b given the relatively long period of the glitch pattern which can be mimicked up to some extent by a third- or higher-order polynomials. An ad hoc phase was added to the glitch pattern from the toy model given that this model is not expected to correctly predict the phase as described above. This allows us to focus on the period and amplitude of the glitch pattern to see how accurate the predictions are.

Fig. 13. Γ1 profile in model Mreal along the island mode orbit. The stellar surface corresponds to τ = 0 s on the right side, and the intersection of the equatorial plane with the orbit to τ = 6612.9 s on the left side.

As can be seen from Fig. 14, a nice agreement is obtained for models M7b, Mreal, and to a lesser extent M7. This confirms that the toy model is able to correctly predict the periodicity of the glitch pattern, at least in some cases. The agreement on the amplitude is satisfactory for M7b but rather poor for M7. For model Mreal, an ad hoc amplitude was used for the predicted glitch pattern. Indeed, the toy model was specifically constructed for discontinuities and is therefore unable to predict the amplitude of the glitch pattern for a smoother transition such as what takes place in an ionisation zone. It is nonetheless interesting to note that the amplitude of this glitch pattern decreases at higher frequencies, as would be expected for such a transition. Model M is not expected to show a glitch pattern since it contains no discontinuities and the Γ₁ profile is very close to 5/3 throughout the star, as a result of the ideal gas equation of state. The plot shows what is likely to be a fourth order polynomial residual as expected when subtracting a third-order polynomial fit, as confirmed by the much smaller scale of the y-axis. In contrast, no agreement is found between the toy model and the glitch pattern for model M6. The reasons for this lack of agreement are not entirely understood, but we do note that most of its island modes are undergoing avoided crossings in contrast to the other models. Avoided crossings typically cause the frequencies to deviate from their asymptotic values and could therefore easily mask a glitch pattern.

Fig. 14. Frequencies of the $(\stackrel{\sim }{\ell },m)=(0,0)$ modes after subtraction of a third-order polynomial fit (or second order polynomial fit in the case M7b) versus the predicted glitch from the toy model. Each panel corresponds to a different model. An ad hoc phase, the value of which is indicated in each panel, has been added to the toy model to improve the agreement.

5.3. Pulsation mode geometry at the discontinuity

We now investigate in a detailed way the local geometric properties of the islands modes in the region where the periodic orbit intersects the discontinuity. Specifically, we check whether the wave amplitudes match the predictions from a local analysis, and whether the angle between the discontinuity and the orbit matches a numerical estimate based on the island mode.

As recalled in Appendix E, the pulsation mode including the reflected and refracted waves can locally be approximated as:

$$\begin{array}{c}\hfill {\left(\frac{\delta p}{P_0}\right)}^{\pm }=[A_1^{\pm }cos({\mathit{k}}_1^{\pm }·\mathit{x})+A_2^{\pm }cos({\mathit{k}}_2^{\pm }·\mathit{x})]exp(i\omega t),\end{array}$$(26)

where the superscripts ‘+” and ‘−” designate the upper and lower domains, respectively, ${\mathit{k}}_1^{\pm }$ and ${\mathit{k}}_2^{\pm }$ wave vectors, and where the amplitudes, $A_1^{\pm }$, $A_2^{\pm }$, are related via the relation:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{A}_1^{+}\\ A_2^{+}\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1+\eta & 1-\eta \\ 1-\eta & 1+\eta \end{array}\right]\left[\begin{array}{c}{A}_1^{-}\\ A_2^{-}\end{array}\right],\end{array}$$(27)

where

$$\begin{array}{c}\hfill \eta =\frac{{\rho }_0^{+}}{{\rho }_0^{-}}\frac{k_{\perp }^{-}}{k_{\perp }^{+}}.\end{array}$$(28)

and $k_{\perp }^{\pm }$ is the wave vector component perpendicular to the surface. When the tangential component, k_∥, is negligible in front of k_⊥, the factor η reduces to:

$$\begin{array}{c}\hfill \eta \simeq \sqrt{\frac{{\rho }_0^{+}}{{\rho }_0^{-}}}=\frac{c_0^{-}}{c_0^{+}}.\end{array}$$(29)

We then investigate several m = 0 island modes in different models to extract the amplitudes of the refracted and reflected waves and verify the above equations. We start by extracting the δp/P₀ profile as well as its horizontal and vertical gradients, ∇_∥(δp/P₀), and (∇_⊥(δp/P₀))^±, just above and below the discontinuity⁵. Since we are focusing on axisymmetric modes, the horizontal gradient is in the meridional plane – there is no component in the e_ϕ direction. Figure 15 shows a zoom on part of an island mode in model M6 and Fig. 16 shows the extracted profiles. The amplitudes of these profiles are estimated thanks to their maximum absolute values. Given that the (∇_⊥(δp/P₀))^± profiles have the opposite sign to the ∇_∥(δp/P₀) profile, these have negative amplitudes. The tangential wave vector component (which is the same above and below) is estimated by calculating the ratio between the amplitudes of ∇_∥(δp/P₀) and δp/P₀. The normal components above and below the discontinuity are obtained via the dispersion relation, thus enforcing Snell-Descartes’ law. The individual wave amplitudes are obtained by calculating appropriate linear combinations of (∇_⊥(δp/P₀))^±/k_⊥ and (∇_∥(δp/P₀))/k_∥.

Fig. 15. Zoom in on the island mode shown in Fig. 11 (in model M6). The discontinuity is shown using the dotted line, and the solid line corresponds to the island mode orbit. As can be seen, the wavelength decreases just above the discontinuity.

Fig. 16. Extracted δp/P0, ∇∥(δp/P0), and ∇⊥(δp/P0)± profiles as a function of θ along the discontinuity shown in Fig. 15. The vertical dashed line shows colatitude where the island mode periodic orbit crosses the discontinuity.

Table 5 gives the wave vector components and amplitudes for island modes in three of the models. Given that the mode amplitude is arbitrary, we normalised the amplitudes by $A_2^{\pm }$. The quantities “$A_1^{+}$(theo)” and “$A_2^{+}$(theo)” correspond to the amplitudes deduced from $A_1^{-}$ and $A_2^{-}$ via Eq. (27). Apart from the $A_1^{+}$(theo) for model M6, these values accurately reproduce the numerically obtained amplitudes, $A_1^{+}$ and $A_2^{+}$, thus showing that the relationship on amplitudes is respected.

Another comparison carried out in Table 5 is between the incidence/departure angles of the wave vectors and the predictions from ray dynamics analysis. The quantities k_∥/$k_{\perp }^{\pm }$ correspond to the numerically determined values of tanϑ^± where ϑ^± are the angles between the surface normal and the wave vector⁶. The quantities “k_∥/$k_{\perp }^{\pm }$(rays)” are determined via ray dynamics. A comparison between the two shows some discrepancies but the values remain of the same order. These differences are likely due to the limited accuracy of our approach for extracting the wave vector components, and the fact that the mode behaviour is more complex than what is predicted by ray dynamics.

Finally, as can be seen from the values of $A_1^{\pm }$, the amplitudes of the secondary waves, although smaller, is not negligible. Hence, these can be expected to affect the phase shift that occurs at the discontinuity in the primary wave and may possibly lead to increased coupling with other modes due to the modified mode geometry, thereby leading to more avoided crossings.

5.4. Frequency patterns

We now briefly address the question of whether discontinuities can adversely affect frequency patterns to the point of hindering their detection in observed stars. A tool frequently used in solar-like stars is the so-called echelle diagram (e.g. Bedding et al. 2020), in which the frequencies are plotted as a function of the frequencies modulo the large separation. Due to the nearly equidistant frequency patterns in such stars, modes with the same spherical degree line up on vertical ridges in echelle diagrams. In Fig. 17, we produce similar echelle diagrams using the pseudo large separation, ${\mathrm{\Delta }}_{\stackrel{\sim }{n}}$, and only plotting m = 0 modes for the sake of clarity. Although the discontinuities lead to a more irregular behaviour, clear ridges remain for the different $\stackrel{\sim }{\ell }$ values. Reese et al. (2014) also reached a similar conclusion using histograms of frequency differences for 3 M_⊙ models with discontinuities. Indeed, they found that the pseudo large separation could still be identified in the discontinuous models.

Fig. 17. Echelle diagrams for the axisymmetric modes in four of the models.

We also carry out a more quantitative comparison between the pulsation frequencies and a fit based on a simplified version of the asymptotic formula for island mode frequencies (e.g. Reese et al. 2009):

$$\begin{array}{c}\hfill \omega =\stackrel{\sim }{n}{\mathrm{\Delta }}_{\stackrel{\sim }{n}}+\stackrel{\sim }{\ell }{\mathrm{\Delta }}_{\stackrel{\sim }{\ell }}+m^2{\mathrm{\Delta }}_{\stackrel{\sim }{m}}-m{\mathrm{\Omega }}_{\mathrm{fit}}+\stackrel{\sim }{\alpha },\end{array}$$(30)

where ${\mathrm{\Delta }}_{\stackrel{\sim }{n}}$, ${\mathrm{\Delta }}_{\stackrel{\sim }{\ell }}$, ${\mathrm{\Delta }}_{\stackrel{\sim }{m}}$, and $\stackrel{\sim }{\alpha }$ are various parameters related to the stellar structure (Lignières & Georgeot 2009; Pasek et al. 2012), and Ω_fit an average value of the rotation rate appropriate for the set of modes under consideration. We therefore fit these parameters to reproduce the pulsation spectra of our models for the same set of modes as described in Sect. 3.6.2. Table 6 provides the root mean square differences and maximal differences between the numerical and asymptotic frequencies, normalised by the pseudo large separation, ${\mathrm{\Delta }}_{\stackrel{\sim }{n}}$. As expected, the frequencies of model M are the closest to the asymptotic formula. The mean difference in the realistic model is intermediate between the best and worst model. In terms of maximal differences, model Mreal is among the best whereas model M6 is the worst model, very likely as a result of the increased number of avoided crossings affecting the modes. In all cases, the differences are a few percent of the pseudo large separation (which itself is half the classical large separation), meaning the frequency pattern is still well-preserved and should be possible to identify with a suitable analysis.

Nonetheless, other factors may hinder finding the above frequency pattern. Indeed, the presence of chaotic modes with their own independent semi-random frequency organisation (Lignières & Georgeot 2009; Evano et al. 2019), or the lack of a clear understanding of the mechanisms responsible for mode selection and pulsation amplitudes both contribute to masking the frequency pattern associated with acoustic island modes.

6. Conclusion

In this work, we calculated, thanks to an adiabatic version of the TOP code, acoustic pulsation modes in rapidly rotating continuous and discontinuous stellar models based on the ESTER code. This allowed us to investigate various topics namely the variational principle for general 2D rotation profiles in discontinuous models, generalised rotational splittings, and acoustic glitches. Some of the important results are:

1. Generalised rotational splittings are well approximated via weighted integrals of the rotation profile using rotation kernels deduced from the variational principle, except for specific cases where avoided crossings lead to discrepancies. This raises the question as to how accurately the rotation profile can be recovered using inverse theory. In a forthcoming article, we plan to investigate this question using a variety of different rotation profiles. In this regard, the automatic mode classification algorithm described in Mirouh et al. (2019) can be used to efficiently identify pairs of prograde-retrograde modes.

2. Discontinuities alter the acoustic frequency patterns, but not to the point of preventing their detection in observed stars (especially taking into account the unrealistic nature of the discontinuities in our models), thus lending credence to recent detections of large frequency separations and ridges in echelle diagrams in δ Scuti stars (e.g. García Hernández et al. 2015; Bedding et al. 2020). The modifications to the frequency spectrum leads to glitch patterns the periodicity of which can be calculated in a simple way. Nonetheless, the presence of avoided crossings and possibly partial wave reflection at the discontinuity cause deviations from theoretical expectations in some cases. Accordingly, it may be possible to determine acoustic depths of sharp transitions using glitch patterns in observed frequencies.

In a forthcoming work, we plan to investigate acoustic pulsations of ESTER models using a non-adiabatic version of TOP. This will allow us to investigate other topics such as mode excitation and mode behaviour near the stellar surface.

Acknowledgments

The authors thank François Lignières, Vincent Prat, and Pierre Houdayer for useful and interesting discussions. DRR acknowledges the support of the French Agence Nationale de la Recherche (ANR) to the ESRR project under grant ANR-16-CE31-0007 as well as financial support from the Programme National de Physique Stellaire (PNPS) of the CNRS/INSU co-funded by the CEA and the CNES. G. M. M. acknowledges funding by the STFC consolidated grant ST/R000603/1. D. R. R., G. M. M., M. R., B. P. benefited from the hospitality of ISSI as part of the SoFAR team early on in 2018 and 2019. This work was granted access to the HPC resources of IDRIS under the allocation 2011-99992 made by GENCI (“Grand Equipement National de Calcul Intensif”) and to HPC resources of Calmip under project P0107.

References

Appendix A: Interface conditions

A.1. Domain continuity

In order to ensure that the domain remains continuous, one needs the boundary on either side to be the same. We will now consider a point on the perturbed boundary. This point will be reached by fluid parcels on either side of the boundary. The spatial coordinates of the fluid parcel just below the boundary will be given by the formula:

$$\begin{array}{c}\hfill {\mathit{r}}_{-}+{\mathit{\xi }}_{-}({\mathit{r}}_{-},t).\end{array}$$(A.1)

An analogous formula applies for the spatial coordinates just outside the boundary. This leads to the following matching condition:

$$\begin{array}{c}\hfill {\mathit{r}}_{-}+{\mathit{\xi }}_{-}({\mathit{r}}_{-},t)={\mathit{r}}_{+}+{\mathit{\xi }}_{+}({\mathit{r}}_{+},t).\end{array}$$(A.2)

This matching relation is illustrated in Fig. A.1. It is important to bear in mind that r₋ is not necessarily equal to r₊, as illustrated in the figure, since the fluid may slip along either side of the boundary. However, one can rearrange this expression as follows:

$$\begin{array}{c}\hfill {\mathit{r}}_{+}-{\mathit{r}}_{-}={\mathit{\xi }}_{-}({\mathit{r}}_{-},t)-{\mathit{\xi }}_{+}({\mathit{r}}_{+},t).\end{array}$$(A.3)

Fig. A.1. Schematic illustration showing the movement of fluid parcels along either side of a boundary that is perturbed by a pulsation mode.

Given that ξ₋ and ξ₊ are arbitrarily small, the difference r₊ − r₋ is a vector tangent to the surface. Hence, calculating the dot product of the above equation and n, the normal to the surface, cancels out the left-hand side and yields:

$$\begin{array}{c}\hfill {\mathit{\xi }}_{-}({\mathit{r}}_{-},t)·\mathit{n}={\mathit{\xi }}_{+}({\mathit{r}}_{+},t)·\mathit{n}.\end{array}$$(A.4)

The difference between ξ₊(r₊, t) and ξ₊(r₋, t) is a second order term. Hence, the above condition reduces to Eq. (7).

A.2. Condition on the pressure perturbation

The following condition ensures that the pressure remains continuous during the oscillatory movements:

$$\begin{array}{c}\hfill P_{\mathrm{Total}}^{-}({\mathit{r}}_{-}+{\mathit{\xi }}_{-},t)=P_{\mathrm{Total}}^{+}({\mathit{r}}_{+}+{\mathit{\xi }}_{+},t),\end{array}$$(A.5)

where we have kept the same notation as above and where P_Total is the total pressure (equilibrium + perturbation). This equation can be developed as follows:

$$\begin{array}{c}\hfill P_0^{-}({\mathit{r}}_{-})+\delta P_{-}({\mathit{r}}_{-},t)=P_0^{+}({\mathit{r}}_{+})+\delta P_{+}({\mathit{r}}_{+},t).\end{array}$$(A.6)

We can then use Eq. (A.2) to develop, say, the right-hand side:

$$\begin{array}{cc}& P_0^{+}({\mathit{r}}_{+})+\delta P_{+}({\mathit{r}}_{+},t)\hfill \\ \hfill & =P_0^{+}({\mathit{r}}_{-}+{\mathit{\xi }}_{-}-{\mathit{\xi }}_{+})+\delta P_{+}({\mathit{r}}_{-}+{\mathit{\xi }}_{-}-{\mathit{\xi }}_{+},t)\hfill \\ \hfill & \simeq P_0^{+}({\mathit{r}}_{-})+({\mathit{\xi }}_{-}-{\mathit{\xi }}_{+})·\mathbf{\nabla }{P}_0^{+}+\delta P_{+}({\mathit{r}}_{-},t),\hfill \end{array}$$

where we have neglected second order terms on the third line. Combining this with the previous equation leads to:

$$\begin{array}{c}\hfill \delta P_{-}({\mathit{r}}_{-},t)-\delta P_{+}({\mathit{r}}_{-},t)=({\mathit{\xi }}_{-}-{\mathit{\xi }}_{+})·\mathbf{\nabla }{P}_0^{+},\end{array}$$(A.7)

where we have used the continuity of the equilibrium pressure. If the boundary coincides with an isobar, the right-hand side of the above equation cancels out because the difference ξ₋ − ξ₊ is within the boundary. If the boundary is not an isobar, then in normal circumstances (that is, when the model is continuous), the difference ξ₋ − ξ₊ will be 0. Either way, this leads to the final condition, Eq. (8). There are, however, cases where the right-hand side may not cancel, for instance at the boundary of a convective core with a different chemical composition than the rest of the star. In such a situation, baroclinic flows are set up within the equilibrium model (Espinosa Lara & Rieutord 2013), and probably require setting a specific condition which takes these flows into account. Nonetheless, it is interesting to note that the above condition is in fact symmetric with respect to either side of the boundary. Indeed, the term $\left({\mathit{\xi }}_{-}-{\mathit{\xi }}_{+}\right)\cdot \nabla P_0^{+}$ could be replaced by $\left({\mathit{\xi }}_{-}-{\mathit{\xi }}_{+}\right)\cdot \nabla P_0^{-}$, since it only involves the gradient of the pressure along the boundary and the pressure is continuous across the boundary.

A.3. Condition on the perturbation to gravitational potential

In much the same way as the pressure, the gravitational potential and its gradient are kept continuous through the following relations:

$$\begin{array}{cc}& {\mathrm{\Psi }}_{\mathrm{Total}}^{-}({\mathit{r}}_{-}+{\mathit{\xi }}_{-},t)={\mathrm{\Psi }}_{\mathrm{Total}}^{+}({\mathit{r}}_{+}+{\mathit{\xi }}_{+},t),\hfill \end{array}$$(A.8)

$$\begin{array}{cc}& \mathbf{\nabla }{\mathrm{\Psi }}_{\mathrm{Total}}^{-}({\mathit{r}}_{-}+{\mathit{\xi }}_{-},t)=\mathbf{\nabla }{\mathrm{\Psi }}_{\mathrm{Total}}^{+}({\mathit{r}}_{+}+{\mathit{\xi }}_{+},t),\hfill \end{array}$$(A.9)

where Ψ_Total is the total gravitational potential (equilibrium + perturbation). At this point, however, we will take a different approach than above since we are dealing with the Eulerian rather than Lagrangian perturbation of the gravitational potential. Firstly, the sums r₊ + ξ₊ can be replaced by r₋ + ξ₋ or vice versa so as to have the same arguments everywhere. Therefore, in what follows we will use the generic notation r + ξ which can be arbitrarily chosen as r₋ + ξ₋ or r₊ + ξ₊. Developing both sides of both equations, making use of the continuity of the equilibrium gravitational potential and its gradient to cancel zeroth order terms, and neglecting second order terms lead to the following equations:

$$\begin{array}{cc}& {\mathrm{\Psi }}_{-}(\mathit{r},t)+\mathit{\xi }(\mathit{r},t)·\mathbf{\nabla }{\mathrm{\Psi }}_0^{-}(\mathit{r},t)\hfill \\ \hfill & ={\mathrm{\Psi }}_{+}(\mathit{r},t)+\mathit{\xi }(\mathit{r},t)·\mathbf{\nabla }{\mathrm{\Psi }}_0^{+}(\mathit{r},t),\hfill \end{array}$$(A.10)

$$\begin{array}{cc}& \mathbf{\nabla }{\mathrm{\Psi }}_{-}(\mathit{r},t)+\mathit{\xi }(\mathit{r},t)·\mathbf{\nabla }\left(\mathbf{\nabla }{\mathrm{\Psi }}_0^{-}(\mathit{r},t)\right)\hfill \\ \hfill & =\mathbf{\nabla }{\mathrm{\Psi }}_{+}(\mathit{r},t)+\mathit{\xi }(\mathit{r},t)·\mathbf{\nabla }\left(\mathbf{\nabla }{\mathrm{\Psi }}_0^{+}(\mathit{r},t)\right).\hfill \end{array}$$(A.11)

Given that ∇Ψ₀ is continuous, the first equation reduces to:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{-}={\mathrm{\Psi }}_{+}.\end{array}$$(A.12)

In tensorial notation, the left-hand side of the second equation becomes:

$$\begin{array}{cc}& {\partial }_i{\mathrm{\Psi }}_{-}{\mathit{E}}_{-}^i+{\stackrel{\sim }{\xi }}^j{\partial }_j\left({\mathit{E}}_{-}^k{\partial }_k{\mathrm{\Psi }}_0^{-}\right)\hfill \\ \hfill & ={\partial }_i{\mathrm{\Psi }}_{-}{\mathit{E}}_{-}^i+{\stackrel{\sim }{\xi }}^j[{\partial }_j\left({\mathit{E}}_{-}^k\right){\partial }_k{\mathrm{\Psi }}_0^{-}+{\mathit{E}}_{-}^k{\partial }_{\mathit{jk}}^2{\mathrm{\Psi }}_0^{-}]\hfill \\ \hfill & ={\partial }_i{\mathrm{\Psi }}_{-}{\mathit{E}}_{-}^i+{\stackrel{\sim }{\xi }}^j{\partial }_j\left({\mathit{E}}_{-}^k\right){\partial }_k{\mathrm{\Psi }}_0^{-}+{\stackrel{\sim }{\xi }}^j{\partial }_{\mathit{ij}}^2{\mathrm{\Psi }}_0^{-}{\mathit{E}}_{-}^i,\hfill \end{array}$$

where E_i is the natural basis, $\stackrel{\sim }{\xi }$ the components of ξ over that basis, and Eⁱ the dual basis. Calculating the dot product of the above equation with ${\mathit{E}}_i^{-}$ yields:

$$\begin{array}{c}\hfill {\partial }_i{\mathrm{\Psi }}_{-}-{\stackrel{\sim }{\xi }}^j{\mathrm{\Gamma }-}_{\mathit{ij}}^k{\partial }_k{\mathrm{\Psi }}_0^{-}+{\stackrel{\sim }{\xi }}^j{\partial }_{\mathit{ij}}^2{\mathrm{\Psi }}_0^{-},\end{array}$$(A.13)

where ${\mathrm{\Gamma }}_{ij}^k={\partial }_i\left({\mathit{E}}_j\right)\cdot {\mathit{E}}^k=-{\partial }_i\left({\mathit{E}}^k\right)\cdot {\mathit{E}}_j$ is the Christoffel symbol of the second kind. With our choice of mapping, only ${\mathrm{\Gamma }}_{\zeta \zeta }^{\zeta }$ is discontinuous across the boundary, hence the notation ${\mathrm{\Gamma }-}_{\zeta \zeta }^{\zeta }$ and ${\mathrm{\Gamma }+}_{\zeta \zeta }^{\zeta }$. All of the other geometric quantities (E_i, Eⁱ, ${\mathrm{\Gamma }}_{ij}^k$ with (i, j, k)≠(ζ, ζ, ζ)) are continuous. Inserting this expression into the left-hand side of Eq. (A.11) and a similar expression in the right-hand side, and simplifying out continuous terms (geometric, Ψ₀, and ∇Ψ₀) yields the following three relations:

$$\begin{array}{cc}& {\partial }_{\zeta }{\mathrm{\Psi }}_{-}+({\partial }_{\zeta \zeta }^2{\mathrm{\Psi }}_0^{-}-{\mathrm{\Gamma }-}_{\zeta \zeta }^{\zeta }{\partial }_{\zeta }{\mathrm{\Psi }}_0){\stackrel{\sim }{\xi }}^{\zeta }\hfill \\ \hfill & ={\partial }_{\zeta }{\mathrm{\Psi }}_{+}+({\partial }_{\zeta \zeta }^2{\mathrm{\Psi }}_0^{+}-{\mathrm{\Gamma }+}_{\zeta \zeta }^{\zeta }{\partial }_{\zeta }{\mathrm{\Psi }}_0){\stackrel{\sim }{\xi }}^{\zeta },\hfill \end{array}$$(A.14)

$${\partial }_{\theta }{\Psi }_{-}={\partial }_{\theta }{\Psi }_{+},$$(A.15)

$${\partial }_{\varphi }{\Psi }_{-}={\partial }_{\varphi }{\Psi }_{+},$$(A.16)

where we have made use of the fact that ${\partial }_{ij}^2{\mathrm{\Psi }}_0$ is continuous if (i, j)≠(ζ, ζ). One will in fact notice that the latter two equations are also a direct consequence of Eq. (A.12).

At this point, it is useful to introduce Poisson’s equation in tensorial notation:

$$\begin{array}{cc}& \mathrm{\Lambda }{\rho }_0=\mathrm{\Delta }{\mathrm{\Psi }}_0=\mathbf{\nabla }·\mathbf{\nabla }{\mathrm{\Psi }}_0={\partial }_i\left({\mathit{E}}^j{\partial }_j{\mathrm{\Psi }}_0\right)·{\mathit{E}}^i\hfill \\ \hfill & =g^{\mathit{ij}}{\partial }_{\mathit{ij}}^2{\mathrm{\Psi }}_0+{\partial }_i\left({\mathit{E}}^j\right)·{\mathit{E}}^i{\partial }_j{\mathrm{\Psi }}_0\hfill \\ \hfill & =g^{\mathit{ij}}{\partial }_{\mathit{ij}}^2{\mathrm{\Psi }}_0+[{\partial }_i\left({\mathit{E}}^j\right)·{\mathit{E}}_k][{\mathit{E}}^i·{\mathit{E}}^k]{\partial }_j{\mathrm{\Psi }}_0\hfill \\ \hfill & =g^{\mathit{ij}}{\partial }_{\mathit{ij}}^2{\mathrm{\Psi }}_0-{\mathrm{\Gamma }}_{\mathit{ik}}^j{g}^{\mathit{ik}}{\partial }_j{\mathrm{\Psi }}_0\hfill \\ \hfill & =g^{\zeta \zeta }({\partial }_{\zeta \zeta }^2{\mathrm{\Psi }}_0-{\mathrm{\Gamma }}_{\zeta \zeta }^{\zeta }{\partial }_{\zeta }{\mathrm{\Psi }}_0)+R,\hfill \end{array}$$(A.17)

where Λ = 4πG or 4π in the dimensional or dimensionless case, respectively, g^ij = Eⁱ ⋅ E^j is the contravariant components of the metric tensor, and R is a sum of terms which are continuous across the boundary. This last expression can then be used to simplify Eq. (A.13):

$$\begin{array}{c}\hfill {\partial }_{\zeta }{\mathrm{\Psi }}_{-}+\frac{\mathrm{\Lambda }{\rho }_0^{-}-R_{-}}{g^{\zeta \zeta }}{\stackrel{\sim }{\xi }}^{\zeta }={\partial }_{\zeta }{\mathrm{\Psi }}_{+}+\frac{\mathrm{\Lambda }{\rho }_0^{+}-R_{+}}{g^{\zeta \zeta }}{\stackrel{\sim }{\xi }}^{\zeta }.\end{array}$$(A.18)

The terms R₋ and R₊ cancel out since R is continuous across the boundary. The remaining equation is then

$$\begin{array}{c}\hfill {\partial }_{\zeta }{\mathrm{\Psi }}_{-}+\frac{\mathrm{\Lambda }{\rho }_0^{-}{\zeta }^2{r}_{\zeta }}{r^2+r_{\theta }^2}{\xi }^{\zeta }={\partial }_{\zeta }{\mathrm{\Psi }}_{+}+\frac{\mathrm{\Lambda }{\rho }_0^{+}{\zeta }^2{r}_{\zeta }}{r^2+r_{\theta }^2}{\xi }^{\zeta },\end{array}$$(A.19)

where we have introduced ξ^ζ, the ζ component of ξ on the alternate basis (see, e.g. Eq. (31) of Reese et al. 2006). At this point, it is useful to recall that ${\stackrel{\sim }{\xi }}^{\zeta }$ and hence ξ^ζ are continuous across the boundary (see Eq. (A.4)). Hence, using r₋ + ξ₋ or r₊ + ξ₊ in Eq. (A.10) leads to the same results.

Appendix B: Variational principle

B.1. General formula

In order to derive the variational formula which relates pulsation frequencies and their associated eigenfunctions, we start by calculating the dot product between Euler’s equation (Eq. (2)) and the product of the equilibrium density and the complex conjugate of a second displacement field, η^*, which at this point can be different from ξ, and integrate the total over the stellar volume, V:

$$\begin{array}{cc}& 0={\displaystyle \underset{V}{\int }\{{(\omega +m\mathrm{\Omega })}^2{\rho }_0\mathit{\xi }·{\mathit{\eta }}^{\ast }\underset{I}{\underbrace{-2i(\omega +m\mathrm{\Omega }){\rho }_0(\mathbf{\Omega }\times \mathit{\xi })·{\mathit{\eta }}^{\ast }}}}\hfill \\ \hfill & \underset{\mathit{II}}{\underbrace{-{\rho }_0[\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathit{\xi })]·{\mathit{\eta }}^{\ast }}}\underset{\mathit{III}}{\underbrace{-{\rho }_0{\mathit{\eta }}^{\ast }·[\mathit{\xi }·\mathbf{\nabla }\left(s{\mathrm{\Omega }}^2{\mathit{e}}_s\right)]}}\hfill \\ \hfill & \underset{\mathit{IV}}{\underbrace{-P_0{\mathit{\eta }}^{\ast }·\mathbf{\nabla }\left(\frac{\delta p}{P_0}\right)}}\underset{V}{\underbrace{+{\mathit{\eta }}^{\ast }·\mathbf{\nabla }{P}_0(\frac{\delta \rho }{{\rho }_0}-\frac{\delta p}{P_0})}}\hfill \\ \hfill & \underset{\mathit{VI}}{\underbrace{-{\rho }_0{\mathit{\eta }}^{\ast }·\mathbf{\nabla }\mathrm{\Psi }}}\underset{\mathit{VII}}{\underbrace{+{\rho }_0{\mathit{\eta }}^{\ast }·\mathbf{\nabla }\left(\frac{\mathit{\xi }·\mathbf{\nabla }{P}_0}{{\rho }_0}\right)}}\hfill \\ \hfill & \underset{\mathit{VIII}}{\underbrace{+{\mathit{\eta }}^{\ast }·\left[\frac{(\mathit{\xi }·\mathbf{\nabla }{P}_0)\mathbf{\nabla }{\rho }_0-(\mathit{\xi }·\mathbf{\nabla }{\rho }_0)\mathbf{\nabla }{P}_0}{{\rho }_0}\right]}}\}\mathrm{d}V.\hfill \end{array}$$(B.1)

At this point, the goal is to reformulate the above integral so that it is manifestly symmetric (in a Hermitian sense) with respect to ξ and η. In what follows, it is very important to bear in mind that Ω and its associated vector Ω = Ωe_z depend on ζ and θ. This is different than the approach taken in Lynden-Bell & Ostriker (1967) where Ω is constant (differential rotation is, instead, taken into account as a background velocity field, v₀). Furthermore, given that the equilibrium model may be discontinuous, it will be important, for some of the terms, to decompose the stellar volume into subdomains, V_i, such that the model is continuous in each subdomain. Obviously, the relation V = ∪_iV_i holds. Finally, when dealing with the gravitational potential, we will introduce the notation V_e to represent an external domain which comprises all of the space outside the star, and V_∞ to represent all of space, including the star.

Terms I and II can easily be rearranged into the following symmetric forms:

$$\begin{array}{cc}& I={\displaystyle {\int }_V-2i(\omega +m\mathrm{\Omega }){\rho }_0\mathbf{\Omega }·(\mathit{\xi }\times {\mathit{\eta }}^{\ast })\mathrm{d}V,}\hfill \end{array}$$(B.2)

$$\begin{array}{cc}& II={\displaystyle {\int }_V\{-{\rho }_0(\mathbf{\Omega }·\mathit{\xi })(\mathbf{\Omega }·{\mathit{\eta }}^{\ast })+{\rho }_0{\mathrm{\Omega }}^2\mathit{\xi }·{\mathit{\eta }}^{\ast }\}\mathrm{d}V.}\hfill \end{array}$$(B.3)

Term III can be rewritten as:

$$\begin{array}{cc}& III=-{\displaystyle {\int }_V{\rho }_0{\mathit{\eta }}^{\ast }·[\mathit{\xi }·\mathbf{\nabla }(\frac{\mathbf{\nabla }{P}_0}{{\rho }_0}+\mathbf{\nabla }{\mathrm{\Psi }}_0)]\mathrm{d}V}\hfill \\ \hfill & ={\displaystyle \sum _i{\int }_{V_i}\{\frac{(\mathit{\xi }·\mathbf{\nabla }{\rho }_0)({\mathit{\eta }}^{\ast }·\mathbf{\nabla }{P}_0)}{{\rho }_0}-{\mathit{\eta }}^{\ast }·[\mathit{\xi }·\mathbf{\nabla }\left(\mathbf{\nabla }{P}_0\right)]}\hfill \\ \hfill & -{\rho }_0{\mathit{\eta }}^{\ast }·[\mathit{\xi }·\mathbf{\nabla }\left(\mathbf{\nabla }{\mathrm{\Psi }}_0\right)]\}\mathrm{d}V,\hfill \end{array}$$(B.4)

where we have made use of the hydrostatic equilibrium equation, in which we have neglected viscosity and meridional circulation. It is important to note that on the second and third line, the integral is carried out over ∑_iV_i. The reason for this is that ∇P₀ may be discontinuous, meaning that ∇(∇P₀) has to be calculated over each separate subdomain, V_i. The last two terms are symmetric. This can be seen, for instance, by expressing them explicitly in terms of their Cartesian coordinates: ${\mathit{\eta }}^{*}\cdot \left[\mathit{\xi }\cdot \nabla \left(\nabla P_0\right)\right]={\left({\eta }^i\right)}^{*}{\xi }^j{\partial }_{ij}^2{P}_0$. The first term cancels out with the last part of term VIII.

When developing term IV, it is important to treat each domain, V_i, separately:

$$\begin{array}{cc}& IV={\displaystyle \sum _i{\int }_{V_i}\{-\mathbf{\nabla }·\left(P_0{\mathit{\eta }}^{\ast }\frac{\delta p}{P_0}\right)+\frac{\delta p}{P_0}\mathbf{\nabla }·\left(P_0{\mathit{\eta }}^{\ast }\right)\}\mathrm{d}V}\hfill \\ \hfill & =-{\displaystyle \sum _i{\int }_{B_i}\delta p{\mathit{\eta }}^{\ast }·\mathbf{dS}+\sum _i{\int }_{V_i}\frac{\delta p}{P_0}\mathbf{\nabla }·\left(P_0{\mathit{\eta }}^{\ast }\right)\mathrm{d}V,}\hfill \end{array}$$(B.5)

where B_i denotes the bounds of subdomain V_i. It turns out that the surface terms cancel out. Indeed, there are two possible cases. In the first case, the surface corresponds to an internal discontinuity. As explained in the previous section, both the normal component of the displacement and the Lagrangian pressure perturbation remain continuous across the discontinuity. Furthermore, there will be two surface terms, one for the domain just below the discontinuity and the other for the domain just above. The vector dS takes opposite signs in both surface terms since it is directed outwards from the domain. As a result, the two terms cancel. In the second case, the surface term corresponds to the stellar surface. As explained in Sect. 3.4, we impose the simple mechanical boundary condition δp = 0, thereby cancelling this surface term.

Terms IV and V may be combined as follows:

$$\begin{array}{cc}& IV+V={\displaystyle \sum _i{\int }_{V_i}[({\mathit{\eta }}^{\ast }·\mathbf{\nabla }{P}_0)(\frac{\delta \rho }{{\rho }_0}-\frac{\delta p}{P_0})}\hfill \\ \hfill & +({\mathit{\eta }}^{\ast }·\mathbf{\nabla }{P}_0+P_0\mathbf{\nabla }·{\mathit{\eta }}^{\ast })\frac{\delta p}{P_0}]\mathrm{d}V\hfill \\ \hfill & ={\displaystyle \sum _i{\int }_{V_i}[({\mathit{\eta }}^{\ast }·\mathbf{\nabla }{P}_0)\frac{\delta \rho }{{\rho }_0}-\frac{\delta {\stackrel{\sim }{\rho }}^{\ast }\delta p}{{\rho }_0}]\mathrm{d}V,}\hfill \end{array}$$(B.6)

where we have made use of the continuity equation and introduced the Lagrangian density perturbation, $\delta \stackrel{\sim }{\rho }$ associated with η.

Term VII is treated as follows:

$$\begin{array}{cc}& VII={\displaystyle \sum _i{\int }_{V_i}[\mathbf{\nabla }·\left({\rho }_0{\mathit{\eta }}^{\ast }\frac{\mathit{\xi }·\mathbf{\nabla }{P}_0}{{\rho }_0}\right)-\frac{\mathit{\xi }·\mathbf{\nabla }{P}_0}{{\rho }_0}\mathbf{\nabla }·\left({\rho }_0{\mathit{\eta }}^{\ast }\right)]\mathrm{d}V}\hfill \\ \hfill & ={\displaystyle \sum _i{\int }_{B_i}(\mathit{\xi }·\mathbf{\nabla }{P}_0){\mathit{\eta }}^{\ast }·\mathbf{dS}+\sum _i{\int }_{V_i}[(\mathit{\xi }·\mathbf{\nabla }{P}_0)\frac{\delta {\stackrel{\sim }{\rho }}^{\ast }}{{\rho }_0}}\hfill \\ \hfill & -\frac{(\mathit{\xi }·\mathbf{\nabla }{P}_0)({\mathit{\eta }}^{\ast }·\mathbf{\nabla }{\rho }_0)}{{\rho }_0}]\mathrm{d}V,\hfill \end{array}$$(B.7)

where we have once more made use of the continuity equation. In the above calculations, the surface terms do not cancel, but they are symmetric since both the discontinuities and the stellar surface follow isobars. The term on the last line cancels out with the first part of term VIII.

Last but not least, we deal with term VI. As has been shown in Unno et al. (1989) and Reese (2006), this term can be rearranged into an integral of the form $\frac{1}{\mathrm{\Lambda }}{\int }_{V_{\infty }}\mathbf{\nabla }\mathrm{\Psi }·\mathbf{\nabla }{\mathrm{\Phi }}^{\ast }\mathrm{d}V$, where Λ = 4πG or 4π in the dimensionless case, and Φ is the gravitational potential associated with the displacement field η. However, surface terms and terms arising from internal discontinuities were not dealt with in the above works. In what follows, we re-derive this expression, while keeping track of such terms:

$$\begin{array}{cc}& VI=-{\displaystyle \sum _i{\int }_{V_i}\mathbf{\nabla }·\left({\rho }_0{\mathit{\eta }}^{\ast }\mathrm{\Psi }\right)\mathrm{d}V+\sum _i{\int }_{V_i}\mathrm{\Psi }\mathbf{\nabla }·\left({\rho }_0{\mathit{\eta }}^{\ast }\right)\mathrm{d}V}\hfill \\ \hfill & =-{\displaystyle \sum _i{\int }_{B_i}{\rho }_0\mathrm{\Psi }{\mathit{\eta }}^{\ast }·\mathbf{dS}-\sum _i{\int }_{V_i}\mathrm{\Psi }{\stackrel{\sim }{\rho }}^{\ast }\mathrm{d}V}\hfill \\ \hfill & =-{\displaystyle \sum _{i+e}{\int }_{B_i}{\rho }_0\mathrm{\Psi }{\mathit{\eta }}^{\ast }·\mathbf{dS}-\frac{1}{\mathrm{\Lambda }}\sum _{i+e}{\int }_{V_i}\mathrm{\Psi }\mathrm{\Delta }{\mathrm{\Phi }}^{\ast }\mathrm{d}V}\hfill \\ \hfill & =-{\displaystyle \sum _{i+e}{\int }_{B_i}{\rho }_0\mathrm{\Psi }{\mathit{\eta }}^{\ast }·\mathbf{dS}-\frac{1}{\mathrm{\Lambda }}\sum _{i+e}{\int }_{V_i}\mathbf{\nabla }·\left(\mathrm{\Psi }\mathbf{\nabla }{\mathrm{\Phi }}^{\ast }\right)\mathrm{d}V}\hfill \\ \hfill & +\frac{1}{\mathrm{\Lambda }}{\displaystyle \sum _{i+e}{\int }_{V_i}\mathbf{\nabla }\mathrm{\Psi }·\mathbf{\nabla }{\mathrm{\Phi }}^{\ast }\mathrm{d}V}\hfill \\ \hfill & =-{\displaystyle \sum _{i+e}{\int }_{B_i}\underset{(a)}{\underbrace{\mathrm{\Psi }({\rho }_0{\mathit{\eta }}^{\ast }+\frac{\mathbf{\nabla }{\mathrm{\Phi }}^{\ast }}{\mathrm{\Lambda }})·\mathbf{dS}}}+\frac{1}{\mathrm{\Lambda }}\sum _{i+e}{\int }_{V_i}\mathbf{\nabla }\mathrm{\Psi }·\mathbf{\nabla }{\mathrm{\Phi }}^{\ast }\mathrm{d}V}\hfill \\ \hfill & =\frac{1}{\mathrm{\Lambda }}{\displaystyle \sum _{i+e}{\int }_{V_i}\mathbf{\nabla }\mathrm{\Psi }·\mathbf{\nabla }{\mathrm{\Phi }}^{\ast }\mathrm{d}V,}\hfill \end{array}$$(B.8)