© ESO 2020

1. Introduction

Over the last decades, asteroseismology has been proven to be a powerful method to probe the interior of a large variety of stars at diverse evolutionary stages (e.g., Di Mauro 2016, and references therein). Among different classes of pulsations are the solar-like oscillations, which are global stable modes that are both excited and damped by turbulence in the outermost convective layers of stars. The frequency pattern of these modes is known to be an important source of constraints on stellar properties (e.g., Samadi et al. 2015). However, because of intrinsically low amplitudes, the observation of such oscillations in distant stars were difficult before the advent of the photometric space-borne missions CoRoT (Baglin et al. 2006a,b) and Kepler (Borucki et al. 2010). Indeed, both missions marked a key turning point in the field. They provided high-quality seismic data that enabled the detection of solar-like oscillations in thousands of low-mass stars from the main sequence to the red giant phases. This large number of constraints subsequently led to a revolution in our understanding of stellar structure and stellar evolution (see Chaplin & Miglio 2013, for a review).

Among the main successes of the aforementioned missions was the detection of mixed modes in red giant and red clump stars (e.g., Bedding et al. 2010; Beck et al. 2011; Mosser et al. 2011). Mixed modes can propagate in two resonant cavities located in the inner radiative zone, where they behave as gravity modes, and in the external envelope, where they behave as pressure modes. Both cavities are coupled by an intermediate region located between the profiles of the Brunt–Väisälä and Lamb frequencies where mixed modes are evanescent (e.g., Hekker & Christensen-Dalsgaard 2017, for a recent review; see also Fig. 1). After the exhaustion of the central hydrogen at the end of the main sequence, the structural changes resulting from the contraction of the helium core and the expansion of the envelope in evolved stars made such a coupling possible. The frequency pattern of mixed modes observed in these stars contains the signature of the innermost layers in which they propagate and thus has the potential to provide stringent constraints on their properties.

First, the detection of mixed modes gave us access to the period spacing of dipolar modes, denoted with ΔΠ₁, which was measured in about 6100 evolved stars observed by the Kepler satellite (e.g., Mosser et al. 2012b; Vrard et al. 2016). This seismic index was shown to be correlated with the mass of the helium core (e.g., Montalbán et al. 2013). Besides, the frequency large separation between two consecutive radial modes, denoted with Δν, is sensitive to the mean envelope density; together, they form a couple of indicators of the evolutionary state during the post-main sequence that permits to unambiguously discriminate between hydrogen-burning shell stars on the red giant branch (hereafter, RGB) and helium-burning core stars on the red clump (e.g., Bedding et al. 2011; Mosser et al. 2012b, 2014). Second, the measurement of rotationally-induced frequency splittings provided an estimate of the core rotation rate in hundreds of evolved stars (e.g., Gehan et al. 2018, and references therein) that turned out to be much lower than predicted in current stellar models (e.g., Marques et al. 2013; Cantiello et al. 2014). These seismic observations represent a goldmine of information to identify and characterize the missing mechanism that is responsible for these slow core rotations among several candidates (e.g., Belkacem et al. 2015a,b; Spada et al. 2016; Pinçon et al. 2017, and references therein).

All the seismic constraints provided by mixed modes resulted not only in a better characterization of the interior of evolved stars themselves, but also of the structural properties of their past and future evolutionary stages (e.g., Montalbán et al. 2013; Cantiello et al. 2014, for such investigations). Moreover, such advancements on a large sample of stars benefited to other connected fields as Galactic archaeology, for which evolved stars represent key targets owing to their high luminosity and their wide mass range (e.g., Miglio et al. 2013; Mosser et al. 2016; Noels et al. 2016). However, all the probing potential of mixed modes has not been explored yet. In particular, the gravity offset and coupling factor of mixed modes, denoted with ε_g and q, respectively, are expected to bring additional constraints on the red giant structure. These seismic parameters, which characterize the frequency pattern of mixed modes, were measured in a large sample of stars by Mosser et al. (2017, 2018). The variations in the gravity offset during the evolution on the RGB were recently shown by Pinçon et al. (2019) to be sensitive to the density contrast between the core and the envelope, and to the location of the base of the convective region. Regarding the coupling factor, the relation between the observed values and the internal properties nonetheless still calls for further investigations.

The coupling factor describes the degree of interaction between the central and the external resonant cavities of mixed modes. The theoretical link between q and the internal structure was originally provided by the asymptotic analysis of mixed modes of Shibahashi (1979) in the limiting case of a thick intermediate evanescent zone (sometimes called the weak coupling hypothesis). This study showed that q depends on the wave transmission coefficient through the evanescent zone, which is located between the hydrogen-burning shell and the neighborhood of the base of the convective zone. This was confirmed later by Takata (2016a,b) from a more general point of view. While Δν and ΔΠ₁ are sensitive to the external layers and the surrounding of the helium core, respectively, the coupling factor thus promises to bring us substantial information on this intermediate region, where most of the star luminosity is produced at this evolutionary stage and where mixing processes at the interface between the convective and radiative zones can occur.

The first measurements of q in Kepler evolved stars were obtained by fitting the asymptotic expression derived by Shibahashi (1979) to the observed frequency pattern of dipolar mixed modes, considering q as a free parameter (e.g., Mosser et al. 2012b; Buysschaert et al. 2016; Vrard et al. 2016). All these works demonstrated that this asymptotic analytical form well reproduces the observed frequency pattern but requires values of q higher than 1/4 in some stars, in disagreement with the analytical predictions obtained in the weak coupling hypothesis. Jiang & Christensen-Dalsgaard (2014) reached the same conclusions by comparing the asymptotic frequencies with the theoretical frequencies computed via an oscillation code in evolved stellar models (see also Hekker et al. 2018). These works underlined the need for going beyond the weak coupling hypothesis, which is a prerequisite for fully grasping the information brought by q on the internal structure. Motivated by such discrepancies, Takata (2016a) then developed an asymptotic analysis of mixed modes in the reverse limiting case of a very thin evanescent zone (sometimes called the strong coupling hypothesis). In this formalism, the expression of mixed modes is similar to that in the weak coupling hypothesis, except that q can take values between zero and unity. This result validated the fitting procedures used in the previous observational studies, which assumed q as a free parameter between zero and unity, as well as their outcomes.

In parallel to the work of Takata (2016a), Mosser et al. (2017) proposed a new automated method to measure q in a large set of stars. This in-depth investigation confirmed previous observational results and provided a coherent view of the coupling factor for dipolar modes during the post-main sequence evolution. The large-scale measurement in more than 5000 evolved stars showed that the value of q changes during evolution with clear signatures at the transition from the subgiant phase to the RGB as well as to the red clump, where high values of q are observed, in qualitative agreement with the recently-proposed strong coupling formalism. Assuming the evanescent region is located in the radiative layers just above the helium core, as suggested in stellar models of young red giant stars, Mosser et al. (2017) showed that the value of q depends both on the variation scale height of the Brunt–Väisälä frequency and on the thickness of the evanescent region; however, they could not directly disentangle the origin of the observed variations in q among these two properties.

In this series of papers, we aim to clarify further the link between the coupling factor of mixed modes and the internal properties of evolved stars from a theoretical point of view. As a stepping stone to the investigation, an analytical approach is considered in this first paper. Using the available asymptotic analyses (i.e., both in the weak and the strong coupling hypotheses), we search for simplified expressions of q and scrutinize in detail the physical information encapsulated inside this parameter and its variations during evolution. The paper is organized as follows. In Sects. 2 and 3, the theoretical background about mixed modes is introduced within the framework of the asymptotic limit, with emphasis on the coupling factor and its physical meaning. Following this introductory material, a simple model is then developed in Sect. 4. It is based on the power-law behavior of the Brunt–Väisälä and Lamb frequencies in the evanescent region, as assumed by Mosser et al. (2017), but in contrast also accounts for the migration of the evanescent region from the radiative to the convective zones during the evolution on the RGB. The relation between q and the properties of the intermediate evanescent region is subsequently examined. The observations are then interpreted in the light of the analytical model in Sect. 5. The results and main assumptions used in this work are discussed in Sect. 6 and the conclusions are formulated in Sect. 7.

2. Mixed mode cavities in evolved stars

We introduce the basic background about mixed modes in evolved stars within the asymptotic framework, which is in the short-wavelength WKB approximation (e.g., Gough 2007). The structures of the resonant cavities and the coupling evanescent region are presented. In the following, we focus on dipolar ℓ = 1-pagination modes, with ℓ the angular degree, because of the peculiar role they play from an observational point of view. Indeed, among mixed modes whose energy is mainly trapped in the core of evolved stars and that are thus mostly sensitive to the properties of the innermost layers, the dipolar ones are the most easily detectable (e.g., Dupret et al. 2009; Benomar et al. 2014; Grosjean et al. 2014). As a result, the measurement of the coupling factor in a large set of stars is currently mainly limited to dipolar modes with Kepler data.

2.1. Dispersion relation

Solar-like oscillation modes result from the constructive interferences of gravito-acoustic waves stochastically excited in the uppermost layers of stars and traveling back-and-forth several times between the center and the surface. The propagation of these waves is physically ruled by a fourth-order linear differential equation under the adiabatic approximation (e.g., Unno et al. 1989). In order to analytically study stellar oscillations, the Cowling approximation (e.g., Cowling 1941), which consists in neglecting the perturbation of the gravitational potential denoted with Φ′ in the following, was usually adopted in previous works. This reduces the adiabatic oscillation equations from the fourth to the second order, and thus makes these tractable with usual asymptotic methods. However, this hypothesis remains questionable for dipolar modes because of their very large horizontal wavelength (see Appendix A.1). Fortunately, owing to the Hamiltonian structure of the adiabatic stellar oscillations and the existence of a general first integral for dipolar modes, Takata (2005, 2006a,b) demonstrated that the dipolar wave equation can be rewritten in the form of a second-order differential equation, which has the advantage of fully considering the perturbation of the gravitational potential and being analytically tractable. To do this, the dependent variables of the governing equations have to be modified from the standard mass displacement to the reduced displacement, which means the displacement of mass elements minus that of the center of mass of the corresponding concentric mass.

In the asymptotic limit, Takata (2006a, 2016a) demonstrated that accounting for the perturbation of the gravitational potential modifies the dispersion relation of dipolar modes; in this case, this takes the form of

$$\begin{array}{c}\hfill k_{\mathrm{r}}^2=\frac{{\sigma }^2}{c^2}(\frac{{\mathcal{N}}^2}{{\sigma }^2}-1)(\frac{{\mathcal{S}}_1^2}{{\sigma }^2}-1),\end{array}$$(1)

where k_r is the local radial wavenumber, σ is the angular mode frequency, c² = (Γ₁p)/ρ is the squared sound speed, Γ₁ is the first adiabatic index, and p and ρ are the equilibrium pressure and density, respectively. In this latter equation, 𝒩 and 𝒮₁ are the modified Brunt–Väisälä and Lamb frequencies (for ℓ = 1), respectively, which account for the perturbation of the gravitational potential, and are defined as

$$\begin{array}{c}\hfill \mathcal{N}=\frac{N}{J}\mathrm{and}{\mathcal{S}}_1=S_{1}J.\end{array}$$(2)

where N and S_ℓ are the usual Brunt–Väisälä and Lamb frequencies defined as (e.g., Unno et al. 1989)

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

$$\begin{array}{cc}\hfill 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}$$(4)

in which r is the radius in the star and g is the gravitational acceleration. The J factor in Eq. (2) is equal to

$$\begin{array}{c}\hfill J=1-\frac{1}{3}\frac{\mathrm{d}ln{M}_r}{\mathrm{d}lnr}=1-\frac{4\pi r^3\rho }{3M_r}=1-\frac{\rho (r)}{{\rho }_{\mathrm{av}}(r)},\end{array}$$(5)

in which M_r and ρ_av(r) represent the mass and mean density inside the spherical shell of radius r, respectively. The continuity equation (dM_r/dr) = 4πr²ρ has been used in the second equality of Eq. (5). Under the Cowling approximation, the dispersion relation has the same form as Eq. (1), except that 𝒩 and 𝒮₁ are replaced by the usual critical frequencies N and S₁, respectively (e.g., Unno et al. 1989). Indeed, the adiabatic oscillation equations for dipolar modes obtained within the Cowling approximation can be retrieved from the set of Eqs. (1)–(9) provided by Takata (2016a) in the non-Cowling case, in the limits of J → 1 and Φ′→0.

As an illustration, the modified Brunt–Väisälä and ℓ = 1 Lamb frequencies are displayed in Fig. 1 as a function of radius for three 1.2 M_⊙ stellar models. These models correspond to stars on the subgiant branch (SG), at the beginning of the ascent of the RGB (YR) and just before the bump of luminosity (ER). The modified Brunt–Väisälä and Lamb frequencies are compared to the usual expressions given in Eqs. (3) and (4) in the same figure. The internal radiative zone and the adiabatic external convective zone, where N² ≈ 0, can be easily distinguished. To be complete, the J factor is also plotted for the three considered models in Fig. 2. The stellar models were computed with the CESTAM evolution code (Marques et al. 2013), considering a solar mixture following Asplund et al. (2009), solar calibrated initial helium and metal abundances Y₀ = 0.25 and Z₀ = 0.013, and a solar calibrated mixing-length parameter α_MLT = 1.65. The input physics is standard, without microscopic diffusion, rotation nor overshooting. The models and their evolutionary tracks are indicated in the Hertzsprung–Russell diagram (HR diagram) in Fig. 3. As shown in Fig. 1, a clear difference can be observed between the usual and the modified critical frequencies below the hydrogen-burning shell, because of very low values of J. In the helium core, the density is nearly uniform and we can show by a second-order Taylor expansion near the center that $J\approx (-7{\rho }_{\mathrm{c}}″/10{\rho }_{\mathrm{c}})r^2$, with ρ_c and ${\rho }_{\text{c}}^{\prime \prime }$ the density and its second-derivative with respect to radius at the center. The effect of the Cowling approximation on the behavior of dipolar mixed modes is thus important in the central layers. In contrast, the usual and modified critical frequencies share comparable behavior above the hydrogen-burning shell (i.e., N ≈ 𝒩 and S₁ ≈ 𝒮₁). However, Jiang & Christensen-Dalsgaard (2014) and Mosser et al. (2017) showed using stellar models that the Cowling approximation can lead to incorrect estimates of the mode coupling through the evanescent region. We therefore keep the influence of the perturbation of the gravitational potential in the following by considering the modified expression of the critical frequencies. The effect of the Cowling approximation on the mode coupling is discussed in more detail in Sect. 3.3.1, after having properly introduced the coupling factor in the asymptotic framework.

2.2. Propagation diagram in evolved stars

As shown by Eq. (1), the propagation and trapping of waves in stars are determined by the profiles of the Brunt–Väisälä and Lamb frequencies, and on the typical oscillation frequency. Excitation models by turbulent convection (e.g., Samadi et al. 2015) and observations (e.g., see Fig. 6 in Mosser et al. 2013a) show that solar-like oscillations are efficiently generated in a narrow frequency range around a certain frequency at maximum oscillation power, which is denoted with ν_max. Empirical and theoretical studies demonstrated that ν_max scales as the acoustic cut-off frequency at the stellar surface (e.g., Brown et al. 1991; Kjeldsen & Bedding 1995; Belkacem et al. 2011, 2013). It results in a scaling relation between ν_max, the mass M, the radius R and the effective temperature T_eff of a given star, that reads

$$\begin{array}{c}\hfill {\nu }_{\max }\approx {\nu }_{\max ,\odot }\left(\frac{M}{M_{\odot }}\right){\left(\frac{R}{R_{\odot }}\right)}^{-2}{\left(\frac{T_{\mathrm{eff}}}{T_{\mathrm{eff},\odot }}\right)}^{-1/2},\end{array}$$(6)

which has been scaled for our purpose to the solar reference value ν_max, ⊙ = 3104 μHz such as prescribed by Mosser et al. (2013b). Typically, mixed modes are visible over a 6Δν-wide frequency range around ν_max (e.g., Grosjean et al. 2014; Mosser et al. 2018). Here, Δν represents the frequency spacing between two pressure radial modes with consecutive radial orders, the asymptotic expression of which is given by (e.g., Unno et al. 1989)

$$\begin{array}{c}\hfill \mathrm{\Delta }\nu ={\left(2{\displaystyle {\int }_0^R}\frac{\mathrm{d}r}{c}\right)}^{-1}·\end{array}$$(7)

Such typical observed frequency ranges are represented in Fig. 1. We find ν_max ≈ 700 260, and 50 μHz and Δν ≈ 43, 20, and 6 μHz for the models SG, YR, and ER, respectively.

According to Eq. (1), propagation diagrams in Fig. 1 show that observed oscillation modes can propagate (i.e., $k_{\text{r}}^2>0$) in a central radiative cavity where ${\sigma }^2<({\mathcal{N}}^2 \text{and} {\mathcal{S}}_1^2)$ and they behave as gravity modes (buoyancy cavity), and in an outer cavity where ${\sigma }^2>({\mathcal{N}}^2 \text{and} {\mathcal{S}}_1^2)$ and they behave as pressure modes (acoustic cavity). In these typical models, we also see that considering an asymptotic description is reasonable since ν_max ≪ N/2π and ν_max ≫ S₁/2π inside both cavities, respectively, so that k_r ≫ 1/r according to Eq. (1) and the number of radial nodes is much larger than unity. Both resonant cavities are coupled by an intermediate evanescent zone (i.e., $k_{\text{r}}^2<0$) located between the profiles of the (modified) Brunt–Väisälä and Lamb frequencies, through which the mode energy is transmitted from a cavity to another cavity. More exactly, the evanescent region is located between the turning points r₁(σ) and r₂(σ) where $k_{\text{r}}^2$ vanishes, that is, where 𝒩(r₁) = σ and 𝒮₁(r₂) = σ. These are the so-called mixed modes with a dual character of pressure and gravity modes. Such a configuration of coupling is not encountered in main sequence stars since ν_max is greater than the maximum value of the Brunt–Väisälä frequency in the radiative zone. Nevertheless, as a star evolves on the post-main sequence, its helium core contracts and becomes denser and denser. As observed in Fig. 1, this leads to an increase in the value of the Brunt–Väisälä frequency near the hydrogen-burning shell, denoted with N_M, that is, just above the helium core. Indeed, since N² ∼ g/r at first order in the radiative zone, with g = GM_r/r², we have $N_{\mathrm{M}}^2\sim G{M}_{\mathrm{He}}/R_{\mathrm{He}}^3\propto G{\overline{\rho }}_{\mathrm{He}}$ where M_He, R_He and ${\overline{\rho }}_{\mathrm{He}}$ are the mass, radius and mean density of the helium core, respectively. Similarly, using the hydrostatic equilibrium and the sound speed definition, we can show that $S_1^2\sim g/r$ at first order in this region (see Eq. (B.5) for instance); the profile of the Lamb frequency thus also increases during evolution in the inner layers. In contrast to the core contraction, the envelope expands (i.e., R increases) while the stellar luminosity L increases and the effective temperature decreases. From the Stefan–Boltzmann relation, which is $L\propto R^2{T}_{\text{eff}}^4$, it is thus obvious that the effective temperature T_eff decreases more slowly than $1/\sqrt{R}$. According to Eq. (6) at a given mass, ν_max thus decreases during evolution, as seen in Fig. 1. As a result, the combined effects of the core contraction and envelope expansion make the coupling between the inner buoyancy cavity and the outer acoustic cavity possible in evolved stars. This happens as soon as the subgiant stage starts (see the transition between the dashed and solid lines in Fig. 3).

Figure 1 also shows that the intermediate evanescent zone that couples the buoyancy and acoustic cavities evolves on the post-main sequence. It even radically changes in nature and in structure at some point during the ascent on the RGB. Indeed, for subgiant stars (e.g., model SG) and young red giant stars (e.g., model YR), the evanescent zone is located inside the radiative zone, between the hydrogen-burning-shell and the base of the convective zone. In the following, we call such a coupling zone Type-a evanescent region. As shown in Fig. 1, the radial extent of Type-a evanescent regions can be very thin on the subgiant branch (e.g., model SG), and becomes slightly thicker in young red giant stars (e.g., model YR). At one point during the evolution, the value of ν_max becomes low enough for the evanescent region to be entirely located in the convective zone and its lower bound to be very close to the radius at the base of the convective zone, denoted with r_b (e.g., model ER). In the following, we call such a coupling zone Type-b evanescent region. As shown in Fig. 1, Type-b evanescent regions are expected to be much thicker than Type-a evanescent regions. The transition phase between both configurations is indicated by a red square in the HR diagram of Fig. 3 and turns out to occur before the luminosity bump. In red clump stars, the structural readjustments resulting from the helium burning ignition lead again to Type-a evanescent regions (e.g., see Fig. 9 in Hekker et al. 2018).

We note that Type-a and Type-b evanescent regions generally match the cases a and b considered in Pinçon et al. (2019), which are associated with stars such as ν_max ≳ 𝒩_b/2π and ν_max ≲ 𝒩_b/2π, respectively, where 𝒩_b is the value of the modified Brunt–Väisälä frequency just below the base of the convective region. Using standard stellar models with the same input physics as used to build the stellar models considered in Fig. 1, Pinçon et al. (2019) estimated that the transition from case a to case b, thus corresponding to the transition from Type-a to Type-b evanescent regions, occurs for values of ν_max equal to ν_max, t ≈ 110, 95, 80, 70, 60, and 50 μHz, for typical stellar masses of 1, 1.2, 1.4, 1.6, 1.8, and 2 M_⊙, respectively. Using the ν_max − Δν scaling law $\mathrm{\Delta }\nu \approx 0.28{\nu }_{\max }^{0.75}$ with Δν and ν_max expressed in μHz (e.g., Mosser et al. 2012a; de Assis Peralta et al. 2018), this transition is expected to occur for values of Δν between 10 μHz and 5 μHz in the same mass range. This change of configuration during evolution impacts the coupling of mixed modes through the intermediate evanescent region and therefore needs to be considered.

2.3. Structure of the evanescent coupling region

The actual shape of the evanescent zone, and thus the coupling between the resonant cavities, depend on the profiles of the (modified) Brunt–Väisälä and Lamb frequencies since, according to Eq. (1), these latter define the decay length of the mode amplitude in this region. Actually, as shown in Appendix B.1, both 𝒩 and 𝒮₁ profiles are intrinsically related to the mass distribution inside the star, represented by the J factor in Eq. (5).

In the radiative zone, between the hydrogen-burning shell and the base of the convective region, Fig. 2 shows that J remains nearly uniform with values higher than 0.7 in the considered models, indicating a high density contrast between the considered layers and the helium core. Using the internal structure equations and simple analytical models, we demonstrate in Appendix B.2 that this high density contrast is the main responsible for the similar power-law behavior of 𝒩 and 𝒮₁ with respect to radius in this region (see Fig. 1). Actually, the parallel profiles of 𝒩 and 𝒮₁ result from their dependence in (g/r)^1/2, as shown through Eqs. (D.7) and (D.8) considering (dlnp/dlnr), (dlnρ/dlnr) and J as nearly uniform. Since by definition 𝒩(r₂) = σ and 𝒮₁(r₁) = σ, we thus can write in a good approximation in the upper radiative layers (e.g., Takata 2016a; Mosser et al. 2017)

$$\begin{array}{c}\hfill \mathcal{N}(r)\approx \sigma {\left(\frac{r_1}{r}\right)}^{\beta }\mathrm{and}{\mathcal{S}}_1(r)\approx \sigma {\left(\frac{r_2}{r}\right)}^{\beta },\end{array}$$(8)

where β is defined as the (minus) logarithmic derivative of 𝒮₁ in the middle of the evanescent region at r_ev, that is,

$$\begin{array}{c}\hfill \beta \equiv -{\left(\frac{\mathrm{d}ln{\mathcal{S}}_1}{\mathrm{d}lnr}\right)}_{r_{\mathrm{ev}}}.\end{array}$$(9)

In this power-law configuration, a simple relation exists between β, J and the local polytropic exponent, γ = (dlnp/dlnρ), which reads (see Appendix B.1)

$$\begin{array}{c}\hfill \beta \approx \frac{3}{2}J\approx \frac{1}{2-\gamma }·\end{array}$$(10)

Considering both limiting cases of an isothermal perfect gas or a relativistic photon gas, we can expect the values of γ to range between about 1 and 4/3 in the radiative layers surrounding the helium core. According to Eq. (10), we thus expect J to range between about 2/3 and 1, and β to range between about 1 and 1.5. We check that such an expected range for β is in agreement with stellar models of typical subgiant and RGB stars. During this phase of the evolution, the density contrast between the core and the envelope increases, and thus the value of J increases too (see Fig. 2). According to Eq. (10), the value of β is therefore also expected to increase. This trend in the slope of 𝒩 and 𝒮₁ is clearly observed in stellar models, as shown in Fig. 1.

In the deep layers of the convective region of evolved red giant stars (e.g., model ER), the equation of state is adiabatic and Fig. 1 shows that the Lamb frequency closely behaves as a power law of radius, so that we can write at first approximation

$$\begin{array}{c}\hfill \mathcal{N}\approx 0\mathrm{and}{\mathcal{S}}_1(r)\approx \sigma {\left(\frac{r_2}{r}\right)}^{\beta }·\end{array}$$(11)

The inner turning point associated with the evanescent region is assumed to coincide with the radius at the base of the convective zone, r_b. Although Eq. (11) represents a first-step description of the evanescent region, we note that the hydrostatic equilibrium and continuity equations do not allow 𝒮₁ to rigorously follow a power law of radius in such an adiabatic region (see Appendix B.1.5 for details). In particular, this implies that the same relation as in Eq. (10) does not apply in this region; instead, the link between β, J and γ is ruled by a more complicated differential equation. Using typical stellar models of evolved red giant stars with Type-b evanescent regions, we find that the value of β lies in a confident range going from 1.2 to 1.5 over the evanescent region, and that the value of J ranges between 0.8 and 1.

In the following, the approximations made in Eqs. (8) and (10) or in Eq. (11) enable us to analytically express q in Type-a or Type-b evanescent regions, respectively, while accounting for the relation between these three variables in a consistent way. Before going further, we first need to introduce in more detail the coupling factor of mixed modes.

3. Theoretical background about the coupling factor

In this section, the physical meaning of the coupling factor of mixed modes is introduced from a general point of view and discussed by means of an instructive toy model. Its general expressions in the asymptotic limit are subsequently presented.

3.1. Physical definition

The eigenfrequencies of mixed modes are selected by the boundary conditions that are imposed at the center and surface, resulting in a discrete spectrum. The asymptotic resonance condition was first provided by the analyses of Shibahashi (1979) and Tassoul (1980) in the hypothesis of a very thick evanescent coupling region. This expression was rewritten in an observational context by Mosser et al. (2012b) in the form of

$$\begin{array}{c}\hfill cot{\mathrm{\Theta }}_{\mathrm{g}}tan{\mathrm{\Theta }}_{\mathrm{p}}=q,\end{array}$$(12)

where Θ_g and Θ_p are frequency-dependent phase terms associated with the propagation of waves in the inner and external cavities, respectively, and q is the so-called coupling factor. The observed frequency pattern of mixed modes in stars is thus partly determined by the value of the coupling factor (e.g., see Mosser et al. 2012b, 2017). Using basic wave principles and the WKB approximation, Takata (2016a,b) later demonstrated that the form of the resonance condition in Eq. (12) actually holds true whatever the thickness of the evanescent region. The author showed in general terms that Θ_g and Θ_p are both equal to the sum of the wavenumber integral over each cavity and of the phase lags introduced at the wave reflection near their boundaries (i.e., the associated turning points), and that q is expressed by

$$\begin{array}{c}\hfill q=\frac{1-\sqrt{1-T^2}}{1+\sqrt{1-T^2}},\end{array}$$(13)

where T² is the transmission coefficient of the wave energy through the intermediate evanescent zone. The value of q must therefore depend on the properties of this region.

From Eqs. (12) and (13), we understand that the coupling factor measures the degree of interaction between both cavities, with values between 0 and 1. For example, in the limit of q → 0 (T² → 0), the cavities are independent: from Eq. (12), either cos Θ_g = 0 or sin Θ_p = 0. These conditions are equivalent to the Borh-Sommerfeld’s quantization condition for a gravity mode trapped in the inner cavity¹ or for a pressure mode trapped in the envelope, respectively. In contrast, in the limit of q → 1 (T² →  1), there is no evanescent zone but a unique resonant cavity. In the intermediate case, Takata (2016b) showed that the value of q rules the squared ratio of the mode amplitude in the acoustic cavity to the mode amplitude in the buoyancy cavity. Actually, Goupil et al. (2013, cf. constants a and c in Appendix A.3 of this paper) demonstrated that this ratio is equivalent to the ratio of the mean time-averaged mode kinetic energy over one wavelength in the acoustic cavity, denoted with ℰ_p, to that in the buoyancy cavity, denoted with ℰ_g. From Eq. (22) of Takata (2016b), it can be rewritten

$$\begin{array}{c}\hfill \frac{{\mathcal{E}}_{\mathrm{p}}}{{\mathcal{E}}_{\mathrm{g}}}=\frac{sin(2{\mathrm{\Theta }}_{\mathrm{g}})}{sin(2{\mathrm{\Theta }}_{\mathrm{p}})}=\frac{1}{q}\frac{{cos}^2({\mathrm{\Theta }}_{\mathrm{g}})}{{cos}^2({\mathrm{\Theta }}_{\mathrm{p}})}=q\frac{{sin}^2({\mathrm{\Theta }}_{\mathrm{g}})}{{sin}^2({\mathrm{\Theta }}_{\mathrm{p}})},\end{array}$$(14)

where both last equations result from the use of Eq. (12). For gravity-dominated mixed mode, which are mixed modes with a gravity-like behavior such as cos Θ_g ≈ 0 and cos Θ_p ≈ 0, ℰ_g in the inner cavity is higher than ℰ_p in the envelope by a factor of 1/q. Inversely, for pressure-dominated mixed modes, which are mixed modes with a pressure-like behavior such as sin Θ_p ≈ 0 and sin Θ_g ≈ 0, ℰ_g is smaller than ℰ_p by a factor of q. In other cases depending on the mode frequency, the ratio vary between these two extrema, that is, q ≤ (ℰ_p/ℰ_g)≤1/q. From this latter inequality, it is obvious that the coupling factor determines how the mode energy is distributed between the inner cavity and the envelope, and thus the character of the mode. Indeed, small values of q preferentially result in eigenmodes with a dominant character (i.e., pressure-dominated or gravity-dominated modes), while values of q close to unity result in eigenmodes with an energy more uniformly distributed between both cavities and thus a pronounced pressure-gravity duality. The value of the coupling factor therefore plays an important role in the detectability of the mixed mode frequency pattern (e.g., Grosjean et al. 2014; Mosser et al. 2018). In order to examine the potential of the coupling factor as a seismic diagnosis, we need to understand how the value of q depends on the properties of the evanescent region.

3.2. An instructive toy model: a parabolic evanescent zone

Equation (13) shows that q depends only on the wave transmission factor through the intermediate evanescent zone. While this formulation is very simple, it is much less easy in practice to grasp how an incident wave behaves as it encounters such a barrier during its travel and how the wave transmission depends on the structure of the barrier. As a first step, it is therefore instructive to study the wave transmission-reflection mechanism in simplified cases that are analytically tractable and can give us qualitative clues on this physical process. The most simple situation is to consider a parabolic evanescent region: one assumes that the squared local radial wavenumber is represented by a second-degree polynomial. This is actually similar to study the propagation of a particle in an inverted harmonic potential in quantum physics (e.g., Heim et al. 2013).

In such a toy model, the wave equation for a given wave function Ψ reads

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathrm{\Psi }}{\mathrm{d}{r}^2}+k_{\mathrm{r}}^2\mathrm{\Psi }=0\mathrm{with}{k}_{\mathrm{r}}^2=\left|k_{\mathrm{r}}{(r_{\mathrm{ev}})}^2\right|[\frac{4{(r-r_{\mathrm{ev}})}^2}{\mathrm{\Delta }{r}^2}-1],\end{array}$$(15)

where r_ev is the radius at the middle of the evanescent zone and Δr is its radial extent. The analytical derivation of the wave energy transmission coefficient through such an evanescent zone is detailed in Appendix C following the work of Phinney (1970). Its final expression reads

$$\begin{array}{c}\hfill T^2=\frac{1}{1+e^{2\mathcal{I}}},\end{array}$$(16)

where ℐ corresponds to the wavenumber integral over the evanescent region, the general definition of which is provided by

$$\begin{array}{c}\hfill \mathcal{I}\equiv {\displaystyle {\int }_{r_{\mathrm{ev}}-\mathrm{\Delta }r/2}^{r_{\mathrm{ev}}+\mathrm{\Delta }r/2}}|k_r|\mathrm{d}r.\end{array}$$(17)

The lower and upper bounds of the integral in Eq. (17) correspond to the turning points r₁ and r₂ defined in Sect. 2.2, respectively. In the case of a parabolic evanescent zone, we find according to Eq. (15) that

$$\begin{array}{c}\hfill \mathcal{I}=\frac{\pi }{4}|k_r(r_{\mathrm{ev}})|\mathrm{\Delta }r.\end{array}$$(18)

Equation (18) shows that ℐ is a proxy for the ratio of the thickness of the evanescent region to the characteristic length scale of decay of the wave amplitude in this region. In this toy model, the wave transmission coefficient T² in Eq. (16) and thus the associated mode coupling factor q computed through Eq. (13) depend only on this ratio (see Fig. C.1 for an illustration).

First, in the limiting case of a very large ratio ℐ ≫ 1, we have T² → 0 and q → 0, so that the resonant cavities on both sides of the evanescent region are independent, as expected from Sect. 3.1. Then, in the intermediate case, the wave transmission increases as ℐ decreases, so does the coupling factor. Finally, in the reverse limiting case of ℐ = 0, Δr = 0 and the incident wave energy is as much reflected as transmitted since T² = 1/2 from Eq. (16) and R² = 1 − T² = 1/2 by energy conservation. This peculiar value of 1/2 is mainly due to the symmetry with respect to r = r_ev of a parabolic evanescent barrier. Such a result can seem counterintuitive at first, since we can naively expect the wave transmission to be total when Δr = 0. It is actually not true since the incident wave still encounters a (double) turning point at r = r_ev and is still reflected. The associated mode coupling factor computed using Eq. (13) is equal to about q (Δr = 0)≈0.17, which remains weak.

The previous example highlights two general points: reflection occurs as soon as the wave encounters a turning point, even if the thickness of the evanescent zone is null; a very thin evanescent zone is not necessarily associated with a strong coupling. It also shows that in general, the coupling factor is likely to depend on the wavenumber integral over the evanescent region and so on its radial extent. By analogy with this example, an evanescent zone is referred in the following to as thick when the wavenumber integral over its extent is much larger than unity, and it is referred to as thin when the wavenumber integral over its extent is equal to or smaller than unity.

3.3. Asymptotic expression of the coupling factor

The general expression of the coupling factor in stellar interiors were obtained by means of asymptotic analyses assuming only two limiting cases: either the case of a thick (i.e., ℐ ≫ 1) or a very thin (i.e., ℐ ≪ 1) evanescent zone. In this paragraph, we introduce both formalisms that will form the theoretical basis of the study performed in this work.

3.3.1. Hypothesis of a thick evanescent zone

Shibahashi (1979) first derived the resonance condition for mixed modes in Eq. (12) under the hypothesis of a thick evanescent zone. In this case, the coupling factor is equal to

$$\begin{array}{c}\hfill q\approx \frac{T^2}{4}\end{array}$$(19)

and

$$\begin{array}{c}\hfill T^2=exp(-2\mathcal{I}),\end{array}$$(20)

with ℐ defined in Eq. (17). Equation (19) is consistent with the general definition of the coupling factor given in Eq. (13) in the limit of a low transmission coefficient (i.e., T² ≪ 1). We also note that Eq. (16) is consistent with Eq. (20) in the limit of a thick evanescent zone (i.e., ℐ ≫ 1).

Shibahashi (1979) originally derived Eqs. (19) and (20) using the Cowling approximation. We demonstrate in Appendix A.1 that the Cowling approximation is valid at leading order either where the short-wavelength WKB approximation is met or where the density contrast compared to the inner bulk is very high, that is, where J ∼ 1. In the buoyancy and acoustic cavities, the WKB approximation is met. As a consequence, the error introduced by the Cowling approximation on the asymptotic period spacing or frequency large separation of mixed modes, which are leading-order terms in the asymptotic expansion, is thus negligible. In contrast, the error introduced on the asymptotic gravity offset of mixed modes, which corresponds to higher-order terms, is on the order of these parameters and thus cannot be neglected (e.g., Takata 2016a; Pinçon et al. 2019). Regarding the coupling factor, the WKB approximation is not valid close to the turning-points r₁ and r₂ associated with the evanescent region; the small deviation of J from unity in this region turns out to be sufficient to modify the thickness of the evanescent region compared to the Cowling approximation, in particular in Type-a evanescent regions of subgiant and young red giant stars (see Fig. 1 and Appendix A.2). This explains the discrepancy observed between the values of the coupling factor computed with and without the Cowling approximation by Jiang & Christensen-Dalsgaard (2014) and Mosser et al. (2017). It appears thus necessary to consider the non-Cowling case when studying the coupling factor of dipolar modes. Actually, to transpose the result of Shibahashi (1979) to the non-Cowling case, we have seen in Sect. 2.1 that we simply need to consider the modified Brunt–Väisälä and Lamb frequencies instead of their usual definitions inside the wavenumber integral ℐ.

Whether the Cowling approximation is used or not, Eq. (19) shows that the value of the coupling factor cannot exceed 1/4 in the hypothesis of a thick evanescent region, in disagreement with observations (e.g., Mosser et al. 2017) and stellar models (Jiang & Christensen-Dalsgaard 2014; Hekker et al. 2018) for subgiant and red clump stars.

3.3.2. Hypothesis of a very thin evanescent zone

When the evanescent zone is very thin, the asymptotic development of the wave functions inside the evanescent zone, used by Shibahashi (1979), is not valid. Given the observational need to go beyond the hypothesis of a thick evanescent zone, Takata (2016b) proposed a new asymptotic description of mixed modes in the reverse limiting case of a very thin evanescent zone. Solving the oscillation equations in this case while accounting for the perturbation of the gravitational potential, he found that the expression of the wave energy transmission coefficient reads

$$\begin{array}{c}\hfill T^2=exp(-2\mathcal{I}-{\mathcal{G}}_0).\end{array}$$(21)

In this formulation, T² depends not only on the wavenumber integral ℐ, but also on a new term denoted with 𝒢₀, which is related to the gradients of the equilibrium quantities inside the evanescent zone. The detailed analytical expression of this extra term is provided in Appendix D.1. The coupling factor obtained through Eq. (13) can therefore vary between zero and unity, depending on the values of ℐ and 𝒢₀. Using typical stellar models, Mosser et al. (2017) showed that this formulation can reproduce the global trend in the observed values of the coupling factor between the subgiant branch and the RGB. However, as said in Sect. 1, their interpretation in terms of internal structure called for a more detailed study. In this work, we propose to further examine the dependence of q with the properties of the evanescent region within the asymptotic formulation in order to help interpret the coupling factor measurements.

4. Relation between q and the properties of the evanescent zone

In this section, we investigate the relation that exists between the coupling factor of mixed modes and the evanescent region in an analytical way within the asymptotic framework. The model assumes that the profiles of the Brunt–Väisälä and Lamb frequencies follow the power-law behavior in Eqs. (8) and (11) for Type-a and Type-b evanescent regions, respectively. This enables us to analytically express q as a function of a few parameters characterizing the evanescent zone. The potential of q to provide information on this region is subsequently analyzed. For sake of simplicity, we focus on stars on the subgiant and red giant branches; the case of red clump stars is addressed later in Sect. 5.

4.1. Analytical expression of q

To clarify the following computations, all the assumptions about the structure of the evanescent region are summarized in Table 1. A schematic view of the considered configuration is provided in Fig. 4.

4.1.1. Type-a evanescent zones

The values of q measured by Mosser et al. (2017) range between 0.15 and 0.60 for subgiant stars, and between about 0.2 and 0.3 for young red giant stars. The model SG in Fig. 1 and the high observed values in the subgiant phase suggest to consider the asymptotic analysis of Takata (2016a) at this stage. In contrast, the observed medium coupling values slightly higher or lower than 1/4, which is the upper limit predicted in the formalism of Shibahashi (1979), makes ambiguous the choice of the proper formalism to use, since the available asymptotic analyses consider only the limiting cases of thick and very thin evanescent zones. Given these theoretical limitations, we thus consider both formalisms and compare them for Type-a evanescent regions. In this framework, the coupling factor of mixed modes is deduced from Eqs. (13) and (21) for very thin evanescent regions, and Eqs. (19) and (20) for thick evanescent regions. To obtain its expression, we thus need to compute the wavenumber integral ℐ and the gradient-related term 𝒢₀.

Firstly, the wavenumber integral in Type-a evanescent regions, denoted with ℐ^r, reads using Eqs. (1)–(3)

$$\begin{array}{cc}\hfill {\mathcal{I}}^r& ={\displaystyle {\int }_{r_1}^{r_2}}|k_r|\mathrm{d}r={\displaystyle {\int }_{r_1}^{r_2}}\frac{\sigma }{c}{(1-\frac{{\mathcal{N}}^2}{{\sigma }^2})}^{1/2}{(\frac{{\mathcal{S}}_1^2}{{\sigma }^2}-1)}^{1/2}\mathrm{d}r\hfill \\ \hfill & ={\displaystyle {\int }_{r_1}^{r_2}}\sqrt{2}J\frac{\sigma }{{\mathcal{S}}_1}{(1-\frac{{\mathcal{N}}^2}{{\mathcal{S}}_1^2}\frac{{\mathcal{S}}_1^2}{{\sigma }^2})}^{1/2}{(\frac{{\mathcal{S}}_1^2}{{\sigma }^2}-1)}^{1/2}\frac{\mathrm{d}r}{r}·\hfill \end{array}$$(22)

For sake of convenience, we define in the following the parameter α such as

$$\begin{array}{c}\hfill \alpha \equiv {\left(\frac{r_1}{r_2}\right)}^{\beta }.\end{array}$$(23)

The power-law assumption in Eq. (8) shows that 𝒩(r)/𝒮₁(r) = α in Type-a evanescent regions. This ratio is thus constant with respect to radius since r₁ and r₂ depend only on the oscillation frequency. Using the change of variable u = σ/𝒮₁ in Eq. (22), the fact that u(r₁) = σ/𝒮₁(r₁) = 𝒩(r₁)/𝒮₁(r₁) = α and u(r₂) = σ/𝒮₁(r₂) = 1 in the bounds of the integral, as well as the equality β = 3J/2, the wavenumber integral becomes (e.g., Gradshteyn et al. 2007)

$$\begin{array}{cc}\hfill {\mathcal{I}}^r& ={\displaystyle {\int }_{\alpha }^1}\frac{\sqrt{2}J}{\beta }\sqrt{1-u^2}\sqrt{u^2-{\alpha }^2}\frac{\mathrm{d}u}{u^2}\hfill \\ \hfill & =\frac{2\sqrt{2}}{3}[({\alpha }^2+1)K(1-{\alpha }^2)-2E(1-{\alpha }^2)],\hfill \end{array}$$(24)

where K are E are the complete elliptic integrals of the first and second kind, respectively, in agreement with the result of Takata (2016a). In the case of a very thin evanescent region, α → 1 and ℐ^r → 0. As α decreases, ℐ^r increases since the evanescent region becomes thicker and thicker, as expected. When the evanescent zone becomes very thick, α ≪ 1 and we note for later purpose that the wavenumber integral in Eq. (24) tends to

$$\begin{array}{c}\hfill {\mathcal{I}}^r\approx \frac{2\sqrt{2}}{3}[-ln\alpha +ln4-2+°(\alpha )],\end{array}$$(25)

where the term °(α) indicates negligible higher-order terms compared to α.

Secondly, the theoretical expression of the gradient-term in very thin evanescent regions, 𝒢₀, which was obtained by Takata (2016a), is detailed in Appendix D.1. Contrary to Mosser et al. (2017), we take the dependence of 𝒢₀ on the thickness of the evanescent region into account. In the simple model represented by Eq. (8), this term is equal for Type-a evanescent regions to (see Appendix D.2)

$$\begin{array}{c}\hfill {\mathcal{G}}_0=\frac{3\pi }{4\sqrt{2}}\frac{ln\alpha }{(\alpha -1)}{(\frac{1}{ln\alpha }+\frac{\alpha }{1-\alpha }-\frac{1}{3}+\sqrt{\frac{9}{4{\beta }^2}-\frac{8{\alpha }^2}{9}})}^2.\end{array}$$(26)

In the limit case of a null thickness, that is, when α → 1, we retrieve the expression derived by Takata (2016a) and considered in Mosser et al. (2017) (see Eq. (D.15) for instance). The predicted value of 𝒢₀ computed with Eq. (26) is shown in Fig. 5 as a function of the parameters α and β. In this figure, we see that in the case of very thin evanescent regions (α ∼ 1), the higher β, the higher 𝒢₀. Indeed, the steeper the variation of the structure, the smaller the wave transmission through the evanescent barrier; therefore, according to Eq. (21) with ℐ ≈ 0, the higher 𝒢₀, as expected. As α slightly decreases from unity, 𝒢₀ decreases, unlike the wavenumber integral. This decrease in 𝒢₀ for α ≲ 1 can be understood as the fact that incident waves feel a smoother evanescent barrier as this latter becomes slightly larger (relatively to the local wavelength).

Although it was derived in the hypothesis of a very thin evanescent region, it is interesting to discuss the extrapolation of the expression of 𝒢₀ toward small values of α (i.e., for thick evanescent regions, which is outside of its derivation domain; see dashed lines in Fig. 5). We see that below a certain value of α where 𝒢₀ vanishes, 𝒢₀ diverges toward +∞ as α tends to zero. As shown by Eq. (D.16), the extrapolation of 𝒢₀ diverges as – lnα when α → 0 over the considered range of values for β. According to Eq. (25), it is on the order of magnitude of the wavenumber integral given in Eq. (25). This prediction differs from the expectations of Takata (2016a) who showed that, in general, the gradient-related term should vanish compared to the wavenumber integral when the evanescent region becomes large, because it represents a higher-order term in the asymptotic expansion. This example therefore emphasizes that the expression of 𝒢₀ in Eq. (D.1), from which is deduced Eq. (26), is an approximation around α ≈ 1 of a more general expression of the gradient-related term that is valid whatever the thickness of the evanescent zone. In other words, the expression of the coupling factor computed using Eq. (D.1) does not converge toward the expression provided by Shibahashi (1979) (i.e., T² in Eq. (21) does not converge toward T² in Eq. (20) as α → 0). This fact reinforces the ambiguity in the formalism to use to interpret medium values of the coupling factor, in the intermediate case between a thick and very thin evanescent region. This point is discussed in Sect. 6.3.

4.1.2. Type-b evanescent zones

In typical evolved red giant stars with Type-b evanescent regions (i.e., for Δν ≲ 10 μHz), the shape of the evanescent region (e.g., model ER in Fig. 1) suggests to consider the weak coupling formalism of Shibahashi (1979); this choice is supported by the small values of q close or lower than about 0.1 inferred observationally by Mosser et al. (2017) and in stellar models by Hekker et al. (2018). We thus need to compute only the wavenumber integral in this case according to Eqs. (19) and (20).

Using Eq. (11) and the change of variable u = σ/𝒮₁, the wavenumber integral, denoted with ℐ^c in this case, reads

$$\begin{array}{cc}\hfill {\mathcal{I}}^c& ={\displaystyle {\int }_{r_{\mathrm{b}}}^{r_2}}\sqrt{2}J{(1-\frac{{\sigma }^2}{{\mathcal{S}}_1^2})}^{1/2}\frac{\mathrm{d}r}{r}\hfill \\ \hfill & ={\displaystyle {\int }_{\alpha }^1}\frac{\sqrt{2}J}{\beta }\sqrt{1-u^2}\frac{\mathrm{d}u}{u},\hfill \end{array}$$(27)

where we used at the boundaries of the integral the fact that u(r₂) = 1 and u(r₁) = σ/𝒮₁(r₁) = α according to Eqs. (11) and (23). Since α = σ/𝒮₁(r_b)≪1 in general for Type-b evanescent regions (e.g., see model ER in Fig. 1), the integrand in Eq. (27) is much larger near its lower bound than near its upper bound. Moreover, since the variation of J over the evanescent region remains relatively small, the error made in the integral remains negligible if we assume that the value of J is about equal to its value at the base of the convective zone. As a result, Eq. (27) can be approximated by (e.g., Gradshteyn et al. 2007)

$$\begin{array}{cc}\hfill {\mathcal{I}}^c& \approx \frac{\sqrt{2}J(r_{\mathrm{b}})}{\beta }{\displaystyle {\int }_{\alpha }^1}\sqrt{1-u^2}\frac{\mathrm{d}u}{u}\hfill \\ \hfill & \approx \frac{\sqrt{2}}{\beta }[ln(\sqrt{1-{\alpha }^2}+1)-ln\alpha -\sqrt{1-{\alpha }^2}],\hfill \end{array}$$(28)

where we also assumed J(r_b)≈1 for evolved red giant stars, as shown for the model ER in Fig. 2. As expected, we find that ℐ_c increases as α decreases. When α ≪ 1 in the case of a thick Type-b evanescent zone, we note that the wavenumber integral in Eq. (28) tends to

$$\begin{array}{c}\hfill {\mathcal{I}}^c\approx \frac{\sqrt{2}}{\beta }[-ln\alpha +ln2-1+°(\alpha )].\end{array}$$(29)

4.1.3. Degeneracy with the structural parameters

According to Eqs. (24), (26) and (28), the coupling factor of mixed modes depends on both parameters α and β. In the framework of the model, the knowledge of this couple of parameters is necessary and sufficient to fully characterize the evanescent region, which is basically represented by a certain thickness and a certain steepness of the profiles of 𝒩 and 𝒮₁.

On the one hand, the α parameter defined in Eq. (23) can actually be rewritten in the form of

$$\begin{array}{c}\hfill \alpha ={\left(\frac{2-\mathrm{\Delta }r/r_{\mathrm{ev}}}{2+\mathrm{\Delta }r/r_{\mathrm{ev}}}\right)}^{\beta },\end{array}$$(30)

where Δr = |r₂ − r₁| and r_ev = (r₁ + r₂)/2. Equivalently, we can also write

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }r}{r_{\mathrm{ev}}}=2\left(\frac{1-{\alpha }^{1/\beta }}{1+{\alpha }^{1/\beta }}\right)·\end{array}$$(31)

The α parameter is thus a proxy for the thickness of the evanescent region.

On the other hand, the β parameter represents by definition the steepness of the profiles of 𝒩 and 𝒮₁ over the evanescent region. First, in Type-a evanescent zones, the high density contrast compared to the helium core implies β ≈ 3J/2 over the region, which, included in Eq. (B.2) with J regarded as locally uniform, leads to

$$\begin{array}{c}\hfill {\beta }_{\mathrm{a}}\approx -\frac{1}{2}{\left(\frac{\mathrm{d}ln\rho }{\mathrm{d}lnr}\right)}_{r_{\mathrm{ev}}}\approx \frac{1}{2}\frac{r_{\mathrm{ev}}}{H_{\rho }(r_{\mathrm{ev}})},\end{array}$$(32)

where H_ρ = −(dlnρ/dr)⁻¹ is the local density scale height. Second, in Type-b evanescent regions, which is located in the convective region, the equation of state is isentropic, so that p ∝ ρ^Γ₁. As a result, β can be expressed over the evanescent region as

$$\begin{array}{cc}\hfill {\beta }_{\mathrm{b}}& \approx -{\left(\frac{\mathrm{d}ln{\mathcal{S}}_1}{\mathrm{d}lnr}\right)}_{r_{\mathrm{ev}}}=-\frac{1}{2}{(\frac{\mathrm{d}lnp}{\mathrm{d}lnr}-\frac{\mathrm{d}ln\rho }{\mathrm{d}lnr})}_{r_{\mathrm{ev}}}+1-{\left(\frac{\mathrm{d}lnJ}{\mathrm{d}lnr}\right)}_{r_{\mathrm{ev}}}\hfill \\ \hfill & \approx \frac{1}{2}({\mathrm{\Gamma }}_1-1)\frac{r_{\mathrm{ev}}}{H_{\rho }(r_{\mathrm{ev}})}+1,\hfill \end{array}$$(33)

where (dlnJ/dlnr)_rev is neglected at first approximation. Therefore, β turns out to be directly related to the local density scale height, H_ρ.

To conclude, the model shows us that the coupling factor of mixed modes depends on the thickness of the evanescent zone represented by Δr/r_ev, and on the local density scale height represented by β. The value of q is therefore in general degenerate with respect to these two internal properties whatever the type of the evanescent region.

4.2. Variations of q with respect to the structure

As shown in the previous section, the link between q and the structural parameters β and Δr/r_ev depends on the type of the evanescent region. Both types of evanescent region are thus taken into account and compared in the present section. Given the ambiguity in the applicability of the weak or the strong coupling formalisms regarding the thickness of the evanescent region (see Sect. 4.1.1), we compute and compare both expressions of q as provided by Eqs. (20) and (21) for a large range of values of Δr/r_ev. Assumptions on the applicability domain of each formalism will be made in a second step for the interpretation of the observations (see Sect. 5).

4.2.1. Type-a evanescent regions

The prediction of the model in Type-a evanescent regions is plotted in Fig. 6. In this figure, the coupling factor was computed for typical values of β using either Eqs. (19), (20) and (24) in the formalism of Shibahashi (1979), or Eqs. (13), (21), (24) and (26) in the formalism of Takata (2016a), as well as Eq. (31). The values of q were taken in the range between 0.02 and 1, which encapsulates the observed range of values (e.g., Mosser et al. 2017).

In the formalism of Shibahashi (1979), Fig. 6 shows that the lower the value of Δr/r_ev or β, the higher the value of q. Indeed, using for instance the approximation in Eq. (25) in the limit of a very thick evanescent region, the coupling factor can be rewritten according to Eqs. (19) and (20) such as

$$\begin{array}{c}\hfill q_{\mathrm{a}}\approx 0.8{\alpha }^{4\sqrt{2}/3}.\end{array}$$(34)

Since α decreases as Δr/r_ev or β increases according to Eq. (30), Eq. (34) confirms that q must also decrease. This trend originates from the fact that the wave transmission factor through the evanescent region is lower either when this region is larger or when the length scale of decay of the wave amplitude is smaller (i.e., when β is higher). We also see in Fig. 6 that at a given value of Δr/r_ev, the influence of β is negligible for q ∼ 0.20, but significant for q ∼ 0.05 (with a change in q of about 60% between β = 1.5 and β = 1.2). The dependence of q on β, which comes from the α parameter in Eq. (34), thus cannot be neglected.

In contrast, in the formalism of Takata (2016a), Fig. 6 shows that q depends on Δr/r_ev and β following two different regimes. At low values of Δr/r_ev, that is, between about 0 and 0.15, the predicted values of q are high, comprised between about 0.25 and 1. The lower β, the higher q, but the lower Δr/r_ev, the lower q. This difference of trend compared to the formalism of Shibahashi (1979) results from the predominant effect of the gradient-related term 𝒢₀ in this case (see Fig. 5). This term also results in a more important sensitivity of q to the parameter β in this regime. Indeed, at a given Δr/r_ev, a change in β between 1.2 and 1.5 results in a large change of q by a factor of between two and four. At values of Δr/r_ev ≳ 0.15, Fig. 6 shows that the higher Δr/r_ev, the lower q, but no substantial dependence on β can be depicted. This results again to the gradient-related term 𝒢₀ whose value is large enough to counterbalance the effect of the wavenumber integral (see second paragraph in Sect. 4.1.1).

In summary, the behavior of q with respect to the properties of the evanescent region is substantially different depending on the formalism we consider. In the formalism of Shibahashi (1979), q decreases as Δr/r_ev or β increases, with a substantial dependence on β. In the formalism of Takata (2016a), the effect of the gradient-related term 𝒢₀ is important and drastically modifies how q varies with respect to Δr/r_ev and β. In particular, for Δr/r_ev ≲ 0.15, q turns out to increase as Δr/r_ev increases and to be very sensitive to β over the typical considered range of values.

4.2.2. Type-b evanescent regions

The predictions of the model in Type-b evanescent regions are shown in Fig. 7. In this figure, q was computed for typical values of β using Eqs. (19), (20) and (28) in the formalism of Shibahashi (1979) as well as Eq. (31). The values of q were chosen between 0.02 and its upper limit equal to 0.25.

As seen in Fig. 7, the higher Δr/r_ev or β, the lower q for the same reasons as in Type-a evanescent regions in the thick limiting case. However, Fig. 7 first shows that at a given value of Δr/r_ev, the value of q obtained in Type-b evanescent zones is about twice lower than in Type-a evanescent regions, and much less sensitive to β. To explain this difference, it is convenient to use Eq. (29) in the limit of a very thick evanescent region; the coupling factor can thus be expressed according to Eqs. (19), (20) and (31) such as

$$\begin{array}{c}\hfill q_{\mathrm{b}}\approx \frac{2.{38}^{2\sqrt{2}/\beta }}{4}{\alpha }^{2\sqrt{2}/\beta }=\frac{2.{38}^{2\sqrt{2}/\beta }}{4}{\left(\frac{2-\mathrm{\Delta }r/r_{\mathrm{ev}}}{2+\mathrm{\Delta }r/r_{\mathrm{ev}}}\right)}^{2\sqrt{2}}.\end{array}$$(35)

The first multiplicative factor in Eq. (35) only slightly varies between 0.51 and 0.45 (i.e., lower than 10%) as β changes from 1.2 and 1.5. Therefore, the comparison between Eqs. (34) and (35) confirms the difference depicted in Fig. 7 between Type-a and Type-b evanescent regions. This is the consequence of two facts. First, N² ≈ 0 in convective regions so that the wavenumber integral over the evanescent region is higher and thus the wave transmission factor is lower; second, the linear relation between β and J is not valid in this case.

4.3. Transition from Type-a to Type-b evanescent regions

As seen in the previous section, the value of q is about twice lower in a Type-a than in a Type-b evanescent region for the same values of β and Δr/r_ev. At first sight, the transition from a Type-a to a Type-b evanescent region during evolution could thus be expected to be associated with a rapid variation in the value of q measured around ν_max. Nevertheless, the transition actually does not occur at a given instant (i.e., a given value of ν_max), but is composed of an intermediate phase where the evanescent region progressively changes from a Type-a to a Type-b regime. Indeed, as stars evolve, the Type-a regime ends as soon as r₂ becomes higher than r_b, while r₁ still remains lower than r_b. This happens when 𝒩_b ≲ 2πν ≲ 𝒮_b, where 𝒮_b = 𝒮₁(r_b) is the value of the Lamb frequency at the base of the convective region. During this transition phase, a Type-a configuration is met for r₁ <  r <  r_b whereas a Type-b configuration is met for r_b <  r <  r₂. The Type-b regime then starts when r₁ ≈ r_b, that is, when 2πν ≲ 𝒩_b (see Table 1 and Fig. 4).

To estimate the effect of this progressive transition on the coupling factor in the framework of the analytical model, the wavenumber integral over the evanescent region when r₁ <  r_b <  r₂, denoted with ℐ^t, has to be computed. It can be easily obtained by combining Eqs. (24) and (28) while adapting the upper bound of Eq. (24), that is, such as

$$\begin{array}{c}\hfill {\mathcal{I}}^{\mathrm{t}}\approx \frac{2\sqrt{2}}{3}{\displaystyle {\int }_{{\alpha }_{\mathrm{b}}}^{u_{\mathrm{b}}}}\sqrt{1-u^2}\sqrt{u^2-{\alpha }_{\mathrm{b}}^2}\frac{\mathrm{d}u}{u^2}+\frac{\sqrt{2}}{\beta }{\displaystyle {\int }_{u_{\mathrm{b}}}^1}\sqrt{1-u^2}\frac{\mathrm{d}u}{u},\end{array}$$(36)

with

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{b}}=\frac{{\mathcal{N}}_{\mathrm{b}}}{{\mathcal{S}}_{\mathrm{b}}}\mathrm{and}{u}_{\mathrm{b}}=\frac{\sigma }{{\mathcal{S}}_{\mathrm{b}}}={\alpha }_{\mathrm{b}}\frac{\sigma }{{\mathcal{N}}_{\mathrm{b}}},\end{array}$$(37)

where we used in the lower bound of the first integral the fact that u(r₁) = σ/𝒮₁(r₁) = 𝒩(r₁)/𝒮₁(r₁) = 𝒩_b/𝒮_b since this ratio is independent of radius for r ≤ r_b (see Sect. 4.1.1). To go further, the integral in Eq. (36) can be expressed as a function of special functions to finally obtain (e.g., Gradshteyn et al. 2007)

$$\begin{array}{cc}& {\mathcal{I}}^{\mathrm{t}}\approx \frac{\sqrt{2}}{\beta }[ln(\sqrt{1-u_{\mathrm{b}}^2}+1)-ln{u}_{\mathrm{b}}-\sqrt{1-u_{\mathrm{b}}^2}]\hfill \\ \hfill & -\frac{2\sqrt{2}}{3}[\sqrt{(1-u_{\mathrm{b}}^2)(1-\frac{{\alpha }_{\mathrm{b}}^2}{u_{\mathrm{b}}^2})}-2E({\gamma }_{\mathrm{b}}|1-{\alpha }_{\mathrm{b}}^2)+2E(1-{\alpha }_{\mathrm{b}}^2)]\hfill \\ \hfill & -\frac{2\sqrt{2}}{3}(1+{\alpha }_{\mathrm{b}}^2)[F({\gamma }_{\mathrm{b}}|1-{\alpha }_{\mathrm{b}}^2)-K(1-{\alpha }_{\mathrm{b}}^2)],\hfill \end{array}$$(38)

where K(x) are E(x) and the complete elliptic integrals of the first and second kind, F(z | x) and E(z | x) are the incomplete elliptic integrals of the first and second kind of parameter k² = x, and

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{b}}=arcsin\left(\sqrt{\frac{1-u_{\mathrm{b}}^2}{1-{\alpha }_{\mathrm{b}}^2}}\right).\end{array}$$(39)

Additional simplifying assumptions can be made to go a step further. Since the considered red giant stars close to the transition are evolved enough, we first assume that the evanescent region is thick and use the weak coupling formalism of Shibahashi (1979). Second, we assume that β is nearly uniform and equal to 1.5 on both sides of r_b. This is a reasonable assumption since, on the one hand, the density contrast in the vicinity of the base of the convective region compared to the helium core is so high that J ≈ 1 in the radiative region, and thus β ≈ 3J/2 ≈ 1.5 here. On the other hand, the sensitivity of q to β in the convective region is small enough to be neglected at first approximation (see Sect. 4.2.2), so that choosing β ≈ 1.5 has only a small impact on the final result.

Within these considerations, the value of the coupling factor can be computed as a function of 2πν/𝒩_b using Eqs. (19) and (20), and considering for the wavenumber integral:

• ℐ^r in Eq. (24) if ${\alpha }_{\mathrm{b}}^{-1}<2\pi \nu /{\mathcal{N}}_{\mathrm{b}}$;

• ℐ^t in Eq. (38) if $1<2\pi \nu /{\mathcal{N}}_{\mathrm{b}}<{\alpha }_{\mathrm{b}}^{-1}$; and

• ℐ^c in Eq. (28) if 2πν/𝒩_b <  1.

The range of values for α_b is chosen between about 0.4 and 0.6, as required by the model to reproduce typical values of q observed between about 0.1 and 0.2 in stars with 10 μHz ≲ Δν ≲ 18 μHz (Mosser et al. 2017; Hekker et al. 2018), which correspond to stars before the expected transition to the Type-b regime (see Sect. 2.2). The result is plotted in Fig. 8. From an observational point of view, the ratio 2πν/𝒩_b for ν ∼ ν_max generally decreases during the evolution on the RGB. Indeed, the decrease in ν_max induced by the envelope expansion is predominant compared to the variation in 𝒩_b that is mainly induced by variations in the location of the base of the convective region, as shown by the rough estimate ${\mathcal{N}}_{\text{b}}^2\sim G{M}_{\text{r}}(r_{\text{b}})/r_{\text{b}}^3$. The value of this ratio is considered between 0.2 and 4, which is representative of the ratio 2πν_max/𝒩_b in typical models of red giant stars close to the transition and just before the luminosity bump (Pinçon et al. 2019, see their Appendix B for details).

We see in Fig. 8 that q decreases as 2πν/𝒩_b or α_b decreases whatever the regime considered. Indeed, first, in the Type-b regime, we can express the parameter α in Eq. (23) as α = 2πνα_b/𝒩_b according to Eq. (11) with r₁ = r_b. The last equality shows that the decrease in q as 2πν/𝒩_b or α_b decreases actually results from a decrease in α, or equivalently, from an increase in Δr/r_ev, as already shown in Sect. 4.2.2. Second, in the Type-a and transition regimes, we have α = α_b since α does not depend on radius (i.e., when r <  r_b). Therefore, at a fixed value of 2πν/𝒩_b, the decrease in q as α_b (i.e., α) decreases results also from an increase in Δr/r_ev (since β remains close to 1.5). However, in the transition regime, an increase in Δr/r_ev is not the only possible reason for the decrease in q observed in Fig. 8. Indeed, at a fixed value of α_b, α and Δr/r_ev are also fixed, but q still decreases as 2πν/𝒩_b decreases. In fact, this decrease in q results, according to Eqs. (36) and (37), from the progressive migration of the evanescent region from the radiative to the convective zones, because of the simultaneous progressive increase in the contribution of the integral between r_b and r₂ to the total wavenumber integral (and since the wavenumber integral over any given region is always higher when 𝒩² ≈ 0).

The analytical predictions can now be used to address the variation of the observed coupling factor q(ν_max) for ν ∼ ν_max during the evolution on the RGB. During this phase, the parameter α_b = 𝒩_b/𝒮_b can be assumed to vary much more slowly than 2πν_max/𝒩_b, because both 𝒩² and ${\mathcal{S}}_1^2$ similarly evolve as g/r in the radiative region whereas ν_max decreases more rapidly than 𝒩_b varies (see Sect. 2.2). Therefore, as 2πν_max/𝒩_b decreases (from large values to unity) at a quasi constant value of α_b, the model first predicts that q(ν_max) starts smoothly decreasing at the beginning of the transition phase because of the progressive migration of the evanescent region from the radiative to the convective zones. Second, when the evanescent region is entirely located in the convective region for 2πν_max/𝒩_b ≲ 1, q(ν_max) keeps progressively decreasing along with the evolution because Δr/r_ev(ν_max) increases as ν_max decreases.

Therefore, we conclude that no sharp variation in q(ν_max) is expected at the transition from a Type-a to a Type-b evanescent region. Instead, our simplified estimate predicts that q(ν_max) smoothly and progressively decreases at this stage, mainly because of the progressive migration of the evanescent region from the radiative to the convective regions (provided that α_b does not vary too rapidly). This is in contrast with the gravity offset of mixed modes, denoted with ε_g, which is another seismic parameter associated with mixed modes appearing inside the term Θ_g of Eq. (12). This parameter was shown by Pinçon et al. (2019) to exhibit a clear drop at the transition to the Type-b regime because of the kink in the Brunt–Väisälä frequency at the base of convective region. The authors did not notice any obvious simultaneous signature in the observed values of q(ν_max) for such evolved stars; the analytical model developed in the present work confirms the absence of such a sharp variation in the global evolution of q(ν_max).

4.4. Lessons learned in a nutshell

To summarize this first part of the study, the structure of the evanescent coupling region changes during evolution. In subgiant, young red giant and red clump stars, the evanescent region is located in the radiative region, where 𝒩 and 𝒮₁ vary at first approximation as power laws of radius with the same exponent β (Type-a regime). In evolved RGB stars, with typically Δν ≲ 10 μHz, the evanescent region is located in the lower part of the convective region, where 𝒩 ≈ 0 (Type-b regime).

In both cases, the analytical model developed shows that the value of q is degenerate with respect to β and Δr/r_ev. In particular, the model highlights the following points:

– In the Type-a regime, the behavior of q with both parameters significantly changes depending on the thickness of the evanescent region. In thick evanescent regions, q monotonically decreases as β or Δr/r_ev increases. In very thin evanescent regions, q turns out to be very sensitive to β (for a typical range of values), in a much more important way than in thick evanescent regions; moreover, q in contrast increases as Δr/r_ev increases for Δr/r_ev ≲ 0.15. This drastic change of trend originates from the effect of the gradient in the equilibrium structure on the wave transmission through the evanescent region, represented by 𝒢₀ in Eq. (21).

– In Type-b evanescent regions, q also decreases as β or Δr/r_ev increases. The value of q is nevertheless about twice lower than in Type-a regions for the same couple of parameters β and Δr/r_ev. Moreover, its sensitivity to β is very small; the value of q measured around ν_max is thus expected to decrease during the evolution, because of the progressive increase in Δr/r_ev(ν_max) as ν_max decreases.

– Unlike the gravity offset of mixed modes, q(ν_max) is not expected to exhibit a sharp variation at the transition from the Type-a to the Type-b regimes on the RGB, which occurs just before the luminosity bump.

Finally, we emphasize that the trend of q with Δr/r_ev predicted by the analytical model in evolved red giant stars with Type-b evanescent regions is comparable to that found in stellar models by Hekker et al. (2018, see their Fig. 7), at least for values of q lower than about 0.15. However, the result of their fit is questionable for the red clump stars they considered in the strong coupling case since the computation of Hekker et al. (2018) does not account for the gradient-term 𝒢₀, which turns out to play an important role in Type-a evanescent regions, as demonstrated in Sect. 4.2.1.

5. Clues on the probing potential of q

In this section, we investigate the ability of the coupling factor to bring us constraints on stellar interiors. We give a qualitative interpretation of the variations in q observed by Mosser et al. (2017) in the light of the analytical model developed in the previous section. The discussion is illustrated by considering typical values of q and their error bars in subgiant (SG), young red giant (YR) and evolved red giant stars (ER). The case of red clump stars (RC) is also considered. The confrontation of the model to observations implicitly assumes that q is taken at the frequency ν_max for a given star. The possible dependence of q with the frequency over the observed oscillation spectrum is addressed in a subsequent step. The present study represents a preparatory stage toward a more quantitative investigation using stellar models that will be presented in the Paper II of the series.

5.1. Qualitative interpretation of the current measurements

5.1.1. Observations

The large-scale measurement of q as a function of Δν (or ν_max) by Mosser et al. (2017, see their Figs. 5 and 6) shows that stars exhibit low values of q around 0.15 at the very beginning of the subgiant branch before rapidly increasing toward a peak value of about 0.5. During this peak phase on the subgiant branch, the spread in the measurements is large, on the order of about 0.15.

In the following stages for typically 18 μHz≲Δν ≲ 30 μHz, the observed value of q decreases from about 0.30 to 0.18. During the subsequent evolution on the RGB, the coupling factor continues to decrease slowly in average, with a spread in the measurements lower than about 0.03, that is, much smaller than in subgiant stars. The most luminous red giant stars of the sample with Δν ∼ 5 μHz have values of q around 0.10, with an observed minimum value close to 0.05.

Regarding red clump stars, the observed values of q remain around a high mean value of about 0.3. This suggests that these stars are in the regime of very thin (Type-a) evanescent regions. The spread in the measurements is low and close to about 0.02 around the mean value².

5.1.2. Coupling regimes vs evolutionary stage

The analytical expressions of q depend both on the type and on the hypothesis about the thickness of the evanescent region (i.e., thick or very thin). To interpret observations of q in the framework of the model along with evolution, we thus have to make additional assumptions about which type and coupling formalism to consider depending on the evolutionary state and the value of q.

First, we estimate through standard stellar models that evolved red giant stars with Type-b evanescent regions have Δν ≲ 5 and 10 μHz for stellar masses between 2 M_⊙ and 1 M_⊙, respectively (see Sect. 2.2). Mosser et al. (2017) observationally determined that q on the RGB follows the fit law q ≈ q_M(Δν/10)^0.1, in which q_M ∈ [0.125, 0.145] in the considered mass range and Δν is expressed in μHz. According to this relation, we thus consider that (thick) Type-b evanescent regions are met when q ≲ 0.13 in RGB stars.

Just above q ≈ 0.13 on the RGB, stars are in the transition regime. In this regime, q depends not only on Δr/r_ev and β, but also on the ratio 2πν_max/𝒩_b, which represents the influence of the migration of the evanescent region from the radiative to the convective region. Using the theoretical picture presented in Sect. 4.3 (assuming thick evanescent regions), we find that the model requires α_b ≈ 0.55 to explain q ≈ 0.13 at the end of the transition phase (i.e., when 2πν_max/𝒩_b ≈ 1, see Fig. 8). Assuming that α_b = 𝒩_b/𝒮_b is fixed at 0.55 during this phase because of the similar evolution of 𝒩² and ${\mathcal{S}}_1^2$ as g/r in the radiative zone, the model predicts q ≈ 0.18 at the beginning of the transition regime. We thus fix that the Type-a regime is met for q ≳ 0.18 on the RGB.

In young red giant stars, corresponding to RGB stars with q ≳ 0.18 following the previous estimate, as well as in subgiant stars with typically q ≈ 0.15 − 0.65, thick and very thin Type-a evanescent regions are considered and compared owing to the ambiguity in the coupling formalism to use. For red clump stars, we consider very thin (Type-a) evanescent regions since this is the only current formalism able to reproduce the high observed values close to 0.3.

All the assumptions considered in the following are summarized in Table 2. Under these considerations, the link between q and the couple of parameters (Δr/r_ev, β) is represented in Fig. 9. Because of its dependence on 2πν_max/𝒩_b, the relation between the value of q and (Δr/r_ev, β) is not unique in the transition regime; it is thus marginally indicated by a blue rectangle in Fig. 9 for a large range of values for Δr/r_ev between 0.2 and 0.6. We emphasize that the uncertainties on the definition of the applicability domains of the Type-a, transition and Type-b regimes as a function of q on the RGB will have no impact on the qualitative interpretation of the observed variations, since the behavior of q with Δr/r_ev and β follows the same global trend for the three regimes.

5.1.3. Subgiant and young red giant stars

As an illustration, we consider the case of a subgiant star with typically q = 0.50. For subgiant stars, the uncertainties on the measurements of q by Mosser et al. (2017) are lower than 20%; we thus adopt a relative uncertainty of 10% on its value. As shown in Fig. 9, q = 0.50 ± 0.05 is compatible with the ranges of parameters β ∈ [1.2, 1.5] and Δr/r_ev ∈ [0, 0.25]. At such a high value of q, the degeneracy with β and Δr/r_ev leads to a large ambivalence on the determination of the structure. Furthermore, even if β is fixed for instance at 1.3, such value of q can either correspond to Δr/r_ev = 0.22 ± 0.03 or Δr/r_ev = 0.03 ± 0.02, because the non-monotonic behavior of 𝒢₀ with Δr/r_ev. Without any other information than the value of q, it is therefore not possible to determine unequivocally the properties of the evanescent region in subgiant stars through a direct confrontation to the asymptotic analytical expression of q.

Just after the end of the subgiant branch, we now consider as an example the case of a young red giant star with typically q = 0.20. For red giant stars, the automated measurement method of Mosser et al. (2017) leads to large uncertainties on q. In contrast, the case-by-case study of Hekker et al. (2018) in few stars showed that the error bars on q can reach at most 0.02 on the RGB. We thus adopt an optimistic uncertainty of 0.02 on the RGB. Because of the ambiguity on the formalism to use in this case, both are compared. In the formalism of Takata (2016a), the model predicts Δr/r_ev = 0.48 ± 0.04 for q = 0.20 ± 0.02, with a small sensitivity regarding β. In contrast, in the formalism of Shibahashi (1979), the model predicts Δr/r_ev = 0.28 ± 0.07 at β = 1.5 and Δr/r_ev = 0.35 ± 0.08 at β = 1.2; in this case, the sensitivity of q to β cannot be neglected with the current uncertainties on q. However, since the considered configuration lies somewhere between the limiting cases of a thick and of a very thin evanescent region, we thus cannot conclude on which prediction is the most relevant. We emphasize the same issue holds for the young subgiant stars with low observed values of q between about 0.15 and 0.25.

Despite these difficulties related to the degeneracy of q and the ambiguity on the coupling formalism to use, the model can bring us some clues on the structural origins of the global observed variations in q during the evolution on the subgiant branch³. In particular, the increase in q from about 0.15 to a peak value around 0.5 at the beginning of the subgiant branch is necessarily associated with a predominant decrease in Δr/r_ev, from about 0.50 to at least 0.15. The subsequent decrease in q toward values lower than 0.2 at the very beginning of the RGB automatically results in a predominant increase in Δr/r_ev, toward about 0.5. Moreover, the large spread in the values of q around the short peak phase even suggests that Δr/r_ev reaches values lower than 0.15. Indeed, the high sensitivity of q to both parameters as well as its non-monotonic behavior with Δr/r_ev in this (very thin coupling region) regime may result, for a given stellar mass, in high amplitude variations over this short evolutionary phase, even if the structure slowly varies. Besides, for two stars with different stellar masses at a given value of ν_max, small differences in their mid-layer structures may also lead to significantly different values of q. The behavior of 𝒢₀ in very thin evanescent regions thus represents a plausible explanation for the large spread in the values of q measured by Mosser et al. (2017) during the subgiant peak phase.

5.1.4. Transition phase and evolved red giant stars

In the subsequent transition phase, we can consider as in Sect. 4.3 that variations in β are small enough to not affect q. In the framework of the model, the decrease in q can thus be associated with the predominance of either an increase in Δr/r_ev, either the progressive migration of the evanescent region from the radiative to the convective region, or the simultaneous effects of these two ingredients.

In later stages, we can consider as an illustration the case of an evolved red giant stars with typically q = 0.10 ± 0.02. As seen in Fig. 9, such a value of q is compatible with Δr/r_ev = 0.48 ± 0.08 for β = 1.2 or with Δr/r_ev = 0.51 ± 0.08 for β = 1.5. Given the current error bars on q, the sensitivity of q to β is thus small enough to be neglected at this evolutionary stage. Therefore, we confirm the prediction obtained in Sect. 4.3; the slight decrease in q observed by Mosser et al. (2017) in the more evolved red giant stars of the sample results from an increase in Δr/r_ev as ν_max decreases during evolution.

5.1.5. Red clump stars

In the red clump, where the helium fusion takes place in the core, we consider as an example a star with typically q = 0.3 and an optimistic error bar close to 0.02 (e.g., Mosser et al. 2017; Hekker et al. 2018). According to Fig. 9, such a value of q is compatible with Δr/r_ev ∈ [0, 0.25] and β ∈ [1.2, 1.5]. As in subgiant stars, the link with the internal properties remains ambivalent. However, unlike subgiant stars, the small spread in the observed values suggests that the structure of the intermediate evanescent region little varies during the evolution in the red clump; otherwise, the high sensitivity of q to the structure in the regime of very thin evanescent regions would lead to larger variations, as in subgiant stars. This expectation seems reasonable since the red clump corresponds to a stable phase of the evolution where the helium-burning core produces enough energy to counterbalance the weight of the above layers, while the subgiant branch corresponds to a short unstable phase where the core collapses and the envelope dilates rapidly.

5.2. Convective boundary signature in evolved RGB stars

As shown in Sect. 5.1.4, q turns out to be mainly sensitive to the thickness of the evanescent region in evolved red giant stars, and hence to the radius of the base of the convective region. Indeed, the lower r_b, the higher Δr and the lower r_ev, so that the higher Δr/r_ev and the lower q. The measurement of the coupling factor in such stars may therefore enable us to probe the boundary between the radiative and convective zones.

To discuss this point, we can consider as an example the case of a 1.2 M_⊙ evolved red giant star with typically ν_max ≈ 60 μHz, Δν ≈ 6 μHz, and q = 0.05. Such a value of q is compatible with the asymptotic value computed for the standard model ER in Fig. 1, or with the minimum value measured by Mosser et al. (2017) at this evolutionary stage. To go further, we follow the picture used in Sect. 4.3; we assume that α_b = 𝒩_b/𝒮_b remains equal to 0.5 on the considered phase and fix the value of β at 1.5. As a result, q only depends on the ratio 2πν_max/𝒩_b, as represented in Fig. 8. We remind that this ratio is expected to decrease during the evolution on the RGB.

Under these considerations, the analytical model requires 2πν_max/𝒩_b ≈ 0.6 to explain q = 0.05; the knowledge of q thus enables us to infer that 𝒩_b/2π ≈ 100 μHz for this star with ν_max ≈ 0 μHz. Using standard stellar models with the same input physics as described in Sect. 2.1, Pinçon et al. (2019) found that the transition to the Type-b regime occurs at ν_max, t ≈ 95 μHz for 1.2 M_⊙ stars. Since 2πν_max, t ≈ 𝒩_b, t at the transition, we have 𝒩_b/𝒩_b, t ≈ 100/95 ≈ 1.05. Therefore, in this example, the analytical model predicts that 𝒩_b increases of about 5% as ν_max decreases from ν_max, t ≈ 95 μHz to about 60 μHz. Since ${\mathcal{N}}_{\text{b}}^2\sim G{M}_{\text{r}}(r_{\text{b}})/r_{\text{b}}^3$ and hence is mainly sensitive to the variations in the factor $r_{\mathrm{b}}^{-3}$, we deduce that this slight increase in 𝒩_b is directly related to a slight decrease in r_b. This seems reasonable since, according to Pinçon et al. (2019), the transition to the Type-b regime happens just before the luminosity bump, that is, before the radius of the base of the convective region attains its minimum during the first dredge-up of the chemical elements.

However, the inward migration of r_b is not the only possible reason for observing such a value of q at ν_max ≈ 60 μHz. This can also be explained considering the value of 𝒩_b fixed (i.e., r_b fixed). Indeed, for 𝒩_b/2π to be remain about equal to 100 μHz, that is, in agreement with the value inferred according to the analytical model in the present example, we can choose a slightly higher value of ν_max, t equal to about 100 μHz (i.e., such as 𝒩_b, t ≈ 100 μHz). As shown by Pinçon et al. (2019), this increase in ν_max, t, compared to the value estimated through standard stellar models such as these considered as examples in the previous sections (e.g., Figs. 1–3), can be physically produced by the inclusion of an adiabatic overshooting region below the base of the convective zone since this process decreases r_b and thus increases 𝒩_b at a given value of ν_max. Moreover, the study of Khan et al. (2018) supports the need for overshooting in such evolved RGB models in order to match the location of the RGB luminosity bump with seismic observations. To estimate the quantity of overshooting needed to increase ν_max, t (i.e., 𝒩_b, t/2π) from 95 μHz to 100 μHz, we can proceed as follows. First, we emphasize that overshooting does not modify the formal shapes of 𝒩 and 𝒮₁, and thus the analytical model can still be applied. Second, at the base of the convective region of the considered evolved stars, it is reasonable to assume J ≈ 1 (e.g., see Fig. 1). As a consequence, this high density contrast compared to helium core implies (see Appendix B.1.3)

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}lnp}{\mathrm{d}lnr}\right)}_{r_{\mathrm{b}}}=\frac{r_{\mathrm{b}}}{H_{\mathrm{p},\mathrm{b}}}\approx 6J-2\approx 4,\end{array}$$(40)

where H_p, b is the pressure scale height at r_b. Then, assuming that the variation in 𝒩_b mainly results from its dependence on $r_{\mathrm{b}}^{-3/2}$, using Eq. (40) and rewriting the amount of decrease in r_b caused by overshooting as δr_b = −α_ovH_p, b, with α_ov the overshooting parameter, we obtain

$$\begin{array}{c}\hfill \frac{\delta {\nu }_{\max ,\mathrm{t}}}{{\nu }_{\max ,\mathrm{t}}}\approx \frac{\delta {\mathcal{N}}_{\mathrm{b},\mathrm{t}}}{{\mathcal{N}}_{\mathrm{b},\mathrm{t}}}\approx -\frac{3}{2}\frac{\delta r_{\mathrm{b}}}{r_{\mathrm{b}}}=\frac{3{\alpha }_{\mathrm{ov}}}{8}·\end{array}$$(41)

Therefore, with δν_max, t/ν_max, t ≈ 0.05, we find that α_ov ≈ 0.15 is sufficient for the analytical expression to predict q = 0.05 at ν_max = 60 μHz, without considering any variation in 𝒩_b along with evolution. Interestingly, this is close to the overshooting value found by Khan et al. (2018) to fit the luminosity bump to observations in such stars.

In reality, the value of the coupling factor in such evolved red giant stars, that is, for Δν ≲ 10 μHz, is certainly the result of the simultaneous influences of the migration rate of r_b during the first dredge-up and the extent of overshooting. This simple estimate clearly emphasizes the potential of q in such stars to provide constraints on the dynamical processes at the interface between the convective and the radiative regions, as for instance the convective mixing of the chemical elements.

5.3. Frequency-dependence of q over the spectrum

When comparing the analytical model to observations in the previous sections, we implicitly assume that the current measurements of q are representative of the predicted values at the frequency ν_max. Indeed, the usual measurement methods assume q as a constant parameter to fit, independent of the oscillation frequency around ν_max (e.g., Buysschaert et al. 2016; Mosser et al. 2017; Hekker et al. 2018). Nevertheless, previous works highlighted the frequency-dependence of q in stellar models (e.g., Jiang & Christensen-Dalsgaard 2014; Mosser et al. 2017; Hekker et al. 2018). From an observational point of view, q can be considered frequency-independent if its variation over the observed frequency range of a given star is lower than the uncertainties on its value. In this section, we question the relevance of this assumption in the framework of the analytical model. In this model, we see in Sect. 4.1.3 that q depends on the α and β parameters. While β is supposed to be fixed over the evanescent region, α, and thus Δr/r_ev, can depend on the oscillation frequency since it depends by definition on the frequency-dependent turning points.

In subgiant, young red giant and red clump stars with Type-a evanescent regions, it is easy to show that α = 𝒩(r)/𝒮₁(r) according to Eq. (8) and (23). This is the consequence of the assumed parallel variations of 𝒩 and 𝒮₁. Unlike r₁ and r₂, α is thus frequency-independent, so that q is frequency-independent too. The analytical model thus supports the assumption of the frequency-independence of q in such stars (e.g., Mosser et al. 2017).

In contrast, for RGB stars in the transition regime, we show in Sect. 4.3 that q explicitly depends on the oscillation frequency (see Fig. 8). In more evolved RGB stars with Type-b evanescent regions, we obtain α = [r₂(σ)/r_b]^β according to Eq. (23), since r₁(σ)≈r_b. Using Eq. (11), it turns out that α = σ/𝒮_b. The model therefore predicts that the coupling factor depends on the oscillation frequency in such RGB stars. The question is then to know whether its variation is significant over a typical observed frequency range around ν_max or not. To answer this question, we can proceed as in Pinçon et al. (2019) and assume that the observed oscillation spectrum of a given star belongs to a 6Δν-wide frequency range centered around ν_max, as expected from theoretical estimates and observations (e.g., Grosjean et al. 2014; Mosser et al. 2018). Using the ν_max − Δν scaling law in RGB stars (e.g., Mosser et al. 2012a), the observed frequency range is contained in [ν_max − δν, ν_max + δν], where $\delta \nu \approx 3\mathrm{\Delta }\nu \approx 0.84{\nu }_{\max }^{0.75}$. Considering a large range of values for ν_max between 40 and 250 μHz in such typical RGB stars (see Table 2), we approximately obtain 0.21 ν_max ≲ δν ≲ 0.33 ν_max. We thus consider a mean value of δν ≈ 0.27 ν_max in the following. Therefore, using the theoretical framework presented in Sect. 4.3, we find that δq = |q(ν_max + δν)−q(ν_max − δν)| ≳ 0.02 whatever the value of α_b between 0.4 and 0.6, for 2πν_max/𝒩_b ∈ [0.4,  2] (see Fig. 8). This range of values for 2πν_max/𝒩_b corresponds to stars between the transition phase and one point just before the luminosity bump in the Type-b regime. It turns out that the difference in q over the observed frequency range reaches a maximum value of about δq ≈ 0.05, 0.07 and 0.08 at 2πν_max/𝒩_b ∼ 1, for α_b = 0.4, 0.5 and 0.6, respectively. Given the current precision on the measurement of q of about 0.02 (e.g., Hekker et al. 2018), this result clearly demonstrates that the variation of q with frequency in evolved red giant stars around the transition to the Type-b regime cannot be neglected. We emphasize that such a large variation over the observed frequency range is consistent with the variations in the asymptotic values of q computed by Hekker et al. (2018) in RGB stellar models. It is also consistent with the frequency-dependence of q that Cunha et al. (2019) had to consider to improve the fit of the asymptotic expression of mixed modes to the frequency pattern of RGB stellar models with typically ν_max ≲ 100 μHz, which correspond to stars with Type-b evanescent regions.

In conclusion, while the developed model supports the frequency-independence hypothesis in subgiant, young red giant and red clump stars with Type-a evanescent regions, it demonstrates in contrast the need to account for the frequency-dependence of q over the observed frequency range of red giant stars in the transition and in the Type-b regimes, that is, with about 5 μHz≲Δν   ≲  18 μHz according to the estimate given in Table 2. In particular, this point may explain the too low asymptotic values of q around 0.05 computed in evolved stellar models compared to the values fit on theoretical and observed frequency spectra assuming q as a frequency-independent parameter (Mosser et al. 2017; Hekker et al. 2018). Further studies beyond this work are needed to understand how the frequency-dependence of q can modify the previous measurements and their physical interpretation.

6. Discussion

6.1. Possible improvements of the model

The analytical model that was developed in this work relies on several simplifying assumptions regarding the internal structure. In this section, we discuss the most critical points and the possible future improvements.

First, in Type-a evanescent regions, the model assumes that 𝒩 and 𝒮₁ vary as power laws of radius with the same exponent β. Under this hypothesis, the predicted value of q is independent of the oscillation frequency in a given star. However, while Fig. 1 shows that 𝒩 and 𝒮₁ vary as power laws of radius at first approximation, it also shows that they can have slightly different power-law exponents. This is actually related to the existing small gradient in the J factor over the evanescent region (see Fig. 2, Appendix B.1.3 and the third paragraph of Appendix B.1.5 for details). As a result, the thickness of the evanescent region becomes dependent on the oscillation frequency; this feature certainly represents the reason for the frequency-dependence of q observed in stellar models of subgiant, young red giant and red clump stars (Jiang & Christensen-Dalsgaard 2014; Mosser et al. 2017; Hekker et al. 2018). In this picture, the coupling factor will depend not only on the thickness of the evanescent region and the power-law exponent of 𝒮₁, but also on that of 𝒩. Besides, in Type-b evanescent regions, we demonstrate in Appendix B.1.5 that 𝒮₁ cannot physically behave as a rigorous power law of radius, unlike what was assumed in the present work. In order to investigate analytically how the expression of q is modified by such considerations in both types of evanescent regions, a first approach may consist in accounting for a first-order perturbation of the power-law structure and subsequently deducing the impact on the profiles of 𝒮₁ and 𝒩 in a consistent way by means of the continuity and hydrostatic equations (e.g., see Appendices B.1.1 and B.1.5).

Second, the model neglects the presence of a sharp μ-gradient observed at the base of the convective zone in evolved red giant stellar models. This sharp μ-gradient is located in the deepest layer attained by the base of the convective region during the first dredge-up of the chemical elements on the RGB. Once the ascending hydrogen-burning shell encounters this layer, the resulting structural readjustment will produce the well-known RGB luminosity bump (e.g., Kippenhahn et al. 2012). The observational detection of the luminosity bump confirms the physical existence of such a feature in evolved red giant stars (e.g., Khan et al. 2018). This sharp μ-gradient creates a spike in the Brunt–Väisälä frequency and a sharp jump in the Lamb frequency near the base of the convective region; this can be seen in the model ER of Fig. 1. These rapid variations can even be sharper in slightly less evolved models (e.g., see Fig. 1 of Mosser et al. 2017). It thus appears obvious that such features can have an impact on the coupling factor of mixed modes during the transition from a Type-a to a Type-b evanescent region. Indeed, close to this sharp μ-gradient, the asymptotic approximation is not valid and the wave transmission factor through the evanescent region can be modified. Moreover, because of the resulting jump in the sound speed and the Lamb frequency, the outer turning point associated with the evanescent region can also be discontinuous as a function of the oscillation frequency, and so does the thickness of the evanescent region. We note that this feature represents, with the frequency-dependence of q, a plausible explanation for the discrepancies between the low asymptotic values computed in evolved stellar models and the observed values (Mosser et al. 2017; Hekker et al. 2018). In order to investigate the effect of the sharp μ-gradient on the coupling factor, future theoretical efforts will have to be made to go beyond the current asymptotic analyses, which are unable to deal properly with these sharp variations.

6.2. Effect of stellar mass

While the typical observed stellar masses range between 1 M_⊙ and 2 M_⊙, only 1.2 M_⊙ stellar models were considered as examples throughout this paper. In particular, the modeling of the evanescent region along with evolution was based on the internal structure of 1.2 M_⊙ stellar models. Since the properties of the stellar interior depends on stellar mass, the extension of the conclusions of the present work to the whole typical mass range needs to be discussed.

On the RGB, low-mass stars all follow the Hayashi line in the Hertzsprung–Russell diagram and thus share similar internal properties. It is easy to check through stellar models that the power-law behavior of 𝒩 and 𝒮₁, the high density contrast compared to the helium core and the transition from a Type-a to a Type-b evanescent region are also observed for typical masses between 1 M_⊙ and 2 M_⊙. The predictions and the interpretation of the variations in q on the RGB thus hold valid for the whole observed mass range.

On the subgiant branch, we checked using stellar models that stars with masses below about 1.4 M_⊙ again share similar properties to those depicted in the considered 1.2 M_⊙ models. Nevertheless, we observed in more massive subgiant stars that the deviation from the parallel variations of 𝒩 and 𝒮₁ discussed in Sect. 6.1 is worsened; furthermore, the power-law assumption can even become questionable (at least, in the early subgiant phase). This appears to be related to a smaller density contrast over the evanescent region of such more massive stars. Indeed, we observed that the minimum value of the J factor between the helium core and the base of the convective region is lower than, for instance, in the model SG (see Fig. 2); as a result, the gradient of J is higher and the departure from the power-law structure is more important (e.g., Appendix B.1.3). Interestingly, stars with masses higher than 1.4 M_⊙ correspond to stars whose mass of the helium core exceeds the Schönberg–Chandrasekhar limit (for a solar chemical composition) at the end of the main sequence. In such stars, the density contrast between the helium core and the envelope needs to rapidly increase on the subgiant branch to reach a new thermal equilibrium (since the core has to contract while the envelope has to expand by mirror effect; e.g., see Kippenhahn et al. 2012). We thus suspect that the observed discrepancy compared to 1.2 M_⊙ stars is associated with the thermal instability of the helium core at the beginning of the subgiant branch in more massive stars. This may lead to a different relation between q and the properties of the evanescent region in subgiant stars with masses higher than about 1.4 M_⊙ compared to that predicted in the present work. This issue deserves to be addressed in more detail in Paper II by means of stellar models.

6.3. Need for a medium coupling formalism

The available asymptotic analyses provide the expression for the coupling factor in two limiting cases only. On the one hand, Shibahashi (1979) assumed that the turning points associated with the evanescent region are far away from each other. The wavenumber integral over the considered region is much higher than unity and the effect of the gradients in the equilibrium structure near the evanescent region can be neglected. This formalism can thus be applied to describe thick evanescent zones. On the other hand, Takata (2016a) assumed that both turning points are very close to each other and accounted for the zeroth-order influence of the gradients in the equilibrium structure. More specifically, he assumed a particular ordering in Eq. (50) of his paper, which reads $\mathrm{\Delta }r/r_{\mathrm{ev}}\approx 2tanh(\stackrel{\sim }{a}{K}^{-1/2})$ in our notations, where K and $\stackrel{\sim }{a}$ are a large real constant and a real function with a modulus lower than about unity, respectively. This formalism can thus be used to deal with very thin evanescent zones.

However, none of these formalisms formally describes the intermediate case of a medium coupling. As pointed out in Sect. 4.1.1, the expression of Takata (2016a) does not even converge to the expression of Shibahashi (1979) when the evanescent zone is artificially made very large. This dichotomy between both limiting cases and their unclear domains of applicability make ambiguous the comparison between asymptotic values of q and observations, and demanded to use simplifying additional assumptions in Sect. 5. New theoretical developments in the hypothesis of a medium coupling, that is, between both previous limiting cases, are therefore needed to overcome these issues. Such improvements are currently investigated (Takata, priv. comm.).

7. Conclusions and prospects

In the present work, we investigated the potential of the coupling factor of dipolar mixed modes, q, in probing the internal properties of post-main sequence stars. In this first paper of the series, an analytical approach, based on the asymptotic theory of mixed modes, was considered.

This seismic parameter is known to be sensitive to the structure of the evanescent coupling region of mixed modes. First, guided by stellar models, we highlighted that the evanescent region of typically-observed mixed modes progressively migrates from the radiative to the convective regions after the main sequence. In subgiant and young red giant stars, the evanescent region is located in the radiative layers between the hydrogen-burning shell and the base of the convective zone. There, the Brunt–Väisälä and Lamb frequencies vary at first approximation as power laws of radius with the same exponent. We clearly showed using the hydrostatic and the continuity equations that this peculiar behavior directly originates from the high local density contrast compared to the helium core. In evolved RGB stars, which correspond to stars with Δν ≲ 10 μHz, the evanescent region is situated in the lowest part of the convective region, where the Lamb frequency also varies at first approximation as a power law of radius but the Brunt–Väisälä frequency vanishes. Assuming such a power-law structure and accounting for the migration of the evanescent region with evolution, we obtained analytical expressions of the coupling factor based on the asymptotic formulations of dipolar mixed modes of Shibahashi (1979) and Takata (2016a). The non-negligible effect of the perturbation of the gravitational potential on the coupling factor of dipolar modes was clarified and taken into account.

The model shows that the value of q depends in general on two structural parameters: first, the radial extent of the evanescent region normalized by its middle radius, denoted with Δr/r_ev; second, the density scale height in the evanescent region, represented by the logarithmic slope of the critical frequencies, denoted with β. The model predictions provide us with some clues about the reasons for the variations of q during evolution observed by Mosser et al. (2017):

– From the subgiant branch to the beginning of the RGB, that is, for Δν ≳ 18 μHz, the peak in the observed value of q is necessarily associated to a decrease in Δr/r_ev to at least about 0.15, followed by a subsequent increase. The large spread in the measurements during the peak phase suggests that the coupling region can even become much thinner.

– In the following transition phase, which we estimated for 10 μHz ≲Δν ≲ 18 μHz, the evanescent region progressively migrates from the radiative to the convective regions. The small observed decrease in q over this phase can result either from a predominant increase in Δr/r_ev, either from the predominance of the progressive migration of the evanescent region, or from both ingredients simultaneously. Unlike the gravity offset, ε_g, we find that the progressive migration of the evanescent region is not associated with a sharp variation in q.

– In more evolved red giant stars, that is, for Δν ≲ 10 μHz, the coupling factor turns out to be mainly sensitive to Δr/r_ev; the slight decrease in q observed in such stars mainly results from an increase in Δr/r_ev as ν_max decreases.

– In red clump stars, the evanescent regions of which are in the radiative region, the negligible variations in q around 0.3 and the small spread in the measurements indicate that the structure of their mid-layers very slowly varies during the helium burning phase.

The model also gives us more specific insights into the probing potential of q. In very thin evanescent regions, as found in subgiant and red clump stars, the degeneracy of q prevents us from unambiguously concluding on the mid-layer properties without any other information than its value. In contrast, in evolved RGB stars, q can provide us with stringent constraints on the location and migration of the base of the convective region, as well as on the local convective extra-mixing. This appears to be promising for instance for the study of the first dredge-up or the RGB luminosity bump. Nevertheless, we find that a relevant comparison between the theoretical predictions and observations in these evolved stars needs to account for the frequency-dependence of the coupling factor over the observed frequency range, which is not considered in the usual measurement methods and could affect the physical interpretation of the observations. Finally, the model also highlights the need for the development of a medium coupling formalism, tackling the intermediate case between very thin and thick evanescent regions. This will enable us to assess properly and unequivocally the relation between q and the internal properties, in particular in young subgiant and young red giant stars; this point will demand to go beyond the current theoretical picture of mixed modes.

All the lessons already learned from the analytical model will help us go further in the investigation. The next step will consist in exploring the link between the asymptotic value of q and the internal properties using stellar models of different masses, and confronting the results to the observations. This will permit to question the results of the present study, as well as the structural origin of the variations in q in subgiant and red clump stars, for which the degeneracy of q with the internal properties prevents us from directly concluding. This will be subject of the forthcoming paper II of this series. However, in order to test stellar models, it is also necessary to develop seismic diagnoses that can independently probe the evanescent region, as for instance the diagnosis developed by Goupil et al. (2013) for the internal rotation. Another seismic index is thus necessary for making such a model-independent probing of the mid-layer structure possible. The simultaneous measurement of the gravity offset of mixed modes with q represents an interesting option. Actually, Takata (2016a) demonstrated that ε_g also depends on Δr/r_ev and β. Moreover, in evolved RGB stars, Pinçon et al. (2019) showed that ε_g is also sensitive to the base of the convective region and overshooting. Therefore, both measurements of q and ε_g could be complementary and provide a model-independent information on the mid-layer structure of red giants. This promising possibility will be studied in a further paper.

---

Acknowledgments

C. P. acknowledges Masao Takata for the fruitful and helpful discussions about mixed modes, for his careful reading of the paper, as well as for his support. We also warmly thank Benoît Mosser for his valuable remarks on the paper. During this work, C. P. was partially funded by a postdoctoral fellowship from F.R.S.-FNRS (Belgium) and from CNRS (France).

References

Appendix A: Impact of the Cowling approximation on dipolar mixed modes

The Cowling approximation consists in neglecting in the stellar oscillation equations all the terms associated with the Eulerian perturbation of the gravitational potential, denoted by Φ′ (e.g., Cowling 1941; Unno et al. 1989). This hypothesis is actually equivalent to assume ρ∇_rΦ′≪ρ′g in the radial linearized momentum equation, where ∇_r is the radial component of the nabla operator, ρ is the density, g is the gravitational acceleration and ρ′ represents the Eulerian perturbation of density. Indeed, it is possible to show that this latter assumption implies ρΦ′≪p′ in the horizontal linearized momentum equation, where p′ represents the Eulerian perturbation of pressure. As a result, the linearized Poisson equation, ∇²Φ′ = 4πGρ′, becomes decoupled from the momentum equation, which greatly simplifies the stellar oscillation equations. In this section, we discuss the validity of the Cowling approximation and in particular, its effect on the physical description of dipolar mixed modes in evolved stars.

A.1. Validity conditions

In the following, we assume ∇_rΦ′∼Φ′/λ_r in order of magnitude. In the regions where the WKB approximation is met, λ_r represents either the local radial wavelength or the length of decay of the wave function, depending on whether we consider oscillating or evanescent regions, respectively. Otherwise, it is on the order of the pressure scale height, H_p. In these considerations, the Poisson equation can be rewritten given the decomposition onto the spherical harmonics such as

$$\begin{array}{c}\hfill (\frac{1}{{\lambda }_r^2}+k_h^2){\mathrm{\Phi }}^{\prime }\sim 4\pi G{\rho }^{\prime },\end{array}$$(A.1)

where $k_h^2=\mathcal{l}(\mathcal{l}+1)/r^2$ is the squared horizontal wavenumber with ℓ the angular degree. To go further, we can also rewrite the gravitational acceleration such as g = GM_r/r² = (4πGr/3) ρ_av, where ρ_av is the mean density inside the shell of radius r, which satisfies the condition ρ_av(r) > ρ(r). Using the latter expression of g, ∇_rΦ′∼Φ′/λ_r and Eq. (A.1), orders of magnitude give

$$\begin{array}{c}\hfill \left(\frac{\rho {\mathrm{\nabla }}_r{\mathrm{\Phi }}^{\prime }}{{\rho }^{\prime }g}\right)\sim 3\frac{\rho }{{\rho }_{\mathrm{av}}}\frac{({\lambda }_r/r)}{1+\ell (\ell +1){({\lambda }_r/r)}^2}·\end{array}$$(A.2)

Equation (A.2) tells us that the Cowling approximation is appropriate when at least one of three following conditions is met: (A) the considered angular degree ℓ is higher than unity; (B) the short-wavelength WKB approximation is met, that is, if λ_r/r ≪ 1; and (C) the density contrast between the considered layer and the inner bulk is high, that is, if ρ/ρ_av ≪ 1 (i.e., J ∼ 1).

The conditions (A) and (B) are related to an averaging effect over the volume of the star. Indeed, if the number of nodes in the horizontal or in the radial directions is large enough, the small-scale density perturbations integrated over the whole volume of the star are expected to have only a small effect on the local gravitational acceleration. The condition (C) is related to the central character of the gravitational interaction. If ρ/ρ_av ≪ 1, which is equivalent to the limit of J → 1 according to Eq. (5), the layers below the shell of radius r can be regarded as a massive central point. The local gravitational acceleration thus remains mostly radial and quasi-unaffected compared to the equilibrium state. Actually, this confirms through orders of magnitude the result of Takata (2016a) who rigorously demonstrated for dipolar mixed modes that the oscillation equations accounting for the effect of Φ′ tends to the equations within the Cowling approximation in the limit of J → 1, as mentioned in Sect. 2.1.

In the case of dipolar modes, the first condition is not met since the horizontal nodes are concentrated on only one equatorial plane. Nevertheless, the second and the third conditions can still be met depending on the regions.

A.2. Effect on dipolar modes in evolved stars

Well inside both resonant cavities of evolved stars, that is, far enough away from the associated turning points, the WKB approximation is met and the Cowling hypothesis is consequently valid. This point can be more thoroughly discussed in the light of the analysis of dipolar mixed modes by Takata (2016a). In the inner buoyancy cavity where σ² ≪ (𝒩² and ${\mathcal{S}}_1^2$), the dispersion relation for dipolar modes in Eq. (1) becomes

$$\begin{array}{c}\hfill k_{\mathrm{r}}^2=\frac{2N^2}{r^2{\sigma }^2}[1+O\left(\frac{{\sigma }^2}{{\mathcal{N}}^2}\right)]\approx \frac{2N^2}{r^2{\sigma }^2},\end{array}$$(A.3)

while in the external acoustic cavity where σ² ≫ (𝒩² and ${\mathcal{S}}_1^2$), it reads

$$\begin{array}{c}\hfill k_{\mathrm{r}}^2=\frac{{\sigma }^2}{c^2}[1+O\left(\frac{{\mathcal{S}}_1^2}{{\sigma }^2}\right)]\approx \frac{{\sigma }^2}{c^2}·\end{array}$$(A.4)

As shown by Eqs. (A.3) and (A.4), the dispersion relation is not modified at leading order by the Cowling approximation in such regions. The correction resulting from the effect of Φ′ (i.e., through the J factor) appears only in higher order terms. This explains why the asymptotic expressions of the period spacing and the frequency large separation of dipolar mixed modes, which correspond to leading-order terms in the asymptotic formulation, are not modified within the Cowling approximation, unlike the gravity offset that corresponds to higher order terms (Takata 2016a; Pinçon et al. 2019). We also note that Eq. (A.3) confirms in the case of dipolar oscillations the result of Dintrans & Rieutord (2001), who demonstrated that the Cowling approximation is valid for low-frequency gravity waves of any angular degree because of their incompressible character (i.e., when σ² ≪ N²).

In contrast, inside the evanescent region and in the vicinity of the associated turning points, the WKB approximation does not apply since λ_r ∼ H_p. First, this is easy to confirm close to the turning points since $k_r^2\to 0$. Second, inside the evanescent zone, Eq. (1) shows that the squared modulus of the radial wavenumber is equal to

$$\begin{array}{c}\hfill |k_r^2|=\frac{2J^2}{r^2}(1-\frac{{\mathcal{N}}^2}{{\sigma }^2})(1-\frac{{\sigma }^2}{{\mathcal{S}}_1^2})\lesssim \frac{2J^2}{r^2}·\end{array}$$(A.5)

As a result, the Cowling approximation is valid if the third condition about the density contrast (i.e., the limit of J → 1) is met in this intermediate region. Figure 1 shows that the value of J can be low enough for the thickness of the intermediate coupling region to be significantly changed compared to the case within the Cowling approximation. Moreover, according to Eq. (A.5), the value of the wavenumber inside the evanescent region is also modified. Therefore, the value of the coupling factor must also be affected. We note that the younger the star, the higher the correction induced by Φ′, as the density contrast between the core and the envelope (and thus the J factor) is lower, as shown in Figs. (1) and (2).

To summarize, our analysis shows that the Cowling approximation can lead to incorrect description of the coupling of mixed modes. The formulation of Takata (2006a) in the non-Cowling case is mandatory for the study of the coupling factor of dipolar mixed modes, in particular for subgiant and young red giant stars.

Appendix B: On the structure of the evanescent regions in evolved stars

As seen in Fig. 1, the (modified) Brunt–Väisälä and Lamb frequencies behave at first approximation as power laws of radius in the evanescent region of evolved stars. In this section, we demonstrate that the origin of such a behavior is directly due to the high density contrast existing between the helium core and the envelope, which is measured by the J quantity defined in Eq. (5).

B.1. Link between J and the critical frequencies

In a first step, we show that the logarithmic derivatives of N and S_ℓ, denoted with β_N and β_S, respectively, as well as the J factor are linked to each other through the hydrostatic equilibrium and the continuity equations. We then demonstrate that an equivalence relation between these quantities exists when β_N and β_S are equal constants. The result is extended to the modified definitions 𝒩 and 𝒮₁ in the case of dipolar modes. The following statements only assume that N and S_ℓ are decreasing functions of radius and that the first adiabatic index Γ₁ is constant. The case of radiative regions, in which N² >  0, is considered before discussing the case of convective regions, in which we assume N² = 0.

B.1.1. Constitutive equations in radiative layers

In this section, we aim to link J, β_N and β_S through the hydrostatic equilibrium and the continuity equations. First, according to Eq. (5), which results from the continuity equation, we obtain

$$\begin{array}{c}\hfill \frac{\mathrm{d}ln{M}_r}{\mathrm{d}x}=3(1-J),\end{array}$$(B.1)

where x = lnr. Then, by differentiating J with respect to x and using again Eq. (5), the derivative of ρ with respect to x, which is renamed as −ℬ(x), can be expressed as a function of J only, that is,

$$\begin{array}{c}\hfill \mathcal{B}(x)\equiv -\frac{\mathrm{d}ln\rho }{\mathrm{d}x}=3J+\frac{1}{1-J}\frac{\mathrm{d}J}{\mathrm{d}x}·\end{array}$$(B.2)

To go further, we consider the hydrostatic equilibrium equation, that is, (dp/dx) = − ρgr, where g = GM_r/r² is the gravitational acceleration. Dividing this equation by the product −Γ₁p, taking the logarithmic derivative, and using Eqs. (B.1) and (B.2), we obtain

$$\begin{array}{cc}\hfill \frac{\mathrm{d}ln\mathcal{V}}{\mathrm{d}x}& =\frac{\mathrm{d}ln\rho }{\mathrm{d}x}+\frac{\mathrm{d}ln{M}_r}{\mathrm{d}x}-1+{\mathrm{\Gamma }}_1\mathcal{V}\hfill \\ \hfill & =-\mathcal{B}+3(1-J)-1+{\mathrm{\Gamma }}_1\mathcal{V},\hfill \end{array}$$(B.3)

where Γ₁ was considered as constant and we defined the positive quantity

$$\begin{array}{c}\hfill \mathcal{V}=-\frac{1}{{\mathrm{\Gamma }}_1}\frac{\mathrm{d}lnp}{\mathrm{d}x}·\end{array}$$(B.4)

The final step consists in rewriting the logarithmic derivatives of the Lamb and Brunt–Väisälä frequencies as a function of J, 𝒱 and ℬ only. First, the squared Lamb frequency (for any angular degree ℓ) can be rewritten from the hydrostatic equilibrium as

$$\begin{array}{c}\hfill S_{\ell }^2=\frac{\ell (\ell +1)c^2}{r^2}=\frac{\ell (\ell +1){\mathrm{\Gamma }}_{1}p}{\rho r^2}\frac{\rho g}{-(\mathrm{d}p/\mathrm{d}r)}=\frac{\ell (\ell +1)}{\mathcal{V}}\frac{g}{r}·\end{array}$$(B.5)

Using the fact dln(g/r)/dx = −3J according to Eq. (B.1), the logarithmic derivative of the Lamb frequency is given by

$$\begin{array}{c}\hfill {\beta }_{\mathrm{S}}(x)\equiv -\frac{\mathrm{d}ln{S}_{\ell }}{\mathrm{d}x}=\frac{1}{2}(3J+\frac{\mathrm{d}ln\mathcal{V}}{\mathrm{d}x}).\end{array}$$(B.6)

In a similar way, the logarithmic derivative of the Brunt–Väisälä frequency is equal, according to Eq. (4) with N = J𝒩, to

$$\begin{array}{c}\hfill {\beta }_{\mathrm{N}}(x)\equiv -\frac{\mathrm{d}lnN}{\mathrm{d}x}=\frac{1}{2}(3J-\frac{\mathrm{d}ln\mathcal{A}}{\mathrm{d}x}),\end{array}$$(B.7)

where

$$\begin{array}{c}\hfill \mathcal{A}=(\frac{1}{{\mathrm{\Gamma }}_1}\frac{\mathrm{d}lnp}{\mathrm{d}x}-\frac{\mathrm{d}ln\rho }{\mathrm{d}x})=-\mathcal{V}+\mathcal{B}.\end{array}$$(B.8)

We emphasize that 𝒜 is a positive quantity in the radiative zone since N² = g𝒜/r ≳ 0.

To derive a differential system linking J, β_S and β_N, we first express 𝒱 from Eq. (B.6) such as

$$\begin{array}{c}\hfill \mathcal{V}(x)={\mathcal{V}}_{0}exp\left({\displaystyle {\int }_{x_0}^x}[2{\beta }_{\mathrm{S}}-3J]\mathrm{d}x\right),\end{array}$$(B.9)

where the subscript 0 denotes the value at r = r₀ (i.e., at x = x₀). Then, by integrating Eq. (B.7), we also obtain

$$\begin{array}{c}\hfill \mathcal{A}(x)=-\mathcal{V}+\mathcal{B}={\mathcal{A}}_{0}exp(-{\displaystyle {\int }_{x_0}^x}[2{\beta }_{\mathrm{N}}-3J]\mathrm{d}x).\end{array}$$(B.10)

By injecting Eq. (B.9) into Eq. (B.3) and (B.10), the ℬ quantity in Eq. (B.2) can be expressed following three different ways as a function of the dependent variables J, β_S and β_N, which are

$$\begin{array}{cc}\hfill \mathcal{B}(x)& =3J+\frac{1}{1-J}\frac{\mathrm{d}J}{\mathrm{d}x}\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill & =-2{\beta }_{\mathrm{S}}+2+{\mathrm{\Gamma }}_1{\mathcal{V}}_{0}exp\left({\displaystyle {\int }_{x_0}^x}[2{\beta }_{\mathrm{S}}-3J]\mathrm{d}x\right)\hfill \end{array}$$(B.12)

$$\begin{array}{cc}\hfill & ={\mathcal{A}}_{0}exp(-{\displaystyle {\int }_{x_0}^x}[2{\beta }_{\mathrm{N}}-3J]\mathrm{d}x)\hfill \\ \hfill & +{\mathcal{V}}_{0}exp\left({\displaystyle {\int }_{x_0}^x}[2{\beta }_{\mathrm{S}}-3J]\mathrm{d}x\right).\hfill \end{array}$$(B.13)

Defining the dependent variables W, X, Y and Z such as

$$\begin{array}{cc}& W={\mathcal{A}}_{0}exp(-2{\displaystyle {\int }_{x_0}^x}{\beta }_{\mathrm{N}}\mathrm{d}x),X={\mathcal{V}}_{0}exp\left(2{\displaystyle {\int }_{x_0}^x}{\beta }_{\mathrm{S}}\mathrm{d}x\right),\hfill \\ \hfill & Y=exp\left({\displaystyle {\int }_{x_0}^x}3J\mathrm{d}x\right),Z=\frac{\mathrm{d}Y}{\mathrm{d}x}=3JY,\hfill \end{array}$$(B.14)

Equations (B.11)–(B.13) can be rewritten in the form of a first-order ordinary nonlinear system of three differential equations, that is,

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}x}\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)=\left(\begin{array}{c}2X+({\mathrm{\Gamma }}_1-1){\displaystyle \frac{X^2}{Y}}-WXY\\ Z\\ 2Z^{2}Y+3Z-3W+WZY-{\displaystyle \frac{X}{Y}}({\displaystyle \frac{3}{Y}}-Z)\end{array}\right),\end{array}$$(B.15)

with W₀ = 𝒜₀, X₀ = 𝒱₀, Y₀ = 1 and Z₀ = 3J₀ as initial values. Therefore, given the definitions of W, X, Y and Z in Eq. (B.14), if at least one of the three variables J, β_S and β_N is specified around r₀, it is possible to deduce the two other ones through Eq. (B.15). Since the coupling factor of mixed modes mostly depend on N, S_ℓ and J (see Sect. 3.3), its expression thus must account for this relation in order to be consistent.

B.1.2. Expression of J with given β_S and β_N

Assuming that both β_S and β_N are known functions of r, the J factor can actually be analytically expressed in a direct way. Indeed, by multiplying Eqs. (B.12) and (B.13) by the positive variable Y, we obtain the following quadratic equation for Y,

$$\begin{array}{c}\hfill A{Y}^2+BY+C=0,\end{array}$$(B.16)

with

$$A={\mathcal{A}}_{0}exp(-2{\displaystyle {\int }_{x_0}^x}{\beta }_{\mathrm{N}}\mathrm{d}x)>0$$(B.17)

$$B=2({\beta }_{\mathrm{S}}-1)$$(B.18)

$$C=-{\mathcal{V}}_0[{\mathrm{\Gamma }}_1-1]exp\left(2{\displaystyle {\int }_{x_0}^x}{\beta }_{\mathrm{S}}\mathrm{d}x\right)<0.$$(B.19)

The sign of C results from the thermodynamical condition Γ₁ >  1. The discriminant of the second-order equation, Δ = B² − 4AC, is strictly positive and greater than B². The unique (positive) solution is thus given by $Y=(-B+\sqrt{\mathrm{\Delta }})/2A$. Computing the logarithmic derivative of Y, we then obtain

$$\begin{array}{c}\hfill J=\frac{2}{3}[{\beta }_{\mathrm{N}}-\frac{1}{\mathit{AY}}(\frac{{\beta }_{\mathrm{S}}^{\prime }}{2}-\frac{{\beta }_{\mathrm{S}}^{\prime }({\beta }_{\mathrm{S}}-1)+AC[{\beta }_{\mathrm{N}}-{\beta }_{\mathrm{S}}]}{\sqrt{\mathrm{\Delta }}})],\end{array}$$(B.20)

where ${\beta }_{\text{S}}^{\prime }(x)=(\text{d}{\beta }_{\text{S}}/\text{d}x)$.

B.1.3. Equivalence relation

Analysis. If β_S and β_N are both fixed to a given constant β, then Eq. (B.20) shows that J is uniform with respect to radius and equal to 2β/3. According to Eq. (B.2), this assumption thus implies that the density varies as a power law with respect to radius, that is, such as ℬ = −(dlnρ/dx) = 3J. Moreover, Eq. (B.6) firstly shows that (dln𝒱/dx) is constant; secondly, Eq. (B.3) shows that this constant is zero and that −(dlnp/dx) = 6J − 2. It is easy to check that this solution also satisfies Eq. (B.11). Therefore, the pressure also varies as a power law with respect to radius. In other words, the equation of state has a polytropic form such as P ∝ ρ^γ, where the polytropic exponent is equal to γ = 2 − 2/(3J). In conclusions, a necessary condition for β_S and β_N to be equal constants is that J is constant and the equation of state is polytropic.

Synthesis. If J is uniform and P ∝ ρ^γ, (dlnρ/dx) and (dlnp/dx) are also uniform according to Eq. (B.2). Therefore, we deduce according to Eq. (B.4) and (B.8) that β_S and β_N are both constants equal to 3J/2, and that γ = 2 − 2/(3J).

Conclusions. The Brunt–Väisälä and Lamb frequency follow decreasing power laws of radius with the same exponent β if and only if J is uniform and the equation of state is polytropic. The equivalence configuration corresponds to a power-law solution of the internal structure, where (dlnp/dx), (dlnρ/dx) and (dlnM_r/dx) are constants. As a corollary, J, β and γ are related to each other through the equations β = 3J/2, as previously found by Takata (2016a), and γ = 2 − 2/(3J).

B.1.4. Extension to 𝒩 and 𝒮₁ for dipolar modes

We now consider the modified Brunt–Väisälä and Lamb frequencies for dipolar modes, the logarithmic derivatives of which are denoted with ${\beta }_{\text{N}}^{\star }$ and ${\beta }_{\text{S}}^{\star }$.

On the one hand, if J is supposed to be uniform and P ∝ ρ^γ, it is then straightforward to show that ${\beta }_{\text{N}}^{\star }={\beta }_{\text{S}}^{\star }=3J/2$ following the same reasoning as in the second paragraph of Appendix B.1.3 and the definitions of 𝒩 and 𝒮₁ in Eq. (2).

On the other hand, using ${\beta }_{\text{N}}={\beta }_{\text{N}}^{\star }+(\text{d}\ln J/\text{d}x)$ and ${\beta }_{\text{S}}={\beta }_{\text{S}}^{\star }-(\text{d}\ln J/\text{d}x)$ in Eqs. (B.11) and (B.13), and defining the new variables W^⋆ and X^⋆ by replacing β_N in Eq. (B.14) with ${\beta }_{\text{N}}^{\star }$ and β_S with ${\beta }_{\text{S}}^{\star }$, it is possible to demonstrate that W^⋆, X^⋆, Y and Z are also linked by a first-order nonlinear system of three differential equations. If we impose that ${\beta }_{\text{N}}^{\star }$ and ${\beta }_{\text{S}}^{\star }$ are equal to a same constant, the differential system reduces to a nonlinear first-order system of two differential equations for Y and Z. The constant $\beta ={\beta }_{\text{N}}^{\star }={\beta }_{\text{S}}^{\star }$ is linked to the derivatives of W^⋆ and X^⋆, so that its choice is not arbitrary but is determined by the initial values of the problem. Actually, using the relation ℬ = 𝒜 + 𝒱, it is easy to show that a trivial solution of the differential system given by Eqs. (B.2), (B.3), (B.6) and (B.7), in which we use ${\beta }_{\text{S}}={\beta }_{\text{S}}^{\star }-(\text{d}\ln J/\text{d}x)$ and ${\beta }_{\text{N}}={\beta }_{\text{N}}^{\star }+(\text{d}\ln J/\text{d}x)$, is also J = J₀. This condition again implies β = 3J₀/2. Assuming that the Cauchy-Lipschitz theorem on the initial value problem is applicable over the evanescent region, which is usually the case, the constant solution J = J₀ is the unique solution. Moreover, since we also have ${\beta }_{\text{N}}={\beta }_{\text{N}}^{\star }$ and ${\beta }_{\text{S}}={\beta }_{\text{S}}^{\star }$ in this case, the equation of state also turns out to be polytropic according to the first paragraph of Appendix B.1.3.

Therefore, the same conclusions are reached for the modified critical frequencies; ${\beta }_{\text{N}}^{\star }$ and ${\beta }_{\text{S}}^{\star }$ are equal constants if and only if J is uniform and the equation of state is polytropic, which implies β = 3J/2, γ = 2 − 2/(3J), and a power-law configuration for the internal structure.

B.1.5. Case of convective regions

In the case of an adiabatic convective zone, we have 𝒜 ≈ 0 since P ∝ ρ^Γ₁. A general analytical expression of J can thus be directly derived. Indeed, we have 𝒱 ≈ ℬ according to Eq. (B.8) with 𝒜 ≈ 0. Using Eqs. (B.3) and (B.6), we can thus express 𝒱 as a function of β_S and inject the result into Eq. (B.3) to finally obtain

$$\begin{array}{c}\hfill J=\frac{2}{3}{\beta }_{\mathrm{S}}-\frac{{\beta }_{\mathrm{S}}^{\prime }}{3({\beta }_{\mathrm{S}}-1)}·\end{array}$$(B.21)

Following the same reasoning as in Appendix B.1.3, it is possible to demonstrate that J is uniform if and only if β_S is constant; in this case, β_S = 3J/2, the polytropic exponent is equal to γ = 2 − 2/(3J) and the internal structure behave as power laws of radius.

However, since the temperature gradient is adiabatic in convective zones, γ = Γ₁ = 5/3 for a totally ionized monatomic ideal gas. If the equivalence was met, we would thus obtain J = 2, which is not possible since J <  1 by definition. Such an equivalence configuration is thus not valid in convective regions. Indeed, while the profile of the Lamb frequency behaves at first sight like a power law of radius in the evanescent region of the model ER in Fig. 1, we can see in Fig. 2 that the value of J decreases of about 20% over the evanescent region, and therefore cannot be regarded as rigorously uniform.

To understand the profile of J, it is useful to consider the derivative of J with respect to radius. According to Eq. (B.2), this reads

$$\begin{array}{c}\hfill \frac{\mathrm{d}J}{\mathrm{d}lnr}=(1-J)(-\frac{\mathrm{d}ln\rho }{\mathrm{d}lnr}-3J).\end{array}$$(B.22)

As the value of J rapidly varies from small values to about unity in the vicinity of the hydrogen-burning shell, the second term in brackets in the right-hand side of Eq. (B.22) becomes negative and J starts decreasing with radius. In contrast, in the upper layers close to the surface, the density vanishes and −(dlnρ/dr)≫1. As a consequence, Eq. (B.22) becomes positive and J starts increasing toward unity. These facts result in the existence of a minimum in the J factor that is located in the external convective zone of the model ER (and in the model YR). We emphasize that, in the model SG, the minimum of J is even located in the radiative region because of a lower density contrast⁴. Because of the necessary variation of J in the convective region, the link between β and the J factor is thus provided by the differential equation in Eq. (B.21).

In the case of the modified expression of the Lamb frequency for dipolar modes, we note that the same conclusions hold true. The general expression linking J and ${\beta }_{\text{S}}^{\star }$ is however more complex than Eq. (B.21) since β_S in this equation must be replaced by ${\beta }_{\text{S}}^{\star }-(\text{d}\ln J/\text{d}x)$ and thus also involves the derivatives of J with respect to radius.

B.2. Structure of the upper radiative layers

In Appendix B.1, we have demonstrated that the Brunt–Väisälä and Lamb frequencies vary as power laws of radius with the same exponents between the hydrogen-burning shell and the base of the convective zone if and only if J is uniform and the equation of state is polytropic. While it is easy to confirm through stellar models that these conditions are met at first approximation in this region of evolved stars (e.g., Maeder 2009; Kippenhahn et al. 2012, see also Figs. 1 and 2), we propose in the following to discuss this point considering simplified analytical solutions of the internal structure equations. Actually, we show that these conditions mainly originate from the high density contrast between these layers and the helium core.

B.2.1. Hypothesis of a high density contrast

In a first step, we assume that the density contrast between the considered layers and the helium core is very high. We then aim to show that such an assumption implies that the equation of state is polytropic and J is uniform in a good approximation. First, the definition of J in Eq. (5) shows that the high density contrast compared to the helium core implies ρ/ρ_av ≪ 1 and J ∼ 1. According to Eq. (B.22), the derivative of J with respect to radius is thus expected to be small so that J remains about uniform, as expected. Second, regarding the polytropic approximation, we need to use the following reasonable considerations to move forward:

1. The medium is non-degenerate and the equation of state can be assumed to follow the ideal gas law. In addition, the mean molecular weight can be supposed to be nearly uniform in this region if we neglect the impact of microscopic diffusion and the existence of local sharp variations resulting, for instance, from the migration of the base of the convective during evolution (e.g., Cunha et al. 2015).

2. The luminosity is constant and equal to the total luminosity of the star since it is fully produced in the underlying hydrogen-burning shell.

3. In the range of temperature expected in this region, the opacity predominantly can originate either from the Bremsstrahlung mechanism or even from bound-free electronic transitions near the base of the convective region (e.g., Kippenhahn et al. 2012). The Rosseland mean opacity, κ_op, thus can be supposed to follow the Kramers law, that is, κ_op = κ_op, 0ρ^rT^−s, where κ_op, 0 is a positive constant, r ≈ 1, and s ≈ 7/2.

4. Owing to the high density contrast between the helium core and the envelope, the mass M_r inside a spherical shell of radius r is assumed to change very slowly with respect to radius, unlike the density ρ. It can thus be supposed to be equal to the mass of the helium core, M_c, in the hydrostatic equilibrium equation.

5. Because of the ideal gas law, the temperature and the pressure are supposed to change rapidly with respect to radius in this region since the density is supposed to vary rapidly.

Then, using the hydrostatic equilibrium and the radiative transport equations, the pressure and the density are related to each other by (e.g., see Eq. (21.16) of Kippenhahn et al. 2012)

$$\begin{array}{c}\hfill p\approx p_{\mathrm{i}}{\left(\frac{\rho }{{\rho }_{\mathrm{i}}}\right)}^{(4+s+r)/(3+s)},\end{array}$$(B.23)

where the assumption 5 is considered in order to neglect the integration constants and the subscript i indicates values taken at radius r_i, just above the hydrogen-burning shell. Therefore, the equation of state is polytropic in a good approximation in the upper layers of the radiative region, with a polytropic exponent that is about equal to γ = (dlnp/dlnρ) = (4 + s + r)/(3 + s) = 1.31.

In consequence, assuming a high density contrast between the upper radiative layers and the underlying helium core is sufficient for the equation of state to be polytropic and the J factor to be nearly uniform with radius in a good approximation. As shown in Appendix B.1, this is equivalent to a power-law behavior of the internal structure.

B.2.2. Power-law solutions of the internal structure

Similarly to Applegate (1988), we now search for power-law solutions of the internal structure equations above the hydrogen-burning shell, that is, in the form of

$$\begin{array}{cc}& M_{\mathrm{r}}=M_{\mathrm{i}}{\left(\frac{r}{r_{\mathrm{i}}}\right)}^a,\rho ={\rho }_{\mathrm{i}}{\left(\frac{r}{r_{\mathrm{i}}}\right)}^{-b}\hfill \\ \hfill & p=p_{\mathrm{i}}{\left(\frac{r}{r_{\mathrm{i}}}\right)}^{-c},T=T_{\mathrm{i}}{\left(\frac{r}{r_{\mathrm{i}}}\right)}^{-d},\hfill \end{array}$$(B.24)

Using the Assumptions 1–3 in Appendix B.2.1 and injecting Eq. (B.24) into the continuity, hydrostatic equilibrium, radiative transport equations and the ideal gas law, Applegate (1988) found

$$\begin{array}{ccc}\hfill a=& {\displaystyle \frac{s-3r}{3+s-r}}\hfill & \hfill b={\displaystyle \frac{9+2s}{3+s-r}}\end{array}$$(B.25)

$$\begin{array}{ccc}\hfill c=& {\displaystyle \frac{12+2s+2r}{3+s-r}}\hfill & \hfill d={\displaystyle \frac{3+2r}{3+s-r}}.\end{array}$$(B.26)

Since the Bremsstrahlung effect and bound-free electronic transitions dominate the opacity, the power-law solutions are characterized by a = 1/11, b = 32/11, c = 42/11, and d = 10/11. As a result, the polytropic exponent is still equal to γ = c/b = 1.31. This result justifies a posteriori the assumption 5 made in Appendix B.2.1 about the behavior of P and T with respect to radius. According to Eq. (B.1), we also find J = b/3 = 32/33, which is very close to unity. Therefore, this model suggests that the high density contrast between the considered layers and the helium core is not only sufficient, but also necessary to explain the power-law behavior of the internal structure.

B.3. Concluding synthesis

In the region located between the hydrogen-burning shell and the base of the convective region, the high density contrast compared to the underlying helium core is the main responsible for the similar power-law behavior of the Brunt–Väisälä and Lamb frequencies. Indeed, under this condition, J is uniform and the equation of state is polytropic, which implies that p, ρ and M_r are power laws with respect to r. As a result, the Brunt–Väisälä and Lamb frequencies are also power laws of r with the same exponent, β. In addition, β is linked to the J factor by the simple relation β ≈ 3J/2. In contrast, in the adiabatic convective envelope, such a configuration is forbidden: J varies with respect to radius and is linked to the exponent β through a more complicated differential equation.

Appendix C: Wave transmission and coupling factor through a parabolic barrier

In this section, we study the problem of the wave transmission through a parabolic evanescent barrier. This case does not require any assumptions about the thickness of the evanescent region to be fully solved, contrary to the asymptotic analyses by Shibahashi (1979) and Takata (2016a) for realistic stellar structures (Sect. 3.3). We propose here to retrieve the expression of the wave transmission coefficient T² for this simplified case using the same method as Phinney (1970) in a geophysical context. In this situation, the wave equation to solve is provided by Eq. (15).

To continue, we use the change of variable

$$\begin{array}{c}\hfill {\zeta }^2=4\frac{|k_r(r_{\mathrm{ev}})|}{\mathrm{\Delta }r}{(r-r_{\mathrm{ev}})}^2,\end{array}$$(C.1)

so that Eq. (15) can be rewritten in the form of the Weber’s equation, that is,

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathrm{\Psi }}{\mathrm{d}{\zeta }^2}+(\frac{{\zeta }^2}{4}-a)\mathrm{\Psi }=0,\end{array}$$(C.2)

where the positive parameter a is defined⁵ as

$$\begin{array}{c}\hfill a=\frac{1}{4}|k_r(r_{\mathrm{ev}})|\mathrm{\Delta }r.\end{array}$$(C.3)

The general solution of Eq. (C.2) is given by a linear combination of the complex independent parabolic cylinder functions, namely E(a, ζ) and its complex conjugate E^⋆(a, ζ) (e.g., Abramowitz & Stegun 1972). When |ζ²|≫|a|, that is, far enough from and on both sides of the region around r = r_ev, the first-order asymptotic development of the parabolic cylinder function is given by

$$\begin{array}{c}\hfill E(a,\zeta )\sim \sqrt{\frac{2}{|\zeta |}}exp[+i(\frac{{\zeta }^2}{4}-aln|\zeta |+\frac{\pi }{4}+\frac{{\mathrm{\Psi }}_2}{2})],\end{array}$$(C.4)

where Ψ₂ = arg[Γ(1/2+ia)], Γ is the gamma-function and i is the imaginary unit. In the following, a time-dependence in e^iσt for the wave function and the positive branch in Eq. (C.1) are considered (i.e., such as ζ and (r − r_ev) have the same sign everywhere). Within this convention, Eq. (C.4) shows us that E(a, ζ) corresponds to a regressive wave (i.e., traveling in the negative r-direction) where ζ >  0; it corresponds to a progressive wave (i.e., traveling in the positive r-direction) where ζ <  0. In contrast, E^⋆(a, ζ) corresponds to a progressive wave where ζ >  0 and a regressive wave where ζ <  0. For later purpose, it is helpful to note that the change of variable ζ′= − ζ does not modify the form of Eq. (C.2), so that E(a, −ζ) and E^⋆(a, −ζ) stand also for the regressive and progressive independent solutions where ζ >  0, respectively, and inversely where ζ <  0. Both sets of independent solutions {E(a, ζ),E^⋆(a, ζ)} and {E(a, −ζ),E^⋆(a, −ζ)} are linked by the following connection formula (e.g., Abramowitz & Stegun 1972)

$$\begin{array}{c}\hfill \sqrt{1+e^{2\pi a}}E(a,\zeta )-e^{\pi a}{E}^{\star }(a,\zeta )=i{E}^{\star }(a,-\zeta ).\end{array}$$(C.5)

The transmission and reflection coefficients through the evanescent region can now be easily determined by means of this last relation. We consider an incident progressive wave coming from ζ = −∞ that is both reflected back and transmitted into the negative and positive ζ-direction in the vicinity of the evanescent barrier. In the domain where ζ <  0, the incident and reflected components are thus proportional to the progressive and regressive independent solutions E(a, ζ) and E^⋆(a, ζ), respectively, that is, such as the wave function is equal to

$$\begin{array}{c}\hfill \mathrm{\Psi }=IE(a,\zeta )+r{E}^{\star }(a,\zeta ),\end{array}$$(C.6)

where the constants I and r are the amplitudes of the incident and reflected waves, respectively. In the domain where ζ >  0, the transmitted wave is proportional to the progressive solution E^⋆(a, −ζ), that is, such as the wavefunction can be written

$$\begin{array}{c}\hfill \mathrm{\Psi }=t{E}^{\star }(a,-\zeta )=-it(\sqrt{1+e^{2\pi a}}E(a,\zeta )-e^{\pi a}{E}^{\star }(a,\zeta )),\end{array}$$(C.7)

where the constant t is the amplitude of the transmitted wave and where Eq. (C.5) was used in the second equality. Since the expressions of the wavefunction in Eqs. (C.6) and (C.7) are valid whatever the value of ζ, and since E(a, ζ) and E^⋆(a, ζ) are linearly independent, the transmission and reflection coefficients of the wave energy flux finally read, by identification,

$$\begin{array}{c}\hfill T^2={\left|\frac{t}{I}\right|}^2=\frac{1}{1+e^{2\pi a}}\mathrm{and}{R}^2={\left|\frac{r}{I}\right|}^2=\frac{e^{2\pi a}}{1+e^{2\pi a}}·\end{array}$$(C.8)

These expressions satisfy the energy conservation since R² + T² = 1. Using the wavenumber expression in Eq. (15), we note that the wavenumber integral over the evanescent region, ℐ, is related to the parameter a by

$$\begin{array}{c}\hfill \mathcal{I}={\displaystyle {\int }_{r_{\mathrm{ev}}-\mathrm{\Delta }r/2}^{r_{\mathrm{ev}}+\mathrm{\Delta }r/2}}|k_r|\mathrm{d}r=\frac{\pi }{4}|k_r(r_{\mathrm{ev}})|\mathrm{\Delta }r=\pi a,\end{array}$$(C.9)

so that the transmission coefficient can be written as a function of ℐ only, that is

$$\begin{array}{c}\hfill T^2=\frac{1}{1+e^{2\mathcal{I}}}·\end{array}$$(C.10)

The wave coefficients and the associated coupling factor q computed with Eq. (13) are displayed in Fig. C.1 as a function of ℐ. As seen in Fig. C.1, the smaller ℐ, the larger both T² and q in the ideal case of a parabolic evanescent zone. The result is discussed in more detail within the framework of mixed mode coupling in Sect. 3.2.

Appendix D: Gradient-related term 𝒢₀ in the formalism of Takata (2016a)

D.1. General expression

In the limiting case of a very thin evanescent zone, the asymptotic analyses of Takata (2016a) resulted in

$$\begin{array}{c}\hfill {\mathcal{G}}_0=\frac{\pi }{\kappa (x_0)}{\left(\frac{\mathrm{d}lnc}{\mathrm{d}x}\right)}_{x_0}^2.\end{array}$$(D.1)

In the latter, we have defined x = lnr, x₀ = (x₂ + x₁)/2, and x_{1, 2} = lnr_{1, 2}, where r₁ and r₂ are the turning points for which 𝒩(r₁) = σ and 𝒮₁(r₂) = σ, respectively. The other quantities are given by

$$\kappa =\sqrt{\mathit{ST}}$$(D.2)

$$c={\left(\frac{S}{T}\right)}^{1/4}$$(D.3)

$$S=\frac{1}{(x_2-x)}\frac{2J}{F^2}(1-\frac{{\sigma }^2}{{\mathcal{S}}_1^2})>0$$(D.4)

$$T=\frac{J{F}^2}{(x-x_1)}(1-\frac{{\mathcal{N}}^2}{{\sigma }^2})>0$$(D.5)

$$F=exp\left[\frac{1}{2}{\int }^x({\mathcal{V}}^{\star }-{\mathcal{A}}^{\star }-J)\mathrm{d}x\right],$$(D.6)

where J is defined in Eq. (5),

$$\begin{array}{c}\hfill {\mathcal{V}}^{\star }=-\frac{1}{{\mathrm{\Gamma }}_{1}J}\frac{\mathrm{d}lnp}{\mathrm{d}x}\mathrm{and}{\mathcal{A}}^{\star }=\frac{1}{J}(\frac{1}{{\mathrm{\Gamma }}_1}\frac{\mathrm{d}lnp}{\mathrm{d}x}-\frac{\mathrm{d}ln\rho }{\mathrm{d}x}).\end{array}$$(D.7)

To go further, we also note that according to Eqs. (3) and (4), 𝒱^⋆ and 𝒜^⋆ are respectively related to the (modified) Lamb and Brunt–Väisälä frequencies through

$$\begin{array}{c}\hfill {\mathcal{V}}^{\star }=\frac{1}{{\mathrm{\Gamma }}_{1}J}\frac{\rho gr}{p}=\frac{2J}{{\mathcal{S}}_1^2}\frac{g}{r}\mathrm{and}{\mathcal{A}}^{\star }=J{\mathcal{N}}^2\frac{r}{g},\end{array}$$(D.8)

where the hydrostatic equilibrium was used in the first equality.

D.2. Expression in a simple model

In this section, we compute 𝒢₀ for dipolar modes assuming that the modified Brunt–Väisälä and Lamb frequencies vary as power laws of radius with the same exponent, β. Following the assumptions in Table 1, we focus only on Type-a evanescent zones. To do so, we follow all the successive derivation steps described in Appendices A.3.1 and A.3.2 of Takata (2016a). However, instead of considering that the thickness of the evanescent zone is null, we explicitly take its effect into account. The demonstration is adapted with our notation, that is, by replacing μ by 2β, and s₀ by −(lnα)/2β, with α = (r₁/r₂)^β.

First, using the profiles of 𝒩 and 𝒮₁ in Eq. (8), the resulting relation J = 2β/3 = const. (see Appendix B.1.3), and the fact that $r(x_0)=\sqrt{r_1{r}_2}$, we obtain from Eq. (D.2), (D.4) and (D.5)

$$\begin{array}{c}\hfill \kappa (x_0)=\frac{4\sqrt{2}}{3}{\beta }^2\frac{(\alpha -1)}{ln\alpha }·\end{array}$$(D.9)

In a similar way, using Eqs. (D.3)–(D.6) and J = 2β/3, the derivative of c gives

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}lnc}{\mathrm{d}x}\right)}_{x_0}=-\beta (\frac{1}{ln\alpha }+\frac{\alpha }{1-\alpha })-\frac{{\mathcal{V}}^{\star }(x_0)-{\mathcal{A}}^{\star }(x_0)}{2}+\frac{\beta }{3}·\end{array}$$(D.10)

For the computation of 𝒱^⋆(x₀) and 𝒜^⋆(x₀), we follow Takata (2016a) and use the assumptions (1) 𝒱^⋆ >  𝒜^⋆ and (2) 𝒱^⋆ + 𝒜^⋆ = 3. For sake of clarity, we justify both assumptions. First, the hypothesis (1) results from the fact that the equation is polytropic in a good approximation in the radiative zones of evolved stars so that the pressure and the density must satisfy p ∝ ρ^γ, with γ ∼ 1 − 4/3 (e.g., see B.2.1). From Eq. (D.7), we can therefore write 𝒜^⋆ ∼ 𝒱^⋆(Γ₁/γ − 1). For a totally ionized monatomic ideal gas, Γ₁ = 5/3 and we retrieve (1). Second, the hypothesis (2) comes from the assumption about the power-law profiles of 𝒩 and 𝒮₁. Indeed, such an hypothesis implies that J must be considered as uniform, so that according to Eq. (D.7) and (B.2), we retrieve (2). Noting also from Eqs. (8), (D.8) and the relation J = 2β/3 that

$$\begin{array}{c}\hfill {\left({\mathcal{V}}^{\star }{\mathcal{A}}^{\star }\right)}_{x_0}={\left(\frac{2J^2{\mathcal{N}}^2}{{\mathcal{S}}_1^2}\right)}_{x_0}=\frac{8}{9}{\beta }^2{\alpha }^2,\end{array}$$(D.11)

the use of the hypotheses (1) and (2) enables us to obtain 𝒜^⋆(x₀) by solving the quadratic equation

$$\begin{array}{c}\hfill {\mathcal{A}}^{\star }{(x_0)}^2-3{\mathcal{A}}^{\star }(x_0)+\frac{8}{9}{\beta }^2{\alpha }^2=0,\end{array}$$(D.12)

which gives

$$\begin{array}{c}\hfill {\mathcal{A}}^{\star }(x_0)=\frac{3}{2}-\beta \sqrt{\frac{9}{4{\beta }^2}-\frac{8{\alpha }^2}{9}}·\end{array}$$(D.13)

Finally, considering Eq. (D.13) and 𝒱^⋆ + 𝒜^⋆ = 3 in Eq. (D.10), and considering Eq. (D.9), we find that the gradient-term 𝒢₀ given in Eq. (D.1) reads in the simple considered model

$$\begin{array}{c}\hfill {\mathcal{G}}_0=\frac{3\pi }{4\sqrt{2}}\frac{ln\alpha }{(\alpha -1)}{(\frac{1}{ln\alpha }+\frac{\alpha }{1-\alpha }-\frac{1}{3}+\sqrt{\frac{9}{4{\beta }^2}-\frac{8{\alpha }^2}{9}})}^2.\end{array}$$(D.14)

In the limit of Δr → 0, that is, around α = 1, it is possible to check using a second-order Taylor expansion that Eq. (D.14) tends to

$$\begin{array}{c}\hfill {\mathcal{G}}_0(\alpha \to 1)=\frac{\pi }{12\sqrt{2}}{(\frac{5}{2}-\sqrt{\frac{81}{4{\beta }^2}-8})}^2,\end{array}$$(D.15)

which corresponds to the expression found by Takata (2016a). Besides, the extrapolation of 𝒢₀ outside of its derivation domain, that is, in the limit of α → 0, is equal according to Eq. (D.14) to

$$\begin{array}{c}\hfill {\mathcal{G}}_0(\alpha \to 0)\approx -\frac{3\pi }{4\sqrt{2}}{(\frac{3}{2\beta }-\frac{1}{3})}^{2}ln\alpha .\end{array}$$(D.16)

All Tables

All Figures