Mode Mixing and Rotational Splittings: II. Reconciling Different Approaches to Mode Coupling

TL;DR

This paper compares two theoretical approaches to mode coupling in mixed-mode asteroseismology, derives relations to reconcile them, and assesses how these impact measurements of stellar interior rotation.

Contribution

It provides the first formal reconciliation of asymptotic and matrix-based mode coupling descriptions, enabling better interpretation of mixed-mode observations.

Findings

Reconciled the asymptotic and matrix formulations of mode coupling.

Derived closed-form expressions for key quantities like the period-stretching function.

Identified limitations in current methods for estimating mode parameters in red giants.

Abstract

In the mixed-mode asteroseismology of subgiants and red giants, the coupling between the p- and g-mode cavities must be understood well in order to derive localised estimates of interior rotation from measurements of mode multiplet rotational splittings. There exist now two different descriptions of this coupling: one based on an asymptotic quantisation condition, and the other arising from coupling matrices associated with "acoustic molecular orbitals". We examine the analytic properties of both, and derive closed-form expressions for various quantities -- such as the period-stretching function $\tau$ -- which previously had to be solved for numerically. Using these, we reconcile both formulations for the first time, deriving relations by which quantities in each formulation may be translated to and interpreted within the other. This yields an information criterion for whether a given configuration of mixed modes meaningfully constrains the parameters of the asymptotic construction, which is likely not satisfied by the majority of first-ascent red giant stars in our observational sample. Since this construction has been already used to make rotational measurements of such red giants, we examine -- through a hare-and-hounds exercise -- whether, and how, such limitations affect existing measurements. While averaged estimates of core rotation seem fairly robust, template-matching using the asymptotic construction has difficulty reliably assigning rotational splittings to individual multiplets, or estimating mixing fractions $\zeta$ of the most p-dominated mixed modes, where such estimates are most needed. We finally discuss implications for extending the two-zone model of radial differential rotation, e.g. via rotational inversions, with these methods.

Tables

Table 1. Table 1: Maximum mixing fraction ζ 𝜁 \zeta measured for 1109 Kepler RGB stars, as a function of frequency, as displayed in red in Fig.   7 .

$\nu (\mu$Hz)	${\zeta }_{\max }$
44	0.995
54	0.997
64	0.997
74	0.996
84	0.995
94	0.994
104	0.993
114	0.992
124	0.989
134	0.990
144	0.988
154	0.988
164	0.986
174	0.984
184	0.982
194	0.981
204	0.979
214	0.975
224	0.972
234	0.970
244	0.962
254	0.948
264	0.766
274	0.876
284	0.862

Equations

$${\nu }_{n_{p}lm}\sim \Delta \nu \left(n_p+\frac{l}{2}+{ϵ}_{p,n_{p}lm}\right),$$

$$\frac{1}{{\nu }_{n_{g}lm}}\sim \Delta {\Pi }_l\left(n_g+{ϵ}_{g,n_{g}lm}\right),$$

$${\zeta }_{nlm}\sim \frac{{\int }_{\text{only radiative core}}|{\mathbf{\xi }}_{nlm}{|}^2\mathrm{d}m}{{\int }_{\text{whole star}}|{\mathbf{\xi }}_{nlm}{|}^2\mathrm{d}m},$$

$$\delta {\omega }_{nl}\sim {\omega }_{\text{rot}}-{\omega }_{\text{nonrot}}\sim m\left[{\beta }_g\zeta {\Omega }_{\text{core}}+{\beta }_p(1-\zeta ){\Omega }_{\text{env}}\right].$$

$$F(\nu )=\tan {\theta }_p(\nu )\cot {\theta }_g(\nu )-q(\nu )\text{ such that }F({\nu }_i)=0,$$

$${\theta }_p\mapsto {\Theta }_p\equiv \pi \frac{\nu -{\nu }_p}{\Delta \nu };{\theta }_g\mapsto {\Theta }_g\equiv \pi \frac{\frac{1}{{\nu }_g}-\frac{1}{\nu }}{\Delta {\Pi }_l}.$$

$$F(\nu )\approx \tan {\Theta }_p(\nu )\tan {\Theta }_g(\nu )-q(\nu )\text{ such that }F({\nu }_i)=0.$$

$${\theta }_p(\nu )=\pi \left(\frac{\nu }{\Delta \nu }-\frac{l}{2}-{ϵ}_p(\nu )\right);{\theta }_g(\nu )=\pi \left(\frac{1}{\nu \Delta {\Pi }_l}-{ϵ}_g(\nu )+\frac{1}{2}\right),$$

$${\Theta }_p(\nu ,{\nu }_p)={\theta }_p(\nu )-{\theta }_p({\nu }_p);{\Theta }_g(\nu ,{\nu }_g)={\theta }_g({\nu }_g)-{\theta }_g(\nu ).$$

$$\frac{\delta P}{\Delta {\Pi }_l}$$

$$\sim \frac{{\nu }_g^2}{{\nu }^2}{\frac{\partial \nu }{\partial {\nu }_g}|}_F=-\frac{{\nu }_g^2}{{\nu }^2}\frac{\partial F/\partial {\nu }_g{|}_{\nu }}{\partial F/\partial \nu {|}_{{\nu }_g}}$$

$$={\left[1+\frac{{\nu }^2\Delta {\Pi }_l}{\Delta \nu }\frac{\tan {\Theta }_g(\nu ,{\nu }_g){\cos }^2{\Theta }_g(\nu ,{\nu }_g)}{\tan {\Theta }_p(\nu ){\cos }^2{\Theta }_p(\nu )}\right]}^{-1}\equiv \zeta (\nu ,{\nu }_g).$$

$$\zeta ({\nu }_i,{\nu }_{g,i})={\zeta }_p({\nu }_i)\equiv {\left[1+\frac{\Delta {\Pi }_l}{\Delta \nu }\frac{{\nu }^2}{q{\cos }^2{\Theta }_p({\nu }_i)+\frac{1}{q}{\sin }^2{\Theta }_p({\nu }_i)}\right]}^{-1},$$

$$\frac{\delta P}{\Delta {\Pi }_l}\sim \zeta (\nu ,{\nu }_g),$$

$$\frac{\mathrm{d}\tau }{\mathrm{d}P}=\frac{1}{\zeta (1/P,{\nu }_g)},$$

$$\mathrm{d}F=0={\frac{\partial F}{\partial {\nu }_g}|}_{\nu }\mathrm{d}{\nu }_g+{\frac{\partial F}{\partial \nu }|}_{{\nu }_g}\mathrm{d}\nu ⟹\frac{\mathrm{d}\nu }{\mathrm{d}{\nu }_g}=-\frac{\partial F/\partial {\nu }_g{|}_{\nu }}{\partial F/\partial \nu {|}_{{\nu }_g}}=\frac{{\nu }^2}{{\nu }_g^2}\zeta (\nu ,{\nu }_g).$$

$$\frac{1}{{\nu }_g}=\tau (\nu )$$

$$\sim \frac{1}{\nu }+\frac{\Delta {\Pi }_l}{\pi }\arctan \left(\frac{q}{\tan {\Theta }_p(\nu )}\right)-\left(\nu -{\nu }_g\right)\frac{\partial {ϵ}_g}{\partial \nu }+\dots$$

$${\zeta }_p(\nu )\sim {\left(1+\frac{{\nu }^2\Delta \Pi }{q\Delta \nu }\frac{1}{1+\frac{{\pi }^2}{q^2(\Delta \nu {)}^2}(\nu -{\nu }_p{)}^2}\right)}^{-1};$$

$$D_{\text{samp}}=\frac{\gamma }{h}\sim \frac{2q\Delta \nu }{\pi {\nu }^2\Delta {\Pi }_l}=\frac{2}{\pi }q{\mathcal{N}}_l\gg 1,$$

$${\nu }_{+}-{\nu }_{-}$$

$$\to |{\nu }_p-{\nu }_{g,0}|\text{ as }q\to 0.$$

$${\nu }_{+}-{\nu }_{-}\ge \delta {\nu }_{\text{res, asy}}=\frac{2{\nu }_c}{\pi }\sqrt{q\Delta \nu \Delta {\Pi }_l}+\mathcal{O}(q^{3/2}).$$

$${\zeta }_{\pm }={\zeta }_p({\nu }_{\pm })={\left(1+\frac{1}{1+q{\mathcal{N}}_l}\right)}^{-1},$$

$$\left({\omega }^2\left[\begin{array}{cc}{\mathbb{I}}_p& \mathbf{D}\\ {\mathbf{D}}^{†}& {\mathbb{I}}_g\end{array}\right]+2\omega \mathbf{R}-\left[\begin{array}{cc}{\mathbf{\Omega }}_p^2& \mathbf{\alpha }\\ {\mathbf{\alpha }}^{†}& {\mathbf{\Omega }}_g^2\end{array}\right]\right)\left[\begin{array}{c}{\mathbf{c}}_p\\ {\mathbf{c}}_g\end{array}\right]=0,$$

$$\left[\begin{array}{cc}{\mathbf{\Omega }}_p^2& \mathbf{\alpha }\\ {\mathbf{\alpha }}^{†}& {\mathbf{\Omega }}_g^2\end{array}\right]\left[\begin{array}{c}{\mathbf{c}}_p\\ {\mathbf{c}}_g\end{array}\right]={\omega }^2\left[\begin{array}{cc}{\mathbb{I}}_p& \mathbf{D}\\ {\mathbf{D}}^{†}& {\mathbb{I}}_g\end{array}\right]\left[\begin{array}{c}{\mathbf{c}}_p\\ {\mathbf{c}}_g\end{array}\right],$$

$${\alpha }_{ij}=\int \mathrm{d}m{\mathbf{\xi }}_{\pi ,i}\cdot \left(N^2{\mathbf{e}}_r\otimes {\mathbf{e}}_r+{\omega }_{\pi ,i}^2\right)\cdot {\mathbf{\xi }}_{\gamma ,j},$$

$$\mathbf{\xi }\sim \sum _i{c}_{p,i}{\mathbf{\xi }}_{\pi ,i}+\sum _j{c}_{g,j}{\mathbf{\xi }}_{\gamma ,j},$$

$$A_{ij}\sim {\alpha }_{ij}-\frac{1}{2}\left({\omega }_i^2+{\omega }_j^2\right)D_{ij}+\mathcal{O}(D^2).$$

$$\left[\begin{array}{cc}{\omega }_p^2& \alpha \\ \alpha & {\omega }_g^2\end{array}\right]\left[\begin{array}{c}1\\ u_{\pm }\end{array}\right]\sim {\omega }_{\pm }^2\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}1\\ u_{\pm }\end{array}\right],$$

$${\omega }_{\pm }^2\sim \frac{{\omega }_p^2+{\omega }_g^2}{2}\pm \frac{1}{2}\sqrt{{\left({\omega }_p^2-{\omega }_g^2\right)}^2+4{\alpha }^2},$$

$$u_{\pm }=\frac{{\omega }_{\pm }^2-{\omega }_p^2}{\alpha }=\frac{\alpha }{{\omega }_{\pm }^2-{\omega }_g^2}.$$

Full text

Asteroseismology (73), Red giant stars (1372), Stellar oscillations (1617), Computational methods (1965), Theoretical techniques (2093)

††journal: ApJ

\makesavenoteenv

longtable

\newunicodecharʻ

‘

1 Introduction and Motivation

The asteroseismology of solar-like oscillators in the age of high-precision space photometry missions — such as CoRoT, Kepler and TESS — increasingly must contend with modes of mixed character, associated with solar-like oscillators which have evolved off the main sequence. In particular, whereas stars like the Sun can be seen to support pressure waves, whose p-mode frequencies approximately satisfy

[TABLE]

(where $\Delta\nu$ is a characteristic frequency spacing called the large frequency separation, and $n_{p}$, $l$, and $m$ are integer quantum numbers), the observed waves of these these mixed-mode oscillators also propagate as gravity waves in their interiors, trapped within their radiative cores. Notionally, these g-mode frequencies should satisfy

[TABLE]

where $\Delta\Pi_{l}$ is a characteristic period spacing. In practice, however, those modes which are observed possess mixed character, propagating buoyantly in the core and acoustically in the envelope. This has proven to be a mixed blessing: on one hand, these mixed modes yield highly sensitive observational probes into otherwise inaccessible stellar interiors. On the other, the coupling between the two mode cavities disrupts the simple comb-like structure of the mode frequencies specified by Eqs. 1 and 2, which may no longer be directly relied upon. The interpretation of these mixed-mode frequencies must thus differ significantly from that of pure p- or g-modes considered in isolation.

One crucial diagnostic quantity for these mixed modes is their mixing fraction, typically denoted as $\zeta$, which describes how much a mode under consideration resembles a pure g-mode (where $\zeta\to 1$) or a pure p-mode ($\zeta\to 0$). When the radiative core is much more compact than the star, this quantity may be found directly from the mixed-mode eigenfunctions as a ratio of partial mode inertiae:

[TABLE]

where $m(r)=\int_{0}^{r}4\pi r^{\prime 2}\rho\mathrm{d}r^{\prime}$ is the mass coordinate, and $\bm{\xi}$ is the Lagrangian displacement eigenfunction. In practice, the structure of these stars is not directly accessible, and so these mixing fractions must instead be inferred. The inferred mixing fractions are in turn used to interpret other seismic measurements. For instance, when these stars rotate, their dipole modes exhibit rotational splitting into triplets. To leading order in perturbation theory, the frequency widths of these triplets may be approximated using a simplified two-zone model of radial differential rotation (with only globally averaged core and envelope rotation rates $\Omega_{\text{core}}$ and $\Omega_{\text{env}}$) as (e.g. Goupil et al., 2013)

[TABLE]

Here $\beta_{p}$ and $\beta_{g}$ are sensitivity constants associated with pure p- and g-modes. This relation is best satisfied when the rotational splitting is small — i.e. in less evolved subgiant or red giant stars, or more generally for the most p-dominated mixed modes. The conversion between the observed rotational splittings and the interpreted core rotation rates in these cases thus depends critically on the inferred values of $\zeta$ (e.g. Deheuvels et al., 2017, 2020; Triana et al., 2017). As such, it is in these situations where accuracy in inferring the mixing fractions is paramount.

However, Eq. 4 requires that the rotational splittings be linear, which is not generally the case; in more evolved red giants a highly nonlinear phenomenology of dense avoided crossings emerges (e.g. Ouazzani et al., 2013), rendering first-order perturbation analysis no longer generally tenable. To treat these avoided crossings, more sophisticated techniques are required. There are now two distinct analytic treatments of this mode mixing by which such inferences may be made. Of older vintage, the asymptotic analysis of mixed-mode oscillators in the JWKB approximation (e.g. Shibahashi, 1979; Unno et al., 1989; Takata, 2016) relates mixed-mode wavefunctions to stellar structure in a local fashion, with the mixed-mode frequencies themselves emerging as roots of some characteristic eigenvalue equation. More recently, an algebraic approach has also been developed (e.g. Deheuvels & Michel, 2010; Ong & Basu, 2020; Ong et al., 2021b, 2022), where the pulsation problem is first modified to yield a nonorthogonal “decoupled” set of pure-p and pure-g-mode basis functions; the mixed-mode frequencies then emerge as the roots of the characteristic equation associated with a generalised Hermitian matrix eigenvalue problem, as represented in this nonorthogonal basis. Given such a set of basis functions, inferences of the mixing fraction $\zeta$ for each mode permit an effective p- and g-mode splitting to be assigned to that mode, potentially then allowing for generalisations to the two-zone model of differential rotation (i.e. rotational inversion, as in e.g. Schunker et al., 2016a, b) to be prosecuted for these evolved red giants (as attempted in e.g. Di Mauro et al., 2016, 2018; Fellay et al., 2021).

Associated with these two analytic constructions are two distinct approaches to operationalising them for use with actual observational data. Thus far, the matrix construction has only been used in the context of direct forward modelling against grids of evolutionary structural models (e.g. Ong et al., 2021a). Such forward modelling is very expensive, but at the end of the day produces guesses at the structure of the star by which Eq. 3 may be evaluated directly, obviating the need for further analytic approximations. By contrast, the analytic properties of the asymptotic approach, along with a generative parameterisation for the pure p- and pure g-mode frequencies, have been put together to yield a more flexible diagnostic construction — that of ‘stretched échelle diagrams’ (Mosser et al., 2015; Gehan et al., 2018; Mosser et al., 2018) — proposed for direct use with actual power spectra, and unencumbered by any requirement for stellar modelling (at the expense of needing to fit free parameters).

This generative property has rendered the asymptotic construction particularly compelling for use on highly evolved red giant stars, where JWKB analysis is assumed to hold good, and where in any case the evolutionary and algebraic calculations associated with the matrix methods become increasingly expensive. At the same time, however, the asymptotic approach is also not free of defects. For one, its free parameters are strictly phenomenological, and no complete interpretation of them has been supplied. Observationally, the deviations between its predictions and actual mixed-mode frequencies are largest in the neighbourhood of p-dominated mixed modes (e.g. Fig. 3 of Mosser et al., 2018) — precisely where its correctness is most necessary, since it is only these modes which are observed with any appreciable amplitude. Conventionally, this is regarded as a consequence of some approximations made in its underlying analytic construction. Since these outstanding issues afflict the p-dominated mixed modes most significantly, existing users of the stretched échelle diagram construction have elected instead to omit the p-dominated mixed modes entirely from their analysis.

In this series of papers, we attempt to extend the two-zone model of differential rotation in the presence of these avoided crossings, with the ultimate goal of performing rotational inversions. As discussed in Ong et al. (2022, hereafter Paper I), the notional pure g-modes in red giants are deep into the asymptotic regime, and thus yield little information beyond a single averaged core rotation rate, even collectively; accordingly, in this work we place a particular focus on the p-dominated mixed modes. In Section 2, we therefore critically reexamine the asymptotic construction, and analytically justify its applicability to these modes specifically. Since the matrix-based analytic approach discussed in Paper I is computationally expensive to operationalise, we also examine in this section the prospects for applying observational results derived from the asymptotic procedure, which is significantly cheaper, within the matrix construction. Even within our fortified asymptotic construction, however, our discussion also indicates that the stretching function of a given configuration of mixed modes may be meaningfully constrained only if a sampling-based information criterion is satisfied. In Section 3, we find that this criterion is not satisfied for much of the existing observational sample. Given that such stretching functions have already extensively been used to characterise stellar properties, we explore the implications of this result, via a hare-and-hounds exercise against one existing observational pipeline for rotational measurements. Finally, we summarise our findings, and discuss future directions, in Section 4.

2 Analytic Developments

We consider mixed modes arising from a dense set of pure g-modes coupling to a sparse set of p-modes. We will limit our attention to the g-modes closest to a particular p-mode with frequency $\nu_{p}$. We label the frequency of the g-mode closest to this p-mode as $\nu_{g,0}$, and those of the adjacent g-modes with integer indices $i$ as $\nu_{g,i}$. The coupling between these pure modes gives rise to a pair of mixed modes $\nu_{\pm}$ close to the on-resonance pair $\nu_{p}$ and $\nu_{g,0}$, as well as further mixed modes $\nu_{i}$. Critically, the frequency of the p-mode $\nu_{p}$ is completely independent of the properties of the g-modes $\nu_{g,i}$. We illustrate this configuration in Fig. 1.

In the presence of rotation, each mode splits into a $(2\ell+1)$-tuplet. The interactions between rotation and mixed-mode coupling are such that for fast enough rotation, the linear expression Eq. 4 may cease to hold for the mixed modes. However, in Paper I, we also demonstrate that the rotational perturbations to the underlying p- and g-modes may still be treated as linear in these cases. Moreover, the interactions between mode coupling and rotation are such that, to a good approximation, the components of the same $m$ from each multiplet may be considered to couple to each other independently of the other $m$. We thus restrict our following discussion to modes of the same $m$, on the understanding that the same calculations must be performed repeatedly for each $m$ to recover the full set of observed modes.

We now derive the basic properties of the mixed-mode frequencies with respect to both of the existing analytical frameworks we have previously described.

2.1 Asymptotic Construction

In the JWKB construction of (Unno et al., 1989), the frequency eigenvalues emerge from an eigenvalue equation associated with a characteristic function

[TABLE]

where $\theta_{p}$, $\theta_{g}$, and $q$ are frequency-dependent phase or tunnelling integrals over the stellar structure with respect to an appropriate JWKB dispersion relation. Note that the characteristic function $F(\nu)$ is not generally 0 away from the eigenvalues. Mosser et al. (2015), and subsequent works, neglect the frequency dependence of $q$, and replace these phase integrals with expressions involving the underlying p- and g-modes as

[TABLE]

These expressions are chosen so that $\Theta_{p}\to n\pi$ as $\nu\to\nu_{p}$ for some integer $n$, and likewise for $\Theta_{g}$ and $\nu_{g}$. Moreover, as multiplication by a constant does not change the roots of the characteristic function, we have chosen an overall sign for our description of it to yield a corresponding positive sign of $q$ in many of our subsequent expressions. Before we proceed, we must now break with the existing analysis in pointing out that this identification of $\theta_{g}$ with $\Theta_{g}$ results in inconsistent behaviour. Specifically, the frequencies of the most g-dominated mixed modes should be best described by those of the underlying g-modes: that is, $\nu_{i}\to\nu_{g,i}$ far from resonance with a p-mode (e.g. Christensen-Dalsgaard, 2012). For these modes, $\tan\Theta_{p}$ takes on very large absolute values, as $\Theta_{p}(\nu_{g})\to\pi\left(n+{1\over 2}\right)$. Accordingly, the eigenvalues of these g-dominated modes must satisfy $\cot\theta_{g}\to 0$ as $\nu\to\nu_{g}$. But $\Theta_{g}\to n\pi$ as $\nu\to\nu_{g}$, so $\cot\Theta_{g}=\tan\left({\pi\over 2}-\Theta_{g}\right)$ is also singular as $\nu\to\nu_{g}$, and therefore cannot satisfy Eq. 5 — a contradiction. This inconsistency is resolved by identifying $\theta_{g}$ with ${\pi\over 2}-\Theta_{g}$ instead, resulting in a phase offset of ${\pi\over 2}$ relative to almost all literature expressions involving $\Theta_{g}$ in our subsequent discussion. This phase offset can be motivated analytically, relative to the formulation of Unno et al. (1989), by either working outside of the Cowling approximation (as in Takata, 2016), or by using the Lagrangian pressure perturbation (rather than displacement) as the dynamical variable (as in Gough, 1993). Lindsay et al. (2022) have demonstrated that such a phase offset is indeed necessary for agreement between values for $\epsilon_{g}$ from the sample of Mosser et al. (2018), and from the most g-dominated mixed modes that emerge from explicit numerical frequency calculations made without the JWKB approximation. Such a correction also resolves the implicit inconsistency in Fig. 12 of Mosser et al. (2018), where the Cowling-approximation asymptotic value of Provost & Berthomieu (1986) — which does not suffer from this issue — appears not to exhibit the predicted offset of 1/2 caused by the Cowling approximation, relative to the values of Takata (2016) — made without the Cowling approximation, but also suffering from this issue. While the expressions of Mosser et al. (2018) may be alternatively (and charitably) interpreted as simply encoding a different convention for $\epsilon_{g}$ in mixed modes vs. pure g-modes, it would still be more desirable for us to, in the first place, choose a convention in which the values of $\epsilon_{g}$ from the pure g-mode frequencies were consistent with those computed directly from the most g-dominated mixed modes, as they should approach each other far from resonance (cf. Fig. 1). We thus rewrite the characteristic function Eq. 5 and the associated eigenvalue equation as

[TABLE]

The definitions of $\Theta_{p}$ and $\Theta_{g}$ given in Eq. 6 have so far been interpreted in a local sense as functions of a single variable $\nu$, where $\nu_{p}$ and $\nu_{g}$ are taken to be specified by the nearest p- and g-mode frequencies to $\nu$, and held fixed. However, we note that these definitions may be generalised to yield bona fide multivariable functions in a self-consistent fashion. In particular, Eqs. 1 and 2 suggest definitions of angular quantities as (Mosser et al., 2012)

[TABLE]

yielding eigenvalues where $\theta_{p}(\nu_{p})=n_{p}\pi$ and $\theta_{g}(\nu_{g})=\left(n_{g}+{1\over 2}\right)\pi$ at integers $n_{p}$ and $n_{g}$. Accordingly, Eqs. 6 and 7 can be seen to be consistent with the definitions

[TABLE]

With these modified definitions, the characteristic function may likewise be treated as a multivariable function in $\nu$, $\nu_{g}$, and $\nu_{p}$. Eq. 7 then specifies mixed-mode eigenvalues when $F(\nu,\nu_{g},\nu_{p})=0$ with $\nu_{p}$ and $\nu_{g}$ set to any combination of p- and g-mode eigenvalues, respectively. In our subsequent discussion, we hold $\nu_{p}$ fixed at just such a p-mode eigenvalue (it does not matter which), but permit $\nu_{g}$ to vary freely.

With this machinery in hand, we may now critically examine the assumptions going into the usual construction of stretched échelle diagrams. We first reframe the standard derivation of the asymptotic parameterisation, with respect to our analytic continuation of $\nu$ and $\nu_{g}$ as independent continuous variables. Far from resonance, the frequency of a given mixed mode $\nu_{i}$ can be well approximated as a function of the nearest associated g-mode frequency $\nu_{g,i}$ (with reference to Fig. 1). Since $F(\nu_{i},\nu_{g,i})=0$ for all mixed modes $\nu_{i}$ and g-modes $\nu_{g,i}$, the local spacing between mixed modes can be shown to approximately satisfy (Christensen-Dalsgaard, 2012)

[TABLE]

where we have used the identity$\left.\partial x\over\partial y\right|_{z}\left.\partial y\over\partial z\right|_{x}\left.\partial z\over\partial x\right|_{y}=-1$, and neglect the frequency dependence of $q$, and of the g-mode phase function $\epsilon_{g}$. The dependences on $\Theta_{g}$, and so $\nu_{g}$, may be further eliminated when this is evaluated at eigenvalues $\nu_{i}$, or more generally on level sets where $F(\nu,\nu_{g})=0$ (via Eq. 7), as

[TABLE]

which, fortuitously, does not suffer from the systematic offset of $\pi\over 2$ which we have discussed above. We are fortunate indeed that it is this expression that has found widespread practical application. Supposing that this function $\zeta_{p}$ were directly accessible, Goupil et al. (2013) find that in the JWKB approximation, its values, evaluated at the mixed mode frequencies $\nu_{i}$, also approximate the mixing fractions specified by inertia ratios, Eq. 3.

We now consider the existing construction of stretched period-échelle diagrams with respect to this analytical framework. Having followed through an analogous derivation of Eq. 10 (also following Christensen-Dalsgaard, 2012), Mosser et al. (2015) then observed that the finite-difference approximation,

[TABLE]

may be replaced by ansatz with an equality,

[TABLE]

such that the coordinate $\tau$ that is so defined accumulates increments of $\Delta\Pi_{l}$ when $P$ is integrated over the intervals between adjacent mixed-mode eigenvalues $1/\nu_{i}$, by construction. In their treatment, the g-mode frequency $\nu_{g}$ is held fixed, yielding highly oscillatory behaviour when the integrand is evaluated as a continuous function of $\nu$ (cf. their fig. 2). We note that this is inconsistent with the definition of $\zeta$ as a directional derivative along connected level sets of $F$, as can be seen from our Eq. 10. In a bid to eliminate this oscillatory behaviour (and better condition their numerical integration), they then noted, as above, that evaluating $\zeta(\nu_{i},\nu_{g})$ only at mixed-mode eigenvalues $\nu_{i}$ (i.e. on a disconnected level set of $F=0$ with constant $\nu_{g}$) yielded values of what appeared to be a far simpler function ($\zeta_{p}$; Eq. 11), which also produced the desired properties upon integration when inserted into Eq. 13. Thus, despite their general inequivalence on lines of constant $\nu_{g}$, $\zeta$ and $\zeta_{p}$ were identified with each other for the purposes of integrating Eq. 13: first by interpolation in Mosser et al. (2015), and then explicitly in Mosser et al. (2018).

We now proceed with an alternative derivation that does not suffer from these issues. First, we note that, upon analytic continuation, there is some redundancy in the definition of this $\tau$ coordinate, via Eqs. 10 and 13, which may not be immediately obvious. Since we are entitled to treat $F(\nu,\nu_{g})$ formally as a function of two variables, we have seen that $\zeta$ and $\zeta_{p}$ are equivalent only on level sets of $F$ on the $\nu$-$\nu_{g}$ plane where $F(\nu,\nu_{g})=0$. This being the case, Eq. 10 defines these level sets near the line of equality (where $\nu=\nu_{g}$) via an ordinary differential equation of the form $\mathrm{d}\nu/\mathrm{d}\nu_{g}=f(\nu,\nu_{g})$. Explicitly, we may parameterise these level sets in differential form as

[TABLE]

Eq. 13 is then recovered upon the change of variables to $\tau\equiv 1/\nu_{g}$ and $P=1/\nu$. The prescription of Mosser et al. (2015) is thus equivalent to producing integral curves of Eq. 14 in the $\nu$-$\nu_{g}$ plane from an initial value problem (IVP), which can be seen to yield connected level sets of $F=0$ only when passing through some fiducial point $(\nu_{i},\nu_{g,i})$ near the line of equality. Only on such an integral curve, we are now correctly entitled to replace $\zeta$ with $\zeta_{p}$ when performing this integration. It is extremely fortunate that Mosser et al. (2015) achieve the correct result through an even number of cancelling errors, so that observational results relying upon this construction remain valid, at least conceptually.

A discrete family of such connected level sets exists on the $\nu$-$\nu_{g}$ plane. To build graphical intuition for this construction, we show in Fig. 2 an example of the values of $F$ one might obtain on the $\nu$-$\nu_{g}$ plane near the line of equality, emphasising the loci of these level sets by saturating our colour map, for an illustrative configuration of these mixed modes (generated with constant values of $q$, $\epsilon_{p}$, and $\epsilon_{g}$). Each of the white loci are a different level set, $\nu=f_{n}(\nu_{g})$. In this restricted scenario, the level sets are periodic as $f_{n}(\nu_{g})=f_{n-k}\left([1/\nu_{g}+k\Delta\Pi_{l}]^{-1}\right)$, although this is not generally true. Morphologically, these level sets resemble diagonal lines parallel to the line of equality, undergoing avoided crossings with horizontal lines $\nu=\nu_{p}$ at the pure p-mode eigenvalues (dotted line).

Mixed-mode frequencies may be read directly off Fig. 2 by searching for intersections between these level sets (white loci) with vertical lines specified by the underlying g-mode frequencies (dashed lines) — i.e. by evaluating $f_{n}(\nu_{g})$ at g-mode eigenvalues $\nu_{g,i}$. Per Eq. 10, $\zeta$ may also be read off the diagram as the gradients of these level sets at such points of intersection which are nearest to the line of equality: that is to say, $\zeta_{i}=\zeta_{p}(\nu_{i})={\partial f_{n}\over\partial\nu_{g}}(\nu_{g}=\nu_{g,i})$, with $n$ chosen to yield the level set nearest the line of equality.

Having demonstrated that the integral curves generated by Eqs. 13 and 14 are indeed level sets of $F=0$, however, we can see immediately that integrating Eq. 13 to recover a level set is actually unnecessary, as Eq. 7 already specifies these level sets in closed form. In particular we have

[TABLE]

This is most useful near the line of equality where $\nu\sim\nu_{g}$, so that any frequency dependence of $\epsilon_{g}$ may be neglected. We note that similar expressions were derived in Mosser et al. (2018); Pinçon & Takata (2022), but were not recognised to be a closed-form expression for the stretching function $\tau$, owing to the previously discussed phase offset of $\pi/2$. Since such an expression exists, however, we recommend that it be used in future applications, in preference over numerical integration of Eq. 13. Because this expression is exact, and does not depend on the finite-difference approximation made by going from Eq. 12 to Eq. 13, it is in principle also applicable to even the most p-dominated mixed modes. As such, the poor performance of the stretched-échelle-diagram construction for these modes in particular (e.g. Fig. 3 of Mosser et al., 2018) cannot be attributed to the poor quality of such a finite-difference approximation, as is customary in existing work using it. Instead, this must be the result of either the stretching function potentially not having been correctly evaluated numerically by integration (which we explore in more detail in Section 4.1), or of an incorrect stretching function, correctly evaluated numerically, being used to describe a given set of mixed modes.

2.1.1 Constraining the stretching function

So far, our discussion has been generative, in the sense that we have described curves on the $\nu$-$\nu_{g}$ plane generated by $\zeta_{p}$, assumed to be known perfectly. However, in practice, neither the pure p- or g-mode eigenvalues, nor $\zeta$, are known: all of them have to be inferred solely from the mixed-moded eigenvalues $\nu_{i}$. In principle, this is an underdetermined problem, as the number of p- and g-mode eigenvalues is unchanged when they are coupled to yield mixed modes, but extra properties of the coupling (i.e. the coupling strength $q$) must additionally be inferred. For most red giants, however, a single effective value of $\Delta\Pi_{l}$ and $\epsilon_{g}$ may be assumed, considerably reducing the number of unknowns. Many prescriptions exist by which these properties may be derived, including autocorrelation-based techniques (e.g. Vrard et al., 2016; Mosser et al., 2017), parametric fits to the mixed-mode frequencies (e.g Jiang et al., 2020, 2022; Li et al., 2022; Deheuvels et al., 2023), or various constructions involving the mixed-mode period differences (e.g. Cunha et al., 2019; Dréau et al., 2020; Farnir et al., 2021). While much of the existing literature has focused on accurate determination of $\Delta\Pi_{l}$, we focus here on constraining the p-mode frequencies $\{\nu_{p,i}\}$ and the coupling strength $q$ instead.

The p-mode frequencies are notionally located at minima of the mixed-mode period differences, although procedures for recovering them built upon this (even incorporating frequency-dependent $q$ and more sophisticated fitting techniques, e.g. Dréau et al., 2020) yield results that systematically deviate from the pure p-modes derived from stellar models (Saunders et al., accepted to ApJ). Our construction provides a further constraint on their values. For the two most p-dominated mixed modes closest to resonance (i.e. $\nu_{\pm}$ of Fig. 1), $\zeta_{\pm}$ can be seen to be the derivatives, evaluated at the same value of $\nu_{g}=\nu_{g,0}$, of two different level sets adjacent to the line of equality. Conversely, perfect knowledge of $\tau$ would yield $\tau(1/\nu_{+})=\tau(1/\nu_{-})=1/\nu_{g,0}$. In a rotationally split multiplet, this criterion must be satisfied separately for each multiplet component indexed by $m$. In principle, Fig. 2 indicates that this is exact, even when the actual period difference between pure g-modes of adjacent order is not a constant $\Delta\Pi_{l}$ (e.g. due to curvature from buoyancy glitches, rotational splitting, or magnetic fields). However, an overall difference of some integer multiple of $\Delta\Pi_{l}$ may be introduced from alternative choices of branch for the arctan function when Eq. 15 is evaluated numerically. We propose this property to be used as a diagnostic check of correctness for whether $\tau$ has been properly constructed to match $\{\nu_{p,i}\}$ for any given set of mixed modes.

We now examine how observational data may also fail to constrain $q$, and thus $\zeta_{p}$. While we will focus on the direct fitting of period differences in our subsequent discussion (as it is more analytically tractable than assessing the properties of a least-squares fit or periodogram), we observe that constraining power on $q$ is always, in these red giants, obtained from the mixed modes closest to the underlying p-mode. Moreover, taking period differences does not substantively decrease the amount of information available from the mixed modes, so these considerations should be generally applicable irrespective of actual observational practice. To simplify matters, we will assume a priori that the pure p-mode frequencies have somehow already been perfectly determined. By inspection, estimating $\zeta(\nu)$ directly using the period differences, as done in Eq. 12, is equivalent to evaluating $f_{n}^{\prime}(\nu_{g})$ using a second-order finite-difference scheme, with step size $h=\nu_{g}^{2}\Delta\Pi_{l}/2$, at values of $\nu_{g}$ located exactly in between the g-mode eigenvalues $\left\{\nu_{g,i}\right\}$. For these to be representative of the derivatives evaluated at $\left\{\nu_{g,i}\right\}$, and thus $\zeta_{p}(\nu_{i})$, the curvature of these level sets must be heuristically “small”, in the sense that these step sizes $h$ must oversample such curvature. In the limit of weak coupling, the avoided crossings between adjacent $f_{n}$ become increasingly sharply curved over vanishingly small zones of avoidance (Fig. 3). As such, the length scale on this diagram for such oversampling must be set by the sizes of these zones of avoidance.

We find this characteristic scale by expanding the trigonometric terms appearing in Eq. 11 in the neighbourhood of $\nu_{p}$ (following Christensen-Dalsgaard, 2012) to obtain

[TABLE]

that is to say, in the neighbourhood of p-mode eigenvalues $\nu_{p}$, the function ${1\over\zeta_{p}(\nu_{p}+\delta\nu)}-1$ is extremely well approximated by a Lorentzian profile of the form ${A\over\delta\nu^{2}+\gamma^{2}}$, with a characteristic resonance width $\gamma=q\Delta\nu/\pi$. If the available mixed-mode period differences are to adequately oversample this width, we therefore require

[TABLE]

where $\mathcal{N}_{l}(\nu)=\Delta\nu/\nu^{2}\Delta\Pi_{l}$ is the relative density of g-modes to p-modes at frequency $\nu$. Correspondingly, when $D_{\text{samp}}\ll 1,$ we claim that no meaningful constraint on the coupling strength $q$ may be obtained from the mode frequencies alone.

To understand the basis for this claim, it is instructive to consider what happens when the sampling criterion Eq. 17 is violated, which happens in the limit of weak coupling as $q\to 0$. We first consider the case where a p-mode is maximally out of resonance with the available g-modes ($\sin\Theta_{g}(\nu_{p})=1$), which we illustrate in Fig. 4a. As $D_{\text{samp}}\to 0$ (and so likewise $q$), the g-dominated mixed mode frequencies tend asymptotically to those of the underlying g-modes, and likewise the p-dominated mixed mode frequencies to the underlying p-modes. We show both the true function $\zeta_{p}$ generating these mixed modes (coloured solid curves), as well as the values obtained using the finite-difference approximation, Eq. 12, with perfect knowledge of $\Delta\Pi_{l}$ (filled circles joined with dashed lines). The finite-difference approximation is such that these are assumed to sample $\zeta_{p}$ exactly in between the mixed-mode frequencies.

At $D_{\text{samp}}\gtrsim 1$, the finite-difference approximation can be seen to hold well. However, as $D_{\text{samp}}\to 0$, the two period differences closest to the p-mode frequency must tend to $\delta P\to\Delta\Pi_{l}/2$, sampled at $\nu_{\pm}\sim\nu_{p}\pm\nu_{p}^{2}\Delta\Pi_{l}/4$. Thus, were Eq. 12 taken at face value, one would assume to have sampled $\zeta={1\over 2}$ at locations which become increasingly discrepant from the true half-width of the negative peak (which instead also shrinks to [math] as $q\to 0$). Consequently, one would obtain a systematic overestimate of $q$, which persists across various methods used to estimate it. For example, were Eq. 11 or Eq. 16 to be fitted against the period differences, a larger value of $q$ than actually present would be required to reproduce the width of the feature implied by the finite-difference approximation. Alternatively, if the ratio $Z=\zeta_{\text{min}}/\zeta_{\text{max}}$ were to be used (as in e.g. Farnir et al., 2021), we see as $q\to 0$ that one would obtain from finite differences $Z\to{1\over 2}$, rather than $\to 0$ as required, again leading to a systematic overestimate.

Let us now also examine the near-resonance configuration shown in Fig. 1. Christensen-Dalsgaard (2012) derives several results for this configuration which we will find useful later. Adapting them to our notation, we have that the difference between the two closest mixed modes to resonance can be found as

[TABLE]

This separation is bounded from below by the case of exact resonance, where $\nu_{p}\to\nu_{g}\to\nu_{c}$, as

[TABLE]

Inserting this into Eq. 16 yields the mixing fraction of both modes,

[TABLE]

which can be seen to tend to ${1\over 2}$ in the limit of weak coupling.

We show in Fig. 4b the same quantities as in Fig. 4a, in this case of exact resonance, where $\nu_{g,0}\to\nu_{p}$. Again, we can see that while the finite-difference approximation holds well for $D_{\text{samp}}\gtrsim 1$, its limitations at low $q$ also conspire to give the appearance of stronger coupling than actually present, in much the same fashion as we have described in the out-of-resonance configuration. Nonetheless, these overestimates can be seen not to be as extreme. For example, Eq. 16 and Eq. 19 together indicate that while $Z$ will still be systematically overestimated from period differences, it does eventually tend to 0 as $q\to 0$.

As such, estimating $\zeta$ from only period differences will yield overestimates of $q$ in configurations of both exact resonance and exact antiresonance with the p-modes, when $D_{\text{samp}}$ is small. Supposing that intermediate proximity to resonance between the two configurations that we have considered above will yield intermediate overestimates of $q$, we make the stronger claim that when $D_{\text{samp}}\ngtr 1$, values of $q$ which may be inferred from only mixed-mode frequencies will be systematically higher than required to actually describe their mode coupling. The amount by which $q$ will be overestimated then depends on the extent to which the observationally available p- and g-modes are variously close to or out of resonance, as well as the value of $D_{\text{samp}}$ itself.

2.2 Algebraic Construction

We now examine the near-resonance configuration using the algebraic construction, in which the rotating mixed-mode frequencies are fully described by a quadratic Hermitian eigenvalue problem

[TABLE]

as written in block matrix form, where $\mathbb{I}_{p}$ and $\mathbb{I}_{g}$ are identity matrices of the same rank as the notional pure p- and g-mode subspaces under consideration; $\mathbf{\Omega}_{p}$ and $\mathbf{\Omega}_{g}$ are diagonal matrices containing the pure p- and g-mode angular frequencies; $\mathbf{R}$ specifies the rotational properties of the star, and $\bm{\alpha}$ and $\mathbf{D}$ are matrices describing the pairwise interactions between each combination of p- and g-modes. In Paper I we show that the rotation matrix $\mathbf{R}$ is well-approximated as being diagonal in the basis of isolated p- and g-modes. This being the case, one may approximate this problem in this basis by solving the simpler Generalised Hermitian Eigenvalue Problem separately for each combination of $m$ and $l$, in the form

[TABLE]

with the diagonal elements of $\mathbf{\Omega}_{i}$ suitably modified to account for the rotational perturbation associated with the corresponding $l$ and $m$.

While this construction is entirely general, numerical prescriptions for evaluating these matrix elements have been provided in Ong & Basu (2020) as overlap integrals between notional “$\pi$” and “$\gamma$” modes (in the sense of Aizenman et al., 1977). In particular, given their prescription for isolating such $\pi$ and $\gamma$ modes, the coupling matrix element for the $i$th p-mode coupling to the $j$th g-mode can be written as

[TABLE]

where $\mathbf{e}_{r}$ is the unit vector in the radial direction, $\otimes$ denotes the tensor product, $\bm{\xi}_{k}$ is the Lagrangian perturbation eigenfunction of the $k$th mode, and $N^{2}$ is the (squared) buoyancy frequency. Within this construction, the mixed-mode frequencies are found as the eigenvalues of Eq. 22. The eigenfunctions of the resulting mixed modes can be written as linear combinations of these $\pi$ and $\gamma$ modes:

[TABLE]

where these coefficients $c$ are specified by the eigenvectors corresponding to each of these eigenvalues.

For these evolved red giants in particular, two simplifying approximations are typically made in this analysis. The first is to reduce Eq. 22 to a conventional Hermitian eigenvalue problem. This may be done (following Lavely & Ritzwoller, 1992) by observing that all of the $|D_{ij}|\ll 1$, and thus expanding the equivalent Hermitian eigenvalue problem in powers of $D$ to yield an effective coupling strength

[TABLE]

We truncate this expansion to zeroth order in our analytic development here (as in e.g. Deheuvels & Michel, 2010) for simplicity of notation, but will use the first-order correction in our later numerical calculations.

The second approximation is that the coupling strengths are sufficiently weak that each g-dominated mixed mode may be assumed to result from the two-mode coupling between a single g-mode and a single p-mode. Following Ong & Basu (2020), Eq. 22 may be written in that case as

[TABLE]

where the frequency eigenvalues are found as

[TABLE]

with corresponding eigenvectors specified by

[TABLE]

These eigenvalue components can be used to define mixing fractions $\zeta$ as

[TABLE]

Ong & Basu (2020) demonstrate that this is equivalent to the original definition of mixing fractions as inertia ratios, in that it reduces to Eq. 3 in the JWKB approximation. However, this construction is also generally applicable to mixed modes outside of the JWKB limit.

Of these two solutions, one of them will be closer in frequency to the g-mode than the other; it will have frequency $\omega_{+}$ if $\omega_{g}>\omega_{p}$, and vice versa. Far from resonance, this frequency, which we denote by $\omega$, again can be thought of as a function of the associated g-mode frequency $\omega_{g}$ (Fig. 1). We may verify, as above, that

[TABLE]

(cf. Appendix B of Ong & Basu, 2020, and Appendix A). Since this property is also shared by the continuous construction resulting from asymptotic analysis, we conclude, conversely, that such an asymptotic construction operates best in the same regime as where three- and higher-mode coupling may be ignored algebraically. We formalise this statement more rigorously in Section 2.3.

Again, we evaluate properties of the resonant pair. We note that

[TABLE]

Moreover,

[TABLE]

Eliminating the coupling strength $\alpha$ in favour of $\zeta$ thus yields

[TABLE]

Note that it does not matter whether $\zeta_{+}$ or $\zeta_{-}$ (associated with the lower or higher root) are used in this expression, since $\zeta_{+}=1-\zeta_{-}$ by orthogonality. Thus,

[TABLE]

again in violation of Eq. 12. Finally, we once again consider the case where $\nu_{p}$ and $\nu_{g}$ are exactly in resonance at frequency $\nu_{c}$. In that case, $\zeta=1/2$ for both modes (compare Eq. 20), and Eq. 34 becomes singular. Rather, we obtain from Eq. 32 that

[TABLE]

2.3 A Weak-Coupling Limit for the Nonasymptotic Mixed-mode Eigenvalue Equation

While both the nonasymptotic algebraic and asymptotic constructions appear to yield good descriptions of the phenomenology of the observed mixed modes, little discussion has taken place regarding how the two are related. A better understanding of how the two formulations are related would both significantly allieviate some known computational and conceptual difficulties, as well as ease the interpretation of existing observational results in comparison to numerical models of stellar structure. In particular, the matrix construction becomes increasingly unwieldy to apply for evolved red giants, where the relevant coupling matrices become very large (Ong et al., 2021b); numerical solutions to the matrix eigenvalue problem, which scale in time complexity with the rank $N$ of these matrices as $\mathcal{O}(N^{3})$, become prohibitively expensive. At the same time, despite extensive observational characterisation, the asymptotic coupling strength is difficult to relate directly to mode frequencies obtained from numerical solutions to the pulsation problem from stellar modelling. Should some equivalence between the two formulations be found, it would allow measurements of the quantities of one formulation to be interpreted and used within the other.

The quantity $q$ carries an interpretation of a product of transmission coefficients in the JWKB limit of the matrix construction (Ong & Basu, 2020, §3.3), but this is not easy to relate to other quantities in the algebraic construction, which primarily depends on global, rather than local, properties of the wavefunctions. However, we note that the mixing fractions and the mode frequencies of both constructions tend to each other in the limit $q\to 0$. If we further demand that their avoided crossings have the same shape at both small $q$ and small $\alpha$, then we may also identify the minimum separations of the avoided crossings, Eqs. 19 and 35, that appear in both constructions when the underlying p- and g-modes are exactly in resonance. Doing so then allows us to effect a translation between the coupling strengths of the algebraic and asymptotic constructions:

[TABLE]

The oversampling criterion, Eq. 17, may therefore be rewritten in terms of $\alpha$ as

[TABLE]

Building on this result, we will now go further and derive a weak-coupling limit for the eigenvalue equation of the algebraic construction, and show that it passes naturally to that of the asymptotic construction. At the heart of this discussion is the asymptotic eigenvalue equation, Eq. 7, which we will rewrite as

[TABLE]

so as to avoid having any singular points, while retaining the same roots. The nonasymptotic analogue to this is the characteristic equation of the Hermitian eigenvalue problem,

[TABLE]

Without loss of generality, we may obtain an equivalent expression with the same roots by taking the determinant of some other matrix $\mathbf{PLQ}$ of equal rank, so long as $\det\mathbf{P}(\omega)$ and $\det\mathbf{Q}(\omega)$ have no roots for $\omega>0$. This being the case, we choose instead as our characteristic equation that

[TABLE]

Let us first restrict our attention to the diagonal blocks of this matrix, which are both also diagonal matrices. The determinants of these blocks are thus products of the diagonal entries: the determinant of the upper left block takes the form

[TABLE]

while that of the lower right block is of the form

[TABLE]

where we have used the asymptotic relations, Eqs. 1 and 2.

We now recall that analytic functions admit Weierstrass factorisations as infinite products over their roots. In particular, the trigonometric functions may be factorised as

[TABLE]

which is clearly very similar in structure to Eqs. 41 and 42. Motivated by this, we define in trigonometric form the characteristic functions

[TABLE]

for the p- and g-mode cavities separately. Note that such trigonometric definitions are in principle correct only up to multiplication by some functions, $f_{p}(\omega)$ and $f_{g}(\omega)$, having no zeros for $\omega>0$. However, we may always choose a different transformation to Eq. 39 (rather than Eq. 40), scaling each block by some appropriate function of frequency, without changing the roots of these characteristic functions. In particular we may choose such a transformation that that sets these functions to 1. In turn, our subsequent expressions turn out not to depend on the precise transformation chosen. We therefore may assume $f_{p}$ and $f_{g}$ both to be 1 without loss of correctness.

Let us now consider Eq. 40 again, to which we apply the identity

[TABLE]

and expand the final multiplicand in powers of the coupling matrices as

[TABLE]

In this fashion we obtain

[TABLE]

where for compactness we write $\alpha_{ij}-\omega^{2}D_{ij}\equiv A_{ij}$. To simplify the sum, we consider the derivatives of Eqs. 41 and 42 with respect to $\omega$. By the product rule, this yields

[TABLE]

Accordingly, we rewrite the sum in Eq. 47 in the fairly suggestive form

[TABLE]

So far we have made no assumptions about the structure of the matrices $\bm{\alpha}$ and $\mathbf{D}$. At this point, however, let suppose that the quantity

[TABLE]

which we have defined by adapting Eq. 36, may be treated as being numerically independent of its matrix indices. In that case it can be brought out of the sum, which collapses to yield

[TABLE]

We can now see that identifying $Q$ with $q$, and substituting the trigonometric expressions Eq. 44 into Eq. 51, does indeed cause the algebraic construction to yield the asymptotic eigenvalue equation, Eq. 38. Moreover, our explicit construction allows us to make more precise statements regarding the conditions under which the two formulations become equivalent. In particular, we require that (1) the elements of the coupling matrices have no parametric dependence on the p-mode frequencies, and depend on the g-mode frequencies as $(\alpha_{ij}-\omega^{2}D_{ij})\propto\omega_{g,j}$; and that (2) we may neglect terms to fourth and higher order in the coupling matrices (i.e. the coupling between mode cavities is weak). In general, however, it is known that condition (1) is not satisfied by stellar models: the coupling matrix elements can be shown at least numerically to exhibit dependence on both the p- and g-mode indices, and their dependence on the g-mode frequency is not necessarily of the required form (as we will examine below). In a local sense, an effective value of $q$ may be obtained using Eq. 36 from the on-resonance matrix element corresponding to each avoided crossing considered in isolation, yielding frequency dependence of the kind considered in e.g. Farnir et al. (2021); Jiang et al. (2022). However, absent some conspiracy of stellar structure, the two formulations are globally inequivalent.

3 Numerical Results

We have now examined how the finite-difference approximation Eq. 12 breaks down when the avoided crossings are not adequately oversampled. We quantify this using an oversampling discriminant, expressed with respect to the coupling strength in both the asymptotic (Eq. 17) and algebraic (Eq. 37) formulations. When this oversampling requirement is not satisfied, the coupling strengths $q$ may be systematically overestimated. The mixing fraction $\zeta$ — inferred values of which depend on $q$ — is itself an observationally significant diagnostic quantity in its own right, and errors in estimating it will in turn contaminate the interpretation of other seismic measurements, such as rotation rates (e.g. via Eq. 4). As such, it is important to understand the extent to which existing estimates of it, and by extension quantities derived from it, remain reliable in the face of these systematic issues.

3.1 Stellar Modelling

We first assess this qualitatively, by examining how well the oversampling criterion is satisfied numerically by stellar models produced with evolutionary calculations. For the purposes of our following discussion, we will use the relative density of g-modes to p-modes at ${\nu_{\text{max}}}$, $\mathcal{N}_{l}({\nu_{\text{max}}})=\Delta\nu/{\nu_{\text{max}}}^{2}\Delta\Pi_{l}$, as an indicator of evolution up the red giant branch. In particular, we focus on the region $\mathcal{N}_{l=1}({\nu_{\text{max}}})<30$, which encompasses essentially the entirety of the current observational sample of first-ascent red giants for which measurements of $\zeta$ using asymptotic techniques have been attempted (e.g. Mosser et al., 2018; Gehan et al., 2018).

Since we seek to assess whether or not the oversampling criterion Eq. 17 holds at all, it would be circular to evaluate it with respect to some phenomenological $q$ as inferred from model mixed modes, as we have earlier established that we may only meaningfully do so if this criterion were to be already satisfied in the first place. Moreover, the asymptotic construction in which $q$ appears is itself strictly speaking not universally applicable to mixed-mode oscillators in general, as the eigenvalue equation from which it is derived, Eq. 5, further relies on the JWKB approximation, which serves as a good description only at high wavenumber (in this case, at high radial order). The JWKB approximation is known not to hold for both subgiants and early red giants at low $\mathcal{N}_{l}$ (where the underlying g-modes are of low degree: cf. Bedding, 2012), or for very evolved red giants at high $\mathcal{N}_{l}$ (where the underlying p-modes are of low degree: cf. Stello et al., 2014). Since well-defined bounds for $\mathcal{N}_{l}$ in which the JWKB approximation may be applied are not generally available, we cannot necessarily guarantee that it holds good within our evolutionary regime of interest. Obversely, our assessment of whether or not the oversampling criterion holds should preferably be independent of the validity of the JWKB approximation per se. For these reasons, we will instead evaluate the oversampling discriminant using the nonasymptotic expression, Eq. 37.

We evaluate our quantities with respect to stellar models in a set of evolutionary tracks produced by the stellar evolution code mesa r15140 (Paxton et al., 2011, 2013, 2015, 2018, 2019), constructed in the same fashion as those considered in Paper I. For each model under consideration, we compute the off-diagonal matrix element corresponding to the resonant pair of dipole $\pi$- and $\gamma$-modes closest in frequency to ${\nu_{\text{max}}}$. We do so using the pulsation code gyre (Townsend & Teitler, 2013) and Eq. 23.

Conventional methods by which $q$ is inferred from mixed modes necessitate either specifying some parametric form of $q$ as a function of frequency, or, more commonly, ignoring its frequency dependence altogether. By evaluating it directly from stellar structure using Eq. 36 instead, we are in a position to examine this frequency dependence without making any such constraining assumptions. We show these effective coupling strengths in Fig. 5 as a function of the normalised frequency $\nu/\Delta\nu$, along an evolutionary track of solar composition and $M=1.4M_{\odot}$, for modes within $\pm 5\Delta\nu$ of ${\nu_{\text{max}}}$. These values are evaluated as a function of $\nu_{c}$ (i.e. the central frequencies of pairs of g- and p-modes closest to resonance); such resonant-pair frequencies are located close to the p-mode frequencies for evolved red giants, and close to the g-mode frequencies for subgiants. Values for each stellar model are coloured by $\mathcal{N}_{1}$, evaluated at ${\nu_{\text{max}}}$ for each stellar model.

Unlike the constant, linear, or quadratic dependences on frequency that are typically assumed (e.g. Cunha et al., 2019; Dréau et al., 2020), the actual variation of the coupling strength with frequency appears quite complicated, and also exhibits further nontrivial morphological variations depending on evolutionary stage. Specifically, we see that it exhibits significant, and seemingly oscillatory, curvature in our subgiant models (which will require more complicated analysis than previously considered, e.g. Pinçon et al., 2020, to explain; although the applicability of the JWKB approximation to subgiants is perhaps questionable). For more evolved red giants it clearly also varies significantly, increasing by more than a factor of 2 between the lowest and highest p-modes within our frequency range of interest — this is consistent with the phenomenological findings of Farnir et al. (2021), performed on a mode-by-mode basis. At the same time, the functional dependence on frequency appears not to adhere well to the power-law parameterisation suggested by the JWKB analysis of Pinçon et al. (2020), which would be incapable of producing inflection points as exhibited by the curve at $\mathcal{N}\sim 20$.

We show in Fig. 6a the value of $q({\nu_{\text{max}}})$ as a representative value for each stellar model, over the course of stellar evolution, for a series of evolutionary tracks ranging from $1M_{\odot}$ to $2M_{\odot}$. We then use it to evaluate the sampling discriminant, Eq. 37, and show the values of this discriminant in Fig. 6b, following the evolution of these quantities up to $\mathcal{N}_{1}({\nu_{\text{max}}})\sim 30$. We find that the various evolutionary tracks exhibit fairly similar morphology in the evolution of this coupling strength. After exhibiting relatively strong mode coupling as subgiants with primarily p-dominated mixed modes, the effective coupling strengths of our model stars level off with increasing $\mathcal{N}_{1}\lesssim 15$, before again falling off precipitously. Correspondingly, the sampling discriminant appears to grow linearly with $\mathcal{N}$ until this turnoff point, after which it decreases sharply. Within the range of $\mathcal{N}_{l}({\nu_{\text{max}}})$ and stellar masses that we have considered, the sampling discriminant never takes values much larger than unity for these tracks. On the contrary, it is significantly smaller than unity for most of our mesa models within this range: the oversampling criteria, Eqs. 17 and 37, are generally speaking not satisfied by mixed modes near ${\nu_{\text{max}}}$, except for higher masses at values of $N_{1}({\nu_{\text{max}}})$ near the upper end of the observable range. We note that this discriminant does exceed unity near the turnoff point in the evolution of the coupling strengths for stars massive enough to host non-degenerate helium cores, and this transition from core electron degeneracy is known to induce other seismic signatures (e.g Deheuvels et al., 2022). However, it is difficult to say for certain from this cursory analysis whether this feature is of genuine astrophysical origin or merely an interesting numerical coincidence in our stellar modelling, and a detailed investigation of this lies beyond the scope of this work.

Since the oversampling criterion $D_{\text{samp}}>1$ is not satisfied by our numerical models, we accordingly conclude that the actual combinations of frequency spacings, period spacings, and coupling strengths associated with first-ascent red giants in this region of parameter space are quite likely insufficiently informative to unambiguously constrain the mixing fractions $\zeta$ of near-resonance mixed modes without further input, e.g. from explicit stellar modelling, or strong priors on $q$. Following our earlier discussion in Section 2.1.1, attempts to infer coupling strengths under these conditions would yield values that are systematic overestimates. To assess the potential impact of this, we also show in both panels of Fig. 6, with gray circles, the values of $q$ (6a), and the implied oversampling discriminant (6b), of first-ascent red giants from the catalogue of coupling strengths published in Mosser et al. (2017). We see that the values of $q$ returned from our numerical mesa models form a lower envelope to most of the reported data points. Even despite being higher than our numerical predictions, the majority of the corresponding values of $D_{\text{samp}}$ can still be seen to be less than 1. These features are entirely consistent with the interpretation supplied in Section 2.1.1, of observationally inferred values of $q$ being systematically overestimated when $D_{\text{samp}}<1$. Supposing that such an interpretation were to be correct, it would imply that the true coupling strengths of these mixed modes would largely be lower than reported, resulting in true values of the oversampling discriminant which would be also lower than those shown here.

3.2 Hare-and-Hounds Exercise

The primary applications of the stretched-échelle-diagram construction have been to make rotational measurements of first-ascent red giants in the regime where the linear expression, Eq. 4, does not hold well. However, measurements made using these diagrams implicitly assume that the values of $q$ used to make them are consistent with those actually generating the mixed modes being stretched. Our prior discussion suggests that this may not be the case for a large fraction of our observational sample. Thus, we now seek to quantify the performance of these existing methods for rotational characterisation that use the stretched-échelle construction, in the face of these systematic issues. In keeping with the overall direction of this series of papers, we also examine prospects for adapting such measurements for recovering quantities required for rotational inversions.

We do this by way of a hare-and-hounds exercise. For our hares, we use rotating mesa evolutionary models, with angular momentum transport following the prescription of Fuller et al. (2019), and the same initial rotational configuration at ZAMS as the evolutionary tracks used in Paper I, generated at initial masses of 1.2 through 1.8 $M_{\odot}$ at a spacing of $0.2M_{\odot}$. From these models we compute the dipole mixed-mode frequencies, first without isolation into the basis of $\pi$ and $\gamma$ modes, at all azimuthal orders. Additionally, we also then evaluate the associated $\pi$ and $\gamma$-mode eigensystems, thereby obtaining a separate set of linearised rotational-splitting estimates. For each mixed mode, $1-\zeta$ is evaluated by taking the inverse of the ratio of the mode inertia of the mixed mode under consideration, and the value obtained by interpolating the mode inertiae of the $\ell=1$ $\pi$-modes to the specified mixed-mode frequency:

[TABLE]

We select models with $\mathcal{N}_{1}$ taking values between 2 and 30 from each of these evolutionary track as our primary set of hares, intended to investigate mass dependences of these methodological systematics; we will refer to these as Set I. In addition, we also generated one evolutionary track at $1.4M_{\odot}$ with the ZAMS rotation rate being half that of the same track in Set I, and another with double the ZAMS rotation rate, as an additional pair of sets of hares, from which we selected models and generated frequencies and rotation estimators in the same fashion. We will refer to these hares as Set II.

As our hound, we use the template-matching rotational measurement techniques described in Gehan et al. (2018), which are built upon the stretched-échelle-diagram construction. These techniques require estimates of the mixed-mode spectrum (in particular, the g-mode period spacing and phase offset, as well as the coupling strength between the mode cavities) to be specified a priori. Rather than simulate an entire fitting and peakbagging pipeline for these quantities, we mimic both internal systematic and statistical errors on them by sampling posterior distributions of $\Delta\Pi_{1}$, $\epsilon_{g}$ and constant $q$. We treat as (simulated) observational data the computed mixed $m=0$ modes (with a nominal frequency measurement uncertainty of 0.005 $\mu$Hz), and evaluate likelihood functions for a generative model for the mixed-mode frequencies (specified via Eq. 7) that accepts these three quantities as input parameters, in conjunction with the $\pi$-mode frequencies at $m=0$ (held fixed with no noise for this procedure). Samples from the posterior distributions with uniform sampling were then constructed using emcee. This procedure was used to supply the hound with $q$, $\Delta\Pi_{1}$, and $\epsilon_{g}$, drawn from the joint posterior distribution, as well as nominal measurement errors, estimated from the $1\sigma$ quantiles of the marginal posterior distributions. This was done independently for each hare. The hound was also supplied with the rotating mixed-mode frequencies, with simulated measurement errors as above, as well as the nonrotating $\pi$-mode frequencies, whose simulated measurement errors were also inflated by the envelope rotation rates of the stellar models, which are ordinarily ignored in this procedure.

We recall that in the first-order approximation, mode frequencies in a rotating star relate to those in a non-rotating configuration through Eq. 4. As such, relating the observed quantities to structural stellar models on a multiplet-by-multiplet basis, as necessary for e.g. rotational inversions, will require accurate determination of $m$, $\zeta$, and $\delta\omega_{\text{rot}}$. While current methods adapted for red giants (including that of Gehan et al., 2018) are able to estimate the first two quantities, they derive only a single aggregate value of $\Omega_{\text{core}}$ from template-matching against stretched échelle diagrams. Nonetheless, we will now assess the robustness of these methods in recovering these quantities when confronted with our hares.

3.2.1 Inference of $\zeta$

The method developed by Gehan et al. (2018) enables the measurement of the mean core rotation rate, which is probed by gravity modes. Hence, further methodological considerations for preparing the hares described here are focused on the gravity-dominated (g-m) dipole mixed modes, per analysis as described in Gehan et al. (2018). Additionally, the operation of the hound’s procedure specifically requires the most p-dominated mixed modes to be excluded from the construction of the stretched échelle diagram, given the conventionally assumed inapplicability of the technique to the most p-dominated modes (Mosser et al., 2015; Gehan et al., 2018). Nonetheless, we include them in the frequency list provided to the hound, and then consider the consistency of their properties as inferred by the hound’s constraints derived from the g-dominated mixed modes.

The amplitude of a mode in the power spectrum depends strongly on $\zeta$: pure g-modes with $\zeta\to 1$ will have no intrinsic photometric amplitude, and would not be observable. Furthermore, the most g-dominated dipole mixed modes would be located in between the dipole p-modes, and therefore would be drowned out by the radial and quadrupole p-modes in the power spectrum. Thus, we must include only those modes in our hare-and-hounds exercise as might be plausible to observe even in ideal conditions.

For this purpose, we used values for the maximal mixing fractions $\zeta$ that are consistent with observations. Fig. 7 represents the measured $\zeta$ values for all mixed modes detected over 1109 Kepler RGB stars analyzed by Gehan et al. (2018), as a function of frequency. We withheld from the hound the dipole mixed modes in our models which have $\zeta$ values above the red envelope displayed in Fig. 7, for each given frequency range (see Table 1 in Appendix B). This serves as a representation of what would be available under ideal observing conditions. However, the most g-dominated mixed modes are also, correspondingly, the most susceptible to degradation under less ideal observing conditions than the best possible ones provide by Kepler. To simulate the effects of progressive degradation of observing conditions, we prepared a third tranche of hares, constructed by taking hares from the evolutionary track with $M=1.4M_{\odot}$ from Set I, applying the frequency-dependent cutoff from Kepler as above, and then further excluding modes with $\zeta$ above some threshold value. For this purpose we adopt threshold values of the maximal mixing fractions as $\zeta_{\mathrm{max}}=[1-10^{-5},1-10^{-4},1-10^{-3},1-10^{-2},1-10^{-1}]$. We will refer to these hares as Set III.

We show in Fig. 8 the relative differences in the values of $\zeta$ inferred by the hound, compared to the ground-truth values available from the hare, for Set III of the hares. The overall qualitative trend is identical for the other two sets of hares, which we omit from discussion as they exhibit no visible dependence on either the mass or rotational configuration. We can see that the relative differences between the inferred and true values of $\zeta$ depend strongly on $\zeta$ itself; these determinations become increasingly unreliable as $\zeta\to 0$ — i.e. for the most p-dominated mixed modes. Unfortunately, it is for these modes where determination of $\zeta$ is most needed, compared to the modes of the most g-like character, which to a good approximation have $\zeta\to 1$. Moreover, we see no significant trend in this pattern as more modes become available to the hypothetical observer. This suggests that this trend is intrinsic to the method by which $\zeta$ is inferred by the hound, rather than being induced by a lack of available g-dominated mixed modes with which to constrain the mixed-mode multiplet configuration. We finally note that while the algorithm employed by the hound automatically rejects modes which it deems too p-dominated to be amenable to reliable analysis, there is no clear difference in how accurately $\zeta$ is inferred between modes that have been accepted (filled circles) vs. those that have been rejected (open triangles) by the hound.

3.2.2 Identification of the azimuthal order of dipole mixed modes

We also check the consistency of the identification of the azimuthal order of dipole g-dominated mixed modes based on the method developed by Gehan et al. (2018), since it can impact the measurement of the mean core rotation rate and, to a larger extent, the measurement of the inclination angle of the rotation axis (Gehan et al., 2021). Fig. 9 (a) presents an example for a 1.2 $M_{\mathrm{\odot}}$ model with a mixed-mode density $\mathcal{N}_{1}=$ 25, for which the identification of the azimuthal order of dipole mixed modes differs significantly between the hare and the hound. Such misidentification essentially occurs near crossings of the rotational components of different azimuthal orders in the stretched échelle diagram, corresponding to local g-mode frequency spacings congruent with integer multiples of the rotational splitting, and are more susceptible to occur both when the mean core rotation rate is large (leading to larger rotational splittings), or for more evolved stars (as the g-mode spacing decreases with evolution along the RGB). This is what can be seen in Fig. 9 (b) for Set I, where the fraction of modes correctly identified by the hound decreases along evolution, for the four masses considered in this exercise.

Gehan et al. (2018, 2021) used this template-matching method to identify the azimuthal order of g-m mixed modes, and from these derive inclination measurements in an automated fashion. Their inclination measurements are consistent with earlier manual measurements, all of them for stars with $\mathcal{N}_{1}\leq$ 15. While Fig. 9(b) suggests that some of these inclination measurements ought to be revisited, we also see that $\sim$ 80 % of the modes, or more, are correctly identified for our hares with $\mathcal{N}_{1}\leq$ 15. 86.2 % — a substantial majority — of the sample of Gehan et al. (2021) lie within this range of evolution. Moreover, planet-hosting red giants in particular, for which such inclination measurements are of crucial astrophysical interest for understanding their dynamical architectures, live low on the red giant branch — i.e. at lower values of $\mathcal{N}_{1}$ — since significant expansion of their host stars over the course of evolution up the RGB results in planetary engulfment.

At the same time, it may be instructive to consider the various reasons behind the increasing misidentification of azimuthal orders for $N_{\mathrm{1}}\gtrsim 15$. One possible reason is structural in origin, in the sense that the mixing physics adopted in generating the hares results in stellar models exhibiting sharp, localised variations of the Brunt-Väisälä frequency. These in turn cause buoyancy glitches: quasiperiodic modulations of the observed frequency of the g-dominated mixed modes (Cunha et al., 2015; Mosser et al., 2015; Cunha et al., 2019). In our models, one of these glitches corresponds to the hydrogen-burning shell, which is located at some depth in the radiative zone. This glitch produces modulations of low amplitudes, and does not significantly interact with the procedure by which the hare’s modes are matched and aligned with the hound’s templates. However, a second glitch is also present in our models, corresponding to the first dredge-up event on the low red giant branch, when the convective envelope reaches regions were nuclear burning previously took place (Salaris et al., 2015). The convective envelope then recedes as the helium core grows, leaving behind a chemical discontinuity in the gravity-mode cavity probed by the Brunt-Väisälä frequency. This peak in the Brunt-Väisälä frequency is localised much closer to the convective boundary — i.e. at a much shallower buoyancy depth — and results in a long-period modulation of the rotational ridges, inducing a gross curvature that is clearly visible for the $m=0$ modes in Fig. 9(a). The template-matching method used by Gehan et al. (2018) assumes that the ridges for each $m$ are straight, which contributes to the misidentification between the hares and the hounds. While curvature of this kind has been detected in observational mode frequencies (Mosser et al., 2018), it may not necessarily be uniquely attributed to buoyancy glitches (Deheuvels et al., 2023), at least for first-ascent red giants (as opposed to stars on the red clump; Vrard et al., 2022). Supposing that these glitches are indeed astrophysically rare (which would suggest the operation of additional mixing processes in actual red giants that are not present in our models; cf. Lindsay et al., 2022), we would expect the fraction of modes correctly identified to be higher for observed red giants than for the models used in this work.

However, we have also found implementation details specific to the template-matching algorithm used in Gehan et al. (2018, 2021) to be responsible for this azimuthal-order misidentification in many other cases. While an exhaustive denumeration of failure modes is neither possible nor productive, we note some of the major ones here. In some cases where (by numerical coincidence) the rotational frequency splitting is close to some integer multiple of $\nu^{2}\Delta\Pi_{l}/2$ (i.e. half the g-mode spacing), the hound may commit an off-by-one error in identifying entire ridges on the stretched échelle diagram, adjusting its adopted value of $\Delta\Pi_{l}$ to compensate. Fortuitously, this appears not to substantially affect the value of $\Omega_{\text{core}}$ that the hound needs to adopt, in order to generate a template that matches the hare, even under this misidentification. In other cases, we found that numerical issues incurred in the hound’s integration of Eq. 13 to yield $\tau$ (which we discuss in more detail in Section 4.1) caused the echelle diagrams obtained by the hound to be differently curved from the ground truth. For example, the curvature seen in Fig. 10, generated by integrating Eq. 13, is markedly different than that exhibited in Fig. 9a when the closed-form expression Eq. 15 is used — in particular, it is curved in the opposite direction. These issues resulted in template-matched misidentifications made preferentially at high and low frequencies, for essentially the same reasons as would be caused by a buoyancy glitch in the absence of these implementation issues.

3.2.3 Measurement of the mean core rotation rate

Finally, we consider the robustness of mean core rotation rates inferred with the method of Gehan et al. (2018). The relative errors in the mean core rotation measurement with respect to the ground-truth value are displayed in Fig. 11 as a function of evolution along the RGB.

In Fig. 11 (a), we show the relative differences between the g-mode rotational splitting for the hares, and the values inferred by the hounds, for Set I (examining dependences on mass). These are expressed as relative differences between the hound’s inferred values, and the JWKB estimator of the core rotation rate from the stellar models, which we evaluate as

[TABLE]

within the radiative core. The rotation rates of both the hound and the JWKB estimator are assumed to match the observed splittings through Eq. 4 with $\zeta=1$, and with the g-mode sensitivity factor, $\beta_{g}$, assumed to take its asymptotic value of ${1\over 2}$. However, this is known to differ from the true sensitivity factor $\beta_{g}$ by a small amount. Since the true value of $\beta_{g}$ is faithfully incorporated into the rotating mode frequencies evaluated by GYRE, but is only accounted for by the hound in such an approximate fashion, we would expect the hound’s values to likewise differ from that of the JWKB estimator by such a small amount even in the event that the hound’s measurements were perfectly accurate. For the $1.4M_{\odot}$ evolutionary track of of Set I, we thus also compare this JWKB estimator with $\beta={1\over 2}$ against a ground-truth g-mode rotational splitting $(\omega_{n,m=1}-\omega_{n,m=-1})/2$, averaged over the three $\gamma$-modes of the stellar model nearest to ${\nu_{\text{max}}}$. We show this comparison with the blue dotted line, up until the point where the g-mode rotational splitting becomes larger than ${\nu_{\text{max}}}^{2}\Delta\Pi_{1}$, the spacing between adjacent g-modes of the same $m$ at ${\nu_{\text{max}}}$.

We find no significant difference between the four masses considered. The relative error in the mean core rotation measurement decreases with evolution as the number of visible mixed modes increases. This is expected since g-dominated mixed modes become more sensitive to the core along evolution, resulting in a more accurate mean core rotation measurement. The relative error in the mean core rotation measurement then tends to suddenly increase above $\mathcal{N}\sim 27$, as a result of the numerous crossings in the stretched échelle diagram that seriously affect the correct identification of the mixed-mode azimuthal order.

In Fig. 11 (b), we show the same relative errors for Set II of hares (examining dependences on the ZAMS rotation rate $\Omega_{0}$). Once again, we obtain results qualitatively similar for the different core rotation rates we have considered. The relative error in the mean core rotation measurement first decreases with evolution, reaches minimum values for $\mathcal{N}\sim[7,13]$, before taking larger values again. At large $\mathcal{N}$ values, we see that the core rotation rate slightly impacts the measured core rotation rate. Indeed, the relative error in the mean core rotation measurement significantly increases for $\mathcal{N}\gtrsim 20$ in the case of $2\,\Omega_{\mathrm{0}}$, while it is the case for $\mathcal{N}\gtrsim 25$ in the case of $\Omega_{\mathrm{0}}$ and we do not see any increase in the case of $\Omega_{\mathrm{0}}/2$. As stated above, this is due to the increase in the number of crossings in the stretched échelle diagram along evolution. For a given mixed-mode density, a larger mean core rotation results in more crossings, leading to less accurate mean core rotation measurements. Hence, less accurate mean core rotation measurements start at lower $\mathcal{N}$ values as the mean core rotation rate increases.

In Fig. 11 (c), we show the same relative differences for Set III of hares. We obtain similar results for the different maximal mixing fractions considered, except for $\zeta_{\mathrm{max}}=1-10^{-1}$ for which we obtain significantly larger relative errors in some cases.

In summary, we find that:

• (1)

the relative systematic errors on the mean core rotation rate tend to be larger than the intrinsic error from the JWKB approximation, in particular for $\mathcal{N}_{\mathrm{1}}\lesssim 10$ where g-dominated mixed modes are less sensitive to the core;

• (2)

the mean core rotation rate as measured from the method of Gehan et al. (2018) tends to be systematically overestimated by at least $\sim 1-2$ %;

• (3)

the minimum relative error in the mean core rotation measurement is on the order of $\sim 4$ % or below for $\mathcal{N}\sim 7-27$;

• (4)

measuring accurate mean core rotation rates with a relative error of $\sim 1-2$ % requires to have some dipole mixed modes with $\zeta>0.90$.

We note that we infer overall quite accurate mean core rotation rates (Fig. 11), even despite mismatches in the mode identification (Fig. 9). The template matching method used by Gehan et al. (2018) counts the peaks that are close enough to the synthetic rotational template, and selects the configuration with the maximum number of peaks aligned with the synthetic ridges. Values of $m$ are assigned to modes post-hoc, based on the nearest synthetic ridge to each mode. As such, measuring the mean core rotation rate with the method developed by Gehan et al. (2018) does not require the azimuthal order of the modes to be identified at all on a mode-by-mode basis. Hence, we do not need to get a very high fraction of modes correctly identified to infer the mean core rotation rate with the method of Gehan et al. (2018), as long as the synthetic stretched échelle diagram matches well the morphology of the observed one.

4 Discussion

We now briefly outline some of the implications of our findings, for various aspects of how rotational characterisation with mixed modes may be performed.

4.1 Numerical Integration of $\tau$

We have derived a closed form expression, Eq. 15, to the solutions of the initial value problem specified by Eq. 13. Since the latter follows from differentiating the former only when the frequency dependence of $q$ is neglected, it follows that Eq. 13 cannot be used to generate diagrams with frequency-dependent $q$. Moreover, we note that the stretched echelle diagrams shown in Fig. 9a and Fig. 10 are constructed with the same modes, same value of $\Delta\Pi$, and same value of $q$. The only difference is that in Fig. 9a, the exact expression Eq. 15 is used, and in Fig. 10, $\tau$ is evaluated by numerical integration. Thus, the differences between the two, which are partially responsible for misidentification of $m$, must necessarily be attributed to systematics caused by such numerical integration.

Our results yield a simple criterion for the convergence of such numerical integration. Specifically, the accuracy of numerical integration of an ODE of the form $y^{\prime}=f(x,y)$ is determined by the relation between the integration step size to the characteristic scales in the $x$-coordinate over which $f$ varies. Our discussion of the oversampling criterion Eq. 17 required to accurately constrain the properties of $\zeta$ from observational data can therefore be adapted to also describe the step sizes required for these numerical integrations of Eq. 13 to accurately reproduce the analytic solutions for the stretching function $\tau$.

For illustrative purposes, we compare stretched échelle diagrams generated using stretching functions $\tau$ computed by various numerical solutions to Eq. 13, integrated using a Runge-Kutta scheme implemented in scipy. For this demonstration, $\zeta_{p}$ is specified by the frequency-dependent $q$ (via Eq. 36), the $\pi$-mode frequencies, and idealised g-modes with the same mean $\Delta\Pi$ and $\epsilon_{g}$ as of a first-ascent red giant stellar model with ${\nu_{\text{max}}}=100\ \mu$Hz, roughly located near the luminosity bump on the HR diagram. The integration of Eq. 13 to obtain $\tau$ is performed with respect to the dimensionless period coordinate $p=P/\Delta\Pi$, with initial condition $\tau=0$ as $p\to p_{0}$, where $p_{0}$ corresponds to the highest-frequency mode of the set of mixed modes for which we wish to generate the stretched échelle diagram. While this solver (RK45) is usually adaptive, we are also able to override the maximum step size permitted for use by the solver. Thus, by analogy with Eq. 17, we may define an oversampling discriminant associated with the maximum stepsizes of the integrator, $h_{0}$, as

[TABLE]

Accordingly, in addition to results from permitting the solver to adapt its stepsizes freely, we also show results obtained from setting the maximal stepsize to achieve various values of $D_{\text{samp}}$. Since the integrator does not accept a variable maximal stepsize, we evaluate $D_{\text{samp}}$ at $p=1/{\nu_{\text{max}}}\Delta\Pi_{1}$. We consider solutions with stepsizes chosen to yield various values of $D_{\text{samp}}$, indicated with different colours in Fig. 12.

We apply these numerical stretching functions to the corresponding idealised mixed modes associated with the specified $q$ and pure modes (generated for consistency with Eq. 7), showing the resulting stretched échelle diagrams in Fig. 12. Since the supplied mixed modes are fully consistent with the $\zeta_{p}$ used to generate the stretching functions, the exact solution should yield an essentially vertical line. The initial condition for integration is specified at high frequencies, and so we see that all of the numerical solutions agree there, and diverge as the integrator passes to lower frequencies. It can be seen that the fully adaptive strategy to selecting stepsizes used by the RK45 integrator yields large excursions from the correct solution in the neighbourhood of the pure p-modes. Even when smaller maximal stepsizes are enforced, it can be seen that integration error is incurred in the vicinity of the p-dominated mixed modes, with some residual variability as the integrator attempts to recover. These are reduced and eventually eliminated as $h_{0}$ is reduced so as to achieve $D_{\text{samp}}=1$. Conversely, in order to achieve $D_{\text{samp}}\geq 1$, we obtain an upper limit on the permissible step sizes to be used in numerical integration of Eq. 13:

[TABLE]

Precise details regarding the numerical integration of Eq. 13 have not generally been made available in most existing works employing the stretched-échelle-diagram construction. The integrand $f(P)\sim{1/\zeta_{p}(P)}$ can be well-approximated as being mostly close to 1 except for narrow, sharp peaks near the p-mode eigenvalues (per Eq. 16). As such, this IVP constitutes a quasiperiodic numerically stiff problem, not unlike that of integrating highly eccentric Keplerian orbits, and the details of this numerical integration are actually quite important. Beyond rotation, these issues may also interfere with the interpretation of other features of these stretched diagrams, such as potentially disrupting observational signatures of buoyancy glitches (as we saw in our discussion of Fig. 10). However, should the step-sizes of the integrators in these existing works be made known, this expression provides a criterion for whether or not their results ought to be reassessed. In such cases, numerical errors relating to the integration of Eq. 13 may be diagnosed by comparison to the analytic expression, as we have done in our hare-and-hounds exercise.

Since this bound is proportional to $q$, we may also surmise at least that the evaluation of the stretching function by numerical integration becomes increasingly computationally expensive in the limit of weak coupling, contraindicating doing so in the regime where we might wish to use it in conjunction with the matrix construction. Ultimately, these numerical issues may be avoided altogether by simply using the closed-form expression for $\tau$, Eq. 15, rather than performing any integration of Eq. 13 in the first place.

4.2 Modifications to existing asymptotic techniques

Template-matching techniques, as used by our hound, ordinarily require a set of mode frequencies to be measured from a power spectrum before they may be applied. However, the actual power spectrum ordinarily contains more information than just the mode frequencies: in particular, the mode linewidths and amplitudes are also known to strongly depend on the mixing fractions $\zeta$ (Mosser et al., 2018). Moreover, since the most g-dominated modes have the lowest amplitudes, they may not be significant enough to be reported by a standard peakbagging routine, despite being the least sensitive to misestimation of $\zeta$, and thus the most informative for the determination of $\Delta\Pi_{l}$ and $\Omega_{\text{core}}$. It would be desirable to include such constraints from these less-significant modes, however scarce, as exist in actual power spectra. Doing so would also ameliorate any ill effects which might otherwise arise when the oversampling criterion Eq. 17 is not satisfied.

Such an inclusion might be effected by template matching against stretched period-échelle power diagrams, as constructed in Gaulme et al. (2022), rather than lists of mode frequencies. This diagram was constructed in that work by analogy with the frequency-échelle power diagrams appearing in p-mode asteroseismology, to provide visual diagnostics for the g-dominated mixed modes directly from the power spectrum without necessitating an initial peakbagging step. Rather than evaluating $\tau$ by ODE-solving for this purpose, we may use the closed-form expression, Eq. 15, in this application as well. For demonstration, we show an illustrative example of such a diagram, constructed in this fashion, in Fig. 13. Aside from avoiding the above-described numerical issues and expense of ODE-solving at low $q$, doing so permits the g-mode phase offset $\epsilon_{g}$ to be read directly off the diagram, completely analogously to $\epsilon_{p}$ in frequency-echelle power diagrams; this information would be otherwise lost as a constant of integration for the ODE. Whereas Gaulme et al. (2022) mask out the regions containing the $\ell=0,2$ p-modes (indicated with the dashed lines), we note that it is precisely these regions that contain the most g-dominated dipole mixed modes. In the residual power in these regions, lower-amplitude g-dominated dipole modes can be seen, which would ordinarily be filtered out by the significance testing used in many peakfinding procedures (e.g. Corsaro & De Ridder, 2014), and which would therefore be inaccessible to template matching against frequency lists. The development of techniques for template-matching directly with respect to such stretched power diagrams, as opposed to with respect to frequency lists, is likely to be a substantial undertaking, which we defer to a future work.

4.3 Implications for rotational inversions

In a linear inversion procedure, it is the departures of the kernel functions from the scaled mean asymptotic kernel (specified by the JWKB wavenumber, as in Eq. 53) which permit localisation of perturbations from a fiducial stellar model. Even in relatively less evolved red giants ($\mathcal{N}_{1}\sim 5$), the kernels associated with the pure g-modes are deep into the asymptotic regime (e.g. fig. 13 of Paper I), and therefore carry very little information about rotation beyond that provided by the JWKB estimator (as further demonstrated by the ability of Ahlborn et al., 2022, to construct combination kernels out of mixed modes that suppress core sensitivity entirely). As such, any localisation provided by the mixed-mode rotational kernels will originate from the most p-dominated modes and their mixing fractions. Unfortunately, it is these mixing fractions which are least well-constrained by the asymptotic procedure. Similar considerations apply to other kinds of inversion problems, such as inversions for stellar structure and composition.

Consequently, seismic constraints on stellar rotation beyond the two-zone model, and specifically within the p-mode cavity, must draw upon additional priors on stellar evolution and structure, such as may be derived from explicit modelling. Fortunately, such modelling has to be performed anyway in order to permit the construction of appropriate inversion kernels in the first place. Given access to a sufficiently good guess to the structure of the star, the coupling strength $q$ may be then supplied a priori (from the coupling matrices via Eq. 36) rather than inferred generatively, which should in principle improve the constraints on $\zeta$, and potentially also on the notional pure p-mode rotational splittings, returned from the asymptotic procedure for the most p-dominated mixed modes. Supplying the coupling strengths in this fashion may also resolve other issues with inferring the pure p-mode frequencies themselves using the asymptotic construction (as encountered in Saunders et al., under review).

At the same time, the closed-form nature of the asymptotic parameterisation makes it far easier to operationalise, and computationally far less expensive, than the matrix decompositions employed in the algebraic approach. Moreover, unlike the matrix procedure, numerical solutions to its eigenvalue problem are not susceptible to truncation effects, and do not require an exhaustive search of eigenvalues. Conversely, the use of isolated p- and g-mode kernels (as we propose in Paper I) precisely suppresses the sensitivity of the inversion procedure to the inferred mixing fractions of the observed modes, thereby avoiding some of the shortcomings of the asymptotic procedure that we have illustrated with our hare-and-hounds exercise.

5 Conclusion

Our ultimate objective in this series of papers is to place constraints on the internal rotation of stars, beyond the presently-used one- or two-zone model of radial differential rotation. In Paper I, we found that the asymptotic formulation underlying the stretched-échelle-diagram construction had to be modified to treat each $m$ separately in order to permit this; in this work we assess various other features of techniques employing this construction. Since certain liberties taken in the $\zeta$-function construction are thought to limit its applicability to the most p-dominated mixed modes, we first produce a derivation that fills some of these conceptual deficiencies, and so show that in fact the current construction of the stretching function ought to be applicable to them as well — i.e. we first produce a "steel man" of the asymptotic construction that is as rigorous as we can make it. Given this exposition, we find that even within this fortified construction, a configuration of mixed modes must still satisfy a certain criterion, Eqs. 17 and 37, to provide enough information to constrain the true function $\tau$ describing it, and that most of this information comes from the p-dominated mixed modes. As a bonus, our geometric formulation yields a closed-form expression for the stretching function $\tau$ underlying stretched échelle diagrams, permitting it to be evaluated in closed form via Eq. 15, rather than by solving the initial value problem specified by Eq. 13, as is usually done.

From stellar modelling we find that red giant models in the range of $\mathcal{N}$ spanned by the observational sample appear not to satisfy this information criterion, indicating that characterisations of the mode coupling strengths of the avoided crossings, and thus mixing fractions for the most p-dominated mixed modes, are potentially misestimated — and systematically so — for this observational sample. We perform a hare-and-hounds exercise to determine how this might affect measurements of not just g-mode core rotation (i.e. the object of astrophysical interest for existing studies), but also of $\zeta$ and of mode identifications (of importance to the next steps in this inversion-oriented plan of work), as recovered by an existing observational pipeline using this asymptotic parameterisation — in particular, the template-matching algorithm developed by and used in Gehan et al. (2018, 2021), which relies on the stretched-échelle construction. The hare-and-hounds exercise shows that robust recovery of the former, but not of the latter (particularly for the most p-dominated modes), is possible when the information criterion is not satisfied. In particular, the identification of multiplet components decreases in accuracy as the coupling strength of the mixed-mode system weakens, and the mixing fraction $\zeta$ is identified least accurately for the most p-dominated mixed modes. However, some of these shortcomings stem from other methodological features of the algorithm that has played the role of the hound in this exercise, rather than being limitations of the asymptotic parameterisation per se, or issues stemming from the mischaracterisation of $q$.

In parallel to this, we compare the algebraic construction (which allows the generation of inversion kernels with desirable localisation properties; see Paper I) to the asymptotic one, so that in cases where the asymptotic construction is usable, we can immediately adapt its results for use with inversions. In Paper I we already show that envelope rotation is not accessible without significantly modifying current applications of the asymptotic construction; as above, we show here that other ingredients required for inversions (namely $\zeta$ for the most p-dominated modes, and the multiplet component identifications in general) are also methodologically inaccessible as well. To complete the comparison, we have derived an explicit translation between quantities appearing in the asymptotic construction, and in the nonasymptotic, matrix-based approach to mixed-mode coupling. Specifically, we find that the asymptotic eigenvalue equation can be obtained as the weak-coupling limit of the nonasymptotic eigenvalue equation, and that strict equivalence between the two formulations requires further constraints to be placed on the frequency dependence of the strength of the coupling between the two mode cavities. We interpret this as indicating that the nonasymptotic description is most well-approximated by the asymptotic one in the limit of weak coupling.

Since we have demonstrated that the two formulations become increasingly interoperable in the limit of weak coupling — particularly after the turnoff in $q$ with increasing $\mathcal{N}$ as seen in Fig. 6a — the use of features from each may remedy limitations of the other, particularly for the most evolved red giants. Such synergy may prove useful in also addressing an outstanding methodological gap in operationalising the algebraic approach for inversions, being that the notional pure p- and g-mode rotational splittings must be inferred from the observed mixed modes. We defer an investigation of this problem, and of the localisation properties of these isolated kernels, to future work.

We thank S. Basu, B. Mosser, J. van Saders, and M. Cunha for constructive feedback on this work. We also thank the anonymous referee for their unusually voluminous criticism and feedback, which have substantially improved the quality of the finished paper. JO acknowledges support from NASA through the NASA Hubble Fellowship grant HST-HF2-51517.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. CG was supported by Max Planck Society (Max Planck Gesellschaft) grant “Preparations for PLATO Science” M.FE.A.Aero 0011. We make available our implementation of the analytic expressions derived in this paper, and a helper routine to construct stretched échelle power diagrams, through a Python package hosted at https://gitlab.com/darthoctopus/zeta.

Appendix A Verification of an identity involving $\zeta$ in the algebraic construction

We proceed by direct differentiation, almost exactly following the derivation of eq. B12 in Ong & Basu (2020). We will first consider both mixed-mode eigenvalues at once, and specialise to the g-dominated mixed mode later. On the right-hand-side of Eq. 30, we have

[TABLE]

For the left-hand-side of Eq. 30, we first note that Eq. 28 implies

[TABLE]

Thus, we obtain

[TABLE]

For the more g-dominated of the two solutions, $\omega_{g}\sim\omega_{\pm}$, thereby yielding the desired result. However, for a pair of modes close to resonance, this is a good approximation for both modes in the limit of weak coupling.

Appendix B Maximum measured mixing fractions

In Table 1 we indicate the maximum mixing fractions $\zeta_{\mathrm{max}}$ that we measure for each given frequency range in Fig. 7. We report the maxima of modes binned in intervals 10 $\mu$Hz wide; the frequency values shown here are the midpoints of these intervals.