Oscillations of 2D ESTER models. I. The adiabatic case

Recent numerical and theoretical considerations have shown that low-degree acoustic modes in rapidly rotating stars follow an asymptotic formula and recent observations of pulsations in rapidly rotating delta Scuti stars seem to match these expectations. However, a key question is whether strong gradients or discontinuities can adversely affect this pattern to the point of hindering its identification. Other important questions are how rotational splittings are affected by the 2D rotation profiles expected from baroclinic effects and whether it is possible to probe the rotation profile using these splittings. Accordingly, we numerically calculate pulsation modes in continuous and discontinuous rapidly rotating models produced by the 2D ESTER (Evolution STEllaire en Rotation) code. This spectral multi-domain code self-consistently calculates the rotation profile based on baroclinic effects and allows us to introduce discontinuities without loss of numerical accuracy. Pulsations are calculated using an adiabatic version of the Two-dimensional Oscillation Program (TOP) code. The variational principle is used to confirm the high accuracy of the pulsation frequencies and to derive an integral formula that closely matches the generalised rotational splittings, except when modes are involved in avoided crossings. This potentially allows us to probe the the rotation profile using inverse theory. Acoustic glitch theory, applied along the island mode orbit deduced from ray dynamics, can correctly predict the periodicity of the glitch frequency pattern produced by a discontinuity or the Gamma1 dip related to the He II ionisation zone in some of the models. The asymptotic frequency pattern remains sufficiently well preserved to potentially allow its detection in observed stars.


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, 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), 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. (2009Reese et al. ( , 2013, Ouazzani et al. (2015), and Ouazzani et al. (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 con-A&A proofs: manuscript no. article text, 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 , 2009) before being justified using ray dynamics (Lignières & Georgeot 2008, 2009, Pasek et al. 2011. 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, 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 lowdegree 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. 2006Reese et al. , 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.

Stellar models based on the ESTER code
The aim of the ESTER project is to produce and evolve selfconsistent 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 1 (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.
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 Γ 1 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 Γ 1 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 Γ 1 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, in-A&A proofs: manuscript no. article Table 1. Characteristics of the models used in this study. X int and X ext are the hydrogen contents below and above the discontinuity, respectively. R disc /R eq gives the equatorial radius at the discontinuity, normalised by the star's equatorial radius. We note that models Mreal and M are smooth.

Pulsation calculations
The pulsation modes are calculated using the Two-dimensional Oscillation Program (TOP, Reese et al. 2006Reese et al. , 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.

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

Non-dimensionalisation
The following reference length, pressure, and density scales are used: where R eq is the equatorial radius and M the mass. As a result of this choice of reference scales, the frequencies are nondimensionalised by the inverse of the dynamic time scale: Using this non-dimensionalisation leads to the same set of pulsation equations as previously (Eqs. (1)-(4)) except that Λ is now equal to 4π.

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: 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 App. A.1. A second condition is that the pressure remains continuous across the perturbed boundary. This condition is simply expressed as follows (see App. A.2): The third condition is the continuity of Ψ and its gradient across the perturbed boundary. This is enforced by the following conditions (see App. A.3):

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 0 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): where r ζ = ∂ ζ r = 1 − ε, r ext = 2, and where we have used a harmonic decomposition of Ψ, ℓ being the spherical harmonic degree.

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

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 to 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. Table 2. Maximum relative differences on pulsation frequencies using various resolutions in three of the models.
7.9 × 10 −4 2.8 × 10 −11 -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.

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 A&A proofs: manuscript no. article 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 App. B, we give a full derivation for baroclinic models with 2D rotation profiles and discontinuities. The final expression is: 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 (∇P + 0 = 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ñ = 19 to 30,l = 0 to 1, m = −3 to 3, was used. We recall thatñ is the number of nodes along an island mode's orbit whereasl 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ñ = 2n+ε andl = ℓ−|m|−ε 2 , where ε ≡ ℓ+m mod 2 ≡ñ 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 −11 to 10 −4 , apart from an outlier in model M6 2 . 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 −5 in the best cases. This is higher than the errors on most of the variational frequencies, except for model M6.
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 3 which lead to island modes which are 'polluted' by contributions from neighbouring modes.

An approximate rotation kernel
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 4 . 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: where we have used the following definitions/approximations: ( 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 two modes, the frequencies do not cross but the modes progressively exchange their geometric characteristics, thereby leading to a mixture of the two geometries when the frequencies are closest. Figure 3 of Espinosa et al. (2004) provides a nice illustration of a rotationally induced avoided crossing. 4 We are using the 'retrograde' convention, that is, retrograde modes have positive azimuthal orders.
Article number, page 6 of 21 D. R. Reese et al.: Oscillations of 2D ESTER models following equation: (18) This can be re-expressed as: Taking the square-root of both sides and assuming C ± ≪ (ω ± ∓ |m|Ω eff ± ) leads to: (ω + − |m|Ω eff . (20) This equation can finally be rearranged to yield: 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. 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 righthand side of Eq. (21)) is excellent thus potentially providing the basis for probing the rotation profile via inversions.

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). Figure 9 shows the pulsation frequencies obtained for the various models for modes withñ = 19 to 30,l = 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 ∆ñ = ω˜n +1,l, m − ω˜n ,l, m . In Fig. 10, we plot averaged large separations, ∆ñ = ω˜n +1,l − ω˜n ,l , where ω˜n ,l 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.

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: 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: . Generalised splittings versus weighted integrals of the rotation profile, and different terms from Eq. (13) for two pairs of prograde and retrograde modes. δω var /ω corresponds to the relative error on the variational frequency, and δω approx. var /ω is the same error when ω var is calculated using the approximations in Eqs. (15) and (17). The modes from the second pair, Solutions 3 and 4, are undergoing avoided crossings, whereas the other two modes are not.

Quantity
Solution 1   where ϑ is the angle between the surface normal and the wavevector, 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 App. 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. Figure 12 then shows the sound velocity, density, and perturbed pressure (δp/ √ 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.
These observations provide the basis for a simple toy model which is described in App. D. According to this model, the fre-quencies are given to first order by: where τ T = eq. surf. dr c is the acoustic travel time from the surface to the equator along the ray path, and τ 1 = disc. surf. dr c , the acoustic travel time from the surface to the discontinuity, as illustrated in Fig. 11. The quantity ǫ is given by the relation: and is treated as a small parameter. As shown in App. 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 τ 1 for the different models in our study. Although model Mreal is continuous, we included the τ 1 value for the He II ionisation zone. Indeed, the Γ 1 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 (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 cross- Table 4. Acoustic travel times in various models. The quantities τ T and τ 1 are illustrated in Fig. 11 ings 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. 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 Γ 1 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.

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 App. E, the pulsation mode including the reflected and refracted waves can locally be approximated as: where the superscripts '+' and '−' designate the upper and lower domains, respectively, k ± 1 and k ± 2 wave vectors, and where the amplitudes, A ± 1 , A ± 2 , are related via the relation: and k ± ⊥ is the wave vector component perpendicular to the surface. When k ≪ k ⊥ , the tangential component, the factor η reduces to: 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 ex-A&A proofs: manuscript no. article tracting the δp/P 0 profile as well as its horizontal and vertical gradients, ∇ (δp/P 0 ), and (∇ ⊥ (δp/P 0 )) ± , just above and below the discontinuity 5 . Since we are focusing on axisymmetric modes, the horizontal gradient is in the meridional planethere 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 0 )) ± profiles have the opposite sign to the ∇ (δp/P 0 ) 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 0 ) and δp/P 0 . 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 0 )) ± /k ⊥ and ∇ (δp/P 0 ) /k . 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 . 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. The quantities k and k ± ⊥ are the horizontal and vertical components of the wave vectors. The quantities k /k + ⊥ and k /k − ⊥ are given a negative sign since the amplitudes A ± 2 (corresponding to the wave vectors k ± 2 = k e − k ± ⊥ e ⊥ ) are larger (in absolute value) than the amplitudes A ± 1 .
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 ± ⊥ correspond to the numerically determined values of tan ϑ ± where ϑ ± are the angles between the surface normal and the wave vector 6 . The quantities 'k /k ± ⊥ (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 , 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.

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 solarlike 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, ∆ñ, 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 differentl 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.
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): where ∆ñ, ∆l, ∆m, andα 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, ∆ñ. 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 wellpreserved 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. D. R. Reese et al.: Oscillations of 2D ESTER models 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. Pressure profiles δp/P 0 ∇ ⟂ (δp/P 0 ) + ∇ ⟂ (δp/P 0 ) − ∇ ∥ (δp/P 0 ) Fig. 16. Extracted δp/P 0 , ∇ (δp/P 0 ), and ∇ ⊥ (δp/P 0 ) ± 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.

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.
A&A proofs: manuscript no. article Appendix A.2: Condition on the pressure perturbation The following condition ensures that the pressure remains continuous during the oscillatory movements: 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: We can then use Eq. (A.2) to develop, say, the right-hand side: where we have neglected second order terms on the third line.
Combining this with the previous equation leads to: 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 ξ − − ξ + · ∇P + 0 could be replaced by ξ − − ξ + · ∇P − 0 , since it only involves the gradient of the pressure along the boundary and the pressure is continuous across the boundary.

Appendix 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: 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 and its gradient to cancel zeroth order terms, and neglecting second order terms lead to the following equations: Ψ − (r, t) + ξ(r, t) · ∇Ψ − 0 (r, t) = Ψ + (r, t) + ξ(r, t) · ∇Ψ + 0 (r, t), (A.10) ∇Ψ − (r, t) + ξ(r, t) · ∇ ∇Ψ − 0 (r, t) = ∇Ψ + (r, t) + ξ(r, t) · ∇ ∇Ψ + 0 (r, t) . Given that ∇Ψ 0 is continuous, the first equation reduces to: In tensorial notation, the left-hand side of the second equation becomes: where E i is the natural basis,ξ the components of ξ over that basis, and E i the dual basis. Calculating the dot product of the above equation with E − i yields: is the Christoffel symbol of the second kind. With our choice of mapping, only Γ ζ ζζ is discontinuous across the boundary, hence the notation Γ − ζ ζζ and Γ + ζ ζζ . All of the other geometric quantities (E i , E i , Γ k i j 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, Ψ 0 , and ∇Ψ 0 ) yields the following three relations: where we have made use of the fact that ∂ 2 i j Ψ 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: (A.17) where Λ = 4πG or 4π in the dimensional or dimensionless case, respectively, g i j = E i · 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.14): The terms R − and R + cancel out since R is continuous across the boundary. The remaining equation is then 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ξ ζ 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
Appendix 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: 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 0 ). 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 = ∪ i V 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: Term III can be rewritten as: 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 i V i . The reason for this A&A proofs: manuscript no. article is that ∇P 0 may be discontinuous, meaning that ∇ (∇P 0 ) 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: η * · ξ · ∇ (∇P 0 ) = η i * ξ j ∂ 2 i j 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: 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: where we have made use of the continuity equation and introduced the Lagrangian density perturbation, δρ associated with η.
Term VII is treated as follows: 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 1 Λ V ∞ ∇Ψ · ∇Φ * dV, 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: whereρ is the Eulerian density perturbation associated with η, the notation 'i + e' stands for internal domains plus the external domain V e . Various steps in the above developments need further explanation. Firstly, on the third line, the external domain was incorporated along with internal domains. This step is justified because the supplementary terms are equal to zero. It must be noted that the surface associated with V e , only includes the lower bound, that is, the stellar surface. Secondly, the divergence (or Ostrogradsky's) theorem was used to transform volume integrals on lines one and four into surface integrals. While straightforward in most cases, it is not as obvious on line four for the external domain V e . Indeed, it is not clear if some external boundary at infinity should be included or not. This problem can be dealt with in a rigorous way by considering an external domain,Ṽ e , which is bounded by a sphere of radius R e and then taking the limit as R e goes to infinity. Such an approach was taken in Reese (2006) who showed that the external surface term goes to zero in such conditions. Finally, it is necessary to show that the surface terms cancel out. We first start by noting that the surface element, dS, is parallel to the vector E ζ , and so may be written asdSE ζ . Hence, the integrand in the surface terms may be written: Now, we recall that each internal boundary (either from a discontinuity or from the stellar surface) gives rise to two surface terms, one from the domain just below the boundary and the other from the domain just above. As can be seen in the above expression, these two surface terms will cancel out. Indeed, based on the interface conditions (Eqs. (9) and (10)) and the continuity of g i j (for our choice of coordinate system), the first part of the above expression is continuous across the boundary. OnlydS changes signs given that the vector dS is always directed outwards from the corresponding domain.
In this section, we consider a 1D toy model representative of a sound wave travelling along an island mode ray path in the presence of a discontinuity. Figure D.1 illustrates half of this model, the other half being deduced by symmetry. For the sake of simplicity, constant sound velocities, denoted c 1 and c 2 , are used over the domains [0, x 1 [ and [x 1 , x T ] (as well as their symmetric counterparts), where x T = x 1 + x 2 . We assume that the density is discontinuous between the two domains whereas the pressure and first adiabatic exponent are continuous, in accordance with our stellar models.
We then consider the following set of simplified pulsation equations: along with the boundary conditions: δp P 0 (x = 0) = 0 (D.4) at the surface and δp P 0 (x = x T ) = 0 or d dx δp P 0 (x = x T ) = 0 (D.5) at the equator (that is, at x T ). The latter conditions come from the fact that pulsation modes are either antisymmetric or symmetric with respect to the equator. Finally, the following interface conditions apply at x = x 1 : The pressure perturbation then takes on the following form: A 2 cos (k 2 (x − x T )) (D.7) where k i = ω/c i . The two options for the solution in the [x 1 , x T ] domain correspond to antisymmetric (or odd) and symmetric (or even) solutions, respectively. Enforcing the interface condition then leads to the following discriminants which define the eigenvalues: sin(ωτ 1 ) cos(ωτ 2 ) k 2 + sin(ωτ 2 ) cos(ωτ 1 ) k 1 = 0 (D.8) for antisymmetric modes, or − sin(ωτ 1 ) sin(ωτ 2 ) k 2 + cos(ωτ 2 ) cos(ωτ 1 ) k 1 = 0 (D.9) for symmetric modes, where τ i = x i /c i . In the simple case where c 1 = c 2 , the solutions are ω k = kπ τ T or ω k = k + 1 2 π τ T (D.10) for odd and even modes respectively, and where τ T = τ 1 + τ 2 . When c 1 and c 2 differ, we can perform a first order perturbative analysis by introducing a small parameter ǫ as follows: k 2 = k 1 (1 + ǫ).
(D.11) This leads to the following corrections on odd and even modes respectively: δω = (−1) k 2τ T ǫ sin(ω 0 (τ 1 − τ 2 )), (D.12) δω = (−1) k 2τ T ǫ cos(ω 0 (τ 1 − τ 2 )), (D.13) where ω 0 corresponds to the unperturbed frequencies given in Eq. (D.10). Combining the even and odd cases and including both the zeroth and first order components yields: where odd values of n correspond to even solutions and vice versa. The period of the frequency perturbation is analogous but somewhat simplified compared to the more general formula given in Monteiro et al. (1994) for non-rotating stars. Figure D.2 compares large frequency separations using the first order expression above (Eq. (D.14)) and those obtained from exact solutions to the discriminants given in Eqs. (D.8) and (D.9) for of the values τ 1 and τ T from model M7, the most extreme case. As can be seen, the first order expression gives an accurate idea of the period of the frequency deviation, and a rough idea of its amplitude.

Appendix E: Wave refraction and reflection
In this section, we recall some of the basic principles behind the Snell-Descartes law including partial wave reflection. A more complete treatment can be found in various textbooks such as Brekhovskikh (1980). We begin with a simple plane-parallel model using Cartesian coordinates. This can also be thought of as a local approximation to a more complex system. A discontinuity in density is located at z = 0. The media below and above this discontinuity is assumed to be uniform. Under these conditions, the fluid dynamic equations take on the following expressions: The interface conditions, as explained in App. A, ensure the continuity of δp P 0 and ξ z . Combining these equations leads to: Because of the partial reflection at the boundary, we cannot consider a plane-parallel wave in isolation but have to include the reflected wave. For the sake of generality, we consider such a combination both above and below the discontinuity. The leads to following generic solution: where the superscripts '+' and '−' designate the upper and lower domains, respectively. The standing wave equivalent to the above solution would take on the expression: δp P 0 ± = A ± 1 cos(k ± 1 · x) + A ± 2 cos(k ± 2 · x) exp(iωt). (E.6) The wave vectors take on the following form: The horizontal wave vector, k , is preserved between the two domains as a result of the continuity of the horizontal gradient of δp P 0 at the discontinuity. When combined with the dispersion relation, this leads to Snell-Descartes' law.
The continuity of δp/P 0 leads to the relation: The continuity of ξ z leads to: Combining these two equations leads the following matrix relation between the amplitudes: When the wave vector is nearly perpendicular to the discontinuity, that is, k ≪ k z , η takes on the following approximate expression: (E.12)