Issue 
A&A
Volume 645, January 2021



Article Number  A46  
Number of page(s)  20  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201935538  
Published online  06 January 2021 
Oscillations of 2D ESTER models
I. The adiabatic case
^{1}
LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France
email: daniel.reese@obspm.fr
^{2}
Astrophysics Research Group, Faculty of Engineering and Physical Sciences, University of Surrey, Guildford GU2 7XH, UK
^{3}
Universidad de Alcalá, Space Research Group, 28805 Alcalá de Henares, Spain
^{4}
Universté de Toulouse, UPSOMP, IRAP, Toulouse, France
^{5}
CNRS, IRAP, 14 avenue Edouard Belin, 31400 Toulouse, France
Received:
25
March
2019
Accepted:
12
October
2020
Context. Recent numerical and theoretical considerations have shown that lowdegree acoustic modes in rapidly rotating stars follow an asymptotic formula. In parallel, recent studies have revealed the presence of regular pulsation frequency patterns in rapidly rotating δ Scuti stars that seem to match theoretical expectations.
Aims. In this context, a key question is whether strong gradients or discontinuities can adversely affect the asymptotic frequency 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.
Methods. In order to address these questions, we numerically calculate stellar pulsation modes in continuous and discontinuous rapidly rotating models produced by the 2D Evolution STEllaire en Rotation (ESTER) code. This code selfconsistently calculates the rotation profile based on baroclinic effects and uses a spectral multidomain approach, thus making it possible to introduce discontinuities at the domain interfaces without loss of numerical accuracy. The pulsation calculations are carried out using an adiabatic version of the Twodimensional Oscillation Program (TOP) code. The variational principle is then used to confirm the high numerical accuracy of the pulsation frequencies and to derive an integral formula for the generalised rotational splittings. Acoustic glitch theory, combined with ray dynamics, is applied to the discontinuous models in order to interpret their pulsation spectra.
Results. Our results show that the generalised rotational splittings are very well approximated by the integral formula, except for modes involved in avoided crossings. This potentially allows the application of inverse theory for probing the rotation profile. We also show that glitch theory applied along the island mode orbit can correctly predict the periodicity of the glitch frequency pattern produced by the discontinuity or Γ_{1} dip related to the He II ionisation zone in some of the models. Furthermore, the asymptotic frequency pattern remains sufficiently well preserved to potentially allow its detection in observed stars.
Key words: stars: oscillations / stars: rotation / stars: interiors
© D. R. Reese et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 selfconsistent way using energy conservation. In contrast, the ESTER code deduces the rotation profile in a selfconsistent 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 firstorder perturbative approach (e.g. Deheuvels et al. 2014; Schou et al. 1998; Thompson et al. 2003). At high rotation rates, higherorder 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 thirdorder 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 lowdegree 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 selfconsistent 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 meandensity. 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 nonrotating 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 nonrotating 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 nonrotating 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 gmodes. In particular, they show the presence of a periodic component in the period spacings of such modes, analogous to what was found in the nonrotating 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 gmodes (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 multidomain spectral approach, ideal for introducing discontinuities while retaining a high numerical accuracy. The pulsation modes are calculated using a multidomain spectral version of the Twodimensional 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 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 surfacefitting 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 breakup 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 wellstudied δ 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 GaussLegendre 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 crosssection 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. 
Characteristics of the models used in this study.
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 nearsurface 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 gravitymode 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 nearsurface layers of stars such as the one shown in Fig. 3 for a 2 M_{⊙} nonrotating 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 SnellDescartes law in different situations.
Fig. 3.
Nearsurface density profile in a 2 M_{⊙} nonrotating main sequence model from grid B of Marques et al. (2008). 
3. Pulsation calculations
The pulsation modes are calculated using the Twodimensional 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 multidomain 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:
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 surfacefitting radial coordinate, and θ, the colatitude), Λ = 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. Nondimensionalisation
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 nondimensionalisation 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:
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):
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):
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 nonzero 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.
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 shiftinvert approach to target frequencies around a given shift, σ, before being solved through the ArnoldiChebyshev approach (e.g. Braconnier 1993; Chatelin 2012).
The multidomain spectral approach used in the radial direction leads to matrices A and B which are block tridiagonal. The matrix A − σB (which intervenes in the shiftinvert approach) can be efficiently factorised using successive factorisations of the diagonal blocks (including a corrective term from the nondiagonal 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.
Maximum relative differences on pulsation frequencies using various resolutions in three of the models.
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. ChristensenDalsgaard 1982). A general formulation of the variational principle in differentially rotating bodies has previously been obtained by LyndenBell & 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, numericallyfriendly 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:
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 ( 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 to 30, to 1, m = −3 to 3, was used. We recall that is the number of nodes along an island mode’s orbit whereas the number of nodes parallel to it. These are related to the usual quantum numbers, (n, ℓ, m), of pulsation modes in the nonrotating case via the relations and , where 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.
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^{3} 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^{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 following equation:
This can be reexpressed as:
Taking the squareroot of both sides and assuming leads to:
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 midlatitudes.
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 righthand 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.
Fig. 6.
Comparison between generalised rotational splittings and the corresponding weighted integrals of the rotation profile for model M (see Eq. (21)). 
Fig. 8.
Meridional crosssections 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 to 30, 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 . In Fig. 10, we plot averaged large separations, , where 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, , 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:
where k is the norm of the wavevector. A simple reflection was used at the stellar surface, rather than a more realistic but complex approach involving the cutoff frequency (e.g. Lignières & Georgeot 2009). Furthermore, we applied the SnellDescartes refraction law at the discontinuity:
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 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 () 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 (, 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:
where is the acoustic travel time from the surface to the equator along the ray path, and , 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 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 τ_{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 frequencies minus a second or thirdorder 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 thirdorder 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 higherorder 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. 
Acoustic travel times in various models.
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 thirdorder polynomial fit, as confirmed by the much smaller scale of the yaxis. 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 modes after subtraction of a thirdorder 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:
where the superscripts ‘+” and ‘−” designate the upper and lower domains, respectively, and wave vectors, and where the amplitudes, , , are related via the relation:
where
and is the wave vector component perpendicular to the surface. When the tangential component, k_{∥}, is negligible in front of k_{⊥}, 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 extracting 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 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_{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 SnellDescartes’ law. The individual wave amplitudes are obtained by calculating appropriate linear combinations of (∇_{⊥}(δp/P_{0}))^{±}/k_{⊥} and (∇_{∥}(δp/P_{0}))/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/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. 
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 . The quantities “(theo)” and “(theo)” correspond to the amplitudes deduced from and via Eq. (27). Apart from the (theo) for model M6, these values accurately reproduce the numerically obtained amplitudes, and , thus showing that the relationship on amplitudes is respected.
Wave vector components and amplitudes at the discontinuity for island modes in three of the models.
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_{∥}/ correspond to the numerically determined values of tanϑ^{±} where ϑ^{±} are the angles between the surface normal and the wave vector^{6}. The quantities “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 , 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 solarlike stars is the socalled 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 different 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):
where , , , 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.
Root mean square and maximal differences between numerical and asymptotic frequencies for the different models.
Nonetheless, other factors may hinder finding the above frequency pattern. Indeed, the presence of chaotic modes with their own independent semirandom 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:

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 prograderetrograde modes.

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 nonadiabatic version of TOP. This will allow us to investigate other topics such as mode excitation and mode behaviour near the stellar surface.
We recall that avoided crossings occur when the frequencies of two coupled modes approach one another as a function of some stellar parameter such as age or rotation rate. Due to the coupling between 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.
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 ANR16CE310007 as well as financial support from the Programme National de Physique Stellaire (PNPS) of the CNRS/INSU cofunded 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 201199992 made by GENCI (“Grand Equipement National de Calcul Intensif”) and to HPC resources of Calmip under project P0107.
References
 Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baglin, A., Auvergne, M., Barge, P., et al. 2009, in IAU Symposium, IAU Symp., 253, 71 [Google Scholar]
 Ballot, J., Lignières, F., Reese, D. R., & Rieutord, M. 2010, A&A, 518, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bedding, T. R., Murphy, S. J., Hey, D. R., et al. 2020, Nature, 581, 147 [CrossRef] [Google Scholar]
 Borucki, W., Koch, D., Batalha, N., et al. 2009, in IAU Symposium, IAU Symp., 253, 289 [Google Scholar]
 Bouabid, M.P., Dupret, M.A., Salmon, S., et al. 2013, MNRAS, 429, 2500 [NASA ADS] [CrossRef] [Google Scholar]
 Bouchaud, K., Domiciano de Souza, A., Rieutord, M., Reese, D. R., & Kervella, P. 2020, A&A, 633, A78 [CrossRef] [EDP Sciences] [Google Scholar]
 Bowman, D. M., & Kurtz, D. W. 2018, MNRAS, 476, 3169 [NASA ADS] [CrossRef] [Google Scholar]
 Braconnier, T. 1993, The ArnoldiTchebycheff Algorithm for Solving Large non Symmetric Eigenproblems, Technical Report TR/PA/93/25, CERFACS, Toulouse, France [Google Scholar]
 Breger, M., Fossati, L., Balona, L., et al. 2012, ApJ, 759, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Brekhovskikh, L. M. 1980, Waves in Layered Media, 2nd edn. (New York: Academic Press, Inc.) [Google Scholar]
 Chatelin, F. 2012, Eigenvalues of Matrices: Revised Edition ( SIAM) [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 1982, MNRAS, 199, 735 [Google Scholar]
 Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27 [Google Scholar]
 Deupree, R. G. 2011, ApJ, 742, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Espinosa, F., Pérez Hernández, F., & Roca Cortés, T. 2004, ESA SP559: SOHO 14 Helio and Asteroseismology: Towards a Golden Future, 424 [Google Scholar]
 Evano, B., Lignières, F., & Georgeot, B. 2019, A&A, 631, A140 [CrossRef] [EDP Sciences] [Google Scholar]
 García Hernández, A., MartínRuiz, S., Monteiro, M. J. P. F. G., et al. 2015, ApJ, 811, L29 [NASA ADS] [CrossRef] [Google Scholar]
 García Hernández, A., Moya, A., Michel, E., et al. 2009, A&A, 506, 79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 García Hernández, A., Moya, A., Michel, E., et al. 2013, A&A, 559, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Irwin, A. W. 2012, FreeEOS: Equation of State for Stellar Interiors Calculations (Astrophysics Source Code Library) [Google Scholar]
 Jackson, S. 1970, ApJ, 161, 579 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, S., MacGregor, K. B., & Skumanich, A. 2005, ApJS, 156, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Lignières, F., & Georgeot, B. 2008, Phys. Rev. E, 78, 016215 [NASA ADS] [CrossRef] [Google Scholar]
 Lignières, F., & Georgeot, B. 2009, A&A, 500, 1173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lignières, F., Rieutord, M., & Reese, D. 2006, A&A, 455, 607 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lovekin, C. C., & Deupree, R. G. 2008, ApJ, 679, 1499 [NASA ADS] [CrossRef] [Google Scholar]
 Lovekin, C. C., Deupree, R. G., & Clement, M. J. 2009, ApJ, 693, 677 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293 [Google Scholar]
 MacGregor, K. B., Jackson, S., Skumanich, A., & Metcalfe, T. S. 2007, ApJ, 663, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars. Astronomy and Astrophysics Library (SpringerVerlag) [CrossRef] [Google Scholar]
 Mantegazza, L., Poretti, E., Michel, E., et al. 2012, A&A, 542, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marques, J. P., Monteiro, M. J. P. F. G., & Fernandes, J. M. 2008, Ap&SS, 316, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michel, E., Dupret, M. A., Reese, D., et al. 2017, in European Physical Journal Web of Conferences, Eur. Phys. J. Web Conf., 160, 03001 [CrossRef] [Google Scholar]
 Miglio, A., Montalbán, J., Noels, A., & Eggenberger, P. 2008, MNRAS, 386, 1487 [NASA ADS] [CrossRef] [Google Scholar]
 Mirouh, G. M., Reese, D. R., Rieutord, M., et al. 2017, in SF2A2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, P. Di Matteo, F. Herpin, et al. [Google Scholar]
 Mirouh, G. M., Angelou, G. C., Reese, D. R., & Costa, G. 2019, MNRAS, 483, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Monteiro, M. J. P. F. G., ChristensenDalsgaard, J., & Thompson, M. J. 1994, A&A, 283, 247 [NASA ADS] [Google Scholar]
 Ostriker, J. P., & Mark, J. W.K. 1968, ApJ, 151, 1075 [NASA ADS] [CrossRef] [Google Scholar]
 Ouazzani, R.M., & Goupil, M.J. 2012, A&A, 542, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ouazzani, R.M., Roxburgh, I. W., & Dupret, M.A. 2015, A&A, 579, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ouazzani, R.M., Salmon, S. J. A. J., Antoci, V., et al. 2017, MNRAS, 465, 2294 [NASA ADS] [CrossRef] [Google Scholar]
 Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paparó, M., Benkő, J. M., Hareter, M., & Guzik, J. A. 2016, ApJS, 224, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, A&A, 598, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pasek, M., Georgeot, B., Lignières, F., & Reese, D. R. 2011, Phys. Rev. Lett., 107, 121101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Pasek, M., Lignières, F., Georgeot, B., & Reese, D. R. 2012, A&A, 546, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., Lignières, F., & Ballot, J. 2016, A&A, 587, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. 2006, PhD thesis, Université Toulouse III  Paul Sabatier, http://tel.archivesouvertes.fr/tel00120334 [Google Scholar]
 Reese, D. 2008, J. Phys. Conf. Ser., 118, 012023 [NASA ADS] [CrossRef] [Google Scholar]
 Reese, D., Lignières, F., & Rieutord, M. 2006, A&A, 455, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D., Lignières, F., & Rieutord, M. 2008, A&A, 481, 449 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. R., MacGregor, K. B., Jackson, S., Skumanich, A., & Metcalfe, T. S. 2009, A&A, 506, 189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. R., Prat, V., Barban, C., van ’t VeerMenneret, C., & MacGregor, K. B. 2013, A&A, 550, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. R., Lara, F. E., & Rieutord, M. 2014, in Precision Asteroseismology, eds. J. A. Guzik, W. J. Chaplin, G. Handler, & A. Pigulski, IAU Symp., 301, 169 [Google Scholar]
 Reese, D. R., Lignières, F., Ballot, J., et al. 2017, A&A, 601, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003 [Google Scholar]
 Rieutord, M., & Espinosa Lara, F. 2009, Commun. Asteroseismol., 158, 99 [NASA ADS] [Google Scholar]
 Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
 Roxburgh, I. W. 2006, A&A, 454, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roxburgh, I. W., Griffith, J. S., & Sweet, P. A. 1965, Z. Astrophys., 61, 203 [Google Scholar]
 Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [Google Scholar]
 Soufi, F., Goupil, M.J., & Dziembowski, W. A. 1998, A&A, 334, 911 [Google Scholar]
 Suárez, J. C., Moya, A., Amado, P. J., et al. 2009, ApJ, 690, 1401 [NASA ADS] [CrossRef] [Google Scholar]
 Suárez, J. C., García Hernández, A., Moya, A., et al. 2014, A&A, 563, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thompson, M. J., ChristensenDalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars, 2nd edn. (Tokyo: University of Tokyo Press) [Google Scholar]
 Van Reeth, T., Tkachenko, A., Aerts, C., et al. 2015, ApJS, 218, 27 [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
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:
An analogous formula applies for the spatial coordinates just outside the boundary. This leads to the following matching condition:
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:
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 lefthand side and yields:
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:
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 righthand 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 righthand 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 righthand 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 could be replaced by , 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:
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:
Given that ∇Ψ_{0} is continuous, the first equation reduces to:
In tensorial notation, the lefthand 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 yields:
where 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}, with (i, j, k)≠(ζ, ζ, ζ)) are continuous. Inserting this expression into the lefthand side of Eq. (A.11) and a similar expression in the righthand side, and simplifying out continuous terms (geometric, Ψ_{0}, and ∇Ψ_{0}) yields the following three relations:
where we have made use of the fact that 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:
where Λ = 4πG or 4π in the dimensional or dimensionless case, respectively, g^{ij} = 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.13):
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
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 LyndenBell & 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 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: . 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 , 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 rederive 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, , 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 as . 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^{ij} (for our choice of coordinate system), the first part of the above expression is continuous across the boundary. Only changes signs given that the vector dS is always directed outwards from the corresponding domain.
We now combine the above formulas to obtain the following expression:
We note that at this stage, we have not yet made use of the adiabatic approximation (apart from cancelling out a surface term, thanks to the boundary condition δp = 0, a condition which is usually applied in adiabatic calculations). We now use the adiabatic relation, , to replace the Lagrangian density variations by Lagrangian pressure perturbations, and then develop these in terms of Eulerian pressure perturbations and displacement fields. The final result is:
where π is the Eulerian pressure perturbation associated with the displacement fields η. This expression is manifestly symmetric in (ξ, P, Ψ) and (η, π, Φ) and consequently leads to the variational principle (LyndenBell & Ostriker 1967). A useful consequence of this is the quadratic convergence of the variational frequencies (obtained by assuming that (ξ, P, Ψ) = (η, π, Φ) and solving the above equation for ω) to the true frequency, as a function of the error on the eigenfunctions (e.g. ChristensenDalsgaard 1982).
B.2. Explicit formulas
In what follows, we provide explicit expressions for the different terms which intervene in Eq. (B.11), a number of which were already given in Reese et al. (2006). Such expressions are needed when evaluating numerically the variational frequency:
The term ξ^{*} ⋅ [ξ⋅∇(∇P_{0})] can be developed through tensor analysis:
where we have expressed the displacement using the alternate components on the last line.
Appendix C: Ray dynamics
Ray trajectories were calculated in the usual spherical coordinate system (r, θ, ϕ) using the system of equations provided in Prat et al. (2016). We used as the Hamiltonian function. This leads to the following system:
Although the ray trajectories are calculated in the spherical coordinate system, the various derivatives of are first calculated in the spheroidal coordinate system before being converted to the spherical system via the following relations:
The system of Eqs. (C.1)–(C.4) is solved numerically for an initial position and wave vector using a fourth order RungeKutta method, except near discontinuities and the stellar surface where Heun’s thirdorder method (also a RungeKutta method) is used instead. Indeed at these locations, the step is adjusted so as to fall precisely on the relevant boundary, thus making it easier to apply a wave reflection or SnellDescartes’ law, and the use of Heun’s method reduces the risk of overstepping this boundary, unlike the fourth order RungeKutta method.
Appendix D: Toy model for glitches
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.
Fig. D.1.
Sound velocity profile in toy model. Only half the model is shown, the other half (beyond the centre, that is, x_{1} + x_{2}) being symmetric. 
We then consider the following set of simplified pulsation equations:
along with the boundary conditions:
at the surface and
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:
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:
for antisymmetric modes, or
for symmetric modes, where τ_{i} = x_{i}/c_{i}.
In the simple case where c_{1} = c_{2}, the solutions are
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:
This leads to the following corrections on odd and even modes respectively:
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 nonrotating 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.
Fig. D.2.
Large frequency separations for various calculations of the frequency: ω_{0} corresponds to the zeroth order expression (i.e. no perturbations to the sound velocity are included), ω to the first order approximation, and ω_{exact} to exact solutions of the discriminant equations. 
Appendix E: Wave refraction and reflection
In this section, we recall some of the basic principles behind the SnellDescartes law including partial wave reflection. A more complete treatment can be found in various textbooks such as Brekhovskikh (1980).
We begin with a simple planeparallel 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 Appendix A, ensure the continuity of and ξ_{z}. Combining these equations leads to:
Because of the partial reflection at the boundary, we cannot consider a planeparallel 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:
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 at the discontinuity. When combined with the dispersion relation, this leads to SnellDescartes’ 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:
where
When the wave vector is nearly perpendicular to the discontinuity, that is, k_{∥} ≪ k_{z}, η takes on the following approximate expression:
All Tables
Maximum relative differences on pulsation frequencies using various resolutions in three of the models.
Wave vector components and amplitudes at the discontinuity for island modes in three of the models.
Root mean square and maximal differences between numerical and asymptotic frequencies for the different models.
All Figures
Fig. 1.
Meridional crosssection of model M6 showing where the discontinuity is located. The other models have a discontinuity closer to the surface. 

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

In the text 
Fig. 3.
Nearsurface density profile in a 2 M_{⊙} nonrotating main sequence model from grid B of Marques et al. (2008). 

In the text 
Fig. 4.
Relative differences between numerical and variational frequencies. 

In the text 
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 the text 
Fig. 6.
Comparison between generalised rotational splittings and the corresponding weighted integrals of the rotation profile for model M (see Eq. (21)). 

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

In the text 
Fig. 8.
Meridional crosssections 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. 

In the text 
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). 

In the text 
Fig. 10.
Averaged pseudo large frequency separations, , for the different models. 

In the text 
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. 

In the text 
Fig. 12.
Sound velocity (first row), density (second row), and perturbed pressure (, 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. 

In the text 
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. 

In the text 
Fig. 14.
Frequencies of the modes after subtraction of a thirdorder 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. 

In the text 
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. 

In the text 
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. 

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

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

In the text 
Fig. D.1.
Sound velocity profile in toy model. Only half the model is shown, the other half (beyond the centre, that is, x_{1} + x_{2}) being symmetric. 

In the text 
Fig. D.2.
Large frequency separations for various calculations of the frequency: ω_{0} corresponds to the zeroth order expression (i.e. no perturbations to the sound velocity are included), ω to the first order approximation, and ω_{exact} to exact solutions of the discriminant equations. 

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