© C. Pinçon and M. Takata 2022

1. Introduction

Mechanical forcing, turbulent motions, or thermal instabilities can perturb the equilibrium state of stars and generate internal waves (e.g., Samadi et al. 2015, and references therein). These waves can propagate back and forth several times between the center and surface of stars and a resulting positive interference can give rise to global resonant modes, the oscillation frequencies of which directly depend on the stellar interior properties. The study of the oscillation power spectrum of these modes (asteroseismology) can, in turn, provide us with precious information on the stellar structure. Over recent decades, the high-quality seismic data first provided for the Sun by the spacecraft SoHO (Domingo et al. 1995) and ground-based telescope networks (e.g., GONG and BiSON projects, Chaplin 1997; Leibacher & GONG Project Team 1998) as well as subsequent data gathered for thousands of distant stars by space-borne missions CoRoT, Kepler, and TESS (Baglin et al. 2006; Borucki et al. 2010; Ricker et al. 2015), have indeed brought stringent constraints on the stellar dynamics and evolution (e.g., Chaplin & Miglio 2013; Grotsch-Noels et al. 2016; Mosser & Miglio 2016; Christensen-Dalsgaard 2021). The exploitation of such an amount of information and its physical interpretation have greatly relied on analytical descriptions of the linear stellar oscillations. In particular, theoretical expressions of the eigenfrequency patterns appeared to be crucial for extracting the prevailing features of the observed oscillation power spectra, defining seismic indicators that are relevant of the stellar structure, and enabling statistical studies on large samples of stars via automated methods (e.g., see Farnir et al. 2019; Appourchaux 2020; Gehan et al. 2021, for recent examples).

Most of the available analytical descriptions of the linear stellar oscillations are based on asymptotic methods within the JWKB approximation. This approximation assumes that the wavelength is much smaller than the variation scale height of the background state almost everywhere in stars (e.g., Olver 1975; Gough 2007). For example, asymptotic analyses predict that the high-frequency acoustic modes, the propagation cavity of which is mostly located in the uppermost layers of stars, are nearly evenly spaced in frequency, with a frequency spacing (or large frequency separation) directly linked to the mean density of stars. In contrast, the low-frequency gravity modes, which mostly propagate in the inner layers of stars, are asymptotically expected to be nearly evenly spaced in period, with a period spacing related to the stellar core density (Vandakurov 1968; Tassoul & Tassoul 1968; Smeyers 1968; Shibahashi 1979; Tassoul 1980). This leading-order description of the mode frequency distribution is nevertheless insufficient for reproducing all the diversity and details observed in the high-quality oscillation spectra of the Sun and thousands of other stars.

First, in low-mass main sequence stars, the deviation of the observed acoustic mode frequencies from the nearly regular asymptotic pattern was shown to take the form of a low-amplitude sinusoidal-like signal (e.g., Hill & Rosenwald 1986). The origin of these small frequency perturbations is attributed to sharp structural variations in some regions of stars where the JWKB approximation fails, the so-called glitches (Vorontsov 1988; Gough & Thompson 1988; Gough 1990). For instance, glitches are expected at the interface between convective and radiative zones where discontinuities in the temperature gradient, in its derivatives, or in the molecular weight can occur, depending on the mixing processes at work (e.g., Gough & Sekii 1993; Roxburgh & Vorontsov 1994; Monteiro et al. 1994, 2000; Audard & Provost 1994), as well as in the helium ionization zone where the first adiabatic index abruptly varies (e.g., Gough & Vorontsov 1995; Gough 2002; Basu & Antia 2004; Houdek & Gough 2007). Most prior works have assumed that the amplitude of these acoustic glitches is small enough so that the induced deviation from the asymptotic frequency pattern can be analytically modeled using first-order perturbation methods based on the variational principle. These developments have permitted in-depth studies of acoustic glitches in the Sun and brought about stringent constraints on the position of the base of the convective envelope, on the extent of the overshooting region beneath (e.g., Basu et al. 1994; Christensen-Dalsgaard et al. 1995; Basu 1997), and on the location of the ionization zone and the surface helium abundance (e.g., Vorontsov et al. 1991; Perez Hernandez & Christensen-Dalsgaard 1994; Monteiro & Thompson 2005). The exquisite store of seismic data collected by space-borne missions CoRoT and Kepler subsequently enabled similar studies in dozens of other distant main-sequence stars (e.g., Mazumdar et al. 2012, 2014; Lebreton & Goupil 2012; Verma et al. 2017; Farnir et al. 2020a). All the seismic information extracted from acoustic glitches currently stand for one of the most important sources of constraints for stellar modeling (e.g., Verma & Silva Aguirre 2019; Farnir et al. 2020b).

Furthermore, in intermediate-mass main sequence stars, the frequency pattern of the observed gravity modes can also be affected by glitches. In these stars, the expansion and recession of the convective core create a gradient in the mean molecular weight at its upper boundary, resulting in a large discontinuity in the Brunt–Väisälä frequency (Berthomieu & Provost 1988; Provost et al. 1990). Such large-amplitude glitches locally induce an important wave reflection and result in an unequal distribution of the mode energy on both sides of the glitch: this is referred to as the mode-trapping phenomenon. Unlike low-amplitude glitches in the Sun, the deviations from the asymptotic frequency pattern of gravity modes that are induced by such large-amplitude glitches cannot be treated as small perturbations. In this situation, the eigenfrequencies of gravity modes can be better analytically expressed using the matching of the asymptotic solutions on each side of the glitch (e.g., Miglio et al. 2008a,b). The period spacing between adjacent gravity modes predicted by such models significantly differ from the uniform asymptotic value predicted in the absence of discontinuity, which offers an interesting potential to constrain the chemical mixing at the edge of convective cores in these stars (e.g., Degroote et al. 2010). The frequency pattern of gravity modes in white dwarfs can similarly be affected by compositional layering and therefore can similarly provide information on the internal structure and transport processes in these advanced evolutionary stages (e.g., Brassard et al. 1992; Kawaler & Bradley 1994).

Finally, in post-main sequence stars, the density contrast between the core and surface is so large that the oscillation modes can propagate both in an inner cavity, where they behave as gravity modes, and in an outer cavity, where they behave as acoustic modes: these are the so-called mixed modes (Scuflaire 1974; Aizenman et al. 1977). Both cavities are coupled by an intermediate region where the modes are evanescent (e.g., Hekker & Christensen-Dalsgaard 2017, for a review). Neglecting the possible effect of internal glitches, the asymptotic frequency pattern of mixed modes presents the characteristics of both the acoustic and gravity modes spectra (Shibahashi 1979; Tassoul 1980; Takata 2016a; Loi & Papaloizou 2020). The asymptotic analyses of mixed modes appeared essential to analyze the large amount of seismic data collected by the satellite CoRoT and Kepler for evolved stars and represent a solid theoretical basis for interpreting these observations in terms of internal structure (e.g., Mosser et al. 2018). The study of mixed modes put stringent constraints for instance on the nuclear-burning state (e.g., Montalbán et al. 2010; Mosser et al. 2014; Vrard et al. 2016), the core rotation (e.g., Goupil et al. 2013; Gehan et al. 2018; Deheuvels et al. 2020), and the mid-layer structure of these stars (e.g., Mosser et al. 2017; Hekker et al. 2018; Khan et al. 2018; Pinçon et al. 2019, 2020), as well as on the amount of core overshooting on the main sequence (e.g., Montalbán et al. 2013; Noll et al. 2021). Delving into more details of the mixed mode oscillation spectra, acoustic glitches associated with the helium ionization zone could also be detected (Vrard et al. 2015; Dréau et al. 2020). Other seismic signatures, still not observed, have also been theoretically predicted. On the one hand, Cunha et al. (2015, 2019) proposed a more complex description of mixed modes including the presence of buoyancy glitches, that is, the sharp gradient in the molecular weight induced by the migration of the base of the convective zone during the first dredge-up. On the other hand, Deheuvels & Belkacem (2018) investigated the properties of mixed modes during the ignition of the helium burning (i.e., helium sub-flashes). During this phase, the temporary helium-burning shell is convective and mixed modes become evanescent inside. The inner propagation cavity is thus split into two parts. Mixed modes can thus propagate in three cavities separated from each other by two evanescent coupling regions. Both previous internal features were shown to produce remarkable seismic signatures in the mixed mode spectra, with promising probing potentials.

All the previous examples demonstrate the diversity of potential resonance configurations throughout the Hertzsprung-Russel diagram and the large amount of associated analytical descriptions. Despite this diversity, the effect of low- and large-amplitude glitches as well as evanescent regions on the mode frequencies results from a similar physical phenomenon. Indeed, either near sharp gradients or close to the boundaries between evanescent regions and resonant cavities (i.e., close to turning points where the radial wave number vanishes), the variation scale height of the background state is smaller than the local oscillation wavelength, the JWKB approximation locally fails, and an incident wave energy ray is partially reflected and transmitted (see, e.g., Appendix C of Pinçon et al. 2020, for a simple example). Based on basic linear wave principles, Takata (2016b) considered such a physical picture and reformulated the mixed mode frequency pattern of red giant stars in a very general way by describing the evanescent region as a simple barrier associated with a wave reflection coefficient. Pinçon (2019) adapted the same picture to model the effect of one glitch on the oscillation frequencies of star, which had been actually already anticipated by Roxburgh & Vorontsov (2001). In this paper, we aim to extend these previous representations to a multitude of resonant cavities and barriers and to obtain a generalized formulation for the resonance condition of oscillation modes in spherical stars that is applicable to any configurations and evolutionary states in a practical way.

The paper is organized as follows. In Sect. 2, we introduce the theoretical background about gravito-acoustic oscillations in stars and the modeling of the wave transmission-reflection problem. This introductory material is subsequently applied to obtain a very general expression of the resonance condition following two complementary pictures in Sects. 3 and 4. In addition, the distribution of the mode energy and the mode amplitudes among the different resonant cavities are addressed in Sect. 5. As a first illustration, the resulting resonance expression is then applied on simple usual oscillation configurations in Sect. 6 and the compatibility with previous formulations is discussed. We present our concluding remarks in Sect. 7.

2. Setting the stage

In this first section, we introduce the theoretical background and set the main physical description of the stellar oscillation modes used in the present paper.

2.1. General theoretical framework

In this work, we focus on the linear asymptotic and adiabatic global standing modes of spherical stars. In other words, we make the three following main assumptions on the oscillations.

First, we assume that the oscillations are small-amplitude perturbations of the stellar equilibrium state (i.e., within the linear approximation) and that there is at least one region where they can propagate as progressive waves and where the short-wavelength JWKB approximation is met (i.e., within the asymptotic limit). Such a region is referred to as a resonant cavity.

Second, we assume that the waves are generated at a time t = −∞, and we neglect non-adiabatic effects during their propagation (i.e., within the adiabatic limit). We therefore focus on stationary modes that persist indefinitely over time with a time dependence in e^−iσt, where i is the imaginary unit and σ is the angular oscillation frequency. This seems reasonable in a first step since, otherwise, excited progressive waves would be rapidly damped and could not constructively interfere to form eigenmodes.

Third, we neglect the effects of stellar rotation and internal magnetic effects on the oscillations; in other words, there is no preferential axis from the point of view of the waves. Owing to the spherical symmetry, the perturbations associated with the oscillations are separable into angular and radial parts, and the angular part is represented by the orthonormal set of scalar and vector spherical harmonics, $Y_{\mathcal{l}}^m$, in the stellar frame, with angular degrees, ℓ, and azimuthal numbers, m (e.g., Ledoux & Walraven 1958; Unno et al. 1989). In the following, we explicitly focus on a (ℓ,m) harmonic; the Eulerian perturbation of pressure, p′, and the oscillation displacement field, ξ, are thus expressed as:

$$\begin{array}{cc}& p^{\prime }(\mathit{r},t)={\stackrel{\sim }{p}}^{\prime }(r)Y_{\ell }^m(\theta ,s)e^{-\mathrm{i}\sigma t},\hfill \end{array}$$(1)

$$\begin{array}{cc}& \mathit{\xi }(\mathit{r},t)=[{\stackrel{\sim }{\xi }}_{\mathrm{r}}(r)Y_{\ell }^m(\theta ,s){\mathit{e}}_{\mathrm{r}}+{\stackrel{\sim }{\xi }}_{\mathrm{h}}(r)r\mathbf{\nabla }{Y}_{\ell }^m(\theta ,s)]e^{-\mathrm{i}\sigma t},\hfill \end{array}$$(2)

where (r, θ, s) are the spherical coordinates in the stellar frame, e_r is the radial unit vector, ∇ is the gradient operator, and ${\stackrel{\sim }{\xi }}_{\mathrm{r}}$ and ${\stackrel{\sim }{\xi }}_{\mathrm{h}}$ are the radial and poloidal components of the mode displacement, respectively.

The present theoretical framework is therefore appropriate to investigate the frequency pattern of eigenmodes in the slow rotator, low magnetic field, and adiabatic limits. We note that the study of the small deviations from this leading-order pattern that can be induced by rotational, magnetic, and non-adiabatic effects is theoretically tractable in a subsequent step using perturbative methods (e.g., Ledoux 1951; Dziembowski 1977; Dziembowski & Goode 1984; Pinçon et al. 2021), which is beyond the scope of this paper. Once the framework is set, we introduce in the following the basics on the propagation of gravito-acoustic waves that are useful for our purpose while distinguishing two kinds of regions: first, the resonant cavities where progressive waves can propagate and where the JWKB approximation is met; and second, the regions where waves are evanescent or where the JWKB hypothesis fails at some point (i.e., in the vicinity of turning points or near sharp structural variations). These regions are referred to as barriers in the following.

2.2. Wave propagation in resonant cavities

2.2.1. Adiabatic linear wave equation

Assuming an adiabatic equation of state, the continuity, momentum and Poisson equations linearized around the hydrostatic equilibrium of stars results in a fourth-order differential system with respect to radius for the radial wave displacement, the Eulerian pressure perturbation, the Eulerian perturbation of the gravitational potential, and the radial derivative of this latter (e.g., Ledoux & Walraven 1958). Nevertheless, within the short-wavelength JWKB assumption, the effect of the perturbation of the gravitational potential on the wave displacement field is known to be negligible at leading order. This is the so-called Cowling approximation (Cowling 1941). It relies on the fact that because of the large number of radial nodes in the vicinity of the considered layer, the effect of the small-scale density perturbations integrated over the whole stellar volume is negligible on the local gravitational acceleration according to the Poisson equation (see, e.g., Appendix A of Pinçon et al. 2020, for a scaling-based justification). The validity of the Cowling approximation for any angular degree, ℓ, inside the resonant cavities thus permits us to reduce the adiabatic linear wave equation to the second order, which has the advantage of being analytically tractable using the usual asymptotic methods (e.g., Olver 1975).

In these considerations, the formulation of Shibahashi (1979) within the Cowling approximation can be used to describe the leading-order behavior of the oscillations inside resonant cavities. According to the work of Shibahashi (1979), the wave equation can be expressed as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathrm{\Psi }}{\mathrm{d}{r}^2}+[{\mathcal{K}}_{\mathrm{r}}^2-\frac{M(r)}{H_{\mathrm{p}}^2}]\mathrm{\Psi }=0,\end{array}$$(3)

where the wave function and the squared local asymptotic radial wave number are respectively defined as

$$\begin{array}{cc}& \mathrm{\Psi }={\rho }^{1/2}cr{|1-\frac{S_{\ell }^2}{{\sigma }^2}|}^{-1/2}{\stackrel{\sim }{\xi }}_{\mathrm{r}}{e}^{-\mathrm{i}\sigma t},\hfill \end{array}$$(4)

$$\begin{array}{cc}& {\mathcal{K}}_{\mathrm{r}}^2=\frac{{\sigma }^2}{c^2}(\frac{N^2}{{\sigma }^2}-1)(\frac{S_{\ell }^2}{{\sigma }^2}-1).\hfill \end{array}$$(5)

In these equations, ρ is the equilibirum density, c² = Γ₁p/ρ is the squared sound speed, with Γ₁ the first adiabatic index and p the equilibrium pressure, H_p is the pressure scale height, and S_ℓ and N are the Lamb and Brunt–Väisälä frequencies, respectively, the expressions of which are provided by

$$\begin{array}{cc}& S_{\ell }^2=\frac{\ell (\ell +1)c^2}{r^2}\hfill \end{array}$$(6)

$$\begin{array}{cc}& N^2=\frac{g}{r}(\frac{1}{{\mathrm{\Gamma }}_1}\frac{\mathrm{d}lnp}{\mathrm{d}lnr}-\frac{\mathrm{d}ln\rho }{\mathrm{d}lnr}),\hfill \end{array}$$(7)

where g is the gravitational acceleration. Finally, the M(r) function in Eq. (3) is a radial function that depends on the variation of the equilibrium structure (see Appendix A for an expression).

As mentioned earlier, the second-order wave equation in Eq. (3) can be considered as a very good approximation inside regions where the JWKB approximation is met. By definition, these regions are far enough away from sharp variations in the stellar structure (i.e., relative to the local wavelength) in such a way that the M function in Eq. (3) remains on the order of unity at most inside these regions¹. We therefore understand that the resonant cavities correspond to regions where ${\mathcal{K}}_{\text{r}}^2>0$ that are located between consecutive turning points (i.e., where ${\mathcal{K}}_{\text{r}}^2=0$ for σ² = N² or ${\sigma }^2=S_{\mathcal{l}}^2$) and sharp gradients, but far enough away from these latter layers for the JWKB approximation to be met (i.e., ${\mathcal{K}}_{\text{r}}^2\gg 1/H_{\text{p}}^2$). According to Eq. (5), this is the case where ${\sigma }^2\ll (S_{\mathcal{l}}^2\text{and}{N}^2)$ or ${\sigma }^2\gg (S_{\mathcal{l}}^2\text{and}{N}^2)$, which corresponds to low-frequency gravity-dominated waves or high-frequency pressure-dominated waves, respectively.

From a general point of view, one particular solution of the wave equation in Eq. (3), denoted by ψ, can be formally written in a plane wave form as

$$\begin{array}{c}\hfill \psi (r)=e^{\mathrm{i}\varphi (r)},\end{array}$$(8)

where ϕ(r) is a complex phase function. It is easy to check that its complex conjugate, ψ^⋆, is also a solution. Moreover, their Wronskian 𝒲_ψ is equal to

$$\begin{array}{c}\hfill {\mathcal{W}}_{\psi }=({\psi }^{\star }\frac{\mathrm{d}\psi }{\mathrm{d}r}-\psi \frac{\mathrm{d}{\psi }^{\star }}{\mathrm{d}r})=2\mathrm{i}{\mathcal{R}}_{\mathrm{e}}\left[\frac{\mathrm{d}\varphi }{\mathrm{d}r}\right]{\left|\psi \right|}^2,\end{array}$$(9)

where ℛ_e[⋅] denotes the real part. Therefore, ψ and ψ^⋆ are linearly independent and thus form a vectorial basis for the solutions of the wave equation if and only if ℛ_e[dϕ/dr]≠0. This condition is met inside resonant cavities, since ${\mathcal{K}}_{\text{r}}^2>0$. It is also worth noting by differentiating Eq. (9) with respect to r and using Eq. (3) that their Wronskian 𝒲_ψ is conserved in resonant cavities, as expected from the Liouville formula.

2.2.2. Asymptotic form of the wave solution

In resonant cavities, we denote ${\mathcal{K}}_{\mathrm{r}}=\sqrt{{\mathcal{K}}_{\mathrm{r}}^2}$. We assume that the phase ϕ in Eq. (8) varies on a length scale on the order of 1/𝒦_r, which is supposed to be much smaller than the variation scale height of the medium on the order of H_p. Under these considerations, injecting Eqs. (8) in (3), we can show that the ϕ function in the cavity is equal at leading order in 𝒦_r, up to addition by a constant, to (e.g., Gough 2007)

$$\begin{array}{c}\hfill \varphi (r;\overline{r})\approx \frac{\mathrm{i}}{2}ln{\mathcal{K}}_{\mathrm{r}}+\phi (r;\overline{r}),\end{array}$$(10)

where

$$\begin{array}{c}\hfill \phi (r;\overline{r})={\displaystyle {\int }_{\overline{r}}^r{\mathcal{K}}_{\mathrm{r}}\mathrm{d}r,}\end{array}$$(11)

with $\overline{r}$ an arbitrary reference point in the considered cavity. In the latter equation, we have chosen as a convention the positive branch in Eq. (10) such as (dφ/dr) = 𝒦_r > 0 and the ψ function is associated with the progressive component of the wave function Ψ whose phase travels upward (given the temporal dependence on e^−iσt). Therefore, in resonant cavities, the ψ function merely reads at leading order

$$\begin{array}{c}\hfill \psi (r;\overline{r})\approx \frac{1}{\sqrt{{\mathcal{K}}_{\mathrm{r}}}}{e}^{\mathrm{i}\phi (r;\overline{r})},\end{array}$$(12)

which corresponds to a plane wave with a slowly varying amplitude. Given Eq. (10), the Wronskian of the two basis solutions in Eq. (9) is equal at leading order to 𝒲_ψ ≈ 2i𝒦_r/𝒦_r = 2i, which is constant as expected.

2.2.3. Wave energy luminosity

Physically speaking, the propagation of gravito-avoustic waves is not only described by the propagation of their phases, but also by the propagation of their energy. To analyze this, we find convenient to consider the radial wave energy luminosity. Within the Cowling approximation (resulting from the short-wavelength hypothesis), it is defined at leading order by the integrated quantity at time t and radius r (e.g., Lighthill 1978; Unno et al. 1989):

$$\begin{array}{c}\hfill \overline{{\mathcal{L}}_{\mathrm{w}}}\approx \frac{1}{T}{\displaystyle {\int }_{t-T/2}^{t+T/2}(\int {\int }_{\theta ,s}{\mathcal{R}}_{\mathrm{e}}\left[p^{\prime }(\mathit{r},\tau )\right]{\mathcal{R}}_{\mathrm{e}}\left[v_{\mathrm{r}}(\mathit{r},\tau )\right]r^{2}sin\theta \mathrm{d}\theta \mathrm{d}s)\mathrm{d}\tau ,}\end{array}$$(13)

where v_r is the Eulerian perturbation of radial velocity and T is the oscillation period. By recalling that the wave velocity field is given by v = (∂ξ/∂t) and by expressing within the considered framework its radial component as

$$\begin{array}{c}\hfill v_{\mathrm{r}}(\mathit{r},t)={\stackrel{\sim }{v}}_{\mathrm{r}}(r)Y_{\ell }^m(\theta ,s)e^{-\mathrm{i}\sigma t},\end{array}$$(14)

we obtain the simple relation

$$\begin{array}{c}\hfill {\stackrel{\sim }{v}}_{\mathrm{r}}(r)=-\mathrm{i}\sigma {\stackrel{\sim }{\xi }}_{\mathrm{r}}(r).\end{array}$$(15)

In order to express Eq. (13) further within the JWKB approximation, we then write the wave function Ψ in a general way as the linear combination of ψ and ψ^⋆, that is,

$$\begin{array}{c}\hfill \mathrm{\Psi }(r,t)\approx e^{-\mathrm{i}\sigma t}[a_{\mathrm{p}}\psi (r,\overline{r})+a_{\mathrm{r}}{\psi }^{\star }(r;\overline{r})],\end{array}$$(16)

where a_p and a_r are two complex constants representing the amplitudes of the progressive and regressive components, respectively. Using Eqs. (4), (5), (12), and (15), the (complex) radial part of the velocity in resonant cavities is therefore equal to:

$$\begin{array}{c}\hfill {\stackrel{\sim }{v}}_{\mathrm{r}}\approx -\mathrm{i}{\left(\frac{\sigma }{\rho r^{2}c}\right)}^{1/2}{\left|\frac{S_{\ell }^2-{\sigma }^2}{N^2-{\sigma }^2}\right|}^{1/4}(a_{\mathrm{p}}{e}^{\mathrm{i}\phi }+a_{\mathrm{r}}{e}^{-\mathrm{i}\phi }).\end{array}$$(17)

To obtain the Eulerian perturbation of pressure, we can then use the relations derived by Shibahashi (1979), and given in Eqs. (A.5) and (A.6), while neglecting at leading order the variations of the structure equilibrium compared to that of the wave phase, which provides

$$\begin{array}{c}\hfill {\stackrel{\sim }{p}}^{\prime }\approx \mathrm{sgn}(S_{\ell }^2-{\sigma }^2)\mathrm{i}{\left(\frac{\sigma \rho c}{r^2}\right)}^{1/2}{\left|\frac{N^2-{\sigma }^2}{S_{\ell }^2-{\sigma }^2}\right|}^{1/4}(a_{\mathrm{p}}{e}^{\mathrm{i}\phi }-a_{\mathrm{r}}{e}^{-\mathrm{i}\phi }),\end{array}$$(18)

where sgn() is the sign function. Therefore, using Eqs. (1), (14), (17), and (18), and taking advantage of the time dependence on e^−iσt and the orthonormality of the spherical harmonics, it is straightforward to show that the radial wave energy luminosity in Eq. (13) finally reads

$$\begin{array}{c}\hfill \overline{{\mathcal{L}}_{\mathrm{w}}}\approx \mathrm{sgn}({\sigma }^2-S_{\ell }^2)\frac{\sigma }{2}(|a_{\mathrm{p}}{|}^2-|a_{\mathrm{r}}{|}^2).\end{array}$$(19)

First, we check that the radial wave energy luminosity does not depend on r and t within the hypothesis of adiabatic oscillations, as expected. Second, we see that it is equal to the sum of two distinct parts that are proportional to |a_p|² and |a_r|², which result, respectively, from the progressive and regressive components². Third, we see that in general, the mean wave energy flux is provided by

$$\begin{array}{c}\hfill \overline{{\mathcal{F}}_{\mathrm{w}}}(r)\equiv \frac{\overline{{\mathcal{L}}_{\mathrm{w}}}}{4\pi r^2}\approx \frac{1}{16\pi }({\stackrel{\sim }{p}}^{\prime }{\stackrel{\sim }{v}}_{\mathrm{r}}^{\star }+{\stackrel{\sim }{p}}^{\prime \star }{\stackrel{\sim }{v}}_{\mathrm{r}}).\end{array}$$(20)

In the case where there is only one progressive or one regressive component (i.e., a_r = 0 or a_p = 0, respectively), we retrieve the well-known formula for plane waves given by ${\overline{\mathcal{F}}}_{\mathrm{w}}=\stackrel{\sim }{p}\prime {\stackrel{\sim }{v}}_{\mathrm{r}}^{\star }/8\pi$. Finally, it is worth mentioning that the Wronskian for the wave function Ψ is equal to

$$\begin{array}{c}\hfill {\mathcal{W}}_{\mathrm{\Psi }}=({\mathrm{\Psi }}^{\star }\frac{\mathrm{d}\mathrm{\Psi }}{\mathrm{d}r}-\mathrm{\Psi }\frac{\mathrm{d}{\mathrm{\Psi }}^{\star }}{\mathrm{d}r})=(|a_{\mathrm{p}}{|}^2-|a_{\mathrm{r}}{|}^2){\mathcal{W}}_{\psi }.\end{array}$$(21)

At leading order in the resonant cavities, 𝒲_ψ ≈ 2i, and we see that the Wronskian of the wave function Ψ is thus directly related to the radial wave energy luminosity, which explains its conservation within the adiabatic limit.

Another important point is the presence of the term $\text{sgn}({\sigma }^2-S_{\mathcal{l}}^2)$ in Eq. (19). This term is actually related to the direction of the propagation of the wave energy. Indeed, the radial group and phase velocities in resonant cavities, v_g and v_φ, are equal, according to Eq. (5), to:

$$\begin{array}{c}\hfill v_{\phi }=\pm \frac{\sigma }{{\mathcal{K}}_{\mathrm{r}}}\mathrm{and}{v}_{\mathrm{g}}=\frac{\partial \sigma }{\partial {\mathcal{K}}_{\mathrm{r}}}=\frac{c^2{\mathcal{K}}_{\mathrm{r}}^2{\sigma }^2}{({\sigma }^4-N^2{S}_{\ell }^2)}{v}_{\phi },\end{array}$$(22)

where the plus and minus signs correspond to the progressive and regressive components, respectively. Therefore, two cases have to be distinguished. On the one hand, in the case of pressure-dominated oscillations, we have ${\sigma }^2\gg (S_{\mathcal{l}}^2\text{and}{N}^2)$ so that v_g and v_φ have the same sign. As a consequence, the wave energy propagates in the same direction as the phase wavefront. Therefore, the progressive (regressive) component is associated with a positive (negative) mean wave energy flux propagating upward (downward). On the other hand, in the case of gravity-dominated oscillations, we have ${\sigma }^2\ll (S_{\mathcal{l}}^2\text{and}{N}^2)$ so that v_g and v_φ have opposite signs. The wave energy propagates in the opposite direction of the phase wavefront. In this case, the progressive (regressive) component is associated with a negative (positive) mean wave energy flux propagating downward (upward). This behavior is in agreement with the term $\text{sgn}({\sigma }^2-S_{\mathcal{l}}^2)$ in Eq. (19).

2.3. Barriers and basic wave reflection-transmission problem

There are some regions in stars where the waves are evanescent (i.e., where ${\mathcal{K}}_{\text{r}}^2<0$) or where the JWKB approximation is not met, as in the vicinity of sharp variations in the medium (e.g., discontinuity in the chemical composition, thin ionization region) or turning-points (i.e., where ${\mathcal{K}}_{\text{r}}^2=0$). These are the so-called barriers. While in the case of one infinite-length cavity a purely progressive or regressive solution can exist (i.e., as in Eq. (16) with a_r = 0 or a_p = 0), this does not hold true in the presence of one single barrier. Indeed, an incident wave coming from a first cavity is partially reflected back near the barrier and partially transmitted toward a second cavity on the other side of the barrier. It is thus necessary to consider the reflected and transmitted waves as well, the amplitudes of which depend both on the amplitude of the incident wave and on the properties of the barrier, and they have to satisfy the conservation of the wave energy flux. In general, this physical problem can be described by reflection and transmission coefficients, as already proposed by Roxburgh & Vorontsov (2001) or Takata (2016b) in the context of stellar pulsations.

To describe the wave reflection-transmission problem, we considered a first cavity C_n underlying a second cavity C_n + 1 with n ∈ ℕ^⋆; both are separated by an intermediate barrier. We chose two reference radii $r_n^{+}$ or $r_{n+1}^{-}$ that are located below and above the considered barrier, at the interfaces with the resonant cavities C_n and C_n + 1, respectively; these layers define the upper and lower boundaries of the cavities C_n and C_n + 1. Within such a framework, we can first express the wave function inside the cavity C_n in a similar way to Eq. (16), but with an explicit representation in terms of the propagation of the energy and with the origin of the wave phase set to the upper boundary, $r_n^{+}$, that is,

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_n(r,t)=a_{n,+}^{(\to )}{\psi }^{(\to )}(r,t;r_n^{+})+a_{n,+}^{(\leftarrow )}{\psi }^{(\leftarrow )}(r,t;r_n^{+}),\end{array}$$(23)

where the superscripts ( → ) and ( ← ) denote the wave components in the cavity C_n whose energy propagates upward and downward with the complex amplitudes $a_{n,+}^{(\to )}$ and $a_{n,+}^{(\leftarrow )}$, respectively. As the phase and group velocities for pressure-dominated (gravity-dominated) modes have the same (opposite) direction, we can directly conclude:

$$\begin{array}{cc}\hfill {\psi }^{(\to )}(r,t;\overline{r})& =\{\begin{array}{cc}\psi (r;\overline{r})e^{-\mathrm{i}\sigma t}\hfill & \mathrm{for}\mathrm{pressure}\mathrm{modes}\hfill \\ {\psi }^{\star }(r;\overline{r})e^{-\mathrm{i}\sigma t}\hfill & \mathrm{for}\mathrm{gravity}\mathrm{modes}\hfill \end{array}\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\psi }^{(\leftarrow )}(r,t;\overline{r})& ={\psi }^{(\to )\star }(r,-t;\overline{r}),\hfill \end{array}$$(25)

where we recall that the $\psi (r;\overline{r})$ function is defined in Eq. (12). The same representation can be used to express the wave function inside the cavity C_n + 1, using the lower boundary $r_{n+1}^{-}$ as the origin of the wave phase, that is,

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{n+1}(r,t)=a_{n+1,-}^{(\to )}{\psi }^{(\to )}(r,t;r_{n+1}^{-})+a_{n+1,-}^{(\leftarrow )}{\psi }^{(\leftarrow )}(r,t;r_{n+1}^{-}),\end{array}$$(26)

where $a_{n+1,-}^{(\to )}$ and $a_{n+1,-}^{(\leftarrow )}$ denote the complex amplitudes of the wave components in the cavity C_n + 1 whose energy propagates upward and downward, respectively. We emphasize that Ψ_n and Ψ_n + 1 are two different representations inside different cavities of the same and unique global wave function Ψ, which denotes the solution of the wave equation throughout the star.

Such a representation of the oscillations in the resonant cavities is appropriate to define properly the reflection and transmission coefficients while accounting for the conservation of the wave energy flux, as already formalized in Takata (2016b). First, in the base wave reflection-transmission problem, an upward incident energy ray propagating in the bottom cavity C_n encounters the barrier and is reflected back downward, while a part of the energy ray is transmitted through the barrier and propagates upward in the overlying resonant cavity C_n + 1. In other words, we assume $a_{n+1,-}^{(\leftarrow )}=0$ (see Fig. 1). In this base configuration, the reflection and transmission coefficients, R_{n, n + 1} and T_{n, n + 1}, are respectively defined as the ratio of the amplitudes of the reflected and transmitted components to the amplitude of the incident component inside the resonant cavities, that is,

$$\begin{array}{cc}\hfill R_{n,n+1}& =\frac{a_{n,+}^{(\leftarrow )}}{a_{n,+}^{(\to )}}=\left|R_{n,n+1}\right|e^{\mathrm{i}{\delta }_{n,n+1}}\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill T_{n,n+1}& =\frac{a_{n+1,-}^{(\to )}}{a_{n,+}^{(\to )}}=\left|T_{n,n+1}\right|e^{\mathrm{i}{\gamma }_{n,n+1}},\hfill \end{array}$$(28)

where δ_{n, n + 1} and γ_{n, n + 1} are the phase lags introduced at reflection and transmission, respectively. In this work, we will assume that the phase lags take values between [ − π, π]. We note that the notion of phase lags is readily defined well inside the cavities where the JWKB is met and where the progressive and regressive wave components are distinctly defined. In addition, the conservation of the mean wave energy flux throughout the star, provided by Eq. (19) inside both cavities, translates into the constraint:

$$\begin{array}{c}\hfill {\left|a_{n,+}^{(\to )}\right|}^2-{\left|a_{n,+}^{(\leftarrow )}\right|}^2={\left|a_{n+1,-}^{(\to )}\right|}^2\Rightarrow {\left|R_{n,n+1}\right|}^2+{\left|T_{n,n+1}\right|}^2=1.\end{array}$$(29)

Second, in the adjoint wave reflection-transmission problem, a downward incident energy ray propagating in the resonant cavity C_n + 1 encounters the barrier and is reflected back upward, while a part of the energy ray is transmitted through the barrier and propagates downward in the underlying resonant cavity C_n. In other words, we assume $a_{n,+}^{(\to )}=0$ (see Fig. 2). In this configuration, the reflection and transmission coefficients, R_{n + 1, n} and T_{n + 1, n}, are defined, respectively, as

$$\begin{array}{cc}\hfill R_{n+1,n}& =\frac{a_{n+1,-}^{(\to )}}{a_{n+1,-}^{(\leftarrow )}}=\left|R_{n+1,n}\right|e^{\mathrm{i}{\delta }_{n+1,n}}\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill T_{n+1,n}& =\frac{a_{n,+}^{(\leftarrow )}}{a_{n+1,-}^{(\leftarrow )}}=\left|T_{n+1,n}\right|e^{\mathrm{i}{\gamma }_{n+1,n}},\hfill \end{array}$$(31)

where δ_{n + 1, n} and γ_{n + 1, n} are the phase lags introduced at reflection and transmission. The modulus of the reflection and transmission coefficients must also satisfy the conservation of the wave energy flux as in Eq. (29).

Using basic wave principles of time-reversal symmetry, linear superposition and energy conservation, it is possible using a reasoning similar to Takata (2016b, namely, their Sect. 2), to connect the wave coefficients of the base problem and those of the adjoint problem. As demonstrated in Appendix B, they follow the relations

$$\begin{array}{cccc}\hfill \left|R_{n+1,n}\right|& =\left|R_{n,n+1}\right|,\hfill & \hfill \left|T_{n+1,n}\right|& =\left|T_{n,n+1}\right|,\hfill \\ \hfill {\gamma }_{n+1,n}& ={\gamma }_{n,n+1},\hfill & \hfill {\delta }_{n+1,n}& =\pi -{\delta }_{n,n+1}+2{\gamma }_{n,n+1},\hfill \end{array}$$(32)

which is similar to the results found by Takata (2016b) when γ_{n, n + 1} = 0. As a result, we finally see that a barrier between two cavities C_n and C_n + 1 is entirely characterized by only three parameters: R_{n, n + 1}, δ_{n, n + 1}, and γ_{n, n + 1}. The modulus of the transmission coefficient can then be retrieved using the conservation of the wave energy flux.

To complete these definitions, we emphasize that the choice of the origin of the phase in the cavities, that is, the so- called cavity boundaries within our framework, is somehow arbitrary and that any change in the latter is compensated by a modification of the values of the phase lags, δ_{n, n + 1} and γ_{n, n + 1}; the final physical solution of the problem in contrast does not depend on this choice. We note that in practice, it is often convenient to take the boundaries as equal to the turning-points or the middle radius of glitches, as we see later in Sect. 6.

2.4. Central and surface stellar boundary conditions

Close to the center and surface of stars, the wave energy luminosity must vanish at some point. This condition is required by the regularity of the oscillation displacement near the center and the fact that the density vanishes beyond the surface. The core and surface of stars can thus be modeled as totally reflective barriers (e.g., Unno et al. 1989).

Inside the first cavity just above the stellar core, which is denoted by C₁, the wave function can be generally written within the considered convention as

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_1(r,t)=a_{1,-}^{(\to )}{\psi }^{(\to )}(r,t;r_1^{-})+a_{1,-}^{(\leftarrow )}{\psi }^{(\leftarrow )}(r,t;r_1^{-}).\end{array}$$(33)

The reflective boundary condition at the center then requires

$$\begin{array}{c}\hfill a_{1,-}^{(\to )}=e^{\mathrm{i}{\delta }_{\mathrm{c}}}{a}_{1,-}^{(\leftarrow )},\end{array}$$(34)

where δ_c is the phase lag introduced during the reflection. This ensures that the modulus of the amplitudes of the upward and downward components are equal and hence that the incident wave energy flux is totally reflected. Similarly, inside the last cavity just below the stellar surface where the stellar density vanishes, which is denoted by C_N with N ∈ ℕ^⋆, the wave function reads

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_N(r,t)=a_{N,+}^{(\to )}{\psi }^{(\to )}(r,t;r_N^{+})+a_{N,+}^{(\leftarrow )}{\psi }^{(\leftarrow )}(r,t;r_N^{+}).\end{array}$$(35)

The reflective boundary condition at the surface then requires

$$\begin{array}{c}\hfill a_{N,+}^{(\leftarrow )}=e^{\mathrm{i}{\delta }_{\mathrm{s}}}{a}_{N,+}^{(\to )},\end{array}$$(36)

where δ_s is the phase lag introduced at the reflection.

At this point, all the basic ingredients have been introduced to formulate in a general way the resonance condition of global modes in stars by taking simultaneously into account the effect of an ensemble of barriers, as we see in the next sections.

3. Infinite-time reflection picture for multi-cavity oscillation modes

In this section, we consider a star composed of N successive resonant cavities {C_i}_{1 ≤ i ≤ N} between the center and surface, which are separated from each other by N − 1 intermediate barriers. In order to make explicit the properties of the oscillation eigenmodes, we follow the propagation of an incident energy ray along its infinite-time travel through the ensemble of cavities and impose a condition for constructive interferences. The problem is graphically represented in Fig. 3. This is actually very similar to the computation of the transmission transfer function of a Fabry-Pérot in optics, but with an ensemble of resonant cavities and totally reflective boundaries.

3.1. Cavities C₁ and C₂

As a first step, we focus on the two lowest cavities C₁ and C₂ while considering the reflective condition at the bottom boundary of C₁ (i.e., close to the stellar center). In the following, $\{{\stackrel{\sim }{r}}_i\}$ denotes a set of given radii in the middle of the cavities {C_i}. We then consider the case of an incident downward energy ray at ${\stackrel{\sim }{r}}_2$ and at a given time t₀ with an amplitude $a_{2,-}^{(\leftarrow )}$. According to Eq. (26), the wave function at ${\stackrel{\sim }{r}}_2$ is provided by

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_2({\stackrel{\sim }{r}}_2,t_0)=a_{2,-}^{(\to )}{\psi }^{(\to )}({\stackrel{\sim }{r}}_2,t_0;r_2^{-})+a_{2,-}^{(\leftarrow )}{\psi }^{(\leftarrow )}({\stackrel{\sim }{r}}_2,t_0;r_2^{-}).\end{array}$$(37)

The goal is then to express $a_{2,-}^{(\to )}$ as a function of $a_{2,-}^{(\leftarrow )}$. Unlike the adjoint wave reflection-transmission configuration presented in Sect. 2.3, the amplitude of the reflected component has to account for the reflective boundary condition at the stellar center and the fact that the wave energy that is transmitted from C₂ to C₁ can come back into C₂. This reflected part can thus be decomposed into two components, namely as:

$$\begin{array}{c}\hfill a_{2,-}^{(\to )}={\alpha }_{2,R}^{(\to )}+{\alpha }_{2,T}^{(\to )}.\end{array}$$(38)

The first term ${\alpha }_{2,R}^{(\to )}$ corresponds to the part of the incident energy ray that is directly reflected on the barrier between C₁ and C₂ as in the associated adjoint wave reflection-transmission problem, that is,

$$\begin{array}{c}\hfill {\alpha }_{2,R}^{(\to )}=R_{2,1}{a}_{2,-}^{(\leftarrow )}.\end{array}$$(39)

The second term ${\alpha }_{2,T}^{(\to )}$ corresponds to the part of the incident energy ray that is transmitted to C₁, that then indefinitely travels back and forth between the center and the upper boundary of C₁, and that is finally transmitted back to C₂. Based on Eq. (23), the wave function at ${\stackrel{\sim }{r}}_1$ in C₁ can be written as

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_1({\stackrel{\sim }{r}}_1,t_0)=a_{1,+}^{(\to )}{\psi }^{(\to )}({\stackrel{\sim }{r}}_1,t_0;r_1^{+})+a_{1,+}^{(\leftarrow )}{\psi }^{(\leftarrow )}({\stackrel{\sim }{r}}_1,t_0;r_1^{+}).\end{array}$$(40)

The wave amplitudes in Eq. (33) can be expressed as the superposition of an infinite number of contributions resulting from the infinite-time multiple reflections of the energy ray transmitted from C₂, that is,

$$\begin{array}{cc}\hfill a_{1,+}^{(\leftarrow )}& ={\alpha }_{1,T}^{(\leftarrow )}+{\displaystyle \sum _{k=1}^{+\infty }{\alpha }_{1,k}^{(\leftarrow )}}\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill a_{1,+}^{(\to )}& ={\alpha }_{1,R}^{(\to )}+{\displaystyle \sum _{k=2}^{+\infty }{\alpha }_{1,k}^{(\to )}.}\hfill \end{array}$$(42)

First, ${\alpha }_{1,T}^{(\leftarrow )}$ corresponds to the part of the incident energy ray that is transmitted from C₂ to C₁, that is,

$$\begin{array}{c}\hfill {\alpha }_{1,T}^{(\leftarrow )}=T_{2,1}{a}_{2,-}^{(\leftarrow )}.\end{array}$$(43)

Second, ${\alpha }_{1,R}^{(\to )}$ corresponds to the transmitted part that is in addition reflected close to the core and returns back to ${\stackrel{\sim }{r}}_1$. To express this amplitude as a function of the wave coefficients within the conventions presented in Sect. 2, we need to change the origin of the phase in Eq. (40) from $r_1^{+}$ to $r_1^{-}$ where the reflection occurs. From Eqs. (12), (24), and (25), we deduce for any integer i that

$$\begin{array}{cc}\hfill {\psi }^{(\to )}({\stackrel{\sim }{r}}_i,t_0;r_i^{+})& =e^{-\mathrm{i}{\mathrm{\Theta }}_i}{\psi }^{(\to )}({\stackrel{\sim }{r}}_i,t_0;r_i^{-}),\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill {\psi }^{(\leftarrow )}({\stackrel{\sim }{r}}_i,t_0;r_i^{+})& =e^{+\mathrm{i}{\mathrm{\Theta }}_i}{\psi }^{(\leftarrow )}({\stackrel{\sim }{r}}_i,t_0;r_i^{-}),\hfill \end{array}$$(45)

with

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_i=\pm {\displaystyle {\int }_{r_i^{-}}^{r_i^{+}}{\mathcal{K}}_{\mathrm{r}}\mathrm{d}r,}\end{array}$$(46)

where the plus (minus) sign has to be chosen when the radial group and phase velocities have the same (opposite) directions, that is, in cases of pressure-dominated (gravity-dominated) modes. In comparing Eq. (33) with the expression of the wave function at ${\stackrel{\sim }{r}}_1$, but with $r_1^{-}$ as the origin of the wave phase, that is,

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_1({\stackrel{\sim }{r}}_1,t_0)=a_{1,-}^{(\to )}{\psi }^{(\to )}({\stackrel{\sim }{r}}_1,t_0;r_1^{-})+a_{1,-}^{(\leftarrow )}{\psi }^{(\leftarrow )}({\stackrel{\sim }{r}}_1,t_0;r_1^{-}),\end{array}$$(47)

we then can understand, by using Eqs. (44) and (45), that a change in the origin of the phase from the upper boundary of the cavity C_i, $r_i^{+}$, to the lower boundary, $r_i^{-}$, is associated with the transformation of the amplitudes

$$\begin{array}{cc}\hfill a_{i,-}^{(\to )}& =e^{-\mathrm{i}{\mathrm{\Theta }}_i}{a}_{i,+}^{(\to )},\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill a_{i,-}^{(\leftarrow )}& =e^{+\mathrm{i}{\mathrm{\Theta }}_i}{a}_{i,+}^{(\leftarrow )}.\hfill \end{array}$$(49)

Starting from the amplitude ${\alpha }_{T,1}^{(\leftarrow )}$, we first change the origin of the phase to $r_1^{-}$, which is equivalent to multiply the amplitude by e^iΘ₁ according to Eq. (49). With this convention, the reflected amplitude is just lagged by a phase δ_c compared to the incident one according to Sect. 2.4. To conclude, we have to change the origin of the phase back to $r_1^{+}$ to retrieve the convention used in Eq. (40), which is equivalent to multiply the amplitude again by e^iΘ₁ according to Eq. (48). As a result, this gives

$$\begin{array}{c}\hfill {\alpha }_{1,R}^{(\to )}=e^{2\mathrm{i}{\mathrm{\Theta }}_1}{e}^{\mathrm{i}{\delta }_{\mathrm{c}}}{\alpha }_{1,T}^{(\leftarrow )}.\end{array}$$(50)

Third, ${\alpha }_{1,k}^{(\leftarrow )}$ and ${\alpha }_{1,k}^{(\to )}$ in Eqs. (41) and (42) correspond to the parts of the downward and upward components having traveled back and forth k times inside C₁. For instance, the downward term having made only one back and forth is deduced from the reflection of the amplitude ${\alpha }_{1,R}^{(\to )}$ on the barrier between C₁ and C₂, that is,

$$\begin{array}{cc}\hfill {\alpha }_{1,1}^{(\leftarrow )}& =R_{1,2}{\alpha }_{1,R}^{(\to )}=R_{1,2}{e}^{2\mathrm{i}{\mathrm{\Theta }}_1}{e}^{\mathrm{i}{\delta }_{\mathrm{c}}}{\alpha }_{1,T}^{(\leftarrow )}.\hfill \end{array}$$(51)

We retrieve a factor e^2iΘ₁ that results from the total “optic” path covered by the initial downward energy ray transmitted from C₂ to come back to its initial position during a back and forth in C₁. Based on the same reasoning as before, it is straightforward to express the other components by recurrence for k > 1 as:

$$\begin{array}{cc}\hfill {\alpha }_{1,k}^{(\leftarrow )}& =R_{1,2}{e}^{2\mathrm{i}{\mathrm{\Theta }}_1}{e}^{\mathrm{i}{\delta }_{\mathrm{c}}}{\alpha }_{1,k-1}^{(\leftarrow )}\hfill \end{array}$$(52)

$$\begin{array}{cc}\hfill {\alpha }_{1,k}^{(\to )}& =e^{2\mathrm{i}{\mathrm{\Theta }}_1}{e}^{\mathrm{i}{\delta }_{\mathrm{c}}}{\alpha }_{1,k-1}^{(\leftarrow )}.\hfill \end{array}$$(53)

As a consequence, the upward and downward amplitudes in C₁ can be linked through Eqs. (41), (42), (50), and (53), which merely results in

$$\begin{array}{c}\hfill a_{1,+}^{(\to )}=e^{2\mathrm{i}{\mathrm{\Theta }}_1}{e}^{\mathrm{i}{\delta }_{\mathrm{c}}}{a}_{1,+}^{(\leftarrow )}.\end{array}$$(54)

In other words, the upward and downward components in C₁ have the same modulus and are just phase lagged during the infinite-time travel of the wave energy ray inside C₁. This means that the downward and upward energy rays carry the same amount of energy (but in the opposite direction) in such a way that the total wave luminosity vanishes in C₁, as expected from Eq. (19) and the conservation of the wave energy flux under the totally reflective core constraint. Finally, the amplitude of the wave transmitted to C₂ is equal to:

$$\begin{array}{c}\hfill {\alpha }_{2,T}^{(\to )}=T_{1,2}{a}_{1,+}^{(\to )},\end{array}$$(55)

so that using Eqs. (42), (43) and (50)–(53), and expressing T_2, 1 as a function of T_1, 2 through Eq. (32), it can be rewritten as:

$$\begin{array}{cc}\hfill {\alpha }_{2,T}^{(\to )}& =a_{2,-}^{(\leftarrow )}{\left|T_{1,2}\right|}^2{e}^{\mathrm{i}({\delta }_{\mathrm{c}}+2{\gamma }_{1,2}+2{\mathrm{\Theta }}_1)}{\displaystyle \sum _{k=0}^{+\infty }{\left[R_{1,2}{e}^{\mathrm{i}({\delta }_{\mathrm{c}}+2{\mathrm{\Theta }}_1)}\right]}^k}\hfill \\ \hfill & =a_{2,-}^{(\leftarrow )}{\left|T_{1,2}\right|}^2\frac{e^{\mathrm{i}({\delta }_{\mathrm{c}}+2{\gamma }_{1,2}+2{\mathrm{\Theta }}_1)}}{1-R_{1,2}{e}^{\mathrm{i}({\delta }_{\mathrm{c}}+2{\mathrm{\Theta }}_1)}}·\hfill \end{array}$$(56)

Therefore, the total amplitude of the upward component in C₂ is equal to, according to Eqs. (32), (38), (39), and (56):

$$\begin{array}{c}\hfill a_{2,-}^{(\to )}=a_{2,-}^{(\leftarrow )}(\left|R_{1,2}\right|e^{\mathrm{i}(\pi -{\delta }_{1,2}+2{\gamma }_{1,2})}+\frac{{\left|T_{1,2}\right|}^2{e}^{\mathrm{i}({\delta }_{\mathrm{c}}+2{\gamma }_{1,2}+2{\mathrm{\Theta }}_1)}}{1-\left|R_{1,2}\right|e^{\mathrm{i}({\delta }_{\mathrm{c}}+{\delta }_{1,2}+2{\mathrm{\Theta }}_1)}}).\end{array}$$(57)

Using the energy constraint |T_1, 2|² = 1 − |R_1, 2|², Eq. (57) leads after some manipulations to

$$\begin{array}{c}\hfill a_{2,-}^{(\to )}=e^{\mathrm{i}{\mathrm{\Delta }}_{2,1}}{a}_{2,-}^{(\leftarrow )},\end{array}$$(58)

where Δ_2, 1 represents the total phase lag introduced by the reflection on the intermediate barrier and the infinite-time travel inside C₁ of the initial downward incident wave. It reads

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{2,1}=\pi +2{\gamma }_{1,2}-{\delta }_{1,2}-2arctan(q_{1,2}cot{\mathrm{\Phi }}_1),\end{array}$$(59)

with

$$\begin{array}{cc}& {\mathrm{\Phi }}_1={\mathrm{\Theta }}_1+\frac{{\delta }_{\mathrm{c}}}{2}+\frac{{\delta }_{1,2}}{2}\hfill \end{array}$$(60)

$$\begin{array}{cc}& q_{1,2}=\frac{1-\left|R_{1,2}\right|}{1+\left|R_{1,2}\right|}·\hfill \end{array}$$(61)

The factor q_1, 2 in the last equation is called the coupling factor of the cavities C₁ and C₂, which can take values between zero and unity. Equation (58) also shows that as in C₁, the upward and downward components in C₂ have the same modulus because of both the totally reflective core constraint and the energy conservation; the total wave luminosity thus also vanishes in C₂.

We thus conclude that the combination of the cavity C₁ and the overlying intermediate barrier is equivalent to one single totally reflective barrier underlying the cavity C₂ and introducing an effective phase lag Δ_2, 1 at the reflection of a downward incident wave.

3.2. Ensemble of N cavities

As the next step, we add a third cavity C₃ and an intermediate barrier above C₂. We then consider the case of an incident downward energy ray at ${\stackrel{\sim }{r}}_3$ with an amplitude $a_{3,-}^{(\leftarrow )}$ at a given time t₀. According to Eq. (26), the wave function at ${\stackrel{\sim }{r}}_3$ is provided by

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_3({\stackrel{\sim }{r}}_3,t_0)=a_{3,-}^{(\to )}{\psi }^{(\to )}({\stackrel{\sim }{r}}_3,t;r_3^{-})+a_{3,-}^{(\leftarrow )}{\psi }^{(\leftarrow )}({\stackrel{\sim }{r}}_3,t_0;r_3^{-}).\end{array}$$(62)

As previously, the goal is to express the total reflected amplitude $a_{3,-}^{(\to )}$ as a function of the incident amplitude $a_{3,-}^{(\leftarrow )}$. Actually, the computation is similar to that performed in Sect. 3.1 when considering two cavities because we have previously shown that the cavity C₁ and the overlying barrier just below C₂ together can be represented as a totally reflective barrier associated with an effective phase lag Δ_2, 1. Therefore, we can conclude in a straightforward way similarly to Eq. (58) that

$$\begin{array}{c}\hfill a_{3,-}^{(\to )}=a_{3,-}^{(\leftarrow )}{e}^{\mathrm{i}{\mathrm{\Delta }}_{3,2}},\end{array}$$(63)

where Δ_3, 2 is the total phase lag introduced during the reflection on the barrier between C₂ and C₃ and the infinite-time travel throughout the cavities C₁ and C₂. Its expression is provided by the set of Eqs. (59)–(61) while replacing δ_c by Δ_2, 1 and the subscripts 1 and 2 by the subscripts 2 and 3, respectively (e.g., |R_1, 2| must be replaced by |R_2, 3|). The region between the stellar core and the barrier underlying the cavity C₃ can thus also be considered as a totally reflective barrier associated with an effective phase lag at reflection Δ_3, 2.

At this point, it is then obvious that the generalization to the case of N cavities can be obtained by adding one by one supplementary overlying cavities and using the same reasoning as before at each step. To do so, we consider a downward incident energy ray in the cavity C_N. Analogously to Eqs. (58) and (63), we understand that the amplitudes of the downward and upward components in each cavity are linked for 1 ≤ i ≤ N by the expression

$$\begin{array}{c}\hfill a_{i,-}^{(\to )}=e^{\mathrm{i}{\mathrm{\Delta }}_{i,i-1}}{a}_{i,-}^{(\leftarrow )},\end{array}$$(64)

where the total effective phase lag, Δ_{i, i − 1}, introduced by the infinite-time reflections and back-and-forth travels through all the underlying cavities, {C_j}_j < i, can be formulated by the recurrence relation for i ≥ 2:

$$\begin{array}{cc}& {\mathrm{\Delta }}_{i,i-1}=\pi +2{\gamma }_{i-1,i}-{\delta }_{i-1,i}-2arctan(q_{i-1,i}cot{\mathrm{\Phi }}_{i-1})\hfill \end{array}$$(65)

$$\begin{array}{cc}& {\mathrm{\Phi }}_{i-1}={\mathrm{\Theta }}_{i-1}+\frac{{\mathrm{\Delta }}_{i-1,i-2}}{2}+\frac{{\delta }_{i-1,i}}{2}·\hfill \end{array}$$(66)

The general definition of the coupling factor between the cavities C_i − 1 and C_i is provided by

$$\begin{array}{c}\hfill q_{i-1,i}=\frac{1-\left|R_{i-1,i}\right|}{1+\left|R_{i-1,i}\right|}·\end{array}$$(67)

Finally, the initialization of the recurrence for i = 1 is ruled by the core boundary condition in Eq. (60), which is reduced to

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{1,0}={\delta }_{\mathrm{c}}.\end{array}$$(68)

We also note that Eq. (64) implies that the modulus of the upward and downward components in all cavities are equal and thus that the total wave energy luminosity is null everywhere. As mentioned before, this is the consequence of both the totally reflective core condition and the mean energy flux conservation.

3.3. Upper surface boundary and resonance condition

In the final step, we add the surface totally reflective barrier above the cavity C_N, so that Eq. (36) must apply. Simultaneously, using Eqs. (48), (49), and (64) for i = N, we also have to impose:

$$\begin{array}{c}\hfill a_{N,+}^{(\to )}=e^{2\mathrm{i}{\mathrm{\Theta }}_N}{e}^{\mathrm{i}{\mathrm{\Delta }}_{N,N-1}}{a}_{N,+}^{(\leftarrow )}.\end{array}$$(69)

This provides the following resonance condition:

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_N+\frac{{\delta }_{\mathrm{s}}}{2}+\frac{{\mathrm{\Delta }}_{N,N-1}}{2}=n\pi ,\end{array}$$(70)

where n is an integer corresponding to the mode radial order. In the adopted stationary configuration, the resonance condition is actually equivalent to consider that a downward incident energy ray in the cavity C_N has to travel throughout the ensemble of cavities and come back to its initial position with the exact same amplitude in order to constructively interfere. In order to give a usual physical meaning to the radial order n in Eq. (70), it may be convenient to choose the branch of the arctangent function in such way that for any real ε, we have:

$$\begin{array}{c}\hfill \underset{q\to 1}{lim}arctan(qtan\epsilon )=\epsilon .\end{array}$$(71)

Using this convention, the radial order n can be interpreted throughout the paper as the difference between the number of radial oscillation nodes in the p-dominated cavities, n_p, and that in the g-dominated cavities, n_g, over all the cavities (i.e., n = n_p − n_g). Using the principal branch for the artangent function, n in Eq. (70) would be instead interpreted as n_N, the number of radial nodes over the cavity C_N.

Hence, the mode eigenfrequency spectrum is obtained by solving the set of equations Eqs. (65), (66), (68), and (70) for any radial order n and accounting for the implicit frequency dependence of the wave number integral and the barrier parameters. To express the resonance condition in a practical form, we can define the phase Υ_i for 1 ≤ i ≤ N such as

$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}_i={\mathrm{\Theta }}_i+\frac{{\delta }_{i,i+1}}{2}+\frac{{\delta }_{i,i-1}}{2}-\frac{\pi }{2},\end{array}$$(72)

and we can set

$$\begin{array}{c}\hfill {\delta }_{1,0}={\delta }_{\mathrm{c}}\mathrm{and}{\delta }_{N,N+1}=\pi +{\delta }_{\mathrm{s}}.\end{array}$$(73)

Using the relation between δ_{i, i − 1} and γ_{i − 1, i} in Eq. (32), the recurrence resonance relation in Eqs. (65) and (66) can be written for 2 ≤ i ≤ N in the form of

$$\begin{array}{c}\hfill tan{\mathrm{\Phi }}_{i-1}=q_{i-1,i}tan({\mathrm{\Phi }}_i-{\mathrm{{\rm Y}}}_i),\end{array}$$(74)

where the core and surface boundary conditions in Eqs. (68) and (70) translate respectively into

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_1={\mathrm{{\rm Y}}}_1+\frac{\pi }{2}\mathrm{and}{\mathrm{\Phi }}_N=(n+\frac{1}{2})\pi .\end{array}$$(75)

For N = 1, Eq. (75) is sufficient alone and leads to Υ₁ = nπ. Depending on the problem, it can also be convenient to express the resonance condition as a series of frequency-dependent sine terms, as we see in Sect. 6.4. Such an alternative formulation is provided in Appendix C.

3.4. On the need for the multi-cavity approach

Before going further, it is worth discussing the need for the multi-cavity approach, since we show in Sects. 3.1 and 3.2 that it is possible to reinterpret the multi-cavity problem as the single-cavity problem. Indeed, the combination of the cavities and barriers below a given cavity C_i > 1 can always be reinterpreted as one single totally reflective barrier with an effective reflection phase lag redefined appropriately. Actually, the same conclusion holds true for the combination of the cavities and barriers above a given cavity, C_i < N. This can be easily shown following the same reasoning as in Sects. 3.1 and 3.2, except that the computation of the mode amplitude has to be made starting from the surface totally reflective condition and going toward deeper cavities. We may therefore wonder whether the multi-cavity picture, which is more complicated, is justified or not. The answer to this question depends on how the wave number integrals and the reflection coefficients vary with frequency over a considered range. On the one hand, if the wave number integral and the reflection coefficients on both sides of the considered cavity are constant over the frequency range of interest, then the single-cavity picture appears sufficient to describe the oscillations. On the other hand, if the mode parameters vary even slowly with frequency, then we can expect that the effective phase lags resulting from the combination of the cavities above and below the considered cavity behave in a complicated way with frequency. In this case, the use of the single-cavity picture is not judicious. In practice, for real stars, the wave number integrals and the reflection coefficients always vary with frequency and choosing a multi-cavity approach to describe the oscillation modes therefore appears necessary to understand their frequency spectra and develop useful seismic diagnoses.

4. Linear boundary value problem picture

In Sect. 3, we describe the eigenmodes using a physical picture that is analogous to ray tracing in optics. In this section, we aim to describe the eigenmodes in a more mathematical way that considers them as the solution to a linear boundary value problem. Besides checking the validity of the previous scenario, such an approach also has the advantage of being more convenient to discuss the distribution of the mode energy and the mode amplitudes ratios throughout stars, which is addressed in Sect. 5.

4.1. Amplitude vector

In each resonant cavity, the general solution for the wave function Ψ takes a similar form to Eqs. (23) and (26) depending on whether the origin of the phase is taken at the lower or upper boundary of the cavity, respectively. In each cavity C_i, we thus define the amplitude vector in both cases as

$$\begin{array}{c}\hfill {\mathit{a}}_{i,-}=\left(\begin{array}{c}{a}_{i,-}^{(\to )}\\ a_{i,-}^{(\leftarrow )}\end{array}\right)\mathrm{and}{\mathit{a}}_{i,+}=\left(\begin{array}{c}{a}_{i,+}^{(\to )}\\ a_{i,+}^{(\leftarrow )}\end{array}\right).\end{array}$$(76)

To find the stationary modes oscillating between the central and surface boundaries of the star, we need to make the link between all the amplitude vectors from the cavities C₁ to C_N.

4.2. Connection through the intermediate barriers

Around each barrier located between the cavities C_i and C_i + 1, a first solution for the wave function is the solution of the base wave reflection-transmission problem presented in Sect. 2.3 and is denoted by Ψ⁽¹⁾. This solution is associated in both cavities with the amplitude vectors up to a given proportionality constant:

$$\begin{array}{c}\hfill {\mathit{a}}_{i+1,-}^{(1)}\propto \left(\begin{array}{c}{T}_{i,i+1}\\ 0\end{array}\right)\mathrm{and}{\mathit{a}}_{i,+}^{(1)}\propto \left(\begin{array}{c}1\\ R_{i,i+1}\end{array}\right).\end{array}$$(77)

As shown in Appendix B, another solution can then be obtained by complex conjugation and time reversal. In other words, Ψ⁽²⁾ = Ψ^(1)⋆(r, −t) is also a solution. According to Eq. (25) (see also the example in Appendix B), this solution is associated with the amplitude vectors, up to the same proportionality constant as in Eq. (77),

$$\begin{array}{c}\hfill {\mathit{a}}_{i+1,-}^{(2)}\propto \left(\begin{array}{c}0\\ T_{i,i+1}^{\star }\end{array}\right)\mathrm{and}{\mathit{a}}_{i,+}^{(2)}\propto \left(\begin{array}{c}{R}_{i,i+1}^{\star }\\ 1\end{array}\right).\end{array}$$(78)

The two solutions Ψ⁽¹⁾ and Ψ⁽²⁾ are linearly independent and constitute a basis for the solution of the wave equation around the considered barrier (indeed, it is straightforward to show that their Wronskian is not null). They are thus sufficient to deduce the general transformation making the connection between a_{i + 1, −} and a_i, + for 1 ≤ i ≤ N − 1 in the following matrix form

$$\begin{array}{c}\hfill {\mathit{a}}_{i,+}={\mathsf{B}}_{i,i+1}{\mathit{a}}_{i+1,-},\end{array}$$(79)

where B_{i, i + 1} is the transformation matrix:

$$\begin{array}{cc}\hfill {\mathsf{B}}_{i,i+1}& =\left(\begin{array}{cc}{T}_{i,i+1}^{-1}& R_{i,i+1}^{\star }{T}_{i,i+1}^{\star -1}\\ R_{i,i+1}{T}_{i,i+1}^{-1}& T_{i,i+1}^{\star -1}\end{array}\right)\hfill \\ \hfill & =\mathsf{C}(-\frac{{\delta }_{i,i+1}}{2}){\mathsf{A}}_{i,i+1}\mathsf{C}(\frac{\pi }{2}-\frac{{\delta }_{i+1,i}}{2}),\hfill \end{array}$$(80)

where we have used Eq. (32) to decompose the matrix and where we have defined

$$\begin{array}{cc}\hfill \mathsf{C}(\epsilon )& =\left(\begin{array}{cc}{e}^{+\mathrm{i}\epsilon }& 0\\ 0& e^{-\mathrm{i}\epsilon }\end{array}\right)\hfill \end{array}$$(81)

$$\begin{array}{cc}\hfill {\mathsf{A}}_{i,i+1}& =\frac{1}{\left|T_{i,i+1}\right|}\left(\begin{array}{cc}1& \left|R_{i,i+1}\right|\\ \left|R_{i,i+1}\right|& 1\end{array}\right).\hfill \end{array}$$(82)

For our purpose, it is useful to specify how the C and A_{i, i + 1} matrices transform a given amplitude vector in the form of:

$$\begin{array}{c}\hfill \mathit{w}(\epsilon )=\left(\begin{array}{c}{e}^{+\mathrm{i}\epsilon }\\ e^{-\mathrm{i}\epsilon }\end{array}\right).\end{array}$$(83)

It is straightforward to find that

$$\begin{array}{cc}& \mathsf{C}({\epsilon }_2)\mathit{w}({\epsilon }_1)=\mathit{w}({\epsilon }_2+{\epsilon }_1)\hfill \end{array}$$(84)

$$\begin{array}{cc}& {\mathsf{A}}_{i,i+1}\mathit{w}(\epsilon )={\mathcal{A}}_{i,i+1}(\epsilon )\mathit{w}(arctan[q_{i,i+1}tan\left(\epsilon \right)]),\hfill \end{array}$$(85)

where, using the convention in Eq. (71) for the branch of the arctangent function³,

$$\begin{array}{c}\hfill {\mathcal{A}}_{i,i+1}(\epsilon )=\sqrt{\frac{q_{i,i+1}^2+(1-q_{i,i+1}^2){cos}^2\left(\epsilon \right)}{q_{i,i+1}}}·\end{array}$$(86)

4.3. Wave function matching in the middle of the cavities

At this point, we need to connect the general solution around the barrier between C_i and C_i + 1 to the general solution around an adjacent barrier, for instance, between C_i − 1 and C_i. It is thus sufficient to make the link between a_i, + and a_i, − inside the cavity C_i for 1 ≤ i ≤ N. This is provided by Eqs. (48) and (49), and can be written in the following matrix form:

$$\begin{array}{c}\hfill {\mathit{a}}_{i,-}={\mathsf{C}}_i{\mathit{a}}_{i,+}\equiv \mathsf{C}(-{\mathrm{\Theta }}_i){\mathit{a}}_{i,+},\end{array}$$(87)

where Θ_i and C are defined in Eqs. (46) and (81). This simple transformation is actually equivalent to match the two JWKB solutions for the wave functions coming from $r_i^{-}$ and $r_i^{+}$, respectively, inside the cavity C_i, as usually done for instance in the usual asymptotic analyses of stellar pulsations (e.g., Shibahashi 1979; Tassoul 1980; Takata 2016a).

4.4. General transformation between adjacent cavities

The two last operations can be composed to describe the general transformation of the amplitude vector from a cavity to an adjacent cavity, that is, the transformation from a_{i + 1, +} to a_i, +. Using Eqs. (79) and (87), we get for 1 ≤ i ≤ N − 1:

$$\begin{array}{c}\hfill {\mathit{a}}_{i,+}={\mathsf{B}}_{i,i+1}{\mathsf{C}}_{i+1}{\mathit{a}}_{i+1,+}\equiv {\mathsf{E}}_{i,i+1}{\mathit{a}}_{i+1,+}.\end{array}$$(88)

It is also useful for the following to deduce the transformation of a vector amplitude w(ε) as defined in Eq. (83) by the matrix E_{i, i + 1}. Using Eqs. (80), (84), (85), and (87), we find for 1 ≤ i ≤ N − 1 that

$$\begin{array}{c}\hfill {\mathsf{E}}_{i,i+1}\mathit{w}(\epsilon )={\mathcal{A}}_{i,i+1}(\epsilon -{\mathrm{{\rm Y}}}_{i+1}+\frac{{\delta }_{i+1,i+2}}{2})\mathit{w}\left[{\mathrm{\Omega }}_{i+1}(\epsilon )\right],\end{array}$$(89)

with

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{i+1}(\epsilon )=arctan[q_{i,i+1}tan(\epsilon -{\mathrm{{\rm Y}}}_{i+1}+\frac{{\delta }_{i+1,i+2}}{2})]-\frac{{\delta }_{i,i+1}}{2},\end{array}$$(90)

where we have used Eqs. (32) and (72) to express the result as a function of the phase Υ_i.

4.5. Surface boundary condition

We first apply the totally reflective boundary condition at the surface of the star, which is represented by Eq. (36). In terms of the amplitude vector in C_N, it reads:

$$\begin{array}{c}\hfill {\mathit{a}}_{N,+}=a_N\mathit{w}(-\frac{{\delta }_{\mathrm{s}}}{2}),\end{array}$$(91)

where a_N is a complex constant. Imposing this surface condition, it is then possible to deduce the amplitude vectors in the underlying cavities by successively applying the linear general transformation E_{i, i + 1} in Eq. (88) from i = N − 1 to i = k. This gives

$$\begin{array}{c}\hfill {\mathit{a}}_{k,+}={\displaystyle \left(\prod _{i=k}^{N-1}{\mathsf{E}}_{i,i+1}\right)a_N\mathit{w}(-\frac{{\delta }_{\mathrm{s}}}{2}).}\end{array}$$(92)

Using Eq. (89) and setting for convenience

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{N+1}(\epsilon )=\epsilon -\frac{{\delta }_{\mathrm{s}}}{2},\end{array}$$(93)

this can be expressed in the simple form

$$\begin{array}{c}\hfill {\mathit{a}}_{k,+}=a_k\mathit{w}[{\mathrm{\Omega }}_{k+1}°\cdots °{\mathrm{\Omega }}_{N+1}(0)],\end{array}$$(94)

where (° ) is the composition operator and a_k is given by

$$\begin{array}{c}\hfill \frac{a_k}{a_N}={\displaystyle \prod _{i=k}^{N-1}{\mathcal{A}}_{i,i+1}[{\mathrm{\Omega }}_{i+2}°\cdots °{\mathrm{\Omega }}_{N+1}(0)-{\mathrm{{\rm Y}}}_{i+1}+\frac{{\delta }_{i+1,i+2}}{2}],}\end{array}$$(95)

We note, according to Eq. (94), that in a given cavity, the amplitudes of the upward and downward components have the same modulus. This is the result of the totally reflective surface boundary condition in Eq. (91) and has already been demonstrated using the infinite-time reflection picture.

4.6. Central boundary and eigenvalue conditions

In the final step, we apply the bottom boundary condition close to the core of the star, which is represented by Eq. (34). In terms of the vector amplitude, this translates into

$$\begin{array}{c}\hfill {\mathit{a}}_{1,-}={\stackrel{\sim }{a}}_1\mathit{w}\left(\frac{{\delta }_{\mathrm{c}}}{2}\right),\end{array}$$(96)

where ${\stackrel{\sim }{a}}_1$ is a complex constant. Another expression of a_1, − can also be obtained but satisfying, this time, the surface boundary condition. Using Eqs. (87) and (94), it is provided by

$$\begin{array}{c}\hfill {\mathit{a}}_{1,-}={\mathsf{C}}_1{\mathit{a}}_{1,+}=a_1\mathit{w}[{\mathrm{\Omega }}_2°···°{\mathrm{\Omega }}_{N+1}\left(0\right)-{\mathrm{\Theta }}_1].\end{array}$$(97)

Equations (96) and (97) have therefore to be met simultaneously to obtain a solution of the boundary value problem. The eigenvalue condition is thus met if the amplitudes are chosen such as ${\stackrel{\sim }{a}}_1={(-1)}^n{a}_1$ and

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_1°\cdots °{\mathrm{\Omega }}_{N+1}(0)=n\pi ,\end{array}$$(98)

where we have defined

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_1(\epsilon )=-\epsilon +{\mathrm{{\rm Y}}}_1-\frac{{\delta }_{1,2}}{2}+\frac{\pi }{2},\end{array}$$(99)

and where n is the radial order. It is straightforward to show that Eq. (98) is equivalent to the recurrence condition found in the infinite-time reflection picture and provided by Eqs. (74) and (75). Both approaches are therefore equivalent. The boundary value problem scenario nevertheless turns out to be more convenient to discuss the mode amplitude, as we show in the next section.

5. Local mean mode energy and amplitudes

In addition to the eigenfrequencies, it is also interesting to know the distribution of the mode energy throughout the ensemble of cavities. This can provide information about the regions that the modes can efficiently probe and this is essential for predicting the surface mode displacement and, thus, the observed oscillation power spectra (e.g., Chaplin et al. 2005). In this section, we briefly address this point based on the present basic formulation of eigenmodes.

5.1. Mean mode energy in each cavity

Owing to the equipartition of the potential and kinetic energy of short-wavelength gravito-acoustic waves (e.g., Lighthill 1978), the mean mode energy averaged over one oscillation period T = 2π/σ in the cavity C_k is defined as:

$$\begin{array}{c}\hfill {\mathcal{E}}_k=\frac{1}{T}{\displaystyle {\int }_{-T/2}^{+T/2}({\int }_{r_k^{-}}^{r_k^{+}}\int {\int }_{\mathrm{\Sigma }}\rho {\left|{\mathcal{R}}_{\mathrm{e}}[\mathrm{i}\sigma \mathit{\xi }(\mathit{r},t)]\right|}^2{r}^2{\mathrm{d}}^2\mathrm{\Sigma }\mathrm{d}r)\mathrm{d}t,}\end{array}$$(100)

where d²Σ = sin θdθds is the solid angle in the direction (θ, s). First, to express Eq. (100), we note that the horizontal mode displacement can be related to the Eulerian pressure perturbation through the momentum conservation in the horizontal direction and the equality (e.g., Unno et al. 1989):

$$\begin{array}{c}\hfill {\stackrel{\sim }{\xi }}_{\mathrm{h}}=\frac{{\stackrel{\sim }{p}}^{\prime }}{\rho r{\sigma }^2}·\end{array}$$(101)

Then, using Eqs. (2), (15), (17), (18), and (101), and taking advantage of the properties of the spherical harmonics, we find at leading order for ${\sigma }^2\gg (S_{\mathcal{l}}^2\text{and}{N}^2)$ or ${\sigma }^2\ll (S_{\mathcal{l}}^2\text{and}{N}^2)$ that (see Appendix D for details)

$$\begin{array}{c}\hfill {\mathcal{E}}_k\approx |a_k{|}^2|{\mathrm{\Theta }}_k|.\end{array}$$(102)

The local mean mode energy is therefore first directly proportional to the wave number integral in the cavity C_k defined in Eq. (46); we can approximately write Θ_k ≈ πn_k, where n_k is an integer whose modulus represents the number of oscillation nodes in the radial direction in the cavity C_k. Second, it is proportional to the squared modulus of the mode amplitude. Knowing the actual value of this amplitude would require studies of the excitation and damping of the modes, which is beyond the scope of this work. The previous analysis can nevertheless provide us with the amplitude ratios between adjacent cavities, as we show in the next section.

5.2. Amplitude ratios

Using Eqs. (86) and (95), the squared amplitude in the cavity C_k relatively to that in the surface cavity C_N is equal to

$$\begin{array}{c}\hfill {\left|\frac{a_k}{a_N}\right|}^2=\prod _{i=k}^{N-1}\frac{q_{i,i+1}^2+(1-q_{i,i+1}^2){cos}^2\left({\mathrm{\Xi }}_{i+1}\right)}{q_{i,i+1}},\end{array}$$(103)

where we define the quantity

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_{i+1}={\mathrm{\Omega }}_{i+2}°···°{\mathrm{\Omega }}_{N+1}(0)-{\mathrm{{\rm Y}}}_{i+1}+\frac{{\delta }_{i+1,i+2}}{2}·\end{array}$$(104)

The squared amplitude ratio between the adjacent cavities C_k and C_k + 1 is thus merely equal to

$$\begin{array}{c}\hfill {\left|\frac{a_k}{a_{k+1}}\right|}^2=\frac{q_{k,k+1}^2+(1-q_{k,k+1}^2){cos}^2\left({\mathrm{\Xi }}_{k+1}\right)}{q_{k,k+1}}·\end{array}$$(105)

We can thus see that the squared amplitude ratio of the wave function between adjacent cavities depends on two ingredients: first, the coupling factor associated with the intermediate barrier, which measures how the energy is transmitted from one cavity to the other cavity; second, the “optical” path during a back and forth travel of a wave inside the ensemble of cavities above the intermediate barrier, which is represented by Ξ_k + 1 and measures the level of constructive interference in the cavity C_k + 1. In the case when q_{k, k + 1} is close to unity, Eq. (105) show that |a_k/a_k + 1| is on the order of unity. Indeed, in this case, both cavities are well coupled and exchange comparable wave energy fluxes. In contrast, when q_{k, k + 1} is much smaller than unity, the result depends on the level of constructive interference in the cavity C_k + 1. In the limiting case where the level of constructive interference is maximum, which is equivalent to Ξ_k + 1 = (n_k + 1 + 1/2)π with n_k + 1 a given integer, Eq. (105) shows that the squared amplitude of the wave function is smaller in the cavity C_k than in the cavity C_k + 1 by a factor of q_{k, k + 1} ≪ 1. Conversely, when the level of constructive interference is minimum, that is Ξ_k + 1 = n_k + 1π, the squared amplitude of the wave function is larger in the cavity C_k than in the cavity C_k + 1 by a factor of 1/q_{k, k + 1} ≫ 1.

Finally, we first recall that the previous discussion address the amplitude ratio of the wave function in both cavities and that, in order to translate these results in terms of mode energy, we also have to take the Θ_k factor in Eq. (102) into account, which represents the local mode inertia in the asymptotic limit. Second, it is worth mentioning that Eq. (105) expresses a_k as a function of the mode amplitude in the upper cavity C_N since it was obtained by imposing only the surface boundary condition; it is thus valid under this sole condition, which is met by definition for eigenmodes (i.e., in the resonance condition). Imposing instead the sole core boundary condition, it is also possible to derive an expression of a_k as a function of the mode amplitude a₁ in the inner cavity C₁, as shown in Appendix E. Such an expression is equivalent to Eq. (105) only and only for eigenmodes.

6. Simple cases

In this section, we apply the present formulation on simple cases and check the compatibility with the eigenfrequency conditions already obtained in previous works.

6.1. Gravity and acoustic modes with one single cavity

First, we consider the case of pure gravity or acoustic modes propagating in one single cavity. The cavity is supposed to be located between two turning-points close to the core and surface and beyond which the modes are evanescent, in agreement with the totally reflective boundary conditions. Using Eq. (75), or equivalently Eq. (98) for N = 1, we generally find that:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{r_1^{-}}^{r_1^{+}}{\mathcal{K}}_{\mathrm{r}}\mathrm{d}r\pm (\frac{{\delta }_{\mathrm{c}}}{2}+\frac{{\delta }_{\mathrm{s}}}{2})=\pm n\pi ,}\end{array}$$(106)

where the plus and minus signs correspond to the case of acoustic and gravity modes, respectively, and n is a positive integer representing the mode radial order. In general, the values of the phase lags δ_s and δ_c can be computed in a second step using a single turning point asymptotic analysis, such as that performed by Shibahashi (1979). Retaining the Cowling approximation and neglecting the gradients in the equilibrium structure in the wave equation everywhere inside the star, Shibahashi (1979) showed that the wave function, Ψ, takes the form of an Airy function of the first kind in the vicinity of a given single turning point. Choosing the origins of the wave phase in Eq. (11) as equal to the turning points, that is, such as ${\mathcal{K}}_{\text{r}}(r_1^{-})={\mathcal{K}}_{\text{r}}(r_1^{+})=0$, we show in Appendix F that δ_c + δ_s = π if the core and surface turning points have the same nature (i.e., both turning points satisfy either σ² = N² or ${\sigma }^2=S_{\mathcal{l}}^2$); and then δ_c + δ_s = 0 otherwise.

As an illustration, we can consider the examples of the low-frequency gravity modes and high-frequency acoustic modes in low-mass main sequence stars. For the former, the core and surface turning points (i.e., $r_1^{-}$ and $r_1^{+}$ in the chosen convention) are such that ${\sigma }^2=N^2(r_1^{-})=N^2(r_1^{+})$ (e.g., Appourchaux et al. 2010). In this case, we thus have:

$$\begin{array}{c}\hfill {\delta }_{\mathrm{c}}+{\delta }_{\mathrm{s}}=\pi .\end{array}$$(107)

For the latter, the core and surface turning points are such that ${\sigma }^2=S_{\mathcal{l}}^2(r_1^{-})$ and ${\sigma }^2=N^2(r_1^{+})$, so that we get:

$$\begin{array}{c}\hfill {\delta }_{\mathrm{c}}+{\delta }_{\mathrm{s}}=0.\end{array}$$(108)

We note that we retrieve the same result as Shibahashi (1979) for high-frequency acoustic modes, but that we find an additional phase lag of π/2 in the quantization condition of low-frequency gravity modes. This results from the fact that Shibahashi (1979) considered in contrast that the upper turning point is such that ${\sigma }^2=S_{\mathcal{l}}^2(r_1^{+})$.

When accounting for the perturbation of the gravitational potential, that is in the non-Cowling case, Takata (2005, 2006) demonstrated for the ℓ = 1 modes that the nature of the core turning point changes. For high-frequency acoustic modes (respectively, low-frequency gravity modes), the situation is actually similar to assume simultaneously the Cowling approximation and ${\sigma }^2=N^2(r_1^{-})$ (respectively, ${\sigma }^2=S_{\mathcal{l}}^2(r_1^{-})$); as a consequence, δ_c + δ_s = 0 (respectively, δ_c + δ_s = π), that is, with a phase shift of π compared to within the Cowling approximation (Takata 2016a; Pinçon et al. 2019). This point emphasizes the importance to take the perturbation of the gravitational potential into account for dipolar modes to study the phase lag at reflection and transmission inside the deep interior of stars.

6.2. Two-cavity mixed modes

In a second example, we consider the case of mixed modes, which can exist not only in red giant stars where they have already been detected, but also in main sequence stars or the Sun. Mixed modes can propagate through an inner buoyancy cavity C₁, where they behave as gravity modes, and an outer pressure cavity C₂, where they behave as acoustic modes. Using Eqs. (74) and (75), or, equivalently, Eq. (98) for N = 2, we can easily get the general resonance condition that states:

$$\begin{array}{c}\hfill cot(-{\mathrm{\Theta }}_1-\frac{{\delta }_{\mathrm{c}}}{2}-\frac{{\delta }_{1,2}}{2}+\frac{\pi }{2})tan({\mathrm{\Theta }}_2+\frac{{\delta }_{\mathrm{s}}}{2}+\frac{{\delta }_{2,1}}{2})=q_{1,2}.\end{array}$$(109)

Equation (109) can be shown to be similar to Eq. (18) of Takata (2016b) using the relation δ_2, 1 = π − δ_1, 2 + 2γ_1, 2, with γ_1, 2 = 0.

At this point, we can also discuss the value of the phase lags in the case of an evolved red giant star within the Cowling approximation, following again the work by Shibahashi (1979). In such considerations, the cavity C₁ is located between two turning points such as ${\sigma }^2=N^2(r_1^{-})=N^2(r_1^{+})$, and the cavity C₂ is located between a lower and upper turning points such as ${\sigma }^2=S_{\mathcal{l}}^2(r_2^{-})$ and ${\sigma }^2=N^2(r_2^{+})$, respectively (e.g., Hekker & Christensen-Dalsgaard 2017). The boundaries of the cavities are chosen equal to the turning points. The main assumption of the analysis of Shibahashi (1979) then consists in considering that the four turning points are far away from each other; this permits us to exploit the results obtained from a single turning point analysis of the stellar oscillation equations around each of them. Within this context, neglecting the gradients of the equilibrium structure in the wave equation everywhere, we can deduce δ_s = −δ_c = −π/2 according to Appendix F and Eq. (F.4). Moreover, as the intermediate evanescent region is thick (i.e., weak coupling between C₁ and C₂), the phase lags associated with the intermediate evanescent barrier δ_1, 2 and δ_2, 1 can also be shown to follow the same rules as in Appendix F (see, e.g., Takata 2016a; Pinçon et al. 2019). We therefore deduce δ_1, 2 ≈ δ_2, 1 ≈ π/2 and γ_1, 2 = 0. In other words, in the Cowling and weak coupling paradigm, Eq. (109) is reduced to cot(−Θ₁)tanΘ₂ = q_1, 2, which is similar to Eq. (31) of Shibahashi (1979), since we recall that Θ₁ is minus the wave number integral in the cavity C₁. To be precise, we also first note that in the non-Cowling case, we must apply a shift of π on δ_c, so that δ_c = −π/2. Moreover, the more complex expressions of δ_1, 2 and δ_2, 1 in the approximation of a very thin evanescent region (i.e., strong coupling hypothesis) are also available (Takata 2016a).

6.3. Three-cavity mixed modes

More recently, Cunha et al. (2015), Deheuvels & Belkacem (2018) tackled the case of mixed modes propagating in three resonant cavities. Using Eqs. (74) and (75), or equivalently Eq. (98) for N = 3, we find for three cavities C₁, C₂, and C₃ that the resonance condition reads in a general way

$$\begin{array}{c}\hfill tan{\mathrm{{\rm Y}}}_3(1-q_{1,2}tan{\mathrm{{\rm Y}}}_{2}tan{\mathrm{{\rm Y}}}_1)+q_{2,3}tan{\mathrm{{\rm Y}}}_2+q_{1,2}{q}_{2,3}tan{\mathrm{{\rm Y}}}_1=0.\end{array}$$(110)

When q_2, 3 = 0, we note that the term in the brackets is equal to zero and we retrieve the resonance condition in Eq. (109) for two-cavity eigenmodes with the substitution δ_2, 3 ← δ_s.

As an additional check, it is interesting to compare in more details this result with the asymptotic analysis of mixed modes by Deheuvels & Belkacem (2018) in the case of helium-core flash red giant stars. In contrast with less evolved red giant stars, the presence of a convective region at the border of the helium core locally creates an evanescent region that splits the radiative core into two cavities. where the modes behave as gravity modes (see, e.g., Fig. 2 of Deheuvels & Belkacem 2018). Assuming the Cowling approximation, the turning points of the two inner cavities are such as σ² = N² while the configuration in the upper cavity is similar to the case of acoustic modes studied in Sect. 6.2. Choosing the turning points as the origin of the wave phase and using the low-coupling hypothesis in both evanescent regions, the same arguments as in Sect. 6.2 thus hold true and we can write

$$\begin{array}{c}\hfill {\delta }_{\mathrm{c}}+{\delta }_{1,2}={\delta }_{2,1}+{\delta }_{2,3}=\pi \mathrm{and}{\delta }_{3,2}+{\delta }_{\mathrm{s}}=0.\end{array}$$(111)

Therefore, using Eqs. (72) and (73) with Eqs. (110) and (111) can be expressed as

$$\begin{array}{cc}\hfill cot& (-{\mathrm{\Theta }}_1)cot(-{\mathrm{\Theta }}_2)tan{\mathrm{\Theta }}_3-q_{1,2}tan{\mathrm{\Theta }}_3-q_{2,3}cot(-{\mathrm{\Theta }}_1)\hfill \\ \hfill & -q_{1,2}{q}_{2,3}cot(-{\mathrm{\Theta }}_2)=0.\hfill \end{array}$$(112)

We retrieve here, based on basic arguments, the same expression as in Eqs. (9) and (10) of Deheuvels & Belkacem (2018), keeping in mind the fact that the minus sign in front of Θ₁ and Θ₂ comes from the definition of Θ_i in Eq. (46) for gravity-dominated modes. We finally note that Cunha et al. (2015, 2019) formulated the resonance condition in an equivalent way to Eq. (110), but, in addition, these authors expressed the q_1, 2 factor by assuming that the barrier corresponds to either a sharp Dirac or Gaussian peak in the Brunt–Väisälä frequency.

6.4. Low-amplitude glitches

As a final illustration, we consider the important case of low-amplitude glitches, that is, low-amplitude sharp and very localized features in the equilibrium structure. In most previous studies, the effect of such sharp gradients on the eigenfrequencies was usually treated as a small perturbation of the “smooth” case, that is, the case where these local sharp features are not taken into account. In order to make a comparison with these previous results, we thus need to express the eigenfrequency deviations induced by low-amplitude glitches that is predicted by the present framework, while additionally using the small perturbations hypothesis. This is actually equivalent to assuming that (a) the barriers are very thin and localized at a frequency-independent radius, namely, $r_{i+1}^{-}(\sigma )=r_i^{+}(\sigma )=r_{\text{g},i}$; and (b) the barriers are weakly reflective, that is, |R_{i, i + 1}|≪1. The goal is then to find the frequency perturbation σ₁ induced by the barriers, that is,

$$\begin{array}{c}\hfill {\sigma }_1=\sigma -{\sigma }_0,\end{array}$$(113)

where σ is the actual mode angular eigenfrequency and σ₀ indicates its value in the smooth case without any barrier.

6.4.1. One single glitch

We first consider the most familiar case of one single glitch located between two resonant cavities C₁ and C₂. Within the perturbation framework, the phase terms in Eq. (72) can be rewritten in any cavity C_i as

$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}_i(\sigma )={\mathrm{{\rm Y}}}_i^{(0)}(\sigma )+{\mathrm{{\rm Y}}}_i^{(1)}(\sigma ),\end{array}$$(114)

where the superscript (1) indicates in the following a perturbation induced by the barriers whereas the superscript (0) indicates a value (or any structural function) in the smooth case when the barrier is not taken into account. Neglecting the perturbation of the wave number integral function because of the hypothesis (a) and making use of the fact that the core and surface phase lags are not perturbed by the barriers, we can deduce for the two-cavity case that

$$\begin{array}{cc}& {\mathrm{{\rm Y}}}_1^{(0)}(\sigma )={\mathrm{\Theta }}_1^{(0)}(\sigma )+\frac{{\delta }_{\mathrm{c}}^{(0)}(\sigma )}{2}\hfill \end{array}$$(115)

$$\begin{array}{cc}& {\mathrm{{\rm Y}}}_2^{(0)}(\sigma )={\mathrm{\Theta }}_2^{(0)}(\sigma )+\frac{{\delta }_{\mathrm{s}}^{(0)}(\sigma )}{2},\hfill \end{array}$$(116)

and

$$\begin{array}{cc}& {\mathrm{{\rm Y}}}_1^{(1)}(\sigma )=\frac{{\delta }_{1,2}^{(1)}(\sigma )}{2}-\frac{\pi }{2}\hfill \end{array}$$(117)

$$\begin{array}{cc}& {\mathrm{{\rm Y}}}_2^{(1)}(\sigma )=\frac{{\delta }_{2,1}^{(1)}(\sigma )}{2},\hfill \end{array}$$(118)

which is consistent with the definitions in Eqs. (72) and (73). First, in the “smooth” configuration without any barrier, the frequency condition is similar to Eq. (106) in the case of one single cavity located between the center and surface of the star, that is:

$$\begin{array}{c}\hfill {\overline{\mathrm{{\rm Y}}}}_2^{(0)}({\sigma }_0)=n\pi ,\end{array}$$(119)

where n is the radial order and where the overline notation for any integers i and k is defined as

$$\begin{array}{c}\hfill {\overline{\mathrm{{\rm Y}}}}_i^{(k)}={\displaystyle \sum _{j=1}^i{\mathrm{{\rm Y}}}_j^{(k)}.}\end{array}$$(120)

Second, adding the intermediate glitch, we choose to write the general resonance condition in Eqs. (74) and (75), or equivalently in Eq. (98), as a series of sine terms as in Eq. (C.12), that is, for N = 2:

$$\begin{array}{c}\hfill sin[{\mathrm{{\rm Y}}}_1(\sigma )+{\mathrm{{\rm Y}}}_2(\sigma )]=-\left|R_{1,2}(\sigma )\right|sin[{\mathrm{{\rm Y}}}_2(\sigma )-{\mathrm{{\rm Y}}}_1(\sigma )].\end{array}$$(121)

We emphasize that this resonance condition is also valid to describe mixed modes, except that in the present case, the modes have the same nature (i.e., gravity or acoustic) on both sides of the intermediate barrier (and thus the same sign in front of the wave number integral in both cavities). Within the small perturbation limit, we assume that the frequency perturbation σ₁ is small enough that the perturbations of the reflection coefficient and of the total “optical” path of the waves during one travel through the ensemble of cavities is also small compared to the smooth case, that is,

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}ln|R_{1,2}|}{\mathrm{d}\sigma }\right)}_{{\sigma }_0}{\sigma }_1\ll 1,{\overline{\mathrm{{\rm Y}}}}_2^{(1)}(\sigma )\ll \pi ,{\left(\frac{\mathrm{d}{\mathrm{{\rm Y}}}_i^{(k)}}{\mathrm{d}\sigma }\right)}_{{\sigma }_0}{\sigma }_1\ll \pi .\end{array}$$(122)

Using Eqs. (114) and (119), a first-order Taylor expansion of Eq. (121) in |R_1, 2| by hypothesis (b) and in the quantities in Eq. (122) leads to

$$\begin{array}{c}\hfill \frac{\pi {\sigma }_1}{\mathrm{\Delta }\sigma }+{\overline{\mathrm{{\rm Y}}}}_2^{(1)}({\sigma }_0)\approx \left|R_{1,2}({\sigma }_0)\right|sin\left(2[{\mathrm{{\rm Y}}}_1^{(0)}({\sigma }_0)+{\mathrm{{\rm Y}}}_1^{(1)}({\sigma }_0)]\right),\end{array}$$(123)

where we have defined the “pseudo” pulsation large separation in the general case of N cavities as

$$\begin{array}{c}\hfill \mathrm{\Delta }\sigma =\pi {\left(\frac{\mathrm{d}{\overline{\mathrm{{\rm Y}}}}_N^{(0)}}{\mathrm{d}\sigma }\right)}_{{\sigma }_0}^{-1}.\end{array}$$(124)

Using Eqs. (32), (117), and (118), Eq. (123) can be rewritten as

$$\begin{array}{c}\hfill \frac{{\sigma }_1}{\mathrm{\Delta }\sigma }\approx -\frac{\left|R_{1,2}({\sigma }_0)\right|}{\pi }sin(2{\mathrm{{\rm Y}}}_1^{(0)}({\sigma }_0)+{\delta }_{1,2}({\sigma }_0))-\frac{{\gamma }_{1,2}^{(1)}({\sigma }_0)}{\pi }·\end{array}$$(125)

Moreover, using Eqs. (46) and (115), we obtain:

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_1^{(0)}({\sigma }_0)& =\pm {\displaystyle {\int }_{r_1^{-}\left({\sigma }_0\right)}^{r_{\mathrm{g},1}}{\mathcal{K}}_{\mathrm{r}}\left({\sigma }_0\right)\mathrm{d}r+\frac{{\delta }_{\mathrm{c}}^{(0)}\left({\sigma }_0\right)}{2}·}\hfill \end{array}$$(126)

The frequency deviation induced by the low-amplitude glitch from the smooth case therefore takes the form of a small offset term plus a sinusoidal term. The magnitude of the sinusoidal perturbation is proportional to the reflection coefficient associated with the glitch. Moreover, its argument depends on the wave number integral from the core to the radius of the glitch, and is related to the so-called acoustic radius. The phase offset results from the phase lag introduced at reflection on the glitch. We note that it is also possible to formulate Eq. (125) as a function of ${{\rm Y}}_2^{(0)}({\sigma }_0)$ instead of ${{\rm Y}}_1^{(0)}({\sigma }_0)$ using Eq. (119), that is as a function of the wave number integral from the glitch radius to the surface, which is related to the so-called acoustic depth.

The description of the first-order effect of low-amplitude glitches on the eigenfrequencies that is provided by Eq. (125) finally appears to be consistent with the computation of Pinçon (2019). It is also compatible with previous works based on the variational principle (see Sect. 1 for references); in addition, it has the advantage of clarifying the general physical meaning of the different terms inside the glitch-induced deviation (e.g., amplitude, offsets). While the approach based on the variational principle is limited to the first-order small-amplitude glitches only, we stress that the formulation developed in this paper is more general and can also be used to treat big glitches (e.g., mode-trapping phenomena).

6.4.2. Multiple glitches

In a second step, it is straightforward to extend the last first-order expansion to the case of N − 1 low-amplitude glitches. Using the hypothesis (a) and the non-perturbation of the core and surface phase lags, the phase terms in Eq. (114) can be expressed for 1 ≤ i ≤ N within the perturbation framework as:

$$\begin{array}{cc}& {\mathrm{{\rm Y}}}_i^{(0)}(\sigma )={\mathrm{\Theta }}_i^{(0)}(\sigma )+{\delta }_{\mathrm{K}}(i-1)\frac{{\delta }_{\mathrm{c}}^{(0)}(\sigma )}{2}+{\delta }_{\mathrm{K}}(i-N)\frac{{\delta }_{\mathrm{s}}^{(0)}(\sigma )}{2},\hfill \end{array}$$(127)

$$\begin{array}{cc}& {\mathrm{{\rm Y}}}_i^{(1)}(\sigma )=\frac{{\delta }_{i,i+1}^{(1)}(\sigma )}{2}+\frac{{\delta }_{i,i-1}^{(1)}(\sigma )}{2}-\frac{\pi }{2},\hfill \end{array}$$(128)

where δ_K(i) is the Kronecker delta function and we set

$$\begin{array}{c}\hfill {\delta }_{1,0}^{(1)}=0\mathrm{and}{\delta }_{N,N+1}^{(1)}=\pi \end{array}$$(129)

in order to be consistent with the definitions in Eqs. (72) and (73). First, in the “smooth” configuration without any barrier, the frequency condition is expected to be similar to Eq. (119) but for N cavities, that is,

$$\begin{array}{c}\hfill {\overline{\mathrm{{\rm Y}}}}_N^{(0)}({\sigma }_0)=n\pi .\end{array}$$(130)

Second, we add the N − 1 low-amplitude glitches. Within the perturbation limit, we assume as previously that the frequency perturbation is small enough that:

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}ln|R_{i,i+1}|}{\mathrm{d}\sigma }\right)}_{{\sigma }_0}{\sigma }_1\ll 1,{\overline{\mathrm{{\rm Y}}}}_N^{(1)}(\sigma )\ll \pi ,{\left(\frac{\mathrm{d}{\mathrm{{\rm Y}}}_i^{(k)}}{\mathrm{d}\sigma }\right)}_{{\sigma }_0}{\sigma }_1\ll \pi .\end{array}$$(131)

Using the resonance condition in Eq. (C.12) for N cavities, a first-order Taylor expansion as a function of |R_{i, i + 1}| by hypothesis (a) and of the small quantities in Eq. (131) leads to:

$$\begin{array}{c}\hfill \frac{{\sigma }_1}{\mathrm{\Delta }\sigma }\approx {\displaystyle \sum _{i=1}^{N-1}\frac{\left|R_{i,i+1}({\sigma }_0)\right|}{\pi }sin(2{\overline{\mathrm{{\rm Y}}}}_i^{(0)}({\sigma }_0)+2{\overline{\mathrm{{\rm Y}}}}_i^{(1)}({\sigma }_0))-\frac{{\overline{\mathrm{{\rm Y}}}}_N^{(1)}({\sigma }_0)}{\pi }·}\end{array}$$(132)

Moreover, using Eqs. (46) and (127), and hypothesis (a) on the one hand, and Eqs. (32), (128), and (129) on the other hand, we can deduce for i ≤ N − 1 that:

$$\begin{array}{cc}& {\overline{\mathrm{{\rm Y}}}}_i^{(0)}({\sigma }_0)=\pm {\displaystyle {\int }_{r_1^{-}({\sigma }_0)}^{r_{\mathrm{g},i}({\sigma }_0)}{\mathcal{K}}_{\mathrm{r}}({\sigma }_0)\mathrm{d}r+\frac{{\delta }_{\mathrm{c}}^{(0)}({\sigma }_0)}{2}}\hfill \end{array}$$(133)

$$\begin{array}{cc}& {\overline{\mathrm{{\rm Y}}}}_i^{(1)}({\sigma }_0)={\displaystyle \sum _{j=1}^{i-1}{\gamma }_{j,j+1}^{(1)}({\sigma }_0)+\frac{{\delta }_{i,i+1}^{(1)}({\sigma }_0)}{2}-\frac{\pi }{2}·}\hfill \end{array}$$(134)

The total frequency deviation from the smooth case therefore reduces to the superposition of the perturbations resulting from each glitch, which are related to the acoustic radius and the phase lags at transmission and reflection.

7. Concluding remarks and discussion

In this work, we derive a general analytical expression for the asymptotic resonance condition of global oscillation modes in spherical stars. While this can seem complex at first glance, we show that this is, in fact, merely analogous to regard a star as a one-dimensional giant Fabry-Pérot interferometer, composed of a multitude of resonant cavities and impermeable surface and core boundaries. In the adopted view, a star is decomposed into an ensemble of resonant cavities where waves can propagate within the short-wavelength JWKB approximation. The cavities are separated by intermediate barriers corresponding to evanescent regions or glitches. Each barrier is associated with a reflection and transmission coefficient. The core and surface boundary conditions are represented by totally reflective barriers. Within this framework, we obtained the resonance condition while considering two different physical pictures. In the infinite-time reflection picture, we follow the back-and-forth travel of a wave energy ray through the ensemble cavities and assume constructive interferences; the resonance condition is provided by Eqs. (74) and (75). In the second picture, the eigenmodes are considered as the solution of a linear boundary value problem and the resonance condition is provided by Eq. (98). Both pictures are equivalent and predict the same resonance condition, which turns out to depend on a number of parameters: the wave number integral over each cavity; the coupling factor associated with each barrier between adjacent cavities; the phase lags at reflection and transmission through the stars. In addition, the amplitude ratio between adjacent cavities is also expressed analytically as a function of these parameters in Sect. 5 and can inform us on the distribution of the mode energy through the star. The present formulation enables us to retrieve in a convenient way simple cases already widely studied in the past, as one single cavity acoustic or gravity mode or two- and three-cavity mixed modes as well as the case of multiple low-amplitude glitches.

This general formulation provides a physically grounded interpretation of the mode resonance condition in stars. It also provides a useful tool for analyzing, interpreting, or predicting the oscillation power spectra in a practical way and across a broad variety of configurations and evolutionary stages. This new diagnosis tool is for instance expected to bring a valuable help for the development of automated seismic analysis pipelines, which are needed in the current era of space missions and especially for the future PLATO mission (e.g., Rauer et al. 2014). The next step in the investigation will consist in expressing the wave reflection and transmission coefficients around the barriers for specific cases using either toy models or more realistic stellar structure models. By means of either analytical or numerical methods to solve the wave equation, such a study can allow for a link to be forged between the mode parameters, the mode frequencies, and the stellar interior properties. This will be further developed in future works and practical applications on observed spectra will be undertaken.

Finally, we also note that the present general formulation can be applied to study the transmission of progressive waves through layered media in stellar interiors or planetary atmospheres since it offers the possibility to easily retain and remove the boundary conditions according to the configuration. This can be useful in addressing the problem of the transport of angular momentum, heat, or chemical elements by waves in these objects and the present results provide a general formalism for tackling such issues (e.g., André et al. 2019; Cai et al. 2021).

---

Acknowledgments

During this work, C.P. was funded by a postdoctoral fellowship of Chargé de Recherche from F.R.S.-FNRS (Belgium). C.P. also acknowledges financial support from Sorbonne Université (EMERGENCE project). M.T. is partially supported by JSPS KAKENHI Grant Number 18K03695.

References

Appendix A: Gravito-acoustic wave equations in the Cowling approximation

In this section, we briefly recall the form of the second-order equations for gravito-acoustic waves as formulated by Shibahashi (1979) under the Cowling approximation.

A first version is provided by Eqs. (3)–(5) where the M function is equal to

$$\begin{array}{c}\hfill M(r)=f[P(r)],\end{array}$$(A.1)

where

$$\begin{array}{cc}& f(y)=H_p^2{\left|y\right|}^{1/2}\frac{{\mathrm{d}}^2{\left|y\right|}^{-1/2}}{\mathrm{d}{r}^2}=-\frac{H_p^2}{2}\frac{{\mathrm{d}}^{2}ln\left|y\right|}{\mathrm{d}{r}^2}+\frac{H_p^2}{4}{\left(\frac{\mathrm{d}ln\left|y\right|}{\mathrm{d}r}\right)}^2\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& P(r)=\frac{r^{2}h(r)}{c^2}(\frac{S_{\ell }^2}{{\sigma }^2}-1)\hfill \end{array}$$(A.3)

$$\begin{array}{cc}& h(r)=exp{\displaystyle \left[{\int }_0^r(\frac{N^2}{g}-\frac{g}{c^2})\mathrm{d}r\right].}\hfill \end{array}$$(A.4)

We also have the following relations linking the radial displacement with the perturbation of pressure denoted by $\stackrel{\sim }{p}\prime (r)$ in Eq. (1) (e.g., Unno et al. 1989),

$$\begin{array}{cc}\hfill \stackrel{\sim }{\mathrm{\Psi }}& ={\rho }^{-1/2}r{|N^2-{\sigma }^2|}^{-1/2}{\stackrel{\sim }{p}}^{\prime }{e}^{-\mathrm{i}\sigma t}\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill \stackrel{\sim }{\mathrm{\Psi }}& =\mathrm{sgn}(P)\frac{1}{|{\mathcal{K}}_r|}(\frac{\mathrm{d}\mathrm{\Psi }}{\mathrm{d}r}+\frac{1}{2}\frac{\mathrm{d}ln|P|}{\mathrm{d}r}\mathrm{\Psi }),\hfill \end{array}$$(A.6)

where sgn() is the sign function. We note that the M function is on the order of unity at most over regions that are far away from the turning points where P = 0 (i.e., where ${\sigma }^2=S_{\mathcal{l}}^2$) or from the sharp gradients in the structure.

In the vicinity of a turning point where P = 0, the first version of the wave equation in Eq. (3) is singular. Around such a turning point, it is therefore more convenient to use the second version of the wave equation, which is provided by

$$\begin{array}{cc}& \frac{{\mathrm{d}}^2\stackrel{\sim }{\mathrm{\Psi }}}{\mathrm{d}{r}^2}+[{\mathcal{K}}_r^2-\frac{L(r)}{H_p^2}]\stackrel{\sim }{\mathrm{\Psi }}=0,\hfill \end{array}$$(A.7)

where

$$\begin{array}{cc}& L(r)=f[Q(r)]\hfill \end{array}$$(A.8)

$$\begin{array}{cc}& Q(r)=\frac{1}{r^{2}h(r)}({\sigma }^2-N^2).\hfill \end{array}$$(A.9)

The link between the perturbation of pressure and the radial displacement can in this case be obtained through the equation

$$\begin{array}{c}\hfill \mathrm{\Psi }=\mathrm{sgn}(Q)\frac{1}{|{\mathcal{K}}_r|}(\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Psi }}}{\mathrm{d}r}+\frac{1}{2}\frac{\mathrm{d}ln|Q|}{\mathrm{d}r}\stackrel{\sim }{\mathrm{\Psi }}).\end{array}$$(A.10)

The L function is on the order of unity at most over regions that are far away from the turning points where Q = 0 (i.e., where σ² = N²) or from the sharp gradients in the structure. Moreover, at turning points such as Q = 0, this second version of the wave equation is conversely singular, and it is therefore more convenient to consider the first version of the wave equation in Eq. (3).

Appendix B: Relation between the base and adjoint wave transmission-reflection coefficients

As shown by Takata (2016b), the reflection and transmission coefficients in a base wave transmission-reflection problem are related to those of the adjoint wave transmission-reflection problem (see Sect. 2.3 for a description of the two configurations). In this section, we recall the result obtained in Sect. 2 of Takata (2016b) and generalize it to the case of a complex transmission coefficient (i.e., accounting for a transmission phase lag).

In the base wave reflection-transmission problem, the global wave function Ψ(r, t) is provided in the cavity C_n, within the convention considered in Sect. 2.3 and up to a proportionality constant, by

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_n(r,t;r_n^{+})={\psi }^{(\to )}(r,t;r_n^{+})+R_{n,n+1}{\psi }^{(\leftarrow )}(r,t;r_n^{+}),\end{array}$$(B.1)

while in the overlying cavity C_n + 1, it is written

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{n+1}(r,t;r_{n+1}^{-})=T_{n,n+1}{\psi }^{(\to )}(r,t;r_{n+1}^{-}).\end{array}$$(B.2)

The goal is then to build a solution of the adjoint problem from the two representations of the global wave function Ψ(r, t) in Eqs. (B.1) and (B.2). To do so, we can use the fact that the oscillation equations are invariant by the operations of (a) time reversal and (b) complex conjugation. As a justification, (a) the wave equation introduced in Eqs. (3)–(5) and Appendix A is independent of the sign of σ, and (b) its coefficients are real. Therefore, if Ψ(r, t) is a solution, then Ψ′(r, t) = Ψ^⋆(r, −t) is also a solution. According to Eq. (25), applying the operations of time reversal and complex conjugation on Eqs. (B.1) and (B.2) provides the form of the new solution Ψ′(r, t), which can be respectively expressed in the cavities C_n and C_n + 1 as

$$\begin{array}{cc}& {\mathrm{\Psi }}_n^{\prime }(r,t;r_n^{+})=R_{n,n+1}^{\star }{\psi }^{(\to )}(r,t;r_n^{+})+{\psi }^{(\leftarrow )}(r,t;r_n^{+})\hfill \end{array}$$(B.3)

$$\begin{array}{cc}& {\mathrm{\Psi }}_{n+1}^{\prime }(r,t;r_{n+1}^{-})=T_{n,n+1}^{\star }{\psi }^{(\leftarrow )}(r,t;r_{n+1}^{-}).\hfill \end{array}$$(B.4)

Another solution can then be built from a linear combination of the solution Ψ represented by Eqs. (B.1)–(B.2) and the solution Ψ′ represented by Eqs. (B.3)–(B.4). This must be done in such a way that the upward component in the cavity C_n vanishes in order to retrieve the adjoint problem. Such a solution Ψ″ can be obtained by the following linear combination

$$\begin{array}{c}\hfill {\mathrm{\Psi }}^{″}(r,t)=e^{2\mathrm{i}{\gamma }_{n,n+1}}\frac{[-R_{n,n+1}^{\star }\mathrm{\Psi }(r,t)+{\mathrm{\Psi }}^{\prime }(r,t)]}{T_{n,n+1}},\end{array}$$(B.5)

which must be simultaneously applied in the cavity C_n, that is,

$$\begin{array}{c}\hfill e^{2\mathrm{i}{\gamma }_{n,n+1}}\frac{[-R_{n,n+1}^{\star }{\mathrm{\Psi }}_n(r,t;r_n^{+})+{\mathrm{\Psi }}_n^{\prime }(r,t;r_n^{+})]}{T_{n,n+1}}=T_{n,n+1}{\psi }^{(\leftarrow )}(r,t;r_n^{+}),\end{array}$$(B.6)

and in the cavity C_n + 1, that is,

$$\begin{array}{cc}\hfill e^{2\mathrm{i}{\gamma }_{n,n+1}}& \frac{[-R_{n,n+1}^{\star }{\mathrm{\Psi }}_{n+1}(r,t;r_{n+1}^{-})+{\mathrm{\Psi }}_{n+1}^{\prime }(r,t;r_{n+1}^{-})]}{T_{n,n+1}}\hfill \\ \hfill & =-R_{n,n+1}^{\star }{e}^{2\mathrm{i}{\gamma }_{n,n+1}}{\psi }^{(\to )}(r,t;r_{n+1}^{-})+{\psi }^{(\leftarrow )}(r,t;r_{n+1}^{-}).\hfill \end{array}$$(B.7)

By comparing Eq. (B.6) with Eq. (23) and Eq. (B.7) with Eq. (26), and using the definitions in Eqs. (30) and (31), we can conclude by identification:

$$\begin{array}{c}\hfill T_{n+1,n}=T_{n,n+1}\mathrm{and}{R}_{n+1,n}=|R_{n,n+1}|e^{\mathrm{i}(\pi -{\delta }_{n,n+1}+2{\gamma }_{n,n+1})}.\end{array}$$(B.8)

Appendix C: Alternative formulation for the resonance condition

Depending on the problem, it can be better to express the resonance condition as a series of sine terms that are functions of the Υ_i phases in Eq. (72). To do so, we first rewrite, for sake of the convenience, the resonance condition provided in Eqs. (74) and (75) in such a way to isolate the recurrence variable on the level i rather than on the level i − 1. This can be done in a simple way by defining τ_i = arctan[q_{i − 1, i}cotΦ_i − 1], so that Eqs. (74) and (75) are equivalent, for i ≥ 2, to the conditions

$$\begin{array}{c}\hfill tan{\tau }_i=-q_{i-1,i}tan({\mathrm{{\rm Y}}}_{i-1}-{\tau }_{i-1}),{\tau }_1=0\mathrm{and}{\tau }_N={\mathrm{{\rm Y}}}_N.\end{array}$$(C.1)

Second, we can proceed step by step. Using Eq. (67) and usual trigonometric formula, Eq. (C.1) for i = N provides

$$\begin{array}{c}\hfill sin({\mathrm{{\rm Y}}}_N+{\mathrm{{\rm Y}}}_{N-1}-{\tau }_{N-1})+\left|R_{N-1,N}\right|sin({\mathrm{{\rm Y}}}_N-{\mathrm{{\rm Y}}}_{N-1}+{\tau }_{N-1})=0.\end{array}$$(C.2)

Then, using Eq. (C.1) to express τ_N − 1 as a function of τ_N − 2 and using the identities cos(arctanx) = 1/(1 + x²)^1/2 and sin(arctanx) = x/(1 + x²)^1/2, we can express Eq. (C.2) as

$$\begin{array}{c}\hfill A_N+\left|R_{N-1,N}\right|B_N=0,\end{array}$$(C.3)

with

$$\begin{array}{cc}\hfill A_N\propto & sin({\mathrm{{\rm Y}}}_N+{\mathrm{{\rm Y}}}_{N-1}+{\mathrm{{\rm Y}}}_{N-2}-{\tau }_{N-2})\hfill \\ \hfill & +\left|R_{N-2,N-1}\right|sin({\mathrm{{\rm Y}}}_N+{\mathrm{{\rm Y}}}_{N-1}-{\mathrm{{\rm Y}}}_{N-2}+{\tau }_{N-2})\hfill \end{array}$$(C.4)

$$\begin{array}{cc}\hfill B_N\propto & sin({\mathrm{{\rm Y}}}_N-{\mathrm{{\rm Y}}}_{N-1}-{\mathrm{{\rm Y}}}_{N-2}+{\tau }_{N-2})\hfill \\ \hfill & +\left|R_{N-2,N-1}\right|sin({\mathrm{{\rm Y}}}_N-{\mathrm{{\rm Y}}}_{N-1}+{\mathrm{{\rm Y}}}_{N-2}-{\tau }_{N-2}),\hfill \end{array}$$(C.5)

where the proportionality constant is the same for A_N and B_N. We note that the expressions of A_N and B_N are very similar to the left-hand side of Eq. (C.2). Successively performing the same operation as before, that is, using Eq. (C.1) to express τ_N − 2 as a function of τ_N − 3, we can express Eqs. (C.4) and (C.5) as

$$\begin{array}{cc}\hfill A_N& \propto A_{N-1}^{(1)}+\left|R_{N-2,N-1}\right|B_{N-1}^{(1)}\hfill \end{array}$$(C.6)

$$\begin{array}{cc}\hfill B_N& \propto A_{N-1}^{(2)}+\left|R_{N-2,N-1}\right|B_{N-1}^{(2)},\hfill \end{array}$$(C.7)

where

$$\begin{array}{cc}\hfill A_{N-1}^{(1)}=& sin({\mathrm{{\rm Y}}}_N+{\mathrm{{\rm Y}}}_{N-1}+{\mathrm{{\rm Y}}}_{N-2}+{\mathrm{{\rm Y}}}_{N-3}-{\tau }_{N-3})\hfill \\ \hfill & +\left|R_{N-3,N-2}\right|sin({\mathrm{{\rm Y}}}_N+{\mathrm{{\rm Y}}}_{N-1}+{\mathrm{{\rm Y}}}_{N-2}-{\mathrm{{\rm Y}}}_{N-3}+{\tau }_{N-3})\hfill \end{array}$$(C.8)

$$\begin{array}{cc}\hfill B_{N-1}^{(1)}=& sin({\mathrm{{\rm Y}}}_N+{\mathrm{{\rm Y}}}_{N-1}-{\mathrm{{\rm Y}}}_{N-2}-{\mathrm{{\rm Y}}}_{N-3}+{\tau }_{N-3})\hfill \\ \hfill & +\left|R_{N-3,N-2}\right|sin({\mathrm{{\rm Y}}}_N+{\mathrm{{\rm Y}}}_{N-1}-{\mathrm{{\rm Y}}}_{N-2}+{\mathrm{{\rm Y}}}_{N-3}-{\tau }_{N-3})\hfill \end{array}$$(C.9)

$$\begin{array}{cc}\hfill A_{N-1}^{(2)}=& sin({\mathrm{{\rm Y}}}_N-{\mathrm{{\rm Y}}}_{N-1}-{\mathrm{{\rm Y}}}_{N-2}-{\mathrm{{\rm Y}}}_{N-3}+{\tau }_{N-3})\hfill \\ \hfill & +\left|R_{N-3,N-2}\right|sin({\mathrm{{\rm Y}}}_N-{\mathrm{{\rm Y}}}_{N-1}-{\mathrm{{\rm Y}}}_{N-2}+{\mathrm{{\rm Y}}}_{N-3}-{\tau }_{N-3})\hfill \end{array}$$(C.10)

$$\begin{array}{cc}\hfill B_{N-1}^{(2)}=& sin({\mathrm{{\rm Y}}}_N-{\mathrm{{\rm Y}}}_{N-1}+{\mathrm{{\rm Y}}}_{N-2}+{\mathrm{{\rm Y}}}_{N-3}-{\tau }_{N-3})\hfill \\ \hfill & +\left|R_{N-3,N-2}\right|sin({\mathrm{{\rm Y}}}_N-{\mathrm{{\rm Y}}}_{N-1}+{\mathrm{{\rm Y}}}_{N-2}-{\mathrm{{\rm Y}}}_{N-3}+{\tau }_{N-3})\hfill \end{array}$$(C.11)

The final expression can therefore be obtained by performing the same operation several times until reaching the variable τ₁ = 0. At the end of the day, the compact general form of the resonance condition reads

$$\begin{array}{c}\hfill sin{\displaystyle \left(\sum _{j=1}^N{\mathrm{{\rm Y}}}_j\right)+\sum _{i=1}^{N-1}\{\sum _{(p_k)\in E_i}\left(\prod _{k=1}^i\left|R_{p_k,p_k+1}\right|\right)sin\left(\sum _{j=1}^N{c}_j{\mathrm{{\rm Y}}}_j\right)\}=0,}\end{array}$$(C.12)

where E_i represents the set of the arrangements of i distinct elements in I_N − 1 = ⟦1, N − 1⟧ listed in descending order, that is, E_i = {(p_k)_{1 ≤ k ≤ i}} such as p_k ∈ I_N − 1 and p_k + 1 < p_k. The coefficients {c_j}_{1 ≤ j ≤ N} in the sine function are equal to either c_j = −1 if the k index of the lowest p_k value such as p_k ≥ j is odd or c_j = +1 otherwise (i.e., c_N is always equal to +1). We note that each set (E_i)_{i ∈ IN − 1} is composed of C(N − 1, i) arrangements, where C(N − 1, i) are the usual binomial coefficients. The number of terms in Eq. (C.12) is thus equal the sum of the number of arrangements in each set (E_i)_{i ∈ IN − 1} plus one for the first term, which results in 2^N − 1 terms.

Appendix D: Local mean mode energy in cavities

We detail in this section the calculation of Eq. (100) within the asymptotic limit. First, using Eqs. (15) and (17), we can get the JWKB form of the radial displacement inside a given cavity, that is,

$$\begin{array}{c}\hfill {\stackrel{\sim }{\xi }}_r={\left(\frac{1}{\rho r^{2}c\sigma }\right)}^{1/2}{\left|\frac{S_{\ell }^2-{\sigma }^2}{N^2-{\sigma }^2}\right|}^{1/4}(a_{\mathrm{p}}{e}^{i\phi }+a_{\mathrm{r}}{e}^{-i\phi }).\end{array}$$(D.1)

Using Eqs. (18) and (101), we can also get the JWKB form of the horizontal displacement function, that is,

$$\begin{array}{c}\hfill {\stackrel{\sim }{\xi }}_h=\mathrm{sgn}(S_{\ell }^2-{\sigma }^2)\mathrm{i}{\left(\frac{c}{\rho r^4{\sigma }^3}\right)}^{1/2}{\left|\frac{N^2-{\sigma }^2}{S_{\ell }^2-{\sigma }^2}\right|}^{1/4}(a_{\mathrm{p}}{e}^{i\phi }-a_{\mathrm{r}}{e}^{-i\phi }).\end{array}$$(D.2)

Second, we recall that the orthonormal spherical harmonics are given by

$$\begin{array}{c}\hfill Y_{\ell }^m(\theta ,s)=N_{\ell m}{P}_{\ell }^m(cos\theta )e^{\mathrm{i}ms},\end{array}$$(D.3)

where $P_{\mathcal{l}}^m$ are the associated Legendre polynomials. The normalization constant is given by

$$\begin{array}{c}\hfill N_{\ell m}=\sqrt{\frac{(2\ell +1)}{4\pi }\frac{(\ell -m)!}{(\ell +m)!}},\end{array}$$(D.4)

and is such that

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\pi }{\left[P_{\ell }^m(cos\theta )\right]}^{2}sin\theta \mathrm{d}\theta =\frac{1}{2\pi N_{\ell m}^2}.}\end{array}$$(D.5)

We also recall that the spherical harmonics are solutions of the eigenvalue equation $r^2{\nabla }^2{Y}_{\mathcal{l}}^m=-\mathcal{l}(\mathcal{l}+1)Y_{\mathcal{l}}^m$. Multiplying this eigenvalue equation by the complex conjugate of $Y_{\mathcal{l}}^m$, integrating over the sphere, using Eq. (D.5) for the right-hand side of the equations and integration by parts for the left-hand side, we can retrieve the well-known relation

$$\begin{array}{cc}& {\displaystyle {\int }_0^{\pi }[{\left(\frac{\mathrm{d}}{\mathrm{d}\theta }\left[P_{\ell }^m(cos\theta )\right]\right)}^2+\frac{m^2}{{sin}^2\theta }{\left[P_{\ell }^m(cos\theta )\right]}^2]sin\theta \mathrm{d}\theta =\frac{\ell (\ell +1)}{2\pi N_{\ell m}^2}·}\hfill \end{array}$$(D.6)

Finally, taking the square modulus of the real part of Eq. (2) while considering the JWKB form of the displacement and the spherical harmonics expression in Eqs. (D.1)–(D.3), integrating first on t and s, and then integrating over θ using Eqs. (D.5) and (D.6), Eq. (100) can be rewritten to a good approximation as

$$\begin{array}{cc}\hfill {\mathcal{E}}_k& \approx \frac{|a_{\mathrm{p}}{|}^2+{|a_{\mathrm{r}}|}^2}{2}{\displaystyle {\int }_{r_k^{-}}^{r_k^{+}}{\mathcal{K}}_r({|\frac{N^2}{{\sigma }^2}-1|}^{-1}+{|1-\frac{{\sigma }^2}{S_{\ell }^2}|}^{-1})\mathrm{d}r,}\hfill \end{array}$$(D.7)

where 𝒦_r is defined in Eq. (5) and where we have also used the fact that the wave number integral is much higher than unity within the asymptotic limit to neglect the residual oscillating radial functions in the integrand. Either in the case of pressure-dominated modes, namely, ${\sigma }^2\gg (S_{\mathcal{l}}^2\text{\hspace{0.33em}and\hspace{0.33em}}{N}^2)$, or gravity-dominated modes, namely, ${\sigma }^2\ll (S_{\mathcal{l}}^2\text{\hspace{0.33em}and\hspace{0.33em}}{N}^2)$, the term in brackets in the integrand of Eq. (D.7) is close to unity. Moreover, as noted in Sect. 4.5, |a_p|=|a_r|=|a_k| inside the cavity C_k in the case of eigenmodes with totally reflective boundary conditions. This justifies the validity of Eq. (102).

Appendix E: Supplementary expressions of the amplitude ratios

In Sect. 4.5, we use the surface boundary condition to express the mode amplitude in the cavity C_k as a function of that in the upper cavity C_N. Imposing instead the core boundary condition, we aim in this section to express a_k as a function of a₁ in the inner cavity C₁. To do so, we first invert the relation in Eq. (88) to obtain for 1 ≤ i ≤ N − 1

$$\begin{array}{c}\hfill {\mathit{a}}_{i+1,+}\equiv {\mathsf{E}}_{i+1,i}{\mathit{a}}_{i,+}\end{array}$$(E.1)

where

$$\begin{array}{cc}\hfill {\mathsf{E}}_{i+1,i}& \equiv {\mathsf{E}}_{i,i+1}^{-1}=\mathsf{C}\left({\mathrm{\Theta }}_{i+1}\right)\mathsf{C}(\frac{{\delta }_{i+1,i}}{2}-\frac{\pi }{2}){\mathsf{A}}_{i+1,i}\mathsf{C}\left(\frac{{\delta }_{i,i+1}}{2}\right)\hfill \end{array}$$(E.2)

$$\begin{array}{cc}\hfill {\mathsf{A}}_{i+1,i}& \equiv {\mathsf{A}}_{i,i+1}^{-1}=\frac{1}{\left|T_{i,i+1}\right|}\left(\begin{array}{cc}1& -\left|R_{i,i+1}\right|\\ -\left|R_{i,i+1}\right|& 1\end{array}\right),\hfill \end{array}$$(E.3)

with C(ε) and A_{i, i + 1} defined in Eqs. (81) and (82). Using Eq. (84), it is then straightforward to show that the application of E_{i + 1, i} on a vector in the form of Eq. (83) leads to

$$\begin{array}{c}\hfill {\mathsf{E}}_{i+1,i}\mathit{w}(\epsilon )={\mathcal{A}}_{i+1,i}(\epsilon +\frac{{\delta }_{i,i+1}}{2})\mathit{w}\left[{\mathrm{\Omega }}_{i+1}^{-1}(\epsilon )\right],\end{array}$$(E.4)

where

$$\begin{array}{c}\hfill {\mathcal{A}}_{i+1,i}(\epsilon )=\sqrt{\frac{1+(q_{i,i+1}^2-1){cos}^2\left(\epsilon \right)}{q_{i,i+1}}},\end{array}$$(E.5)

and in which we have used the convention for the branch of the arctan function in Eq. (71). Within this convention, ${\mathrm{\Omega }}_{i+1}^{-1}$ is the inverse function of Ω_i + 1 in Eq. (90) and is defined for 1 ≤ i ≤ N − 1 as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{i+1}^{-1}(\epsilon )=arctan[\frac{1}{q_{i,i+1}}tan(\epsilon +\frac{{\delta }_{i,i+1}}{2})]+{\mathrm{{\rm Y}}}_{i+1}-\frac{{\delta }_{i+1,i+2}}{2}·\end{array}$$(E.6)

Second, using the inverse of Eq. (87), the core boundary condition implies that

$$\begin{array}{c}\hfill {\mathit{a}}_{1,+}=\mathsf{C}({\mathrm{\Theta }}_1)a_1\mathit{w}\left(\frac{{\delta }_{\mathrm{c}}}{2}\right),\end{array}$$(E.7)

Finally, applying k − 1 times the transformation in Eq. (E.1) on the upper core amplitude vector in Eq. (E.7), we can obtain the amplitude vector in the cavity C_k for k > 1, that is,

$$\begin{array}{c}\hfill {\mathit{a}}_{k,+}={\displaystyle \left(\prod _{i=k-1}^1{\mathsf{E}}_{i+1,i}\right)\mathsf{C}({\mathrm{\Theta }}_1)a_1\mathit{w}\left(\frac{{\delta }_{\mathrm{c}}}{2}\right).}\end{array}$$(E.8)

Using the properties in Eqs. (84) and (E.4), Eq. (E.8) can be written for k > 1 as

$$\begin{array}{c}\hfill {\mathit{a}}_{k,+}=a_k\mathit{w}[{\mathrm{\Omega }}_k^{-1}°\cdots °{\mathrm{\Omega }}_1^{-1}\left(0\right)],\end{array}$$(E.9)

where (° ) is the composition operator, ${\mathrm{\Omega }}_1^{-1}$ is the inverse of the function in Eq. (99), and a_k is given by

$$\begin{array}{c}\hfill \frac{a_k}{a_1}={\displaystyle \prod _{i=k-1}^1{\mathcal{A}}_{i+1,i}[{\mathrm{\Omega }}_i^{-1}°···°{\mathrm{\Omega }}_1^{-1}\left(0\right)+\frac{{\delta }_{i,i+1}}{2}].}\end{array}$$(E.10)

According to Eqs. (E.5) and (E.10), we can therefore deduce that the squared amplitude ratio between the adjacent cavities C_k and C_k + 1 is equal to

$$\begin{array}{c}\hfill {\left|\frac{a_k}{a_{k+1}}\right|}^2={\displaystyle \frac{q_{k,k+1}}{q_{k,k+1}^2+(1-q_{k,k+1}^2){cos}^2\left({\mathrm{\Xi }}_k^{\prime }\right)}},\end{array}$$(E.11)

where we have defined

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_k^{\prime }={\mathrm{\Omega }}_k^{-1}°\cdots °{\mathrm{\Omega }}_1^{-1}\left(0\right)+\frac{{\delta }_{k,k+1}}{2}+\frac{\pi }{2},\end{array}$$(E.12)

which measures the level of constructive interference in the cavity C_k. The same discussion as that based on Eq. (105) in Sect. 5.2 holds true here when only the core boundary condition is met, except that Ξ_k + 1 and C_k + 1 have to be replaced by ${\Xi }_k^{\prime }$ and C_k. We insist on the fact that Eqs. (105) and (E.11) are strictly equivalent only within the resonance condition, that is, when the core and surface boundary conditions are simultaneously met. Under this condition, we note that we can also obtain another equivalent relation by multiplying both equations, which is reduced to:

$$\begin{array}{c}\hfill {\left|\frac{a_k}{a_{k+1}}\right|}^4={\displaystyle \frac{q_{k,k+1}^2+(1-q_{k,k+1}^2){cos}^2\left({\mathrm{\Xi }}_{k+1}\right)}{q_{k,k+1}^2+(1-q_{k,k+1}^2){cos}^2\left({\mathrm{\Xi }}_k^{\prime }\right)}}.\end{array}$$(E.13)

We notice that Eqs. (105), (E.11) and (E.13) generalize the expressions of the amplitude ratios derived in Sect. 3.2 of Takata (2016b) in the case of mixed modes.

Appendix F: Usual boundary reflection phase lags

In this section, we estimate in a simple way the reflection phase lags introduced near the center and surface, δ_s and δ_c. The surface and core boundaries are supposed to be turning points denoted by r_t and such as ${\mathcal{K}}_r^2(r_{\text{t}})=0$, at the interface between a resonant cavity and an evanescent region. Assuming the JWKB approximation is met in the adjacent resonant cavities, it is then possible to use a single turning point analysis of the stellar oscillations equations such as performed by Shibahashi (1979) to deduce the wave function in the vicinity of each of these turning points. The analysis uses in addition the Cowling approximation, neglects the gradients in the equilibrium structure everywhere in the star and considers only regular solutions in the evanescent region (i.e., with finite values of the wave function). Within this framework, the asymptotic solution of the wave function inside the resonant cavities can be expressed and the associated reflection phase lag turns out to depend on the nature of the considered turning point.

F.1. In the case of a turning point such as σ² = N²(r_t),

Eq. (3) can be solved using a Green-Liouville transformation while neglecting the variations of the equilibrium structure (e.g., the M function). The regular solution of the transformed equation thus takes the form of an Airy function of the first kind. The asymptotic expansion of the wave function in the resonant cavity is finally given using a complex notation by (e.g., Shibahashi 1979; Unno et al. 1989)

$$\begin{array}{c}\hfill \mathrm{\Psi }\approx \frac{a}{\sqrt{{\mathcal{K}}_r}}(\underset{A}{\underbrace{e^{-\mathrm{i}\pi /4}}}{e}^{\mathrm{i}|\phi (r;r_{\mathrm{t}})|}+\underset{B}{\underbrace{e^{\mathrm{i}\pi /4}}}{e}^{-\mathrm{i}|\phi (r;r_{\mathrm{t}})|})e^{-\mathrm{i}\sigma t},\end{array}$$(F.1)

where a is a complex constant, φ(r; r_t) is defined in Eq. (11), and we choose in this case the origin of the wave phase at the turning point r_t. Since the progressive and regressive wave components have the same amplitude modulus, we first note that such a regular solution satisfies the totally reflective condition at the turning point, as assumed close to the core and surface in this work. Second, the phase lag at reflection δ_t can be computed as the argument of the amplitude ratio of the upward to the downward energy ray when the resonant cavity is above the turning point (i.e., r > r_t). When the resonant cavity is below the turning point (i.e., r < r_t), the phase lag is equal to the argument of the amplitude ratio of the downward to the upward energy ray. In both cases, the phase lag at reflection is reduced to

$$\begin{array}{c}\hfill {\delta }_{\mathrm{t}}=\mathrm{sgn}({\sigma }^2-S_{\ell }^2)arg\left(\frac{A}{B}\right)=\mathrm{sgn}(S_{\ell }^2-{\sigma }^2)\frac{\pi }{2},\end{array}$$(F.2)

where the first sign factor estimated well inside the cavity accounts for the case of gravity-dominated waves (i.e., ${\sigma }^2\ll S_{\mathcal{l}}^2$) for which the group velocity is in the opposite direction of the phase velocity.

F.2. In the case of a turning point such as ${\sigma }^2=S_{\mathcal{l}}^2(r_{\text{t}})$,

Eq. (3) is singular at r_t so that we first have to solve Eq. (A.7) for the dependent variable $\stackrel{\sim }{\mathrm{\Psi }}$. According to Shibahashi (1979), the asymptotic expansion of $\stackrel{\sim }{\mathrm{\Psi }}$ in the resonant cavity is also given by Eq. (F.1). To compute the wave function Ψ (which is the reference dependent variable we consider in this paper), we can use Eq. (A.10) while neglecting the derivative of the Q function related to the variations of the equilibrium structure. Considering the fact that (d|φ|/dr) = sgn(r − r_t) |𝒦_r|, we obtain:

$$\begin{array}{c}\hfill \mathrm{\Psi }\approx \frac{\mathrm{i}a}{\sqrt{{\mathcal{K}}_r}}\mathrm{sgn}\left[\frac{{\sigma }^2-N^2}{r-r_{\mathrm{t}}}\right]\left(\underset{A^{\prime }}{\underbrace{e^{-\mathrm{i}\pi /4}}}{e}^{\mathrm{i}|\phi (r;r_{\mathrm{t}})|}\underset{B^{\prime }}{\underbrace{-e^{\mathrm{i}\pi /4}}}{e}^{-\mathrm{i}|\phi (r;r_{\mathrm{t}})|}\right)e^{-\mathrm{i}\sigma t},\end{array}$$(F.3)

Similarly to the previous case we describe above, we thus find that the phase lag at reflection is reduced to:

$$\begin{array}{c}\hfill {\delta }_{\mathrm{t}}=\mathrm{sgn}({\sigma }^2-S_{\ell }^2)arg\left(\frac{A^{\prime }}{B^{\prime }}\right)=\mathrm{sgn}({\sigma }^2-S_{\ell }^2)\frac{\pi }{2}·\end{array}$$(F.4)

F.3. In terms of synthesis,

the comparison between Eqs. (F.2) and (F.4) shows that the change of nature in the turning point is associated with a shift in the phase lag at reflection of π (modulo 2π). Moreover, considering a mode with one single cavity located between two turning points $r_1^{-}$ and $r_1^{+}$ at the core and the surface, we see that δ_c + δ_s = π (modulo 2π) if the turning points have the same nature, namely, if ${\sigma }^2=N^2(r_1^{-})=N^2(r_1^{+})$ or ${\sigma }^2=S_{\mathcal{l}}^2(r_1^{-})=S_{\mathcal{l}}^2(r_1^{+})$; whereas it is δ_c + δ_s = 0 otherwise.

All Figures