Issue 
A&A
Volume 687, July 2024



Article Number  A265  
Number of page(s)  26  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202348369  
Published online  30 July 2024 
Nonlinear threemode coupling of gravity modes in rotating slowly pulsating B stars
Stationary solutions and modeling potential
^{1}
Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
email: jordan.vanbeeck@protonmail.com
^{2}
TAPIR, Mailcode 35017, California Institute of Technology, Pasadena, CA 91125, USA
^{3}
Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA
^{4}
Royal Observatory of Belgium, Ringlaan 3, Brussels, Belgium
^{5}
Dept. of Astrophysics, IMAPP, Radboud University Nijmegen, 6500 GL Nijmegen, The Netherlands
^{6}
Max Planck Institute for Astronomy, Koenigstuhl 17, 69117 Heidelberg, Germany
^{7}
Guest Researcher, Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY 10010, USA
Received:
23
October
2023
Accepted:
21
December
2023
Context. Slowly pulsating B (SPB) stars display multiperiodic variability in the gravitoinertial mode regime with indications of nonlinear resonances between modes. Several have undergone asteroseismic modeling in the past few years to infer their internal properties, but only in a linear setting. These stars rotate fast, so that rotation is typically included in the modeling by means of the traditional approximation of rotation (TAR).
Aims. We aim to extend the set of tools available for asteroseismology, by describing timeindependent (stationary) resonant nonlinear coupling among three gravitoinertial modes within the TAR. Such coupling offers the opportunity to use mode amplitude ratios in the asteroseismic modeling process, instead of only relying on frequencies of linear eigenmodes, as has been done so far.
Methods. Following observational detections, we derive expressions for the resonant stationary nonlinear coupling between three gravitoinertial modes in rotating stars. We assess selection rules and stability domains for stationary solutions. We also predict nonlinear frequencies and amplitude ratio observables that can be compared with their observed counterparts.
Results. The nonlinear frequency shifts of stationary couplings are negligible compared to typical frequency errors derived from observations. The theoretically predicted amplitude ratios of combination frequencies match with some of their observational counterparts in the SPB targets. Other, unexplained observed ratios could be linked to other saturation mechanisms, to interactions between different modes, or to different opacity gradients in the driving zone.
Conclusions. For the purpose of asteroseismic modeling, our nonlinear mode coupling formalism can explain some of the stationary amplitude ratios of observed resonant mode couplings in single SPB stars monitored during 4 years by the Kepler space telescope.
Key words: asteroseismology / stars: evolution / stars: interiors / stars: oscillations / stars: rotation / stars: variables: general
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Slowly pulsating B (SPB) stars are midtolateB variable stars on the main sequence that display a variety of lowfrequency oscillations, and have masses ranging from ∼3 to ∼9 M_{⊙} (e.g., Waelkens 1991; De Cat & Aerts 2002; Pedersen et al. 2021; Szewczuk et al. 2021). Much of their multiperiodic variability is attributed to gravity modes excited by the κ mechanism associated with the Fe opacity bump in their envelope (Gautschy & Saio 1993; Dziembowski et al. 1993; Pamyatnykh 1999). Their rotation rates vary from ∼1% of critical to nearly critical velocity. The Coriolis force is a significant restoring force for most SPB oscillations, in addition to buoyancy (e.g., Lee 2012; Pedersen et al. 2021). We refer to such gravitoinertial modes whenever we mention g modes in this work. The large number of oscillations identified in space photometry of SPB stars have made them the subject of many asteroseismic modeling studies in the past few years (e.g., Degroote et al. 2010; Szewczuk & DaszyńskaDaszkiewicz 2018; Walczak et al. 2019; Wu & Li 2019; Wu et al. 2020; Szewczuk et al. 2021; Pedersen et al. 2021; Szewczuk et al. 2022).
The broad general review of asteroseismology by Aerts (2021, and references therein) makes it clear that stellar modeling is currently mainly done in a linear framework. Signals with frequencies approximately equal to linear combinations of frequencies of other detected signals, termed combination frequencies, are often detected in frequency lists generated by the harmonic analyses of stellar variability commonly used in asteroseismology (Aerts et al. 2010). Some of these combination frequencies can be explained by a nonlinear response of the stellar medium to the pulsation wave (see e.g., Bowman et al. 2016), which is referred to as nonlinear distortion by Degroote et al. (2009). In this work, however, we focus on combination frequencies that are explained by nonlinear coupling among oscillation modes (e.g., Buchler & Goupil 1984 and Van Hoolst 1994b). The amplitudes of heatdriven nonradial oscillations in SPB stars cannot be explained in the linear approximation. Observed amplitudes are therefore currently not used in asteroseismic inference. Nonlinear mode interactions that exchange energy among (coupled) modes must be taken into account to describe g mode amplitude limitation. Such interactions also change the mode frequencies. Instead of resorting to resourcecostly numerical integration of the nonlinear hydrodynamical equations that govern the oscillation dynamics, we consider weakly nonlinear effects of the lowest order. We therefore consider isolated weak nonlinear mode coupling among three g modes, and followed and extended the approach of Lee (2012), hereafter referred to as L12. Our approach is guided by the detected properties of SPB stars in Kepler observations, which are summarized in the sample studies by Pedersen et al. (2021) and Szewczuk et al. (2021).
Weak nonlinear mode coupling among three modes has been a topic of interest for stellar pulsation modes since the 1970s, with several seminal papers written long before space photometry was available (limiting ourselves to mode coupling among nonradial pulsations, see e.g., Dziembowski 1982, 1993; Buchler & Regev 1983; Aikawa 1984; Moskalik 1985; Dziembowski & Krolikowska 1985; Dappen & Perdang 1985; Dziembowski et al. 1988; Van Hoolst & Smeyers 1993; Buchler 1993; Takeuti & Buchler 1993; Van Hoolst 1994a,b, 1995; Van Hoolst 1996; Goupil & Buchler 1994; Buchler et al. 1995, 1997; Goupil et al. 1998; Wu & Goldreich 2001). Most of these formalisms focused on the description of mode coupling in nonrotating stars. A notable exception is the formalism developed by Friedman and Schutz (Friedman & Schutz 1978a,b; Schutz 1979), on which Schenk et al. (2001), hereafter referred to as S01, based their treatment of nonlinear threemode coupling in rotating stars. S01 included the effects of the Coriolis force perturbatively, but their framework is generic, allowing for the derivation of formalisms that do not treat the Coriolis force as a perturbation. Several studies are based on the S01 formalism, modeling nonlinear tides in multiplestar systems (e.g., Fuller & Lai 2012; Burkart et al. 2012, 2013, 2014; Fuller et al. 2013; O’Leary & Burkart 2014; Weinberg 2016; Fuller 2017; Guo 2020; Yu et al. 2020, 2021a; Zanazzi & Wu 2021) and in starexoplanet systems (Essick & Weinberg 2016; Vick et al. 2019; Yu et al. 2021b, 2022), nonlinear interactions among modes in neutron stars (e.g., Morsink 2002; Arras et al. 2003; Lai & Wu 2006; Weinberg & Quataert 2008; Weinberg et al. 2013), tidal migration in the moon systems of Jupiter and Saturn (Fuller et al. 2016), nonlinear interactions among mixed modes in red giant stars (Weinberg & Arras 2019; Weinberg et al. 2021), or resonant mode coupling in δ Sct stars (Mourabit & Weinberg 2023).
L12 used the S01 formalism as a basis for an extension that described rapidly rotating stars in which the Coriolis force cannot be treated as a perturbation. To do so, they adopted the socalled traditional approximation for rotation (TAR), in which the latitudinal component of the rotation vector in a spherical coordinate system was neglected, assuming spherical symmetry. This decoupled the radial and horizontal components of the pulsation equations (e.g., LonguetHiggins 1968; Lee & Saio 1997; Townsend 2003; Mathis 2013). For highradial order g modes it is justified to ignore the centrifugal force within the TAR, because the dominant contribution to the mode energy occurs deep inside the star, close to the convective core, where rotational deformation is small and spherical symmetry is a good approximation (Mathis & Prat 2019; Henneco et al. 2021; Dhouib et al. 2021a,b). There is a marked difference in scale of the rotation frequency Ω and the BruntVäisälä frequency N in those deep nearcore regions. The assumption that the Coriolis force is weaker than the buoyancy force in the direction of stable entropy or chemical stratification is therefore fulfilled near the core for lowfrequency (Poincaré) modes^{1}. Their horizontal velocities are also greater than the vertical velocities (see Mathis & Prat 2019), justifying the neglect of the latitudinal rotation vector component within the TAR.
L12 used their quadratic nonlinear mode coupling formalism within the TAR (based on the S01 formalism) to provide a numerical mode coupling example for a specific SPB star model near the zeroage main sequence. The computed mode properties of that example were not compared to observed SPB mode properties. In this work we extend and correct the L12 formalism with the aim of creating a modeling framework that can be used to model nonlinear threemode coupling of g modes in some of the 38 SPB stars considered by Van Beeck et al. (2021), hereafter referred to as V21. We specifically focus on deriving the conditions for which the amplitudes of modes in coupled mode triads, and their combination phase, do not vary over time. The mode parameters inferred by V21 allowed for the discovery of many such potentially ‘locked’ mode couplings. We therefore contrast the theoretically predicted observables computed by our formalism for mode couplings obtained from models typical for the ensemble of SPB stars analyzed by V21 with their detected observational counterparts. We first provide a rigorous overview of our theoretical oscillation model in Sect. 2, followed by Sect. 3, which describes the nonlinear theoretical observables. In Sect. 4 we show the numerical results for resonant mode couplings typical for SPB stars, while Sect. 5 discusses the potential of our theoretical framework for future asteroseismic modeling. Finally, Sect. 6 outlines our conclusions and prospects.
2. Theoretical oscillation model
2.1. Linear free oscillations within the TAR
The linearized momentum equation governing linear stellar oscillations in uniformly rotating stars is expressed in a corotating reference frame as (e.g., Frieman & Rotenberg 1960, LyndenBell & Ostriker 1967 or S01)
$$\begin{array}{c}\hfill \ddot{\mathit{\xi}}+\mathit{B}(\dot{\mathit{\xi}})+\mathit{C}(\mathit{\xi})={\mathit{a}}_{\mathrm{ext}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(1)
where ξ denotes the Lagrangian displacement, the superscripted dot indicates a partial time derivative, $\mathit{B}(\dot{\mathit{\xi}})\equiv 2\phantom{\rule{0.166667em}{0ex}}\mathbf{\Omega}\times \dot{\mathit{\xi}}$ is the Coriolis term, with Ω = Ω e_{z} the (uniform) rotation vector and e_{z} the unit vector along the rotation axis, C(ξ) is the term that describes forces not depending on the oscillation frequency, and a_{ext} is any acceleration due to external forces. For the free oscillations used in linear asteroseismology, a_{ext} = 0. The operators B and C are antiHermitian and Hermitian, respectively (LyndenBell & Ostriker 1967). An equivalent tensor representation of the linearized momentum equation is described in Appendix A and used in Sect. 2.3.
With time dependence Ansatz
$$\begin{array}{c}\hfill \mathit{\xi}(\mathit{x},t)=\mathit{\xi}(\mathit{x})\phantom{\rule{0.166667em}{0ex}}{e}^{i\omega t}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(2)
we can rewrite Eq. (1) for free oscillations as
$$\begin{array}{c}\hfill \phantom{\rule{0.166667em}{0ex}}{\omega}^{2}\mathit{\xi}i\phantom{\rule{0.166667em}{0ex}}\omega \phantom{\rule{0.166667em}{0ex}}\mathit{B}(\mathit{\xi})+\mathit{C}(\mathit{\xi})=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(3)
where ω is the (real or complexvalued) angular frequency in the corotating frame (e.g., Friedman & Schutz 1978a, S01 and Prat et al. 2019). The Hermiticity of the operators iB and C allow one to define a generic orthogonality relation valid for two distinct (ordinary) eigenmodes ξ_{φ} and ξ_{β} of Eq. (3) if ${\omega}_{\beta}\ne {\omega}_{\phi}^{\ast}$ and Ω is timeindependent, where we use the superscript ‘*’ to indicate a complex conjugated quantity,
$$\begin{array}{c}\hfill ({\omega}_{\beta}+{\omega}_{\phi}^{\ast})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}\rangle +\langle {\mathit{\xi}}_{\phi},i\mathit{B}({\mathit{\xi}}_{\beta})\rangle =0\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(4)
In Eq. (4), the inner product of the Hilbert space ℋ spanned by the complex eigenvectors ξ and ξ′ of Eq. (3) is defined as
$$\begin{array}{c}\hfill \langle \mathit{\xi},{\mathit{\xi}}^{\prime}\rangle ={\displaystyle \int \rho}\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}^{\ast}\xb7{\mathit{\xi}}^{\prime}\mathrm{d}V\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(5)
A proof of the orthogonality condition implied in Eq. (4) for the two distinct modes ξ_{φ} and ξ_{β} is given in Appendix B. For a mode φ with complexvalued ω_{φ} described by Eq. (3), the linear heatdriven growth (damping) rate γ_{φ} is defined as γ_{φ} ≡ Im[ω_{φ}]. Because of Ansatz (2), linear growth (damping) occurs when γ_{φ} is positive (negative).
We describe the coupling among nondegenerate modes using their adiabatic eigenfunctions (in Sect. 2.3), for which the generic orthogonality relation implied in Eq. (4) for the modes φ and β can be written as (see Appendix B and S01)
$$\begin{array}{c}\hfill ({\mathrm{\Omega}}_{\phi}+{\mathrm{\Omega}}_{\beta})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}\rangle +\langle {\mathit{\xi}}_{\phi},i\mathit{B}\left({\mathit{\xi}}_{\beta}\right)\rangle ={\delta}_{\beta}^{\phi}\phantom{\rule{0.166667em}{0ex}}{b}_{\phi}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(6)
where ${\delta}_{\beta}^{\phi}$ denotes a Kronecker delta, Re[ω_{φ}] = Ω_{φ} (and similar for mode β), and the realvalued constant b_{φ} is given by
$$\begin{array}{c}\hfill {b}_{\phi}=\langle {\mathit{\xi}}_{\phi},i\mathit{B}\left({\mathit{\xi}}_{\phi}\right)\rangle +2\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\phi}\rangle \phantom{\rule{0.166667em}{0ex}}.\end{array}$$(7)
The constant b_{φ} defined in Eq. (7) is related to the rotatingframe mode energy ϵ_{φ} at unit complex amplitude, if Ω_{φ} ≠ 0 (S01):
$$\begin{array}{c}\hfill {\u03f5}_{\phi}={\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{b}_{\phi}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(8)
Because we describe uniformly rotating nonmagnetic stars and use the Cowling approximation in which the Eulerian perturbation of the gravitational potential is set to zero,
$$\begin{array}{c}\hfill \mathit{C}\left(\mathit{\xi}\right)=\frac{\mathbf{\nabla}\delta P}{\rho}\frac{\mathbf{\nabla}P}{{\rho}^{2}}\delta \rho \phantom{\rule{0.166667em}{0ex}},\end{array}$$(9)
with δP and δρ being the Eulerian perturbations of the pressure and density around their equilibrium values P and ρ (e.g., P19). The explicit dependence of C(ξ) on ξ is defined in, for example, LyndenBell & Ostriker (1967).
We use the TAR, which ignores the latitudinal component − Ω sinθ e_{θ} of the rotation vector (e.g., Lee & Saio 1997). The Lagrangian displacement of a g mode φ can then be expressed in spherical coordinates (r, θ, ϕ) as (Prat et al. 2019)
$$\begin{array}{c}\hfill {\mathit{\xi}}_{\phi}=[{\xi}_{r}^{\phi}(r)\phantom{\rule{0.166667em}{0ex}}{H}_{r}(\theta ),\phantom{\rule{0.166667em}{0ex}}{\xi}_{h}^{\phi}(r)\phantom{\rule{0.166667em}{0ex}}{H}_{\theta}(\theta ),\phantom{\rule{0.166667em}{0ex}}i\phantom{\rule{0.166667em}{0ex}}{\xi}_{h}^{\phi}(r)\phantom{\rule{0.166667em}{0ex}}{H}_{\varphi}(\theta )]\phantom{\rule{0.166667em}{0ex}}{e}^{i\phantom{\rule{0.166667em}{0ex}}({m}_{\phi}\phantom{\rule{0.166667em}{0ex}}{\varphi}_{\phi}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}t)}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(10)
where m_{φ} is the mode azimuthal order, ϕ_{φ} is the mode phase, ${\xi}_{r}^{\phi}(r)$ and ${\xi}_{h}^{\phi}(r)$ are the radial and horizontal mode displacement components, and H_{r}(θ), H_{θ}(θ) and H_{ϕ}(θ) are the realvalued radial, latitudinal and azimuthal Hough functions of a mode φ (see Appendix A of Prat et al. 2019 for their definitions). Following L12, we normalize radial Hough function H_{r}(θ) as
$$\begin{array}{c}\hfill 2\phantom{\rule{0.166667em}{0ex}}\pi \phantom{\rule{0.166667em}{0ex}}{\displaystyle {\int}_{1}^{1}}{{H}_{r}\phantom{\rule{0.166667em}{0ex}}(\theta )\phantom{\rule{0.166667em}{0ex}}}^{2}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\mu =1\phantom{\rule{0.166667em}{0ex}},\end{array}$$(11)
where μ = cosθ. The normalization of H_{r}(θ) also determines the normalization of the latitudinal and azimuthal Hough functions H_{θ}(θ) and H_{ϕ}(θ).
2.2. Nonlinear oscillations within the TAR
Nonlinear mode interactions add acceleration (i.e., forcing) terms to the governing equation of motion, so that the higherorder equation of motion becomes the following quadratic eigenfunction problem (see e.g., Frieman & Rotenberg 1960)
$$\begin{array}{c}\hfill \phantom{\rule{0.166667em}{0ex}}{\omega}^{2}\mathit{\xi}i\phantom{\rule{0.166667em}{0ex}}\omega \phantom{\rule{0.166667em}{0ex}}\mathit{B}\left(\mathit{\xi}\right)+\mathit{C}\left(\mathit{\xi}\right)=\mathit{a}[\mathit{\xi}]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(12)
under the Ansatz given by Eq. (2): ξ(x, t)=ξ(x) e^{−iωt}. In Eq. (12), the first term describes the acceleration. That acceleration is changed by the nonlinear coupling term a[ξ], which can be split up into contributions from pressure gradients (a_{P}[ξ]) and gravity (a_{G}[ξ]).
Because we limit ourselves to nonlinear threemode interactions, it is sufficient to expand a[ξ] to second order (e.g., S01),
$$\begin{array}{c}\hfill \mathit{a}[\mathit{\xi}]={\mathit{a}}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\prime}]+O\left({\mathit{\xi}}^{3}\right)={\mathit{a}}_{P}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\prime}]+{\mathit{a}}_{G}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\prime}]+O\left({\mathit{\xi}}^{3}\right),\end{array}$$(13)
where a^{(2)}[ξ, ξ′], the second order pressure term ${\mathit{a}}_{P}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\mathrm{\prime}}]$, and the second order gravitational acceleration term ${\mathit{a}}_{G}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\mathrm{\prime}}]$ are symmetric bilinear functions of the Lagrangian displacements associated with the interacting modes. The nth components of ${\mathit{a}}_{P}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\mathrm{\prime}}]$ and ${\mathit{a}}_{G}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\mathrm{\prime}}]$ are expressed as (e.g., S01)
$$\begin{array}{cc}\hfill {\mathit{a}}_{G,\phantom{\rule{0.166667em}{0ex}}n}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\prime}]=& \phantom{\rule{0.166667em}{0ex}}\frac{1}{2}\phantom{\rule{0.166667em}{0ex}}{\xi}^{k}\phantom{\rule{0.166667em}{0ex}}{\xi}^{l}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\nabla}}_{k}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\nabla}}_{l}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\nabla}}_{n}\phantom{\rule{0.166667em}{0ex}}\mathrm{\Phi}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(14a)
$$\begin{array}{cc}\hfill {\mathit{a}}_{P,\phantom{\rule{0.166667em}{0ex}}n}^{(2)}[\mathit{\xi},{\mathit{\xi}}^{\prime}]=& \phantom{\rule{0.166667em}{0ex}}\frac{1}{\rho}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\nabla}}_{j}[p\phantom{\rule{0.166667em}{0ex}}({\mathrm{\Gamma}}_{1}1){\mathrm{\Theta}}_{n}^{j}+p\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Xi}}_{n}^{j}+\mathrm{\Psi}\phantom{\rule{0.166667em}{0ex}}{\delta}_{n}^{j}]\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(14b)
where we use the Cowling approximation and the Einstein summation convention, ∇_{φ} denotes the covariant derivative with respect to coordinate x^{φ}, the adiabatic exponent Γ_{1} = (∂lnp/∂lnρ)_{S} with subscript S denoting adiabatic conditions, and Φ is the gravitational potential. We also use the following definitions of the tensors
$$\begin{array}{cc}\hfill {\mathrm{\Theta}}_{i}^{j}=& \phantom{\rule{0.166667em}{0ex}}\left({\mathrm{\nabla}}_{i}\phantom{\rule{0.166667em}{0ex}}{\xi}^{j}\right)\left({\mathrm{\nabla}}_{k}\phantom{\rule{0.166667em}{0ex}}{\xi}^{k}\right)=\left({\mathrm{\nabla}}_{i}\phantom{\rule{0.166667em}{0ex}}{\xi}^{j}\right)(\mathrm{\nabla}\xb7\mathit{\xi})\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(15a)
$$\begin{array}{cc}\hfill {\mathrm{\Xi}}_{i}^{j}=& \phantom{\rule{0.166667em}{0ex}}\left({\mathrm{\nabla}}_{k}\phantom{\rule{0.166667em}{0ex}}{\xi}^{j}\right)\left({\mathrm{\nabla}}_{i}\phantom{\rule{0.166667em}{0ex}}{\xi}^{k}\right)\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(15b)
and of the quantities
$$\begin{array}{cc}& \mathrm{\Psi}=\phantom{\rule{0.166667em}{0ex}}\frac{p}{2}\{\mathrm{\Theta}[{({\mathrm{\Gamma}}_{1}1)}^{2}+{\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}]+\mathrm{\Xi}({\mathrm{\Gamma}}_{1}1)\}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(16a)
$$\begin{array}{cc}& \mathrm{\Theta}=\phantom{\rule{0.166667em}{0ex}}{\delta}_{j}^{i}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Theta}}_{i}^{j}={(\mathrm{\nabla}\xb7\mathit{\xi})}^{2}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(16b)
$$\begin{array}{cc}& \mathrm{\Xi}=\phantom{\rule{0.166667em}{0ex}}{\delta}_{j}^{i}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Xi}}_{i}^{j}=\left({\mathrm{\nabla}}_{k}\phantom{\rule{0.166667em}{0ex}}{\xi}^{i}\right)\left({\mathrm{\nabla}}_{i}\phantom{\rule{0.166667em}{0ex}}{\xi}^{k}\right)\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}$$(16c)
We provide additional information on the expressions for the quantities in Eqs. (15) and (16) in Appendix C when discussing the explicit terms of the mode coupling coefficient defined in Eq. (22) of Sect. 2.3. Appendix C also contains corrections and simplifications for some of the expressions in L12.
2.3. Coupledmode equations
Eigenvalue Eq. (12) needs to be solved to retrieve the eigenfunctions ξ(x) and eigenfrequencies ω for interacting modes in rotating stars. Each pulsation mode φ of the star is characterized by a pair (ξ_{φ}, ω_{φ}). The linear eigenvalue Eq. (3) is used to calculate the linear mode eigenfunctions, which are subsequently employed in the mode coupling computations.
We follow the procedure of S01 to derive a set of equations that describe the couplings among all eigenmodes, which we refer to as the coupledmode equations. This procedure solves the phase space equivalent of the eigenvalue Eq. (12), defined in Eq. (A.6), by employing the phase space mode decomposition for a physically relevant real Lagrangrian displacement ξ(x,t) of the set of coupled modes in a rotating star, which is given by (S01)
$$\begin{array}{c}\hfill \left[\begin{array}{c}\mathit{\xi}(t)\\ \dot{\mathit{\xi}}(t)\end{array}\right]={\displaystyle \sum _{\phi}}{c}_{\phi}(t)\left[\begin{array}{c}{\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\\ i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\end{array}\right]+{c}_{\phi}^{\ast}(t)\left[\begin{array}{c}{\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right)\\ i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right)\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(17)
under the assumption that mode φ is not degenerate and has a realvalued frequency Ω_{φ}. This is justified because the linear mode eigenfrequencies associated with the (linear) adiabatic eigenfunctions and used to compute the coefficient relevant for mode coupling (see Eq. (17)) are real. The realvaluedness of the eigenfrequencies further justifies the use of eigenmode orthogonality relation (6). In principle, the expansion (17) should also include modes with nonordinary eigenvectors (i.e., modes with Jordan chains of length greater than zero, as was explained in e.g., Schutz 1979, Schutz 1980a,b and S01). Similar to the assumption made for rmodes in S01, we assume that such modes with nonordinary eigenvectors are not important dynamically for the saturation of g modes in SPB stars and therefore do not include these modes in the sum in Eq. (17). A formalism that includes these nonordinary eigenvectors can for example be found in Appendix A of S01.
The complex mode amplitudes c_{φ}(t) in expansion (17) and their complex conjugates ${c}_{\phi}^{\ast}(t)$ are given by (see Appendix D)
$$\begin{array}{cc}\hfill {c}_{\phi}(t)=& \phantom{\rule{0.166667em}{0ex}}\frac{1}{{b}_{\phi}}\langle {\mathit{\xi}}_{\phi},{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}\mathit{\xi}(t)+i\dot{\mathit{\xi}}(t)+i\mathit{B}(\mathit{\xi}(t))\rangle \phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(18a)
$$\begin{array}{cc}\hfill {c}_{\phi}^{\ast}(t)=& \phantom{\rule{0.166667em}{0ex}}\frac{1}{{b}_{\phi}}\langle {\mathit{\xi}}_{\phi}^{\ast},{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}\mathit{\xi}(t)i\dot{\mathit{\xi}}(t)i\mathit{B}(\mathit{\xi}(t))\rangle \phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(18b)
for nonzero b_{φ} and Ω_{φ} ≠ −Ω_{β}, where the subscript β refers to a different nondegenerate mode in expansion (17). The rotatingframe mode energy E_{φ} for mode φ with rotatingframe frequency Ω_{φ} and (complex) amplitude c_{φ} is (see Appendix K in S01)
$$\begin{array}{c}\hfill {E}_{\phi}={\u03f5}_{\phi}\phantom{\rule{0.166667em}{0ex}}{c}_{\phi}{}^{2}={\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{b}_{\phi}\phantom{\rule{0.166667em}{0ex}}{{c}_{\phi}}^{2}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(19)
We disregard nonlinear coupling of g modes with toroidal r modes (i.e., r − g mode coupling) in this work because r modes have not been detected in the SPB stars analyzed in V21. Instead we focus on interg mode coupling. For fasterrotating SPB stars, such r − g mode couplings can become significant (see e.g., L12), but only if their coupling initiation threshold is similar to or smaller than that of the interg mode couplings. As discussed by Buchler et al. (1995), stable (constant), periodic or irregular timedependent behavior of the oscillation amplitudes (and frequencies) may be observed. In this work, we focus on the stable (timeindependent) amplitude and frequency solutions of the amplitude equations (see Sect. 2.6) because the observations of mode couplings indicate that this is the dominant observed behavior in many of the SPB stars (see Sect. 2.7). We therefore do not discuss periodic modulation or irregular behavior of amplitudes or frequencies^{2}.
By substituting the phase space mode expansion (17) into the phase space equivalent of the equation of motion, defined in Eq. (A.1), and accounting for mode orthogonality, we derive the coupledmode equations in the quadratic approximation (see Appendix A of S01 for the technical details),
$$\begin{array}{c}\hfill {\dot{c}}_{\phi}+i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{c}_{\phi}=\frac{i}{{b}_{\phi}}({\kappa}_{\phi}^{\beta \gamma}{c}_{\beta}{c}_{\gamma}+{\kappa}_{\phi}^{\overline{\beta}\gamma}{c}_{\beta}^{\ast}{c}_{\gamma}+{\kappa}_{\phi}^{\beta \overline{\gamma}}{c}_{\beta}{c}_{\gamma}^{\ast}+{\kappa}_{\phi}^{\overline{\beta}\overline{\gamma}}{c}_{\beta}^{\ast}{c}_{\gamma}^{\ast})\phantom{\rule{0.166667em}{0ex}},\end{array}$$(20)
in which we employ the Einstein summation convention. If we combine Eqs. (8) and (20) and assume nonzero Ω_{φ}, we get
$$\begin{array}{c}\hfill \frac{{\dot{c}}_{\phi}}{i}+{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{c}_{\phi}={\mathrm{\Omega}}_{\phi}({\eta}_{\phi}^{\beta \gamma}{c}_{\beta}{c}_{\gamma}+{\eta}_{\phi}^{\overline{\beta}\gamma}{c}_{\beta}^{\ast}{c}_{\gamma}+{\eta}_{\phi}^{\beta \overline{\gamma}}{c}_{\beta}{c}_{\gamma}^{\ast}+{\eta}_{\phi}^{\overline{\beta}\overline{\gamma}}{c}_{\beta}^{\ast}{c}_{\gamma}^{\ast})\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(21)
The implicit sums over indices β and γ in Eqs. (20) and (21) run over all modes with ordinary eigenvectors. Coupling coefficient ${\kappa}_{\phi}^{\beta \gamma}$ and energyscaled coupling coefficient ${\eta}_{\phi}^{\beta \gamma}$ involved in the coupling of the three modes φ, β and γ are defined by
$$\begin{array}{c}\hfill {\kappa}_{\phi}^{\beta \gamma}={\u03f5}_{\phi}\phantom{\rule{0.166667em}{0ex}}{\eta}_{\phi}^{\beta \gamma}=\langle {\mathit{\xi}}_{\phi},\phantom{\rule{0.166667em}{0ex}}{\mathit{a}}^{(2)}[{\mathit{\xi}}_{\beta},\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\gamma}]\rangle \phantom{\rule{0.166667em}{0ex}},\end{array}$$(22)
which deviates from the definition in S01: ξ_{φ} is used instead of its complex conjugate. The (realvalued) total displacement is thus expanded in terms of nonconjugated products of mode amplitudes and eigenvectors, instead of their complex conjugated analogues used in S01. The resulting expression (22) is equivalent, and directly yields the specific coupling coefficient selection rules derived in Sect. 2.4. A bar over a subscript or superscript of ${\kappa}_{\phi}^{\beta \gamma}$ or ${\eta}_{\phi}^{\beta \gamma}$ indicates that a complex conjugate of the corresponding mode eigenfunction is used in the definition in Eq. (22). The expressions for the coupling coefficients ${\kappa}_{\phi}^{\beta \gamma}$ and ${\eta}_{\phi}^{\beta \gamma}$ are symmetric under permutation of indices, due to the presence of the symmetric bilinear function a^{(2)}[ξ_{β}, ξ_{γ}]. We normalize the eigenfunctions ξ_{φ}, ξ_{β} and ξ_{γ} such that ϵ_{φ} = ϵ_{β} = ϵ_{γ} = GM^{2}/R, following L12. Additional information on the explicit expressions used in Eq. (22) can be found in Appendix C.
The coupledmode Eqs. (20) and (21) thus describe the temporal dynamics of the amplitudes of all (coupled) modes in a coupling network. To study specific threemode nonlinear interactions among g modes in SPB stars, we derive amplitude equations (AEs) in Sect. 2.6, based on the coupledmode equations. These AEs simplify the study of mode resonances.
2.4. Coupling coefficient selection rules
For the coupling coefficients defined in Eq. (22) to be nonzero, selection rules need to be satisfied (see e.g., S01). Considering the integration of coupling coefficient ${\kappa}_{\phi}^{\beta \gamma}$ over longitude ϕ,
$$\begin{array}{c}\hfill {\kappa}_{\phi}^{\beta \gamma}\propto {\displaystyle {\int}_{0}^{2\phantom{\rule{0.166667em}{0ex}}\pi}}exp\left(i\phantom{\rule{0.166667em}{0ex}}({H}_{\phi}\phantom{\rule{0.166667em}{0ex}}{m}_{\phi}+{H}_{\beta}\phantom{\rule{0.166667em}{0ex}}{m}_{\beta}+{H}_{\gamma}\phantom{\rule{0.166667em}{0ex}}{m}_{\gamma})\phantom{\rule{0.166667em}{0ex}}\varphi \right)d\varphi \phantom{\rule{0.166667em}{0ex}}.\end{array}$$(23)
The same longitudinal dependence is obtained for ${\eta}_{\phi}^{\beta \gamma}$, which yields the generic azimuthal selection rule
$$\begin{array}{c}\hfill {H}_{\phi}\phantom{\rule{0.166667em}{0ex}}{m}_{\phi}+{H}_{\beta}\phantom{\rule{0.166667em}{0ex}}{m}_{\beta}+{H}_{\gamma}\phantom{\rule{0.166667em}{0ex}}{m}_{\gamma}=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(24)
in which we introduce the sign factor H_{φ} that has the value −1 for complex conjugated modes in the coupling coefficient expression (i.e., modes that have a bar over their index) and 1 otherwise. Equation (24) for example leads to the azimuthal selection rule m_{φ} = m_{β} + m_{γ} for coupling coefficients ${\kappa}_{\phi}^{\beta \gamma}$ and ${\eta}_{\phi}^{\beta \gamma}$.
The other selection rule takes into account the symmetry of the modes across the stellar equator. Whether g mode φ is odd or even is determined by (−1)^{kφ}. If mod(k_{φ} , 2) = 1, mode φ is odd, if it is zero, the mode is even (e.g., Lee & Saio 1997). This selection rule requires that two or no odd modes partake in the mode coupling (see Appendix E), which can also be written as a selection rule based on the mode ordering number k,
$$\begin{array}{c}\hfill \mathrm{mod}\left(\phantom{\rule{0.166667em}{0ex}}{k}_{\phi}+{k}_{\beta}+{k}_{\gamma}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}},\phantom{\rule{0.166667em}{0ex}}2\phantom{\rule{0.166667em}{0ex}}\right)=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(25)
where, for a g mode φ, k_{φ} ≡ l_{φ} − m_{φ} (Lee & Saio 1997). We refer to Eq. (25) as the meridional selection rule because Van Reeth et al. (2022) called k the meridional degree.
2.5. Resonant coupling networks
The three potential scenarios for resonant coupling among three modes, hereafter referred to as the isolated resonant threemode coupling networks, are pictured in Fig. 1. A direct resonant coupling is an interaction between a linearly damped daughter mode φ and the linearly driven parent modes β and γ, and is depicted as scenario (i) in Fig. 1. In this scenario, γ_{φ} < 0, and γ_{β}, γ_{γ} > 0. A parametric resonant coupling is an interaction between a linearly driven parent mode φ and the linearly damped daughter modes β and γ, and is depicted as scenario (ii) in Fig. 1. Hence, in this case, γ_{φ} > 0, and γ_{β}, γ_{γ} < 0. We call the resonant interaction between three linearly driven modes φ, β and γ a driven resonant coupling, which is depicted as scenario (iii) in Fig. 1, for which γ_{φ}, γ_{β}, γ_{γ} > 0.
Fig. 1. Representations of isolated resonant threemode coupling networks for modes φ, β and γ. They represent direct resonances (i), parametric resonances (ii), and driven resonances (iii). A linearly unstable (i.e., excited) mode is pictured as an orange circle, whereas a linearly stable (i.e., damped) mode is pictured as a blue square. 
In addition to these three isolated resonant coupling networks, nonisolated resonant coupling networks exist, in which a daughter or parent mode is shared among different resonantly coupled mode triads. O’Leary & Burkart (2014) called this multiplemode or multimode coupling (see also e.g., Guo 2021). The selection rules derived in Sect. 2.4 remain valid in such networks, because coupling occurs at the same quadratic order. Networks with granddaughter and greatgranddaughter modes have also been constructed. The reader is referred to Kumar & Goodman (1996) and Weinberg et al. (2021) for additional information on such networks.
Higherorder coupling terms, such as those that appeared in the cubic fourwave coupling networks developed to describe four interacting modes in nonrotating stars (by for example Van Hoolst & Smeyers 1993 and Van Hoolst 1994a) do not adhere to the selection rules in Sect. 2.4. The fourwave selection rules can however be determined from similar arguments. For example, the fourwave azimuthal selection rule was derived in Weinberg (2016). They however used the complex conjugate of mode φ in their definition of the mode coupling coefficient and did not expand in pairs of complex conjugates (as in Eq. (17)). This yielded a slightly different formulation of the azimuthal selection rule than what would be obtained for a coupling coefficient defined in a similar way as in Eq. (22).
A priori, it is not clear whether multimode or higherorder coupling terms are necessary to explain the (stationary) amplitudes and phases of the threemode resonance candidate couplings detected among the modes of SPB stars analyzed in V21. The most parsimonious models for these identified candidate couplings would be isolated threemode sumfrequency coupling networks that produce stable stationary solutions, if those solutions effectively describe observed (timeindependent) amplitudes, phases and frequencies. We limit ourselves in the remainder of this work to describing such coupling networks.
2.6. Amplitude equations for Ω_{1} ≃ Ω_{2} + Ω_{3}
Under the assumption that the linear driving or linear damping rate γ_{φ} of a mode φ fulfills γ_{φ}/Ω_{φ}≪1 for the modes that are weakly nonlinearly coupled, we can derive the AEs. In this section we use the coupledmode Eqs. (20) and (21) valid for quadratically nonlinear oscillations in rotating stars to derive the specific AEs for the sumfrequency resonant interaction Ω_{1} ≃ Ω_{2} + Ω_{3} between three distinct modes. Similar AEs were derived for the resonant harmonic interaction Ω_{1h} ≃ Ω_{2h} (with subscript h denoting the harmonic nature of the resonance) in Appendix B.2 of Van Beeck (2023). The AEs derived in this Section also describe the temporal evolution of the coupled mode’s properties of differencefrequency resonances Ω_{1d} = Ω_{2d} − Ω_{3d}, as we show in Appendix F.
Following Van Hoolst (1995), we apply the multiple time scales perturbation method (e.g., Nayfeh 1973, 1981; Nayfeh & Mook 1979) to compute the mode amplitudes due to nonlinear coupling at the quadratic level. We introduce time variables
$$\begin{array}{c}\hfill {t}_{0}=t,{t}_{1}=\mathfrak{J}\phantom{\rule{0.166667em}{0ex}}t,\dots \phantom{\rule{0.166667em}{0ex}},\end{array}$$(26)
where 𝔍 is a small, dimensionless ordering parameter. The time derivative in Eq. (20) becomes
$$\begin{array}{c}\hfill \frac{\partial}{\partial t}=\frac{\partial}{\partial {t}_{0}}+\mathfrak{J}\phantom{\rule{0.166667em}{0ex}}\frac{\partial}{\partial {t}_{1}}+O({\mathfrak{J}}^{2})\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(27)
The multiple time scales method then searches for an approximate solution for the complex amplitudes c_{φ} of mode φ as
$$\begin{array}{c}\hfill {c}_{\phi}=\mathfrak{J}\phantom{\rule{0.166667em}{0ex}}{c}_{\phi 1}({t}_{0},{t}_{1},\dots )+{\mathfrak{J}}^{2}\phantom{\rule{0.166667em}{0ex}}{c}_{\phi 2}({t}_{0},{t}_{1},\dots )+O({\mathfrak{J}}^{3})\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(28)
Perturbation series (28) needs to be uniform: each of the higherorder terms should be a small correction to lowerorder terms.
Substituting expansions (27) and (28) into the coupledmode equations (20), subsequently dividing by the nonzero parameter 𝔍 and keeping terms up to first order in 𝔍, yields
$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\partial {c}_{\phi 1}}{\partial {t}_{0}}+\mathfrak{J}\frac{\partial {c}_{\phi 1}}{\partial {t}_{1}}+\mathfrak{J}\frac{\partial {c}_{\phi 2}}{\partial {t}_{0}}+i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}({c}_{\phi 1}+\mathfrak{J}{c}_{{\phi}_{2}})\hfill \\ \hfill & =\mathfrak{J}\frac{i}{{b}_{\phi}}({\kappa}_{\phi}^{\beta \gamma}{c}_{\beta 1}{c}_{\gamma 1}+{\kappa}_{\phi}^{\overline{\beta}\gamma}{c}_{\beta 1}^{\ast}{c}_{\gamma 1}+{\kappa}_{\phi}^{\beta \overline{\gamma}}{c}_{\beta 1}{c}_{\gamma 1}^{\ast}+{\kappa}_{\phi}^{\overline{\beta}\overline{\gamma}}{c}_{\beta 1}^{\ast}{c}_{\gamma 1}^{\ast})\hfill \\ \hfill & +\phantom{\rule{0.166667em}{0ex}}O({\mathfrak{J}}^{2})\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}\end{array}$$(29)
We examine this equation order by order in 𝔍 to determine the complex amplitude variation of the eigenmode φ in the form of AEs. The (linear) equation at order 𝔍^{0} is
$$\begin{array}{c}\hfill \frac{\partial {c}_{\phi 1}}{\partial {t}_{0}}+i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{c}_{\phi 1}\equiv {L}_{\phi}({c}_{\phi 1})=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(30)
in which we introduce the linear operator L_{φ} as the mapping L_{φ}: ℂ → ℂ defined as ${L}_{\phi}(c)=(\frac{\partial}{\partial {t}_{0}}+i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi})[c]$ for a complex amplitude c. The solution of Eq. (30) is
$$\begin{array}{c}\hfill {c}_{\phi 1}={a}_{\phi}({t}_{1},\dots ){e}^{i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{t}_{0}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(31)
where a_{φ} is a complex amplitude factor that varies only on time scales slower than the time scale t_{0} for each mode φ.
The equation at order 𝔍 is
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {L}_{\phi}({c}_{\phi 2})=\frac{\partial {c}_{\phi 1}}{\partial {t}_{1}}+\frac{i}{{b}_{\phi}}& ({\kappa}_{\phi}^{\beta \gamma}{c}_{\beta 1}{c}_{\gamma 1}+{\kappa}_{\phi}^{\overline{\beta}\gamma}{c}_{\beta 1}^{\ast}{c}_{\gamma 1}\hfill \\ \hfill & +\phantom{\rule{0.166667em}{0ex}}{\kappa}_{\phi}^{\beta \overline{\gamma}}{c}_{\beta 1}{c}_{\gamma 1}^{\ast}+{\kappa}_{\phi}^{\overline{\beta}\overline{\gamma}}{c}_{\beta 1}^{\ast}{c}_{\gamma 1}^{\ast})\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}\end{array}$$(32)
Substituting the linear solution (31) into Eq. (32) yields
$$\begin{array}{cc}& {L}_{\phi}({c}_{\phi 2})=\frac{\partial {a}_{\phi}}{\partial {t}_{1}}\phantom{\rule{0.166667em}{0ex}}{e}^{i\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{t}_{0}}\hfill \\ \hfill & +\frac{i}{{b}_{\phi}}({a}_{\beta}\phantom{\rule{0.166667em}{0ex}}{a}_{\gamma}\phantom{\rule{0.166667em}{0ex}}{\kappa}_{\phi}^{\beta \gamma}{e}^{i({\mathrm{\Omega}}_{\beta}+{\mathrm{\Omega}}_{\gamma}){t}_{0}}+{a}_{\beta}^{\ast}\phantom{\rule{0.166667em}{0ex}}{a}_{\gamma}\phantom{\rule{0.166667em}{0ex}}{\kappa}_{\phi}^{\overline{\beta}\gamma}{e}^{i({\mathrm{\Omega}}_{\beta}+{\mathrm{\Omega}}_{\gamma}){t}_{0}}\hfill \\ \hfill & +\phantom{\rule{0.166667em}{0ex}}{a}_{\beta}\phantom{\rule{0.166667em}{0ex}}{a}_{\gamma}^{\ast}\phantom{\rule{0.166667em}{0ex}}{\kappa}_{\phi}^{\beta \overline{\gamma}}{e}^{i({\mathrm{\Omega}}_{\beta}{\mathrm{\Omega}}_{\gamma}){t}_{0}}+{a}_{\beta}^{\ast}\phantom{\rule{0.166667em}{0ex}}{a}_{\gamma}^{\ast}\phantom{\rule{0.166667em}{0ex}}{\kappa}_{\phi}^{\overline{\beta}\overline{\gamma}}{e}^{i({\mathrm{\Omega}}_{\beta}+{\mathrm{\Omega}}_{\gamma}){t}_{0}})\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}$$(33)
Some of the terms in the solution of Eq. (33) increase linearly in time. Such terms are called secular terms (e.g., Nayfeh 1973, 1981 and Nayfeh & Mook 1979) and must vanish to ensure we obtain a uniformly valid perturbation series.
Setting the terms that generate the secular terms in the solution of Eq. (33) equal to zero yields the AEs for the mode amplitudes and phases. These seculartermgenerating terms have a multiplying factor exp(−iΩ_{φ}t_{0}) in Eq. (33), as shown in Appendix B.1 of Van Beeck (2023). The first term on the right side of Eq. (33) always generates a secular term. We obtain additional terms that generate secular terms by substituting the resonance condition Ω_{1} ≃ Ω_{2} + Ω_{3}, defined in terms of the detuning parameters δω and ΔΩ^{l},
$$\begin{array}{c}\hfill \mathfrak{J}\phantom{\rule{0.166667em}{0ex}}\delta \omega ={\mathrm{\Omega}}_{1}{\mathrm{\Omega}}_{2}{\mathrm{\Omega}}_{3}\equiv \mathrm{\Delta}{\mathrm{\Omega}}^{l}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(34)
in Eq. (33). The condition that ensures that the sum of these seculartermgenerating terms vanishes leads to the (extended complex) AEs for the sumfrequency resonance Ω_{1} ≃ Ω_{2} + Ω_{3}:
$$\begin{array}{c}\hfill \frac{\partial \mathit{a}}{\partial {t}_{1}}=\mathit{\gamma}\xb0\mathit{a}+2\phantom{\rule{0.166667em}{0ex}}i\phantom{\rule{0.166667em}{0ex}}(\mathit{\eta}\xb0{\mathit{a}}_{M}\xb0{\mathbf{\Omega}}_{M}\xb0\mathit{e})\phantom{\rule{0.166667em}{0ex}},\end{array}$$(35)
in which we introduce the linear growth or linear damping rates γ_{φ} for mode φ (e.g., Van Hoolst 1995) and use Eq. (22) to convert ${\kappa}_{\phi}^{\beta \gamma}$ into ${\eta}_{\phi}^{\beta \gamma}$ (and similar). In Eq. (35), we also introduce the symbol ‘°’, which denotes a HadamardSchur (i.e., elementwise, Bhatia 2007; Bernstein 2009) product. We use the vectors
$$\begin{array}{c}\hfill \mathit{a}=\left[\begin{array}{c}{a}_{1}\\ {a}_{2}\\ {a}_{3}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{a}}_{M}=\left[\begin{array}{c}{a}_{2}\phantom{\rule{0.166667em}{0ex}}{a}_{3}\\ {a}_{1}\phantom{\rule{0.166667em}{0ex}}{a}_{3}^{\ast}\\ {a}_{1}\phantom{\rule{0.166667em}{0ex}}{a}_{2}^{\ast}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\mathit{e}=\left[\begin{array}{c}exp(i\phantom{\rule{0.166667em}{0ex}}\delta \omega \phantom{\rule{0.166667em}{0ex}}{t}_{1})\\ exp(i\phantom{\rule{0.166667em}{0ex}}\delta \omega \phantom{\rule{0.166667em}{0ex}}{t}_{1})\\ exp(i\phantom{\rule{0.166667em}{0ex}}\delta \omega \phantom{\rule{0.166667em}{0ex}}{t}_{1})\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(36)
and
$$\begin{array}{c}\hfill \mathit{\gamma}=\left[\begin{array}{c}{\gamma}_{1}\\ {\gamma}_{2}\\ {\gamma}_{3}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\mathit{\eta}=\left[\begin{array}{c}{\eta}_{1}\\ {\eta}_{1}^{\ast}\\ {\eta}_{1}^{\ast}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathbf{\Omega}}_{M}=\left[\begin{array}{c}{\mathrm{\Omega}}_{1}\\ {\mathrm{\Omega}}_{2}\\ {\mathrm{\Omega}}_{3}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(37)
in Eq. (35). In Eq. (37), we define the isolated threemode coupling coefficient η_{1}:
$$\begin{array}{c}\hfill {\eta}_{1}^{23}={\eta}_{1}^{32}={\left({\eta}_{2}^{1\overline{3}}\right)}^{\ast}={\left({\eta}_{2}^{\overline{3}1}\right)}^{\ast}={\left({\eta}_{3}^{1\overline{2}}\right)}^{\ast}={\left({\eta}_{3}^{\overline{2}1}\right)}^{\ast}\equiv {\eta}_{1}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(38)
which respects the symmetry of coupling coefficient (22) under permutation of its indices. For exact resonances (i.e., δω = 0), the AEs reduce to a set of trivial equations. In such situations, efficient nonlinear energy transfer is expected, and one enters the regime of intermediate to strong nonlinear coupling, whereas the AEs are valid for weak coupling only.
By introducing real amplitudes A_{φ} and phases ϕ_{φ} as a_{φ} = A_{φ}exp(i ϕ_{φ}) with φ ∈ ⟦3⟧, where ⟦u⟧ denotes the set of integers { x∈ℕ_{0}  x≤u }, and separating the real and imaginary parts of the extended complex AEs (35), we obtain
$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}}{\partial {t}_{1}}& =\mathit{\gamma}\xb0\mathit{A}+2\phantom{\rule{0.166667em}{0ex}}{\eta}_{1}\phantom{\rule{0.166667em}{0ex}}sin\left(\mathrm{{\rm Y}}\right)\phantom{\rule{0.166667em}{0ex}}\left({\mathit{A}}_{N}^{\xb0}{\mathbf{\Omega}}_{M}\right)\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(39a)
$$\begin{array}{cc}\hfill \mathit{A}\xb0\frac{\partial \mathit{\varphi}}{\partial {t}_{1}}& =2\phantom{\rule{0.166667em}{0ex}}{\eta}_{1}\phantom{\rule{0.166667em}{0ex}}cos\left(\mathrm{{\rm Y}}\right)\phantom{\rule{0.166667em}{0ex}}({\mathit{A}}_{P}\xb0{\mathbf{\Omega}}_{M})\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(39b)
where η_{1} = η_{1} e^{−i δ1}, and
$$\begin{array}{c}\hfill \mathit{A}=\left[\begin{array}{c}{A}_{1}\\ {A}_{2}\\ {A}_{3}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\mathit{\varphi}=\left[\begin{array}{c}{\varphi}_{1}\\ {\varphi}_{2}\\ {\varphi}_{3}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{A}}_{N}=\left[\begin{array}{c}{A}_{2}\phantom{\rule{0.166667em}{0ex}}{A}_{3}\\ {A}_{1}\phantom{\rule{0.166667em}{0ex}}{A}_{3}\\ {A}_{1}\phantom{\rule{0.166667em}{0ex}}{A}_{2}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{A}}_{P}=\left[\begin{array}{c}{A}_{2}\phantom{\rule{0.166667em}{0ex}}{A}_{3}\\ {A}_{1}\phantom{\rule{0.166667em}{0ex}}{A}_{3}\\ {A}_{1}\phantom{\rule{0.166667em}{0ex}}{A}_{2}\end{array}\right]\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(40)
In the AEs (39) we also introduce the generic phase coordinate
$$\begin{array}{c}\hfill \mathrm{{\rm Y}}\equiv \delta \omega \phantom{\rule{0.166667em}{0ex}}{t}_{1}+{\varphi}_{1}{\varphi}_{2}{\varphi}_{3}+{\delta}_{1}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(41)
which contains the combination phase Φ = ϕ_{1} − ϕ_{2} − ϕ_{3}. Because the coupling coefficients (22) do not depend on time, we have
$$\begin{array}{c}\hfill \frac{\partial \mathrm{{\rm Y}}}{\partial {t}_{1}}=\delta \omega +cot\left(\mathrm{{\rm Y}}\right)({\gamma}_{\u229e}+\frac{\partial ln{A}_{1}}{\partial {t}_{1}}+\frac{\partial ln{A}_{2}}{\partial {t}_{1}}+\frac{\partial ln{A}_{3}}{\partial {t}_{1}})\phantom{\rule{0.166667em}{0ex}},\end{array}$$(42)
for nonzero resonant mode amplitudes, because of Eq. (39a) and the definition γ_{⊞} ≡ γ_{1} + γ_{2} + γ_{3}. Equations (39a) and (42) thus form an autonomous fourdimensional system equivalent to the sixdimensional system in Eqs. (39a) and (39b), in which the individual mode phases ϕ_{1}, ϕ_{2} and ϕ_{3} do not explicitly appear.
Nonlinear interactions lead to nonlinear frequency shifts, which can be derived from Eq. (39b) for the individual mode phases. By integrating these equations, the first order solution (31) can be expressed as
$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathit{c}}_{1}=\mathit{A}({t}_{1},\dots )\hfill \\ \hfill & exp\left\{i({\mathit{\varphi}}_{0}{\mathbf{\Omega}}_{M}\phantom{\rule{0.166667em}{0ex}}{t}_{0}+2{\displaystyle \int ({\mathit{A}}_{P}\oslash \mathit{A}\xb0{\mathbf{\Omega}}_{M})}{\eta}_{1}cos\left(\mathrm{{\rm Y}}\right)d{t}_{1})\right\}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}\end{array}$$(43)
in which the symbol ‘⊘’ denotes a HadamardSchur (elementwise) division, and where
$$\begin{array}{c}\hfill {\mathit{\varphi}}_{0}=\left[\begin{array}{c}{\left({\varphi}_{1}\right)}_{0}\\ {\left({\varphi}_{2}\right)}_{0}\\ {\left({\varphi}_{3}\right)}_{0}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{c}}_{1}=\left[\begin{array}{c}{c}_{11}\\ {c}_{21}\\ {c}_{31}\\ \end{array}\right]\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(44)
The third term in the exponential factor of Eq. (43) defines the quadratic nonlinear frequency shift multiplied with the elapsed time. If the linear angular mode frequency in the corotating frame – the second term in that factor – is added to the nonlinear frequency shift, we obtain a shifted frequency, hereafter referred to as the nonlinear frequency, in the same reference frame. The expression for the (secondorder) nonlinear correction to the complex amplitude, c_{φ2}, can be found in Appendix B.1 of Van Beeck (2023). The harmonic resonance equivalents of the nonlinear frequency shifts and secondorder nonlinear corrections to the complex amplitudes can be found in Appendix B.2 of that work.
2.7. Stationary solutions of the amplitude equations
The candidate resonant threemode couplings identified by V21 are timeindependent, with the modes in a nonlinear frequency (and phase) lock (i.e., they are synchronized). We therefore derive the timeindependent or stationary solutions of the AEs (i.e., the fixed points of Eqs. (39a) and (42)) in this section.
The trivial stationary solution A^{s} = 0, in which the superscript ‘s’ denotes the stationarity of a quantity, is not oscillatory, and therefore cannot explain the observed couplings. From Eq. (39a), we obtain a oscillatory stationary solution,
$$\begin{array}{c}\hfill \mathit{\gamma}\xb0{\mathit{A}}^{s}=2\phantom{\rule{0.166667em}{0ex}}{\eta}_{1}\phantom{\rule{0.166667em}{0ex}}sin\left({\mathrm{{\rm Y}}}^{s}\right)\phantom{\rule{0.166667em}{0ex}}({\mathit{A}}_{N}^{s}\xb0{\mathbf{\Omega}}_{M})\phantom{\rule{0.166667em}{0ex}},\end{array}$$(45)
in which we define the stationary phase coordinate Υ^{s} as
$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}^{s}=\delta \omega \phantom{\rule{0.166667em}{0ex}}{t}_{1}+{\varphi}_{1}^{s}{\varphi}_{2}^{s}{\varphi}_{3}^{s}+{\delta}_{1}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(46)
with $\frac{\partial {\mathrm{{\rm Y}}}^{s}}{\partial {t}_{1}}=0$. The stationary equivalent of Eq. (42) yields an expression for the detuning parameter of stationary solutions,
$$\begin{array}{c}\hfill \delta \omega =cot\left({\mathrm{{\rm Y}}}^{s}\right)\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\u229e}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(47)
by setting $\frac{\partial {\mathit{A}}^{s}}{\partial {t}_{1}}=\mathbf{0}$ and $\frac{\partial {\mathrm{{\rm Y}}}^{s}}{\partial {t}_{1}}=0$. To derive this equation, we assume that Υ^{s} ≠ p π (p ∈ ℤ). If Υ^{s} = p π, we retrieve the trivial solution, due to Eq. (45).
By considering the products of two equations of the three in Eq. (45), and using Eq. (47), we write the squared stationary amplitudes as
$$\begin{array}{c}\hfill {\left({\mathit{A}}^{s}\right)}^{\xb02}=\frac{\mathit{Q}}{4\phantom{\rule{0.166667em}{0ex}}{\eta}_{1}{}^{2}}[1+{q}^{2}]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(48)
where the superscript ‘°n’ indicates a HadamardSchur (elementwise) nth power. In Eq. (48), we use the detuningdamping ratio q and quality factor vector Q, which we define as
$$\begin{array}{c}\hfill \mathit{Q}=\left[\begin{array}{c}{\displaystyle {\gamma}_{2}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{3}/{\mathrm{\Omega}}_{2}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{3}}\\ {\displaystyle {\gamma}_{1}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{3}/{\mathrm{\Omega}}_{1}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{3}}\\ {\displaystyle {\gamma}_{1}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{2}/{\mathrm{\Omega}}_{1}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{2}}\\ \end{array}\right]=\left[\begin{array}{c}{\displaystyle 1/{Q}_{2}\phantom{\rule{0.166667em}{0ex}}{Q}_{3}}\\ {\displaystyle 1/{Q}_{1}\phantom{\rule{0.166667em}{0ex}}{Q}_{3}}\\ {\displaystyle 1/{Q}_{1}\phantom{\rule{0.166667em}{0ex}}{Q}_{2}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\phantom{\rule{0.166667em}{0ex}}q\equiv \frac{\delta \omega}{{\gamma}_{\u229e}}=cot{\mathrm{{\rm Y}}}^{s}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(49)
with quality factor Q_{φ} ≡ Ω_{φ}/γ_{φ} for φ ∈ ⟦3⟧. The resonant stationary amplitudes are real if
$$\begin{array}{c}\hfill {\mathbf{\Omega}}_{N}{>}_{\xb0}0\phantom{\rule{0.166667em}{0ex}}\vee \phantom{\rule{0.166667em}{0ex}}{\mathbf{\Omega}}_{N}{<}_{\xb0}0\Rightarrow {\mathit{A}}^{s}{\in}_{\xb0}\mathbb{R}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(50)
for a driven resonance scenario, or if
$$\begin{array}{c}\hfill {\mathbf{\Omega}}_{P}{>}_{\xb0}0\phantom{\rule{0.166667em}{0ex}}\vee \phantom{\rule{0.166667em}{0ex}}{\mathbf{\Omega}}_{P}{<}_{\xb0}0\Rightarrow {\mathit{A}}^{s}{\in}_{\xb0}\mathbb{R}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(51)
for direct and parametric resonance scenarios. In Eqs. (50) and (51), Ω_{N} and Ω_{P} are the mode frequency variants of the vectors A_{N} and A_{P} defined in Eq. (40). The operators > _{°} and ∈_{°} are the elementwise equivalents of the operators > and ∈. Equivalent expressions describing stationary solutions for harmonic resonances can be found in Appendix B.2 of Van Beeck (2023).
Three quantities determine the stationary mode amplitudes: the nonlinear interaction described by coupledmode Eqs. (20) and (21), the linear eigenmode properties given by Q, and the detuning of the resonance measured by q. For decreasing values of the nonlinear coupling coefficients η_{1}, the efficiency for nonlinear energy transfer between modes is lower, hence, the stationary amplitudes (48) become larger, because a larger mode energy (and thus mode amplitude) is required to transfer enough energy in order to have a significant nonlinear effect leading to amplitude saturation. Quality factor Q_{φ} expresses the ratio of the mode φ’s efolding damping or driving time (${\gamma}_{\phi}^{1}$) relative to its angular period (${\mathrm{\Omega}}_{\phi}^{1}$). With increasing Q_{φ}, more mode periods are needed for the mode energy (normalized to GM^{2}/R) to damp or grow by a factor of e due to linear heatdriven damping or driving. The stationary amplitude of a certain mode therefore decreases if the other modes involved in the triad have larger values of Q_{φ}, because then either a linearly (heat)driven mode φ gains less energy per cycle or a linearly (heat)damped mode φ loses energy more slowly per cycle. For example, in the parametric resonance scenario, the stationary amplitude of the daughter mode 2 or 3 is increased when the parent mode is driven faster (by the κ mechanism) and when the other daughter mode 3 or 2 is damped faster. The stationary amplitude of the parent mode 1 is larger when it is coupled to daughter modes that are more difficult to nonlinearly excite (i.e., a larger parent mode energy is required to have a measurable nonlinear effect).
Increasing values of q lead to increasing values of the stationary amplitudes. Two factors affect the value of the detuningdamping ratio q: the detuning δω and γ_{⊞} (≡γ_{1} + γ_{2} + γ_{3}). For larger detunings the stationary amplitudes increase, because the increased detuning reduces the efficiency of the nonlinear mode coupling (similar to decreasing η_{1}), requiring a larger amplitude for a nonlinear effect. For smaller values of γ_{⊞}, that is, for the case where the parent and daughter mode linear heatdriven growth and damping rates almost balance out, there is very weak overall linear excitation or damping, requiring a very close resonance (i.e., small detuning) and large amplitudes to have a nonlinear effect.
The relative stationary energies of the modes in the triad are related by the ratios of their quality factors, assuming conditions (50) or (51) hold:
$$\begin{array}{c}\hfill {\left(\frac{{A}_{2}^{s}}{{A}_{1}^{s}}\right)}^{2}=\frac{{Q}_{2}}{{Q}_{1}}\phantom{\rule{0.166667em}{0ex}},{\left(\frac{{A}_{3}^{s}}{{A}_{1}^{s}}\right)}^{2}=\frac{{Q}_{3}}{{Q}_{1}}\phantom{\rule{0.166667em}{0ex}},{\left(\frac{{A}_{3}^{s}}{{A}_{2}^{s}}\right)}^{2}=\frac{{Q}_{3}}{{Q}_{2}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(52)
The energy ratios (52) do not depend on the nonlinear coupling coefficients because of the symmetry of these coefficients, and are determined by linear properties of the modes only. Therefore, of two linearly damped modes, the one with the longest damping time per cycle (largest Q_{φ}) reaches the largest stationary amplitude through nonlinear energy transfer from the linearly excited mode, because it loses energy more slowly.
Nonlinear stationary frequencies Ω^{nl, s} for modes φ ∈ ⟦3⟧ are obtained from Eq. (43) in the form
$$\begin{array}{c}\hfill {\mathrm{\Omega}}_{\phi}^{nl,s}={\mathrm{\Omega}}_{\phi}+\delta {\mathrm{\Omega}}_{\phi}^{s}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(53)
with the nonlinear stationary frequency shifts $\delta {\mathrm{\Omega}}_{\phi}^{s}$ given by
$$\begin{array}{c}\hfill \left[\begin{array}{c}\delta {\mathrm{\Omega}}_{1}^{s}\\ \delta {\mathrm{\Omega}}_{2}^{s}\\ \delta {\mathrm{\Omega}}_{3}^{s}\\ \end{array}\right]=\mathfrak{J}\phantom{\rule{0.166667em}{0ex}}q\phantom{\rule{0.166667em}{0ex}}\left[\begin{array}{c}{\gamma}_{1}\\ {\gamma}_{2}\\ {\gamma}_{3}\\ \end{array}\right]=\mathrm{\Delta}{\mathrm{\Omega}}^{l}\phantom{\rule{0.166667em}{0ex}}\left[\begin{array}{c}{\gamma}_{1}/{\gamma}_{\u229e}\\ {\gamma}_{2}/{\gamma}_{\u229e}\\ {\gamma}_{3}/{\gamma}_{\u229e}\\ \end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(54)
to first order. For parametric and direct resonance scenarios, the sign of the frequency shift is thus determined by the sign of ΔΩ^{l}. The nonlinear stationary frequencies ${\mathrm{\Omega}}_{\phi}^{nl,s}$ are frequencylocked, meaning that they satisfy the resonance condition (34) exactly:
$$\begin{array}{c}\hfill {\mathrm{\Omega}}_{1}^{nl,s}{\mathrm{\Omega}}_{2}^{nl,s}{\mathrm{\Omega}}_{3}^{nl,s}=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(55)
which follows from Eq. (46) and $\frac{\partial {\mathrm{{\rm Y}}}^{s}}{\partial {t}_{1}}=0$. Explicit expressions for harmonic resonance frequency shifts can be found in Appendix B.2 of Van Beeck (2023). As stated in Sect. 2.6, the derived stationary quantities in this Section also describe the stationary properties of modes in differencefrequency resonances.
2.8. Stability of the stationary amplitude equation solution
From a mathematical perspective, the system of amplitude equations defined by Eqs. (39a) and (42) is an autonomous dynamical system, and its stationary solutions are fixed points of that dynamical system. One of the fundamental results of dynamical systems theory is the HartmanGrobman or linearization theorem, which states that if a fixed point is hyperbolic, a linearization of the dynamical system can be used to trace the asymptotic behavior of dynamical system solutions near the fixed point (Guckenheimer & Holmes 1983; Betounes 2010). The stability of a hyperbolic fixed point can thus be determined from the linearization of the system around that point. Two conditions need to be fulfilled for a fixed point $\overline{\mathsf{X}}$ of a dynamical system $\frac{\partial \mathsf{X}}{\partial t}=f(\mathsf{X})$ to be hyperbolic: both the Jacobian and the real parts of the eigenvalues of the Jacobian matrix of $f(\overline{\mathsf{X}})$ cannot be equal to zero (Guckenheimer & Holmes 1983; Betounes 2010). Hereafter, we collectively refer to these two conditions as the hyperbolicity condition. The Jacobian matrix of a dynamical system of amplitude equations depends on the values of the respective nonlinear coupling coefficients because the hyperbolicity condition considers the geometry of the nonlinear solution around the fixed point (see Appendices B.1 and B.2 of Van Beeck 2023 for the expressions).
In the remainder of this section we assume that the hyperbolicity condition is fulfilled for the stationary solution, so that a linearization can be used to guarantee the stability of that solution. Let us now examine the stability of the stationary solutions for each of the three isolated resonant coupling scenarios discussed in Sect. 2.5. Physically, the stability of the (hyperbolic) stationary solutions of the AEs is governed by the reaction of the system (i.e., the pulsating star) to small disturbances away from the stationary state (in this case, away from the stationary pulsation solution). Such disturbances can be damped or get amplified over time. If they are damped, the stationary solution is stable. The solution is unstable if they get amplified.
To derive the linearized dynamical system, we introduce small perturbations of the stationary amplitudes (δA_{φ} for mode φ) and phases (δϕ_{φ} for mode φ), so that the complex amplitudes take the form (Van Hoolst 1995)
$$\begin{array}{c}\hfill {a}_{\phi}=({A}_{\phi}^{s}\phantom{\rule{0.166667em}{0ex}}+\phantom{\rule{0.166667em}{0ex}}\delta {A}_{\phi})exp\left[i({\varphi}_{\phi}^{s}\phantom{\rule{0.166667em}{0ex}}+\phantom{\rule{0.166667em}{0ex}}\delta {\varphi}_{\phi})\right],\forall \phantom{\rule{0.166667em}{0ex}}\phi \in [\phantom{\rule{0.166667em}{0ex}}[3]\phantom{\rule{0.166667em}{0ex}}]\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(56)
Substituting the perturbed complex amplitude factor (56) into the complexvalued AEs (35), subsequently linearizing and separating the realvalued and complexvalued parts, yields the tensor equation that describes the time evolution of the perturbations,
$$\begin{array}{c}\hfill \frac{\partial \mathsf{Z}}{\partial t}=\mathsf{M}\phantom{\rule{0.166667em}{0ex}}\mathsf{Z}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(57)
Here, M is the stability matrix, defined as
$$\begin{array}{c}\hfill \mathsf{M}\equiv \left(\begin{array}{cccc}{\gamma}_{1}& {\gamma}_{1}& {\gamma}_{1}& q\phantom{\rule{0.166667em}{0ex}}{\gamma}_{1}\\ {\gamma}_{2}& {\gamma}_{2}& {\gamma}_{2}& q\phantom{\rule{0.166667em}{0ex}}{\gamma}_{2}\\ {\gamma}_{3}& {\gamma}_{3}& {\gamma}_{3}& q\phantom{\rule{0.166667em}{0ex}}{\gamma}_{3}\\ q\phantom{\rule{0.166667em}{0ex}}{\gamma}_{{\u229e}_{1}}& q\phantom{\rule{0.166667em}{0ex}}{\gamma}_{{\u229e}_{2}}& q\phantom{\rule{0.166667em}{0ex}}{\gamma}_{{\u229e}_{3}}& {\gamma}_{\u229e}\\ \end{array}\right)\phantom{\rule{0.166667em}{0ex}},\end{array}$$(58)
and Z is the perturbation tensor defined as
$$\begin{array}{c}\hfill \mathsf{Z}={(\delta {A}_{1}/{A}_{1}^{s}\phantom{\rule{0.166667em}{0ex}},\phantom{\rule{0.166667em}{0ex}}\delta {A}_{2}/{A}_{2}^{s}\phantom{\rule{0.166667em}{0ex}},\phantom{\rule{0.166667em}{0ex}}\delta {A}_{3}/{A}_{3}^{s}\phantom{\rule{0.166667em}{0ex}},\phantom{\rule{0.166667em}{0ex}}\delta \mathrm{{\rm Y}})}^{\mathrm{T}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(59)
for which
$$\begin{array}{c}\hfill \delta \mathrm{{\rm Y}}=\delta {\varphi}_{1}\delta {\varphi}_{2}\delta {\varphi}_{3}\phantom{\rule{0.166667em}{0ex}},{\gamma}_{{\u229e}_{\phi}}\equiv {\gamma}_{\u229e}2\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\phi},\forall \phi \in [\phantom{\rule{0.166667em}{0ex}}[3]\phantom{\rule{0.166667em}{0ex}}]\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(60)
We infer an analytical stability domain of the stationary solutions of the AEs based on the RouthHurwitz stability criterion (e.g., Hahn 1967; Cesari 1971; Kubicek & Marek 1983). Necessary but not sufficient conditions for stability are that the Hurwitz determinants H_{u} are positive (see Appendix B.4 in Van Beeck 2023):
$$\begin{array}{cc}\hfill {H}_{1}=& \phantom{\rule{0.166667em}{0ex}}{w}_{1}=2\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\u229e}>0\iff {\gamma}_{\u229e}<0\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(61a)
$$\begin{array}{cc}\hfill {H}_{2}=& \phantom{\rule{0.166667em}{0ex}}{w}_{1}\phantom{\rule{0.166667em}{0ex}}{w}_{2}{w}_{3}>0\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(61b)
$$\begin{array}{cc}\hfill {H}_{3}=& \phantom{\rule{0.166667em}{0ex}}{w}_{3}\phantom{\rule{0.166667em}{0ex}}{H}_{2}{w}_{1}^{2}\phantom{\rule{0.166667em}{0ex}}{w}_{4}>0\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(61c)
$$\begin{array}{cc}\hfill {H}_{4}=& \phantom{\rule{0.166667em}{0ex}}{w}_{4}\phantom{\rule{0.166667em}{0ex}}{H}_{3}>0\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(61d)
in which the coefficients w_{u} (u ∈ ⟦4⟧_{0}, with ⟦n⟧_{0} ≡ {0} ∪ ⟦n⟧ = { x∈ℕ  x≤u }) are the coefficients of the characteristic polynomial p_{Z}(λ) of the stability matrix (58),
$$\begin{array}{c}\hfill {p}_{\mathsf{Z}}(\lambda )\equiv {\displaystyle \sum _{u=0}^{4}}{w}_{4u}\phantom{\rule{0.166667em}{0ex}}{\lambda}^{u}={w}_{4}+{w}_{3}\phantom{\rule{0.166667em}{0ex}}\lambda +{w}_{2}\phantom{\rule{0.166667em}{0ex}}{\lambda}^{2}+{w}_{1}\phantom{\rule{0.166667em}{0ex}}{\lambda}^{3}+{w}_{0}\phantom{\rule{0.166667em}{0ex}}{\lambda}^{4}=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(62)
In Eq. (62), λ denotes an eigenvalue of the tensor Eq. (57), and the coefficients w_{u} are
$$\begin{array}{cc}\hfill {w}_{0}=& \phantom{\rule{0.166667em}{0ex}}1\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(63a)
$$\begin{array}{cc}\hfill {w}_{1}=& \phantom{\rule{0.166667em}{0ex}}2\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\u229e}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(63b)
$$\begin{array}{cc}\hfill {w}_{2}=& \phantom{\rule{0.166667em}{0ex}}{\gamma}_{\u229e}^{2}[1+{q}^{2}]4\phantom{\rule{0.166667em}{0ex}}{q}^{2}\phantom{\rule{0.166667em}{0ex}}[{\gamma}_{1}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{2}+{\gamma}_{1}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{3}+{\gamma}_{2}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{3}]\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(63c)
$$\begin{array}{cc}\hfill {w}_{3}=& \phantom{\rule{0.166667em}{0ex}}4\phantom{\rule{0.166667em}{0ex}}{\gamma}_{1}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{2}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{3}\phantom{\rule{0.166667em}{0ex}}[1+3\phantom{\rule{0.166667em}{0ex}}{q}^{2}]\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(63d)
$$\begin{array}{cc}\hfill {w}_{4}=& \phantom{\rule{0.166667em}{0ex}}4\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\u229e}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{1}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{2}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{3}\phantom{\rule{0.166667em}{0ex}}[1+{q}^{2}]\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}$$(63e)
These are the same coefficients as the ones that were determined by Dziembowski (1982) in his formalism for quadratic nonlinear mode coupling of three distinct modes in nonrotating stars. Because H_{3} > 0 due to the stability condition (61c), stability condition (61d) becomes w_{4} > 0.
The necessary but not sufficient stability conditions in Eq. (61) show that a stationary solution at the quadratic coupling level is always unstable for the driven resonant coupling scenario, because Eq. (61a) and γ> _{°}0 are contradictory. This is a logical consequence of only considering modes that are linearly excited and can exchange energy between them, and not including thirdorder couplings that can limit the amplitude. The direct resonant threemode coupling scenario also always yields unstable stationary solutions for distinct modes 2 and 3, similar to what was found by Dziembowski (1982), because Eqs. (61a) and (61d) contradict each other if the necessary stability condition (61c) is to be satisfied. Hence, driven and direct resonance scenarios lead to timevariable amplitudes for their coupled modes. Modes in such coupling scenarios can for example lead to limit cycle behavior of the amplitudes (see e.g., Seydel 2009). For the parametric resonant threemode coupling scenario, Eq. (61d) is trivially fulfilled if Eqs. (61a) and (61c) are fulfilled. The stability of the fixed points of harmonic resonance analogues is discussed in Appendix B.2 of Van Beeck (2023).
Fig. 2. Stability domains as a function of ϑ_{2} and ϑ_{3}, indicated by the dark, shaded area, and evaluated using Eqs. (61a) and (64) for different values of q_{1}. Specifically, q_{1} is set equal to 0.1 (left top panel), 1 (middle top panel), 10 (right top panel), 10^{2} (left bottom panel), 10^{3} (middle bottom panel) and 10^{4} (right bottom panel). The stability domain determined by the stability conditions defined in Dziembowski (1982) is larger and is also indicated in these panels by the hatched areas. Unhatched white areas of the figure panels indicate unstable domains. Note the different axis scales for the different panels, which indicates the importance of the value of q_{1} in determining the stability of stationary solutions. 
In addition to Eq. (61), all w_{u} also need to be positive to obtain stable parametric stationary solutions (see e.g., Hahn 1967 and Appendix B.4 in Van Beeck 2023). Combining these additional conditions with Eq. (61), the stability domain of a parametric threemode resonant coupling with real amplitudes can be described by only three conditions: Eq. (61a), the hyperbolicity check, and the quartic condition
$$\begin{array}{c}\hfill \frac{{\gamma}_{1}^{3}}{{\gamma}_{\u229e}}[{\mathfrak{d}}_{4}+{\mathfrak{d}}_{2}\phantom{\rule{0.166667em}{0ex}}{q}^{2}+{\mathfrak{d}}_{0}\phantom{\rule{0.166667em}{0ex}}{q}^{4}]>0\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\iff \phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}{\mathfrak{d}}_{4}+{\mathfrak{d}}_{2}\phantom{\rule{0.166667em}{0ex}}{q}^{2}+{\mathfrak{d}}_{0}\phantom{\rule{0.166667em}{0ex}}{q}^{4}>0\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(64)
The factor ${\gamma}_{1}^{3}/{\gamma}_{\u229e}$ in Eq. (64) is positive for parametric resonances that fulfill Eq. (61a). The dimensionless quartic coefficients 𝔡_{u} (u = 0, 2, 4) in that equation are given by
$$\begin{array}{cc}& \begin{array}{cc}\hfill {\mathfrak{d}}_{0}=& \phantom{\rule{0.166667em}{0ex}}18\phantom{\rule{0.166667em}{0ex}}{\vartheta}_{2}{\vartheta}_{3}3(1{\vartheta}_{2}{\vartheta}_{3})\hfill \\ \hfill & [{(1{\vartheta}_{2}{\vartheta}_{3})}^{2}+4({\vartheta}_{2}+{\vartheta}_{3}{\vartheta}_{2}{\vartheta}_{3})]\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}\hfill \end{array}$$(65a)
$$\begin{array}{cc}& \begin{array}{cc}\hfill {\mathfrak{d}}_{2}=& \phantom{\rule{0.166667em}{0ex}}12\phantom{\rule{0.166667em}{0ex}}{\vartheta}_{2}{\vartheta}_{3}(1{\vartheta}_{2}{\vartheta}_{3})\hfill \\ \hfill & [2+{({\vartheta}_{2}{\vartheta}_{3})}^{2}+{({\vartheta}_{2}+{\vartheta}_{3})}^{2}]\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}\hfill \end{array}$$(65b)
$$\begin{array}{cc}& {\mathfrak{d}}_{4}=\phantom{\rule{0.166667em}{0ex}}2\phantom{\rule{0.166667em}{0ex}}{\vartheta}_{2}{\vartheta}_{3}+{(1{\vartheta}_{2}{\vartheta}_{3})}^{3}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(65c)
where we define ϑ_{2, 3} ≡ −γ_{2, 3} / γ_{1}, which are positive for parametric resonances. The physical meaning of these dimensionless ratios is similar to that of the quality factor ratios discussed in Sect. 2.7, but inverse: they compare the damping and/or driving time scales.
The quartic stability condition (64) is symmetric in ϑ_{2} and ϑ_{3}, in accordance with the symmetry of the coupling coefficient. Coefficients 𝔡_{2} and 𝔡_{4} are the same as in Eq. (6.14) of Dziembowski (1982). Although the coefficient 𝔡_{0} differs from the one given in that equation, likely due to an error in Dziembowski (1982), the stability condition is derived from the same characteristic polynomial coefficients defined in Eq. (63). The coefficients 𝔡_{u} (u = 0, 2, 4) are a function of ϑ_{2} and ϑ_{3} only, whereas q can be expanded as a function of ϑ_{2}, ϑ_{3}, and the ratio of the linear frequency detuning to the parent’s linear driving rate q_{1} ≡ δω/γ_{1}. We therefore can explore the stability domain of the (hyperbolic) stationary solutions by varying only three dimensionless ratios of linear variables: q_{1} = δω/γ_{1}, ϑ_{2} and ϑ_{3}. The dimensionless ratio q_{1} can be written as
$$\begin{array}{c}\hfill {q}_{1}\equiv {Q}_{1}\phantom{\rule{0.166667em}{0ex}}\left(\frac{\delta \omega}{{\mathrm{\Omega}}_{1}}\right)=\frac{{Q}_{1}\phantom{\rule{0.166667em}{0ex}}\delta \omega \phantom{\rule{0.166667em}{0ex}}}{2\phantom{\rule{0.166667em}{0ex}}\pi}{P}_{1}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(66)
where the cyclic corotating frame period P_{1} = 2 π / Ω_{1}. It is therefore a periodweighted combined measure of the driving time scale of the linearly excited parent mode (Q_{1}) and the efficiency with which nonlinear energy transfer occurs (δω). Hence, when comparing values of q_{1} for triads with linearly excited modes of similar period, one might expect that the triad with larger q_{1} has larger stationary mode amplitudes because of decreased efficiency of nonlinear energy transfer (larger δω) and/or smaller stationary mode energy ratios $({A}_{2,\phantom{\rule{thinmathspace}{0ex}}3}^{s}/{A}_{1}^{s}{)}^{2}$ (larger Q_{1}, when values of Q_{2, 3} are comparable).
Figure 2 displays the domains of stability of stationary solutions in the threedimensional phase space of the parameters q_{1}, ϑ_{2} and ϑ_{3}, covering commonly encountered parameter values when modeling mode interactions in SPB stars (see Sect. 4.3). The necessary condition (61a) for stability of the fixed point, which requires that ϑ_{2} + ϑ_{3} > 1, is clearly recognizable on the left and middle panels in the top row of Fig. 2. Physically, this expresses that daughter modes must be sufficiently damped compared to the linear excitation of the parent mode for stability, otherwise, resonant energy transfer will increase all amplitudes without bound.
The quartic stability condition (64) is more difficult to interpret. Term (1 − ϑ_{2} − ϑ_{3}) in that stability condition is always negative due to Eq. (61a) and can be interpreted as an effective totaldampingtolineardriving ratio γ_{⊞}/γ_{1}. If the absolute value of the ratio is large, the overall damping per unit of driving is large as well. Other terms are less straightforward to interpret. Hence, we base our interpretation of this stability condition primarily on the stability domains pictured in Fig. 2. The quartic stability condition derived in this work is stricter than that of Dziembowski (1982) (displayed as hatched areas in Fig. 2) due to the difference in value of the coefficient 𝔡_{0}. Condition (64) can also be expressed in terms of q_{1}, ϑ_{2} and ϑ_{3}. Hence, conclusions drawn from the visualized stability domains in the chosen phase space (shown in Fig. 2) can be related to the equivalent stability condition that is expressed in terms of these three variables.
For large values of ϑ_{2} and ϑ_{3} a regime of strong damping (relative to the linear excitation of the parent mode) is reached. In such strong damping regimes the fixed point solutions are unstable because the transfer of energy (per cycle) is not enough to overcome linear damping. The domain of stability at the strong damping end of the ϑ_{2} − ϑ_{3} plane moves towards larger values of ϑ_{2} and ϑ_{3} for larger values of q_{1}, as shown in the different panels of Fig. 2. Conversely, the stability domain decreases considerably in size for smaller values of q_{1}. This can be explained by an increase (decrease) in energy transfer efficiency and/or an increase (decrease) in energy available for transfer, leading to faster (slower) rates of energy transfer and corresponding smaller (larger) domains of stability, based on the physical explanation of Eq. (66). Specifically, for faster (slower) rates of energy transfer per cycle, the stationary solutions can endure smaller (larger) perturbations around the stationary solutions before energy transfer renders the fixed points unstable. The stability of the stationary solutions is thus primarily determined by the speed of energy transfer.
2.9. The onset of parametric instability
In this section we derive the minimum amplitude conditions for the onset of the parametric resonance instability for the nonlinear interaction among three distinct modes in a sum or difference frequency. If this instability does not occur for these three distinct modes, the amplitudes of the modes in the triad can only be limited to stable stationary values by higherorder nonlinear mode coupling terms, such as the cubic selfcoupling terms that were described in Van Hoolst (1996). Alternatively, limit cycles characterized by nonstationary amplitudes may occur. Similar conditions for the onset of parametric and direct resonant instability derived for harmonic resonances can be found in Appendix B.2 of Van Beeck (2023).
The initial growth of the daughter modes can be described using the complex AEs (35) for one of the daughter modes and its complex conjugate for the other daughter mode. If we then set a_{k} = S_{k} exp(−i δω t_{1}/2) (for k ∈ {2 , 3}, following Dziembowski 1982) the explicit timedependence disappears. Further assuming that the complex amplitude factor a_{1} stays constant, a plausible assumption in the initial phase of energy transfer, yields
$$\begin{array}{cc}\hfill \frac{\partial {S}_{2}}{\partial {t}_{1}}& =({\gamma}_{2}+i\phantom{\rule{0.166667em}{0ex}}\frac{\delta \omega}{2})\phantom{\rule{0.166667em}{0ex}}{S}_{2}+2\phantom{\rule{0.166667em}{0ex}}i\phantom{\rule{0.166667em}{0ex}}{\eta}^{\ast}\phantom{\rule{0.166667em}{0ex}}{a}_{1}\phantom{\rule{0.166667em}{0ex}}{S}_{3}^{\ast}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{2}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(67a)
$$\begin{array}{cc}\hfill \frac{\partial {S}_{3}^{\ast}}{\partial {t}_{1}}& =({\gamma}_{3}i\phantom{\rule{0.166667em}{0ex}}\frac{\delta \omega}{2})\phantom{\rule{0.166667em}{0ex}}{S}_{3}^{\ast}2\phantom{\rule{0.166667em}{0ex}}i\phantom{\rule{0.166667em}{0ex}}\eta \phantom{\rule{0.166667em}{0ex}}{a}_{1}^{\ast}\phantom{\rule{0.166667em}{0ex}}{S}_{2}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{3}\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}$$(67b)
Under the assumption that S_{k} ∼ exp(σ t_{1}) (for k ∈ {2,3}), Eq. (67) can be solved for the growth parameter σ, yielding
$$\begin{array}{c}\hfill \sigma =\frac{{\gamma}_{2}+{\gamma}_{3}\pm \sqrt{{({\gamma}_{3}{\gamma}_{2}i\phantom{\rule{0.166667em}{0ex}}\delta \omega )}^{2}+16{{\eta}_{1}}^{2}{A}_{1}^{2}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{2}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{3}}}{2}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(68)
equivalent to expressions given by Vandakurov (1981) and Dziembowski (1982). Parametric instability will occur if Re[σ] > 0, because this ensures growth of the daughter mode amplitudes A_{2} and A_{3}. At the onset of parametric instability, Re[σ] = 0. We therefore define the instability threshold amplitude for the parent mode A_{t} as the value of A_{1} for Re[σ] = 0. The growth parameter σ is then imaginary and we can set σ = p i, with p determined by solving Eq. (68):
$$\begin{array}{c}\hfill p=\frac{\delta \omega}{2}\frac{({\gamma}_{3}{\gamma}_{2})}{({\gamma}_{3}+{\gamma}_{2})}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(69)
Using this expression in the real part of the growth parameter (68) yields the parametric instability threshold amplitude
$$\begin{array}{c}\hfill {A}_{t}=\frac{1}{2\phantom{\rule{0.166667em}{0ex}}{\eta}_{1}\sqrt{{Q}_{2}\phantom{\rule{0.166667em}{0ex}}{Q}_{3}}}\sqrt{1+{\left(\frac{\delta \omega}{{\gamma}_{2}+{\gamma}_{3}}\right)}^{2}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(70)
Instability threshold amplitude (70) is equivalent to the ones derived in Dziembowski (1982), Wu & Goldreich (2001) and Arras et al. (2003).
Parametric resonant mode triad interactions require that the parent mode amplitude A_{1} ≥ A_{t}, and ${A}_{1}^{s}$ is always larger than A_{t}. In the limit of very small (nonzero) detuning δω, the threshold amplitude is solely dependent on the coupling coefficient and the quality factors. In that case, A_{t} increases with decreasing η_{1} and faster damping of the daughter modes (expressed by the quality factors) because both terms limit the amplitude growth of the daughter modes due to nonlinear energy transfer, thus requiring a larger parent mode energy for a visible nonlinear effect. A larger detuning increases the threshold amplitude because of the less efficient energy transfer.
3. Theoretically predicted observables
An important observable in linear gmode asteroseismic modeling is a gmode period spacing pattern (which are extensively described in the literature; see e.g., Aerts et al. 2018; Michielsen et al. 2021; Bowman & Michielsen 2021 for some recent examples of how they can be used to probe internal mixing). In this section we derive additional observables based on the theoretical AE formalism described in Sect. 2 and outline how to compare them to observed quantities.
An inherent assumption of our models is that the modes are coherent. That assumption is justified because the detected frequencies of variability in SPB stars have been observed to be stable with a frequency precision of order 10^{−7} d^{−1}, based on longterm groundbased photometric monitoring (De Cat & Aerts 2002). This is not necessarily the case for other pulsators: δ Sct stars, for example, show frequency and amplitude modulation in the majority of detected signals (Bowman et al. 2016). Amplitude and frequency modulation also occurs among g mode triplets in oscillating white dwarfs (see e.g., the pioneering study of the oscillating DB white dwarf star KIC08626021 by Zong et al. 2016b). Moreover, the amplitudes of the oscillating hot B subdwarf star KIC10139564 reveal that the modulation of its observed p modes is larger than that of its g modes (Zong et al. 2016a). Whether this trend is generic among oscillating stars remains to be verified.
3.1. Amplitudes: modelgenerated luminosity fluctuations
We cannot compare the theoretical stationary amplitudes derived in Sect. 2.7 with the surface luminosity amplitudes $\stackrel{\sim}{A}$ determined from observations. These theoretical stationary amplitudes must first be converted to the corresponding observables, the theoretical luminosity fluctuations at the stellar surface, 𝔏.
In this section we follow the approach of Fuller (2017) to compute a conversion factor o_{A} used to convert a theoretical amplitude into a theoretical luminosity fluctuation 𝔏 caused by an oscillation mode of a resonant triad. Analogous to Eq. (82) in Fuller (2017), we estimate the discaveraged luminosity fluctuation for a single SPB star due to g mode φ as
$$\begin{array}{c}\hfill \frac{\mathrm{\Delta}{L}_{\phi}}{L}={\left(\frac{\mathrm{\Delta}{L}_{\phi}}{L}\right)}_{\mathrm{mode}}\phantom{\rule{0.166667em}{0ex}}{H}_{r}({i}_{s}\u037ek)\phantom{\rule{0.166667em}{0ex}},\end{array}$$(71)
within the TAR, where
$$\begin{array}{c}\hfill {\left(\frac{\mathrm{\Delta}{L}_{\phi}}{L}\right)}_{\mathrm{mode}}={t}_{k}\phantom{\rule{0.166667em}{0ex}}\frac{\mathrm{\Delta}{L}_{\phi ,\phantom{\rule{0.166667em}{0ex}}R}(R)}{L(R)}{e}_{k}\phantom{\rule{0.166667em}{0ex}}\frac{{\xi}_{\phi ,\phantom{\rule{0.166667em}{0ex}}r}(R)}{R}\end{array}$$(72)
in which ΔL_{φ, R}(R) is the Lagrangian surface luminosity perturbation due to mode φ, ξ_{φ, r} is the radial part of the adiabatic eigenfunction of that mode, i_{s} is the angle between the rotation axis and the line of sight at time t = 0 s (also called the spin inclination angle, see e.g., Fuller & Lai 2012), and t_{k} and e_{k} are limbdarkening coefficients given by the overlap integrals (in analogy to e.g., Burkart et al. 2012)
$$\begin{array}{cc}\hfill {t}_{k}& ={\displaystyle {\int}_{0}^{1}}\phantom{\rule{0.166667em}{0ex}}\mu \phantom{\rule{0.166667em}{0ex}}{H}_{r}(\mu \u037ek)\phantom{\rule{0.166667em}{0ex}}h(\mu )\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\mu \phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(73a)
$$\begin{array}{cc}\hfill {e}_{k}& ={\displaystyle {\int}_{0}^{1}}[2\phantom{\rule{0.166667em}{0ex}}{\mu}^{2}\phantom{\rule{0.166667em}{0ex}}\frac{\mathrm{d}{H}_{r}(\mu \u037ek)}{\mathrm{d}\mu}(\mu {\mu}^{3})\frac{{\mathrm{d}}^{2}{H}_{r}(\mu \u037ek)}{\mathrm{d}{\mu}^{2}}]\phantom{\rule{0.166667em}{0ex}}h(\mu )\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\mu \phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(73b)
where μ ≡ cosθ and h(μ) is a limbdarkening function. To derive Eq. (71), we assume that the Lagrangian flux perturbation of mode φ, which causes the luminosity perturbation, is equal to the radiative flux perturbation (see e.g., Unno et al. 1989). We use a linear limbdarkening function h(μ)=1 + (3 μ/2) in our computations (as has been customary for decades; see e.g., Osaki 1971 and Aerts et al. 1992, who set μ = 0.36 for B stars). One can however easily account for other, more sophisticated limbdarkening laws by changing the limbdarkening function. Numerical evaluations of t_{k} and e_{k} are necessary, because no analytic closed forms of the classical Hough functions exist.
Equation (71) thus separates contributions to the discaveraged luminosity fluctuation (ΔL_{φ} / L) into two multiplicative factors: a factor attributed to the properties of the mode φ under consideration, (ΔL_{φ} / L)_{mode}, and an angular factor that describes the observer’s orientation in the rotating frame. The theoretical flux fluctuation of that mode at the surface, 𝔏_{φ}, can then be computed by multiplying the (complex) amplitude c_{φ} obtained from solving the AEs with the factor (ΔL_{φ} / L). Its modulus 𝔏_{φ} can directly be compared with observed luminosity amplitudes $\stackrel{\sim}{A}$. The amplitude conversion factor o_{A, φ} for a mode φ is then defined as
$$\begin{array}{c}\hfill {o}_{A,\phantom{\rule{0.166667em}{0ex}}\phi}=\left{\left(\frac{\mathrm{\Delta}{L}_{\phi}}{L}\right)}_{\mathrm{mode}}\right\phantom{\rule{0.166667em}{0ex}}{H}_{r}({i}_{s}\u037ek)\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(74)
We compute the expected theoretical threshold surface luminosity fluctuations 𝔏_{t}, which is the minimum observed luminosity fluctuation that a parent mode in a mode triads needs for (parametric) resonant mode coupling to occur, and the stationary surface luminosity fluctuation $\mid {\mathfrak{L}}_{\phi}^{s}\mid $ of mode φ, as
$$\begin{array}{c}\hfill {\mathfrak{L}}_{t}={o}_{A,1}\phantom{\rule{0.166667em}{0ex}}{A}_{t}\phantom{\rule{0.166667em}{0ex}},{\mathfrak{L}}_{\phi}^{s}={o}_{A,\phantom{\rule{0.166667em}{0ex}}\phi}\phantom{\rule{0.166667em}{0ex}}{A}_{\phi}^{s}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(75)
The theoretical g mode stationary amplitudes ${A}_{\phi}^{s}$ and the theoretical g mode (parametric) threshold amplitudes A_{t} in Eq. (75) are computed by setting the bookkeeping parameter 𝔍 = 1 (similar to e.g., Buchler & Goupil 1984), so that δω = ΔΩ^{l}.
The conversion factors (74) are sensitive to the choice of a limbdarkening function (because this affects t_{k} and e_{k}), as well as the normalization factors for the Hough functions and the radial parts of the mode eigenfunctions defined in Eq. (11) and Sect. 2.3, respectively. The stationary daughterparent surface luminosity ratios $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $, $\mid {\mathfrak{L}}_{3}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ minimize the influence of the choice of limbdarkening function and normalization factors. We determine these ratios as
$$\begin{array}{c}\hfill \frac{{\mathfrak{L}}_{2}^{s}}{{\mathfrak{L}}_{1}^{s}}=\left(\frac{{o}_{A,\phantom{\rule{0.166667em}{0ex}}2}}{{o}_{A,\phantom{\rule{0.166667em}{0ex}}1}}\right)\sqrt{\frac{{Q}_{2}}{{Q}_{1}}}\phantom{\rule{0.166667em}{0ex}},\frac{{\mathfrak{L}}_{3}^{s}}{{\mathfrak{L}}_{1}^{s}}=\left(\frac{{o}_{A,\phantom{\rule{0.166667em}{0ex}}3}}{{o}_{A,\phantom{\rule{0.166667em}{0ex}}1}}\right)\sqrt{\frac{{Q}_{3}}{{Q}_{1}}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(76)
The daughterparent surface luminosity fluctuation ratios (76) are the most robust amplitudebased theoretically predicted observables that can be used in resonant nonlinear asteroseismic modeling. To derive the expressions for these ratios, we assume a parametric resonant mode triad (for a threemode sumfrequency coupling or its difference frequency analogue), and use the definition of the quality factor Q_{φ}, in addition to Eq. (52). Stationary surface luminosity ratios can thus be computed in terms of linear nonadiabatic parameters. The equivalent expression for the daughterparent surface luminosity fluctuation ratio of a harmonic dyad is given in Appendix B.2 of Van Beeck (2023).
In this work, we only envision a rough comparison between theoretical predictions and observables by limiting ourselves to monochromatic predictions. As highlighted by Aerts & Tkachenko (2023), future studies of measured amplitude ratios from multicolor space photometry by combining Gaia (Gaia Collaboration 2016), Kepler (Koch et al. 2010) or PLATO (Rauer et al. 2014) data offer additional opportunities to characterize stellar atmospheric properties. Concrete applications of our theory require integrations over particular passbands instead of the monochromatic predictions for the daughterparent ratios considered here. Such integrations will not be considered in this work but will be considered in followup application papers.
3.2. Frequencies and phases: frequency detuning and combination phase for Ω_{1} ≈ Ω_{2} + Ω_{3}
The inherent assumption made in linear asteroseismic inference is that any nonlinear frequency shifts (e.g., those determined by Eq. (54)) are negligible, so that theoretical frequencies computed within a linear formalism can directly be compared to their observed counterparts. A nonlinear formalism, such as the one we derive in Sect. 2, allows one to verify that assumption.
We define the observed frequency detuning $\mathrm{\Delta}\stackrel{\sim}{\mathrm{\Omega}}$ as
$$\begin{array}{c}\hfill \mathrm{\Delta}\stackrel{\sim}{\mathrm{\Omega}}={\stackrel{\sim}{\mathrm{\Omega}}}_{1,\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}{\stackrel{\sim}{\mathrm{\Omega}}}_{2,\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}{\stackrel{\sim}{\mathrm{\Omega}}}_{3,\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(77)
where the observed frequencies of modes φ are defined as ${\stackrel{\sim}{\mathrm{\Omega}}}_{\phi ,\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}$ with subscript 𝔦 indicating that observed frequencies are measured in the inertial frame. Because of the azimuthal selection rule (24), $\mathrm{\Delta}\stackrel{\sim}{\mathrm{\Omega}}$ is also equal to its corotating frame equivalent $\mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}_{\mathfrak{c}}$, that is, $\mathrm{\Delta}\stackrel{\sim}{\mathrm{\Omega}}=\mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}_{\mathfrak{c}}\equiv {\stackrel{\sim}{\mathrm{\Omega}}}_{1}{\stackrel{\sim}{\mathrm{\Omega}}}_{2}{\stackrel{\sim}{\mathrm{\Omega}}}_{3}$. This equivalence, along with Eq. (55) and the stationary equivalent of Eq. (77), then determine that
$$\begin{array}{c}\hfill \mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}^{s}=\mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}^{\mathit{nl}}=0.0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(78)
needs to be fulfilled for a isolated (resonantly locked) mode triad. In identifying such couplings observationally, one should therefore search for combinations of observed modes with stationary amplitudes, for which
$$\begin{array}{c}\hfill \mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}^{s}+{\sigma}_{\mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}^{s}}\lesssim 2\phantom{\rule{0.166667em}{0ex}}\pi \phantom{\rule{0.166667em}{0ex}}{\mathfrak{R}}_{\nu}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(79)
where ${\sigma}_{\mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}^{s}}$ is the propagated uncertainty of $\mathrm{\Delta}{\stackrel{\sim}{\mathrm{\Omega}}}^{s}$, and ${\mathfrak{R}}_{\nu}\equiv {\textstyle \frac{1}{T}}$ denotes the Rayleigh limit, with T being the total time span of the time series of the SPB star that needs to be modeled.
The observable that can directly be compared with observed frequencies of modes in an inferred candidate resonance is
$$\begin{array}{c}\hfill {\mathrm{\Omega}}_{\phi ,\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}^{\mathrm{NL},\mathrm{s}}\equiv {\mathrm{\Omega}}_{\phi}^{s}+{m}_{\phi}\phantom{\rule{0.166667em}{0ex}}\mathrm{\Omega}+\delta {\mathrm{\Omega}}_{\phi}^{s}={\mathrm{\Omega}}_{\phi ,\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}^{s}+\delta {\mathrm{\Omega}}_{\phi}^{s}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(80)
which includes the quadratic stationary nonlinear frequency shift (54). By summing 𝔏_{φ} = o_{A, φ} A_{φ} with its complex conjugate, we determine the individual stationary phase observables (see Sect. 4.3.2 in Van Beeck 2023 for the explicit manipulations)
$$\begin{array}{c}\hfill \left[\begin{array}{c}{\stackrel{\sim}{\varphi}}_{1,s}^{s}\\ {\stackrel{\sim}{\varphi}}_{2,s}^{s}\\ {\stackrel{\sim}{\varphi}}_{3,s}^{s}\\ \end{array}\right]=\frac{\pi}{2}\left[\begin{array}{c}{\left({\varphi}_{1}^{s}\right)}_{0}\\ {\left({\varphi}_{2}^{s}\right)}_{0}\\ {\left({\varphi}_{3}^{s}\right)}_{0}\\ \end{array}\right]\phantom{\rule{4pt}{0ex}}\left[\begin{array}{c}{\varphi}_{L\phantom{\rule{0.166667em}{0ex}}1}\\ {\varphi}_{L\phantom{\rule{0.166667em}{0ex}}2}\\ {\varphi}_{L\phantom{\rule{0.166667em}{0ex}}3}\\ \end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(81)
where we use the stationary equivalent of Eq. (39b), as well as Eqs. (34), (45), and (49). In Eq. (81), ϕ_{L φ} is the phase of the complex quantity (ΔL_{φ} / L)_{mode} defined in Eq. (72). The observed individual mode phases (81) therefore are independent of time if the modes are part of a locked mode triad.
In analogy with the definition of the stationary equivalent of the theoretical generic phase coordinate Υ in Eq. (41), we define the stationary combination phase observable ${\stackrel{\sim}{\mathrm{\Phi}}}_{s}^{s}$ as
$$\begin{array}{c}\hfill {\stackrel{\sim}{\mathrm{\Phi}}}_{s}^{s}={\stackrel{\sim}{\varphi}}_{1,s}^{s}{\stackrel{\sim}{\varphi}}_{2,s}^{s}{\stackrel{\sim}{\varphi}}_{3,s}^{s}=\frac{\pi}{2}{\left(\mathrm{\Delta}{\varphi}^{s}\right)}_{0}\mathrm{\Delta}{\varphi}_{L}^{s}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(82)
where the last equality holds because of Eq. (81), and in which ${\left(\mathrm{\Delta}{\varphi}^{s}\right)}_{0}={\left({\varphi}_{1}^{s}\right)}_{0}{\left({\varphi}_{2}^{s}\right)}_{0}{\left({\varphi}_{3}^{s}\right)}_{0}$ and $\mathrm{\Delta}{\varphi}_{L}^{s}={\varphi}_{L\phantom{\rule{thinmathspace}{0ex}}1}{\varphi}_{L\phantom{\rule{thinmathspace}{0ex}}2}{\varphi}_{L\phantom{\rule{thinmathspace}{0ex}}3}$. The coupling coefficient η_{1} is realvalued and therefore has no contribution to Eq. (82). Stationary combination phase observable (82) is to be compared with a relative phase computed from the phases of observed candidate resonance signals.
4. Numerical results for SPB models
We simulate stationary resonant parametric threegmode coupling processes by computing a grid of models representative for the SPB oscillators. Numerical stellar evolution models are generated by the stellar evolution code MESA (version 15140; Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019). Linear stellar oscillation models are generated by the stellar oscillation code GYRE (version 6.0.1; Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020), and use the MESA models as input. Numerical mode coupling models use both MESA and GYRE models as input.
We discuss the MESA model grid setup in Sect. 4.1, the GYRE model grid setup in Sect. 4.2, and display and discuss the numerical results for triad ensembles in Sect. 4.3. The link to our inlists for these codes, as well as the link to our mode coupling code repository, can be found in Appendix G.
4.1. MESA model grid setup for SPB stars
We compute MESA stellar evolution models with the parameters given in Table 1. These parameters cover the ranges of inferred initial mass (M_{ini}) and core hydrogen mass fraction (X_{c}) values for SPB stars given in Table 2 of Pedersen (2022). The MESA models have the ‘standard’ initial chemical mixture of nearby Btype stars derived by Nieva & Przybilla (2012) and Przybilla et al. (2013), and an Eddington gray atmosphere. Following Pedersen et al. (2021), we adjust the initial hydrogen and helium mass fractions X_{ini} and Y_{ini} so that the ratio $\frac{{X}_{\mathrm{ini}}}{{Y}_{\mathrm{ini}}}}={\textstyle \frac{{X}_{\ast}}{{Y}_{\ast}}$, with X_{*} and Y_{*} equal to the Galactic standard values for Btype stars in the solar neighborhood (Przybilla et al. 2013). We use Opacity Project (OP) opacity tables (Seaton 2005) that were computed by Moravveji et al. (2015) for this elemental mixture. The full protonproton chain and CNO cycle nuclear reaction networks are used to describe core hydrogen fusion on the main sequence. Beyond the zero age main sequence (ZAMS), we use the Vink et al. (2001) hot wind scheme with a wind scaling factor fixed to a value of 0.3 (see Björklund et al. 2021).
Model parameters of the MESA model grid.
Diffusive isotope mixing processes within the stellar interior are assumed to be described by the simplified transport Eq. (1) of Michielsen et al. (2021) and Pedersen et al. (2021). Mixing processes in the radiative envelope are described by a diffusive mixing profile for internal gravity wave mixing deduced by Rogers & McElwaine (2017) and used by Pedersen et al. (2018) in the context of asteroseismic modeling. This profile scales the radiative envelope mixing level D_{env} with a factor inversely proportional to the (local) mass density. To model coreboundary mixing (CBM) processes, we employ an approach similar to the diffusive exponential overshooting model with efficiency parameters f_{CBM} and f_{0} fixed to 0.02 and 0.005 (see e.g., Michielsen et al. 2019, 2021). The inner boundary of the overshooting zone (i.e., the convective core mass) is determined by the Ledoux criterion for mixing length parameter α_{MLT} = 2.0 within the Cox & Giuli (1968) formalism for mixing length theory. In the implementation, we set the minimal level of diffusive isotope mixing D_{env, min} equal to 100 cm^{2} s^{−1}. If the mixing level drops below that boundary at a certain location within the model, diffusive mixing is halted locally.
Our baseline (fiducial) model is a near terminal age main sequence (TAMS) model with X_{c} = 0.09, M_{ini} = 4 M_{⊙}, solar initial metallicity (Z_{ini} = 0.014), and f_{CBM} equal to 0.02. Models ΔX_{c, 1} and ΔX_{c, 2} have core hydrogen mass fractions X_{c} equal to 0.29 and 0.59, respectively, and are representative for a midMS and nearZAMS SPB star, because X_{c} is a proxy for the main sequence age of the star (hydrogen is depleted in the core during the main sequence). Model ΔX_{c}M_{ini} is representative of a M_{ini} = 6 M_{⊙} midMS SPB star. The fixed values of the parameters D_{env, min}, f_{CBM} and Z_{ini} are representative of the median values of these parameters in Table 2 of Pedersen (2022).
To compute the nonlinear quadratic mode coupling coefficients (22), we need the adiabatic derivative of Γ_{1} with respect to the stellar density, ${\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}$. We compute this quantity in the MESA models with a custom function defined in Appendix C.2.
4.2. GYRE model grid setup for SPB stars
The most commonly observed g modes in SPB stars are prograde (k, m)=(0, 1) and (k, m)=(0, 2) modes (see e.g., Pedersen et al. 2021; Szewczuk et al. 2021). For each of the computed MESA models, we therefore compute the GYRE linear adiabatic and nonadiabatic eigenfrequencies and eigenfunctions of these commonly detected g modes, within the Cowling approximation (Cowling 1941) and the TAR. We assume uniform rotation at 20% of the Roche critical rotation rate Ω_{Roche} for most models, except for the ΔX_{c, 2} model, for which we consider rotation at rates equal to 20% and (the excessively high) 60% of Ω_{Roche}. The positive (negative) linear driving (damping) rates of the modes are estimated as the imaginary part of the nonadiabatic linear eigenfrequencies computed by GYRE.
Radial order ranges {n_{u}} (u ∈ ⟦3⟧) of the (k, m)=(0, 2) parent and (k, m)=(0, 1) daughter modes computed with GYRE and based on the MESA models listed in Table 1, angular rotation rates Ω expressed in percentages of the Roche critical rotation rate Ω_{Roche} and in d^{−1}, model radius R, and the order of magnitude of the bounds for the ranges of ϑ_{2} and ϑ_{3} (denoted by ϑ_{2, 3}), q_{1} and the quality factors Q_{φ} (φ ∈ ⟦3⟧).
The radial orders of the g modes computed for the different MESA stellar evolution models are listed in Table 2. These are similar to those of Pedersen et al. (2021). We ensure that the linear nonadiabatic GYRE computations – which use trial frequencies based on the adiabatic GYRE computation results – do not show evidence of missing nonadiabatic mode radial orders (see e.g., Goldstein & Townsend 2020), so that there is a onetoone mapping between the adiabatic and nonadiabatic mode radial orders. The ranges for the quality factors in Table 2 show that the condition γ_{φ}/ω_{φ}≡1/Q_{φ}≪1 mentioned in Sect. 2.6 is fulfilled for all modes (i.e., all modes are on the slow manifold).
4.3. Search for isolated mode couplings: triad parameter ensembles for SPB stars
Following the stability analysis in Sect. 2.8, only parametric resonance scenarios are able to produce stable stationary solutions for resonances Ω_{1} ≈ Ω_{2} + Ω_{3} (or their difference frequency analogues). We therefore restrict the GYRE computations of potential parent modes to those (k, m)=(0, 2) modes with linear driving rates γ_{1} > 0, and restrict the computations of potential daughter modes to those (k, m)=(0, 1) modes with linear damping rates γ_{2} < 0 and γ_{3} < 0. Gaps in the daughter mode radial order ranges {n_{2, 3}} listed for some of the models in Table 2 indicate the presence of linearly excited (k, m)=(0, 1) modes that cannot partake in a parametric resonance. For the 6 M_{⊙} nearZAMS and 8 M_{⊙} midMS and nearZAMS models with a rotation rate of 20% of Ω_{Roche}, we obtain no linearly excited potential parent modes.
We compute mode coupling coefficients η_{1} for mode triads whose stationary AE solutions satisfy the (linear) stability conditions derived in Sect. 2.8. To guarantee nonlinear stability of the stationary solutions, these mode triads also need to fulfill the hyperbolicity check. This hyperbolicity check can only be done after performing the nonlinear mode coupling computations, because the computed stationary Jacobian matrix contains the coupling coefficient η_{1} (see Eq. (B.39) in Van Beeck 2023). The additional model validity constraints discussed in the next subsection are implemented to ensure that the solutions are away from the edge of the domains of applicability.
4.3.1. Model validity criteria
Buchler et al. (1997) state that resonant AEs are valid not only near a resonance center, but also far from resonance, with the solutions from resonant AEs slowly approaching the nonresonant AE solutions when moving away from that resonance center. This point however only holds under the assumption that all modes participating in the modal coupling are linearly unstable, so that a stable multimode fixed point is reached with the same modes for resonant and nonresonant AEs (Buchler et al. 1997). If that were not the case, the nonresonant solutions would predict negligibly small amplitudes for linearly damped modes. For the cubic couplings considered by Buchler et al. (1997), stable stationary points might exist for the driven resonance scenarios required for this type of behavior. However, for the quadratic nonlinear couplings considered in this work, stable multimode fixed points can only be found for parametric or direct resonant scenarios (with the latter only yielding stable stationary solutions for harmonic resonances; see Sect. 2.8). Linearly damped modes participate in such coupling scenarios, invalidating the Buchler et al. (1997) assumption.
The resonant AEs derived in Sect. 2.6 are therefore valid only near the resonance center, that is, if the absolute value of the linear triad frequency detuning ΔΩ^{l} defined in Eq. (34), is small in comparison to the absolute values of the individual corotatingframe mode frequencies: ΔΩ^{l}≪Ω_{φ} ∀φ ∈ ⟦3⟧. If it is not small, there are no comparatively slow amplitude variations and additional mode interactions need to be considered to describe the dynamics. We therefore estimate the validity domain of the AEs in terms of the parameter:
$$\begin{array}{c}\hfill {\mathrm{\Psi}}_{\mathrm{AE}}\equiv \frac{\left\mathrm{\Delta}{\mathrm{\Omega}}^{l}\right}{min({\mathrm{\Omega}}_{1},\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{2},\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{3})}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(83)
which compares ΔΩ^{l} with the minimum absolute value of the corotatingframe angular mode frequencies of the triad. We require that mode triads considered for modeling have Ψ_{AE} ≤ 0.1, so that ΔΩ^{l} ≪ min(Ω_{1}, Ω_{2}, Ω_{3}). Many mode triads satisfy this model validity criterion, as is shown in Table 3.
Overview of the number of mode triads that satisfy different model validity and stability criteria.
It is furthermore important to consider the threshold surface luminosity fluctuations 𝔏_{t} of all possible mode triads, because these determine the onset of mode interaction. To obtain isolated mode couplings, we determine the lowest and secondlowest 𝔏_{t} for each of the possible mode triads in which a given parent mode participates. Doing so allows us to label the individual modes of those triads with the computed value of 𝔏_{t}, which can be compared for each mode individually. Identified isolated mode couplings must be (1) triads for which the lowest 𝔏_{t} labels of all three modes (across all mode triads in which each mode is involved) are the same; and (2) triads for which $\mid {\mathfrak{L}}_{1}^{s}\mid $ of the lowest𝔏_{t} mode triad is smaller than the secondlowest 𝔏_{t} associated with the parent mode. If at least one of those conditions is not fulfilled for a specific mode triad, the amplitude of at least one of the members of that mode triad is set by parametric energy transfer in a different mode triad. The AEs derived in Sect. 2.6 would then be invalid to describe the multimode coupling that occurs among such mode triads. We collectively refer to these two conditions as the isolation criterion, and call mode triads that satisfy this criterion isolated mode triads.
Fig. 3. Empirical probability distributions of δω / (γ_{2} + γ_{3}) for the 21 identified (stable and AE valid) isolated mode triads (blue and dotted) and their AE valid and stable counterparts (yellow) that do not fulfill the isolation criterion. We choose the bin width in a similar way to what is done in Hogg (2008). 
Arras et al. (2003) have a different implicit ‘isolation criterion’: they stated that if the resonance is sharp, it is plausible that only the parent and two daughter modes are relevant for the mode amplitude dynamics. For sharp resonances, one expects δω / (γ_{2}+γ_{3}) to be small. Figure 3 shows the empirical probability distribution of δω / (γ_{2}+γ_{3}) for the 21 identified (stable, AEvalid) isolated mode triads, and their AEvalid and stable counterparts that do not fulfill the isolation criterion. While it is certainly true that δω / (γ_{2}+γ_{3}) is small for the isolated mode triads, there are nonisolated mode triads that have similarly low values of δω / (γ_{2}+γ_{3}), which reiterates the important role threshold surface luminosity fluctuations play in identifying such isolated mode triads.
Table 3 gives an overview of how restrictive the individual criteria are, and shows the effect of solely enforcing the validity or stability criteria. Many stable mode triads are identified, and fixed point stability is essentially defined by the hyperbolicity check. Both validity criteria are more restrictive than the stability criteria. It does not seem obvious for a g mode to participate in an isolated mode coupling scenario, as is reflected by the small number of identified isolated mode triads in Table 3. Moreover, strict enforcement of the TAR validity frequency hierarchies listed in for example Dhouib et al. (2021a) would further reduce the number of mode triads considered for asteroseismic modeling. This would however restrict the quadratic mode couplings identified for the 4 M_{⊙} nearZAMS models from being used, contradicting observational evidence (e.g., V21). We therefore do not strictly enforce these frequency hierarchies.
4.3.2. Radial coupling contributions
For isolated mode triads that fulfill the stability and validity conditions, we compute a radial profile of the coupling coefficient η_{1}(r) by integrating the overlap integrals in Eq. (22) up to an internal radial coordinate r instead of the stellar model surface radius R (i.e., r ≤ R and η_{1} ≡ η_{1}(R)). The profile η_{1}(r) then indicates the zones in the stellar interior for which contributions to the coupling coefficient are significant for these mode triads.
In the nearcore regions the squared BruntVäisälä frequency (N^{2}) profile has a (local) maximum due to the presence of a μgradient left behind by the receding convective core during the main sequence. This maximum introduces a sharp transition in N^{2} which changes the eigenfunctions: their modal inertia becomes more confined to the nearcore region (Miglio et al. 2008; Michielsen et al. 2021), at the location of the (local) maximum. This can be rationalized by the increasing confinement of modal inertia to the TARvalid nearcore mixing zone with model age. The width of that maximum depends on the evolutionary stage: more evolved SPB models have wider N^{2} maxima, and some modes become trapped in the nearcore regions (Moravveji et al. 2016; Michielsen et al. 2021). For these more evolved models, the instability strip for (k, m)=(0, 2) parent g modes also moves towards higher radial orders, further confining modal inertia to the oscillatory nearcore regions. Both these effects likely explain why the contributions to coupling coefficients are more concentrated in the nearcore regions for more evolved models.
Examples of such radial coupling coefficient profiles are given in Figs. 4 and 5, which display the N^{2} profile along with the coupling coefficient profile η_{1}(r). The profile in Fig. 4 describes the coupling contributions for the mode triad characterized by g mode radial orders (n_{1}, n_{2}, n_{3}) = (−22, −15, −53) in the midMS ΔX_{c, 1} model in Table 4. The nearZAMS example displayed in Fig. 5 describes the coupling contributions for the mode triad characterized by g mode radial orders (n_{1}, n_{2}, n_{3}) = (−14, −10, −26) in the ΔX_{c, 2}_{(a)} model of Table 4. These figures confirm that contributions to the coupling coefficient of the mode coupling in the midMS model are typically more confined to the nearcore region in comparison to those contributions in the nearZAMS model.
Fig. 4. Coupling coefficient and squared BruntVäisälä frequency (N^{2}) profiles as a function of the fractional radius for the isolated mode triad with g mode radial orders (n_{1}, n_{2}, n_{3}) = (−22, −15, −53) in the midMS ΔX_{c, 1} model of Table 4. 
Fig. 5. Same as Fig. 4, but for the isolated mode triad with g mode radial orders (n_{1}, n_{2}, n_{3}) = (−14, −26, −10) in the ΔX_{c, 2}_{(a)} model of Table 4. 
Linear corotatingframe frequencies Ω_{φ}, largest absolute spin parameters s_{φ} (of the triad modes), radial orders n_{φ} and linear driving or linear damping rates γ_{φ} of g modes φ in all identified isolated (g) mode triads with ordering numbers (k_{1}, k_{2}, k_{3}) = (0, 0, 0) and azimuthal orders (m_{1}, m_{2}, m_{3}) = (2, 1, 1) for the models in Table 2 that yield stable and valid stationary solutions.
4.3.3. Properties of isolated g mode triads in SPB stars
Isolated mode triads can be identified for each model listed in Table 2. Their characteristic mode and triad properties are listed in Table 4, and their dimensionless model validity and stability estimators, including the order of magnitude of the coupling coefficient η_{1}, are displayed in Table 5. It is clear from the largest absolute values of the spin parameters that at least one of the modes in the identified isolated mode triads is subinertial (i.e., s_{φ}> 1), requiring the nonperturbative (TAR) description of the Coriolis force.
Dimensionless quantities for the mode triads listed in Table 4: dimensionless AE validity parameter Ψ_{AE}, absolute detuningdamping ratio q and absolute detuningdriving ratio q_{1}, the drivingdamping rate ratios ϑ_{2} and ϑ_{3}, as well as the order of magnitude of the energynormalized coupling coefficient, 𝒪_{η1}.
The Doppler shift causes a logical trend to appear when correlating the average values of Ω_{φ} with a change in the rotation frequency (i.e., upon comparison of mode frequencies for the models ΔX_{c, 2}_{(a)} and ΔX_{c, 2}_{(b)}): Ω_{φ} decreases for increasing rotation rates. This trend is only weakly influenced by the effect of the Coriolis force for these sectoral prograde g or equatorial Kelvin waves, which is not the case for nonsectoral modes (due to geostrophic balance; see e.g., Gill 1982; Townsend 2005; DaszyńskaDaszkiewicz et al. 2015; Szewczuk & DaszyńskaDaszkiewicz 2017).
The corotatingframe mode frequencies also (on average) decrease for highermass models, and decrease with model age. This change with mass can be attributed to the more massive SPB models having a larger model radius, which enlarges the pulsation cavity, although it is to a degree counteracted by the decrease in Ω.
Mode frequency in the corotating frame is instrumental for determining linear, heatdriven mode excitation: it determines the geometrically optimal driving region for a mode, which needs to coincide with the Z bump for the mode to be linearly excited in an SPB star. In general, higher radial order g modes^{3} – which have lower frequencies and are therefore more susceptible to rotational effects – are excited in more evolved models (on the main sequence). The expected higher mode densities for more evolved stars might explain why only one isolated mode triad was identified for the fiducial model. The linear (vibrational) stability sensitively depends on the opacity profile of any particular stellar model (e.g., DaszyńskaDaszkiewicz et al. 2017). A change in opacity profile can therefore lead to the identification of different isolated mode triads.
Stabilization of the (k, m)=(0, 1) g modes in more rapidly rotating stars (e.g., Townsend 2005) can lead to opportunities for nonlinear excitation by a parametric resonance, as it enlarges the pool of potential daughter modes with which combination can be made. A physical explanation of the effect of rotation within the TAR on mode excitation by the κ mechanism can be found in Ushomirsky & Bildsten (1998) and Townsend (2005). That increasing potential for mode combination with larger rotation rate can be observed in the ‘Total’ column of Table 3, where more combinations are considered for the more rapidly rotating model ΔX_{c, 2}_{(b)}, if compared to the number of combinations considered for model ΔX_{c, 2}_{(a)}. The stabilizing effect for the parent (k, m)=(0, 2) g modes comparatively seems weaker (as was also observed by Townsend 2005): only one of the parent modes is stabilized, in comparison to seven of the daughter modes, when we compare the ΔX_{c, 2}_{(a)} and ΔX_{c, 2}_{(b)} models. In this case, we identify four isolated triads for model ΔX_{c, 2}_{(b)}, and two isolated triads for the slowerrotating model ΔX_{c, 2}_{(a)}. However, the stabilizing effect does not necessarily mean that we can identify more isolated mode triads for rapid rotators, because those triads need to fulfill the stability and validity criteria.
The range of q_{1} for the identified isolated mode triads is [10^{4}, 10^{0}] (using the notation of Table 2). This indicates that the two rightmost panels on the top of Fig. 2, as well as the panels on the bottom row of Fig. 2, describe the linear stability of the isolated stationary solutions. These solutions’ resonant nature is important, as is indicated by the values of q being of the order of unity to ten, affecting the stationary amplitude (48).
We can identify only one isolated harmonic resonant triad (for the ΔX_{c, 1} model), which is a small number in comparison to the 20 identified resonances of the combination type Ω_{1} ≈ Ω_{2} + Ω_{3}. It therefore seems comparatively harder to satisfy all stability and validity conditions for harmonic resonances, which provides a reason for why only a few of such resonances are identified in the observational data of V21.
We only consider combinations of a small number of modes in this proofofconcept study and do not study the dependence of nonlinear parameters – such as the coupling coefficient – on stellar parameters in detail.
5. Impact on asteroseismic modeling
In this section, we discuss how the nonlinear theoretically predicted observables defined in Sect. 3 can aid or complement current frequencybased asteroseismic modeling of SPB stars. The various nonlinear observables of interest, or derived quantities thereof, are listed in Table 6.
Maximal absolute nonlinear frequency shifts as a fraction of the linear inertialframe angular frequency $\mid \delta {\mathrm{\Omega}}_{\phi ,\mathrm{m}}^{s}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\Omega}}_{\phi ,\phantom{\rule{thinmathspace}{0ex}}i}\mid $ (in partsperthousand), zeropointcorrected nonlinear combination phases ${\stackrel{\sim}{\mathrm{\Phi}}}_{s,0}^{s}$, and expected stationary surface luminosity fluctuation ratios $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ and $\mid {\mathfrak{L}}_{3}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ for the mode triads listed in Table 4.
5.1. Oscillation frequency and combination phase observables
We find that the nonlinear frequency shifts of the isolated mode triads are small in comparison to the (linear) inertial mode frequencies, as is illustrated by the values in the column displaying values of $\mid \delta {\mathrm{\Omega}}_{\phi ,\mathrm{m}}^{s}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\Omega}}_{\phi ,\phantom{\rule{thinmathspace}{0ex}}i}\mid $ in Table 6, which shows the largest value of the frequency shift for any of the three modes in the isolated triad, expressed in units of the inertial mode frequency. These dimensionless frequency shifts are typically of order 10^{−3} to 10^{−4}, and we find no clear correlations with any of the model parameters (potentially due to the small number of modes considered).
Translated to units of d^{−1}, this means that there are several frequencies whose (nonlinear) frequency shifts are of order 10^{−5} d^{−1} or smaller, rendering them of similar order as the typical errors for the frequencies derived from 4year Kepler light curves (e.g., Bowman & Michielsen 2021). The largest shifts are obtained for the modes in the fasterrotating models, and are of order 10^{−2} or 10^{−3} d^{−1}, which is approximately two orders of magnitude larger than typical uncertainties. In practice, the resonant nonlinear frequency shifts are hard to distill from frequency lists deduced from prewhitening analyses similar to what was done in V21. The nonresonant frequency shifts are even smaller, justifying the approximation of using linear frequencies in asteroseismic modeling.
We show the zeropointcorrected nonlinear combination phases ${\stackrel{\sim}{\mathrm{\Phi}}}_{s,0}^{s}$ for all identified isolated mode triads in Table 6, where ${\stackrel{\sim}{\mathrm{\Phi}}}_{s,0}^{s}={\stackrel{\sim}{\mathrm{\Phi}}}_{s}^{s}+{\left(\mathrm{\Delta}{\varphi}^{s}\right)}_{0}$. We find no specific trends for these combination phases when correlated with other model parameters (potentially due to the small number of modes considered). In principle, such phases can be compared with observed relative phases, if the relative phase at the zero point is accounted for. The latter is however dependent on the initial values of the individual mode phases, ϕ_{0}, which are unconstrained.
5.2. Amplitude ratio observables
The stationary surface luminosity fluctuation ratios $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ and $\mid {\mathfrak{L}}_{3}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ offer the best constraints for asteroseismic modeling. We find that most of these ratios are smaller than unity, indicating that most of the stationarystate (or saturation) energies of the linearly driven (k, m)=(0, 2) parent modes are larger than those of the (k, m)=(0, 1) daughter modes that are linearly damped and parametrically excited. Daughter modes have larger stationary mode amplitude ratios $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ only for certain mode triads in the 4 M_{⊙} nearZAMS models ΔX_{c, 2}_{(a)} and ΔX_{c, 2}_{(b)}, as well as the 6 M_{⊙} nearTAMS models ΔM_{ini, 1}. For these specific mode triads, the other daughter mode’s amplitude is significantly smaller than that of the parent mode.
The stationary amplitude ratios (76) depend on the ratios of the mode quality factors that express the number of periods necessary to linearly increase or decrease modal energy by a factor of e. Such ratios are inversely proportional to ϑ_{2} or ϑ_{3} and proportional to the frequency ratio Ω_{2} / Ω_{1} or Ω_{3} / Ω_{1}. We find that the frequency ratios for these specific isolated mode triads that have largeramplitude daughter modes are not particularly different from those in other identified isolated mode triads, as can be rationalized from the corotating mode frequency values in Table 4. The observed difference in amplitude ratio is thus is mostly due to different values of ϑ_{2}. The smaller ϑ_{2}, the longer the linear damping timescale of daughter mode 2 is, in comparison to the linear driving timescale of parent mode 1. Based on the listed values of ϑ_{2} in Table 5, this effect seems to become important when ϑ_{2} ≲ 0.13 ≡ ϑ_{2, b}, that is, when the linear damping rate of one of the daughter modes is less than 13% of the linear driving rate of the parent mode. We also estimate the value of ϑ_{2} at which this effect becomes important by performing a linear regression between the ϑ_{2} and $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ values of four isolated mode triads in Table 5. Specifically, we construct a regression model using the ϑ_{2} and $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ values obtained from the two mode triads with parentdaughter amplitude ratios greater than unity but approaching the unity limit the closest, and the values of these quantities obtained from those mode triads with parentdaughter amplitude ratios smaller than unity that approach the unity limit the closest. The linear regression yields a linear slope of −0.08 ± 0.04 (ϑ_{2} per unit of $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $) and an intercept 0.22 ± 0.06, when we assume normally distributed residuals. This leads to an estimated boundary ${\widehat{\vartheta}}_{2,\mathrm{b}}=0.14\pm 0.07$ that is consistent with the earlier defined boundary ϑ_{2, b}.
In general, energy transfer occurs from the parent to the daughter mode only if the latter’s energy is smaller than the former’s (e.g., Arras et al. 2003). However, if one of the daughter modes participating in the mode triad is strongly damped (linearly), mode β, when the other daughter mode, mode α, is only weakly damped (linearly), we find that daughter mode α can saturate at larger energies than the parent mode. In such situations the parent mode needs to transfer a lot of energy to the daughter modes per cycle in order to overcome the strong linear damping of daughter mode β and reach a stationary state.
Isolated parametric couplings in which both daughter modes attain amplitudes larger than the parent mode are not observed in our calculations. Of the 276 835 potential triads considered in this work, only 32 nonAE valid triads have both daughterparent luminosity fluctuation energy ratios (76) greater than unity, which all originate from the ΔM_{ini, 1} model. Of those 32 solutions, 26 satisfy q > Q_{1}, the necessary but not sufficient condition for both stationary daughterparent rotatingframe mode energy ratios (52) to be greater than unity, because of Eqs. (49) and (78). The other 6 solutions are relatively close to satisfying that criterion (q and Q_{1} are of the same order of magnitude) and likely obtain daughterparent mode luminosity fluctuation energy ratios (76) greater than unity because of the differing ratios of amplitude conversion factors o_{A, φ}. Given the typical values of Ω_{φ}, γ_{φ} and q listed in Tables 4 and 5, as well as the quality factor ranges listed in Table 2, we deem it unlikely to encounter isolated solutions of this kind.
5.3. Comparison with observations of SPB stars
To assess the impact of mode amplitude ratios on period spacing pattern modeling, we compare daughterparent surface luminosity fluctuation ratios predicted by our formalism with observed equivalents derived for the ensemble of 38 SPB stars of V21. We limit ourselves to comparison among theoretically predicted amplitude ratios and observed equivalents of sum frequencies, because these make up the majority of the isolated mode triads listed in Table 5. Our comparison is also valid for difference frequencies, because the AEs and stationary solutions for these difference frequencies are of the same form as those of sum frequencies (see Appendix F).
Fig. 6. Observed daughterparent amplitude ratios ${A}_{d}^{\mathrm{\prime}}/{A}_{p}^{\mathrm{\prime}}$ of candidate resonances of SPB stars detected in models of light curves measured by Kepler, compared to their theoretically predicted equivalent observables, the monochromatic stationary daughterparent surface luminosity perturbation ratios 𝔏_{d} / 𝔏_{p}, which are computed with the theoretical oscillation model described in Sect. 2. These quantities are shown as a function of the minimal inertialframe frequency ν_{min} of the three signals in the considered candidate resonance, or, the minimal mode frequency observable ${\mathrm{\Omega}}_{\mathrm{min},\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}^{\mathrm{NL},\phantom{\rule{0.166667em}{0ex}}s}\phantom{\rule{0.166667em}{0ex}}/\phantom{\rule{0.166667em}{0ex}}2\pi $ among the three mode frequency observables of the considered triad, respectively. Left panel: ${A}_{d}^{\mathrm{\prime}}/{A}_{p}^{\mathrm{\prime}}$ as a function of ν_{min}. The light curve models were generated using strategy 3 of V21, which uses the signaltonoise ratio of the detected signals to determine their significance. Colors indicate the pseudoclass defined in V21: highf_{sv} (‘hfsv’), highmodedensity (‘hmd’), or outbursting (‘outb’). Right panel: 𝔏_{d} / 𝔏_{p} as a function of ${\mathrm{\Omega}}_{\mathrm{min},\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}^{\mathrm{NL},\phantom{\rule{0.166667em}{0ex}}s}/2\pi $. Colors indicate whether all stability criteria are fulfilled (‘Stab.’), or if additionally the AE validity criterion (‘Stab. + AE’) or all validity criteria (‘Stab. + Val.’; i.e., the isolated triads listed in Table 6) are fulfilled. We show alias signals in the frequency range of 0 to 24.4512 d^{−1} as fainter orangered (‘Stab.’) symbols, because we do not account for the splitting of the alias signals (which will affect the amplitude ratios). 
V21 generated variability models for the light curves, hereafter referred to as (harmonic) light curve models, of 38 SPB stars using five different harmonic analysis strategies. In this work, we use the models generated by their strategy 3, which uses the signaltonoise ratios to determine the significance of detected variability signals, because this resembles commonly used strategies in literature (see Table 1 and the corresponding discussion in Sect. 2.2 of V21 for additional details). They assigned pseudoclasses to the different members of the ensemble of analyzed SPB stars to denote the difficulties encountered during the analysis process: highf_{sv} stars had light curves that were described adequately by the harmonic light curve models; high mode density stars had many closespaced frequencies in their Fourier transforms or LombScargle periodograms; and outbursting stars were found to have distinct variability ‘outbursts’ in their light curves, resulting in many closespaced and highamplitude frequency groups (see V21). To select the relevant observed signals for our comparison, we look for sumfrequency resonance combinations of three signals that are part of the light curve models and fulfill the resonance criterion ν_{1} − ν_{2} − ν_{3}≤ℜ_{ν} + σ_{ν, prop} (with σ_{ν, prop} the propagated uncertainty of the observed frequency difference ν_{1} − ν_{2} − ν_{3} and ν_{1}, ν_{2} and ν_{3} the observed frequencies). These combinations also need to contain (i) at least one of the two highestamplitude signals in the light curve model, and (ii) loweramplitude components (i.e., not one of the two highestamplitude signals) that are not present in other considered combinations. The lowoverhighest frequency amplitude ratios of these combinations, ${A}_{d}^{\mathrm{\prime}}/{A}_{p}^{\mathrm{\prime}}$ (i.e., two ratios per combination), are displayed in the left panel of Fig. 6 as a function of the minimal frequency ν_{min} of the three modes in the combination, as measured within the inertial reference frame. Under the assumption that the nonlinear parametric resonant coupling process dominates the energy transfer among the considered modes, we hereafter refer to ${A}_{d}^{\mathrm{\prime}}/{A}_{p}^{\mathrm{\prime}}$ as the observed daughterparent amplitude ratios; under that assumption, ν_{min} is the minimal daughter mode frequency. Because the stationary combination phase observable ${\stackrel{\sim}{\mathrm{\Phi}}}_{s}^{s}$ defined in Eq. (82) depends on the unconstrained initial zero points of the individual mode phases and because no trends are observed in the zeropointcorrected nonlinear combination phases ${\stackrel{\sim}{\mathrm{\Phi}}}_{s,0}^{s}$ of isolated mode triads listed in Table 6, we do not enforce conditions on the relative phase of these observed combinations of signals, unlike V21.
One compares the observed daughterparent amplitude ratios with the relevant computed daughterparent (stationary) surface luminosity fluctuation ratios 𝔏_{d} / 𝔏_{p} of triads found among the (k, m)=(0, 1) daughter g modes and of the (k, m)=(0, 2) parent g modes computed by GYRE for the models listed in Table 2 (shown in the right panel of Fig. 6). These mode triads have small frequency detunings (34) and satisfy all stability criteria, and we further distinguish between triads with modes that satisfy (i) only the stability criteria, or (ii) the stability criteria and the AE validity criterion, or (iii) all stability and validity criteria (i.e., these are the isolated mode triads). The minimal mode frequency observable ${\mathrm{\Omega}}_{\mathrm{min},\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}^{\mathrm{NL},\phantom{\rule{0.166667em}{0ex}}s}$ (used as the xaxis variable in the right panel of Fig. 6) is smaller than the nonlinear frequency shift (54) for some of the considered modes in mode triads with linearly stable solutions (i.e., triads with modes that only satisfy the stability criteria and none of the validity criteria). In that case, we compute their alias frequencies in the frequency range [24.4512, 0] d^{−1} (using the notation in Table 2; i.e., up to the notional Kepler longcadence Nyquist frequency; Chaplin et al. 2014). Such aliased signals would however have their amplitude ratios affected by their splitting, which we do not account for in this work. The predicted surface luminosity fluctuation ratios of these alias frequencies are therefore shown on the right panel of Fig. 6 as faint (orangered) symbols.
The ensemble of stationary surface luminosity fluctuation ratios of computed triads seems to be able to explain the lowamplituderatio part of the ensemble of observed daughterparent amplitude ratios, as can be derived from the panels in Fig. 6. Some SPB modeling targets, for example KIC008714886, show promise for matching a part of their observed amplitude ratios with their theoretically predicted equivalent. Other targets (notably the slow rotator KIC0010526994) have observed daughterparent amplitude ratios that are all much greater than unity, which we do not have in the calculated ratios.
The absence of high theoretically predicted surface luminosity fluctuation ratios (52) might be caused by our choice of opacity tables, because they determine the linear growth or linear damping rates of the modes. For example, one of the analyzed SPB stars, KIC0010536147, is on the massive end of the SPB instability region (according to the stellar parameters determined in Pedersen 2022), where, according to our simulations that account for the classical κ mechanism driving mechanism, no (k, m)=(0, 2) parent modes are excited (see Table 2). The required changes in opacity to reach the highest observed amplitude ratios are, however, unlikely to be realistic. The modified opacities of DaszyńskaDaszkiewicz et al. (2017) used in modeling the linear stability of the hybrid mainsequence early Btype pulsator ν Eridani for example lead to normalized driving rates that were only two to four times larger than the unmodified driving rates. Asteroseismic modeling of dedicated targets is needed to further assess the effect of opacity on theoretical surface luminosity fluctuation ratios.
Additionally, some of the observed modes might have different nonradial geometries (i.e., different values of m, k) than the ones considered in the limited number of simulations carried out for this work. For example, the period spacing patterns considered by Pedersen et al. (2021) are made up of zonal (k, m)=(1, 0) dipole g modes for 9 of their considered 26 SPB stars, and one star, KIC008714886, has an identified retrograde (k, m)=(0, −1) dipole g mode period spacing pattern (see Supplementary Table 1 of Pedersen et al. 2021). When couplings between g modes other than the ones considered in this work have linear driving or linear damping rates, frequencies, as well as amplitude conversion factors o_{A, φ} of the same order as those obtained for the g modes considered in this work, we expect to obtain similar theoretically predicted luminosity fluctuation ratios. We note, however, that for such different combinations the parent mode can be a lower frequency mode. An example of such a combination that satisfies the coupling coefficient selection rules is the triad that consists of a driven (k, m)=(0, 1) parent mode α, a damped (k, m)=(0, 2) daughter mode β, and a damped (k, m)=(0, −1) daughter mode γ. The damped (k, m)=(0, 2) daughter mode β can have a lower frequency than the parent mode α in the corotating frame, while having a higher frequency in the inertial frame. Such a damped (k, m)=(0, 2) mode β would therefore be labeled as the parent when following the rules used to generate the observed daughterparent amplitude ratios shown in Fig. 6. For comparison with the theoretically predicted daughterparent surface luminosity fluctuation ratios, the relevant observed amplitude ratios then need to either be computed (as the ratio of the current amplitude ratios) or inverted. In our example, we need to compute the amplitude ratio of modes β and γ ${A}_{\gamma}^{\mathrm{\prime}}/{A}_{\beta}^{\mathrm{\prime}}=({A}_{\gamma}^{\mathrm{\prime}}/{A}_{\alpha}^{\mathrm{\prime}})/({A}_{\beta}^{\mathrm{\prime}}/{A}_{\alpha}^{\mathrm{\prime}})$), and we need to invert the currently computed amplitude ratio of modes α and β.
We also deem it likely that some of the high observed amplitude ratios can be attributed to other amplitude saturation mechanisms, such as multimode coupling or higherorder couplings. Correlation of the mean values of the 10 largest mode amplitudes (determined from the light curve models) with the mean values of the 10 largest identified amplitude ratios for each of the different SPB stars considered by V21 for example yields a positive Spearman correlation coefficient. This suggests that higherorder couplings might be needed to explain the high amplitude ratios of highamplitude SPB pulsators.
6. Conclusions
We derive a theoretical oscillation modeling framework that describes nonlinear threemode coupling among gravitoinertial (g) modes of Slowly Pulsating B (SPB) stars within the Traditional Approximation for Rotation (TAR; e.g., LonguetHiggins 1968; Lee & Saio 1997; Townsend 2003; Mathis 2013), extending and correcting terms in the formalism of L12 (see Appendix C). This framework relates the g mode adiabatic eigenfunctions to potential nonlinear energy exchange between the modes, by computing threemode (energyscaled) coupling coefficients η_{1} based on a phase space mode decomposition inserted in the relevant coupled equations of motion. To describe threemode resonant couplings, we derive amplitude equations (AEs) from the coupled equations of motion for a sumfrequency resonance Ω_{1} ≃ Ω_{2} + Ω_{3}, using the multiple time scales perturbation method (e.g., Nayfeh 1973, 1981; Nayfeh & Mook 1979). These coupling coefficients need to satisfy angular selection rules (24) and (25) if energy is to be exchanged among the modes. The isolated stable stationary solutions of these AEs (i.e., their stable fixed points) then describe locked threemode resonant couplings among g modes in rapidly rotating gmode pulsating stars, such as SPB stars.
We use this framework to compute examples of isolated mode couplings in stellar structure and pulsation models that represent SPB stars analyzed in V21. We limit ourselves to computing couplings among g modes with ordering numbers (k_{1}, k_{2}, k_{3}) = (0, 0, 0) and azimuthal orders (m_{1}, m_{2}, m_{3}) = (2, 1, 1); the most frequently observed modes in SPB stars (see e.g., Pedersen et al. 2021; Szewczuk et al. 2021; Pedersen 2022). The locked mode solutions have to fulfill multiple stability and validity criteria to describe physical threemode coupling scenarios in such stars. To ensure stability of the threemode stationary solution, it needs to be hyperbolic, the resonance must be parametric, and the solution must fulfill the condition in Eq. (61a) and the linear quartic stability criterion (64). The validity of the solution is guaranteed if it fulfills the AE validity condition Ψ_{AE} ≤ 0.1 (with Ψ_{AE} defined in Eq. (83)), and the isolation criterion. The most restricting of all these conditions is the isolation criterion, which ensures that no multimode coupling scenarios take place. We find that the restriction to sharp resonances as done by Arras et al. (2003) overestimates the number of isolated mode triads.
By performing coupling computations up to a certain inner radius, we map the regions that contribute significantly to the mode coupling. Typically, we obtain strong contributions to the resonant coupling among the three g modes in the nearcore zones of the SPB models, where N^{2} peaks and modes can become (partially) trapped (Miglio et al. 2008; Michielsen et al. 2021). For more evolved models, this N^{2} maximum is wider and the κ mechanism typically excites higher radial order modes, which are increasingly confined to the nearcore zones. Conversely, for less evolved models, significant contributions to η_{1} can be found outside of the nearcore zones.
Linear heatdriven vibrational instability defines which modes are available for the nonlinear (parametric) couplings we study in this work. The individual mode frequencies depend on a variety of factors, including rotation rate, evolutionary stage, and mass, among others. These frequencies are important, because they define the optimal linear driving regions of the corresponding modes. The linear driving of g modes in SPB stars, however, is ultimately reliant upon the opacity profiles inside the Fe bump. Opacity tables used during stellar modeling define the opacity gradients in the driving zones, and therefore influence (i) the number of possible mode triads that can be formed for a specific model, but also (ii) the linear driving and linear damping rates themselves. These linear driving and linear damping rates affect the nonlinear observables, such as the computed stationary mode luminosity fluctuations and their ratios.
For the few models considered in this work, we find no obvious correlations for changes in the coupling coefficients (22) with changes in model parameters. It is not clear whether such trends can be detected if the number of models and the number of considered modes is increased while varying model parameters such as the opacity tables. Because the absolute values of the detuningdamping ratio, q, are larger than unity, it is clear however that the resonant nature of the coupling strongly determines the stationary state mode amplitudes (48).
The nonlinear frequency shifts (54) for isolated triads are small enough to neglect when compared with the typical errors for the frequencies derived from 4year Kepler light curves, and the combination phases (82) contain an unconstrained initial mode phase parameter. The best constraints for asteroseismic modeling are therefore obtained from the predicted luminosity fluctuation ratios (76) of coupled modes. Most of these ratios indicate that the (k, m) = (0, 2) parent g modes considered in this work have higher apparent mode amplitudes than their coupled (k, m) = (0, 1) daughter g modes. We obtain disparate linear daughter mode damping rates for each of the isolated mode triads. If one of the linearly damped daughter modes in the isolated mode triad is weakly damped (linearly; at a rate ≲13% that of the parent linear driving rate), we find that that daughter mode can saturate at a higher amplitude than its coupled parent mode. We deem it unlikely that isolated solutions exist in which both daughterparent amplitude ratios are greater than unity, because of the typical values for the quality factors and detuningdamping ratios associated with g modes in SPB star models.
The monochromatic stationary luminosity fluctuation ratios are consistent with some of the lowest amplitude ratios ${A}_{d}^{\mathrm{\prime}}/{A}_{p}^{\mathrm{\prime}}$ observed in the Kepler space photometric time series of some of the target SPB stars (see Sect. 5.3). Even the lowest order model of resonant nonlinear mode coupling developed in this work can thus aid future asteroseismic modeling of rapidly rotating SPB stars based on g modes, because it offers additional constraints on some of the observed combination frequencies based on stationary amplitude ratios of specific, resonantly coupled modes. This can in principle also serve as an additional source of mode identification. The logical next step in understanding nonlinear amplitude saturation is to apply this framework to the modeling of some of the SPB stars analyzed in V21 by integrating our monochromatic predictions over the Kepler passband.
It should be investigated whether larger resonantly coupled networks of modes or higherorder and r − g mode coupling can saturate some of the unexplained observed candidate couplings in SPB stars. Higherorder coupling will definitely already include selfsaturation effects that are neglected in the current formalism (see e.g., Van Hoolst et al. 1998 for an example of a study that included these effects for radial modes while neglecting rotation, and Gastine & Dintrans 2008a,b for direct numerical simulations of the κ mechanism and nonlinear saturation for radial modes in Cepheids). Mode triads containing g modes not considered in this work can also explain the presence of some of the highest observed amplitude ratios obtained in this work.
Other authors, such as Lee (2001) and Aprilia & Saio (2011), used a mode expansion formalism that did not invoke the TAR and accounted for linear mode couplings to describe the linear eigenmodes. They found that their formalism linearly stabilizes some of the g modes that were excited when the TAR was assumed (i.e., constituting a systematical error). That method, which truncates the series expansion at some manageable expansion order, is computationally more intensive, and might not be accurate (Dziembowski et al. 2007). Perhaps the hybrid expansion method discussed in Chapter 7 of Goupil et al. (2013) that describes linear mode coupling using expansions of modes whose angular eigenfunctions are described by Hough functions is a good compromise that softens the computational load with increasing expansion order, compared to the original expansion formalism. Because we find that the nonlinear saturation is crucially dependent on the linear excitation properties, indepth investigations of the linear vibrational instability of g modes in SPB stars using the latest improvements in opacity computations are crucial for future nonlinear excitation studies (e.g., DaszyńskaDaszkiewicz et al. 2017). Such investigations will allow us to estimate the systematical error we make in the theoretical predictions of surface luminosity fluctuation ratios due to opacity.
The strength of the Coriolis force is determined by the Coriolis frequency 2Ω. The TAR frequency hierarchies 2Ω≪N and Ω_{φ}≪N, where Ω_{φ} is the realvalued frequency of pulsation mode φ in the corotating reference frame, impose strong stratification in specific regions of the stellar interior, through which the lowfrequency waves described within the TAR propagate (see e.g., Dhouib et al. 2021a).
The repository can directly be accessed using the following link: https://doi.org/10.5281/zenodo.10814654.
Acknowledgments
The research leading to these results received funding from the KU Leuven Research Council (grant C16/18/005: PARADISE, with PIs CA and TVH). The research leading to part of these results has also received funding from the European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement N°101071505: 4DSTAR). Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. JVB and CA acknowledge support from the Research Foundation Flanders (FWO) under grant agreements N°V421221N (Travel Grant) and N°K802922N (Sabbatical leave). This research was also supported in part by the National Science Foundation under Grant No. PHY1748958. JVB is grateful to Dominic Bowman for the useful discussion and comments on the manuscript. JVB is also grateful to Michel Rieutord, Rony Keppens and Timothy Van Reeth for the useful comments on his PhD manuscript, which improved this work. CA is grateful for the kind hospitality offered by the staff of the Center for Computational Astrophysics at the Flatiron Institute of the Simons Foundation in New York during her work visits in the fall of 2022 and spring of 2023. Throughout this work we have made use of the following Python packages: LMFIT (Newville et al. 2020), SCIPY (Virtanen et al. 2020), NUMPY (Harris et al. 2020), PANDAS (McKinney et al. 2010), and MATPLOTLIB (Hunter 2007). We thank their authors for making these great software packages open source.
References
 Aerts, C. 2021, Rev. Mod. Phys., 93, 015001 [Google Scholar]
 Aerts, C., & Tkachenko, A. 2023, arXiv eprints [arXiv: 2311.08453] [Google Scholar]
 Aerts, C., de Pauw, M., & Waelkens, C. 1992, A&A, 266, 294 [NASA ADS] [Google Scholar]
 Aerts, C., ChristensenDalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Springer) [Google Scholar]
 Aerts, C., Molenberghs, G., Michielsen, M., et al. 2018, ApJS, 237, 15 [Google Scholar]
 Aikawa, T. 1984, MNRAS, 206, 833 [NASA ADS] [CrossRef] [Google Scholar]
 Aprilia, L. U., & Saio, H. 2011, MNRAS, 412, 2265 [NASA ADS] [CrossRef] [Google Scholar]
 Arras, P., Flanagan, E. E., Morsink, S. M., et al. 2003, ApJ, 591, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Bernstein, D. S. 2009, Matrix Mathematics: Theory, Facts, and Formulas, Second Edition (Princeton: Princeton University Press) [CrossRef] [Google Scholar]
 Betounes, D. 2010, Differential Equations: Theory and Applications (New York: SpringerVerlag) [CrossRef] [Google Scholar]
 Bhatia, R. 2007, Positive Definite Matrices (Princeton: Princeton University Press) [Google Scholar]
 Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, 648, A36 [EDP Sciences] [Google Scholar]
 Bowman, D. M., & Michielsen, M. 2021, A&A, 656, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bowman, D. M., Kurtz, D. W., Breger, M., Murphy, S. J., & Holdsworth, D. L. 2016, MNRAS, 460, 1970 [NASA ADS] [CrossRef] [Google Scholar]
 Buchler, J. R. 1993, Ap&SS, 210, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Buchler, J. R., & Goupil, M. J. 1984, ApJ, 279, 394 [CrossRef] [Google Scholar]
 Buchler, J. R., & Regev, O. 1983, A&A, 123, 331 [NASA ADS] [Google Scholar]
 Buchler, J. R., Goupil, M. J., & Serre, T. 1995, A&A, 296, 405 [NASA ADS] [Google Scholar]
 Buchler, J. R., Goupil, M. J., & Hansen, C. J. 1997, A&A, 321, 159 [NASA ADS] [Google Scholar]
 Burkart, J., Quataert, E., Arras, P., & Weinberg, N. N. 2012, MNRAS, 421, 983 [NASA ADS] [CrossRef] [Google Scholar]
 Burkart, J., Quataert, E., Arras, P., & Weinberg, N. N. 2013, MNRAS, 433, 332 [NASA ADS] [CrossRef] [Google Scholar]
 Burkart, J., Quataert, E., & Arras, P. 2014, MNRAS, 443, 2957 [NASA ADS] [CrossRef] [Google Scholar]
 Cesari, L. 1971, Asymptotic Behavior and Stability Problems in Ordinary Differential Equations (Berlin, Heidelberg: Springer), 1 [Google Scholar]
 Chaplin, W. J., Elsworth, Y., Davies, G. R., et al. 2014, MNRAS, 445, 946 [NASA ADS] [CrossRef] [Google Scholar]
 Cowling, T. G. 1941, MNRAS, 101, 367 [NASA ADS] [Google Scholar]
 Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (New York: Gordon and Breach) [Google Scholar]
 Dappen, W., & Perdang, J. 1985, A&A, 151, 174 [NASA ADS] [Google Scholar]
 DaszyńskaDaszkiewicz, J., Dziembowski, W. A., Jerzykiewicz, M., & Handler, G. 2015, MNRAS, 446, 1438 [CrossRef] [Google Scholar]
 DaszyńskaDaszkiewicz, J., Pamyatnykh, A. A., Walczak, P., et al. 2017, MNRAS, 466, 2284 [CrossRef] [Google Scholar]
 De Cat, P., & Aerts, C. 2002, A&A, 393, 965 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Degroote, P., Briquet, M., Catala, C., et al. 2009, A&A, 506, 111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Degroote, P., Aerts, C., Baglin, A., et al. 2010, Nature, 464, 259 [Google Scholar]
 Dhouib, H., Prat, V., Van Reeth, T., & Mathis, S. 2021a, A&A, 652, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dhouib, H., Prat, V., Van Reeth, T., & Mathis, S. 2021b, A&A, 656, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dziembowski, W. 1982, Acta Astron., 32, 147 [NASA ADS] [Google Scholar]
 Dziembowski, W. A. 1993, ASP Conf. Ser., 40, 521 [NASA ADS] [Google Scholar]
 Dziembowski, W., & Krolikowska, M. 1985, Acta Astron., 35, 5 [NASA ADS] [Google Scholar]
 Dziembowski, W., Krolikowska, M., & Kosovichev, A. 1988, Acta Astron., 38, 61 [NASA ADS] [Google Scholar]
 Dziembowski, W. A., Moskalik, P., & Pamyatnykh, A. A. 1993, MNRAS, 265, 588 [NASA ADS] [Google Scholar]
 Dziembowski, W. A., DaszyńskaDaszkiewicz, J., & Pamyatnykh, A. A. 2007, MNRAS, 374, 248 [NASA ADS] [CrossRef] [Google Scholar]
 Essick, R., & Weinberg, N. N. 2016, ApJ, 816, 18 [Google Scholar]
 Frieman, E., & Rotenberg, M. 1960, Rev. Mod. Phys., 32, 898 [Google Scholar]
 Friedman, J. L., & Schutz, B. F. 1978a, ApJ, 221, 937 [NASA ADS] [CrossRef] [Google Scholar]
 Friedman, J. L., & Schutz, B. F. 1978b, ApJ, 222, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J. 2017, MNRAS, 472, 1538 [Google Scholar]
 Fuller, J., & Lai, D. 2012, MNRAS, 420, 3126 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J., Derekas, A., Borkovits, T., et al. 2013, MNRAS, 429, 2425 [CrossRef] [Google Scholar]
 Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867 [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gastine, T., & Dintrans, B. 2008a, A&A, 484, 29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gastine, T., & Dintrans, B. 2008b, A&A, 490, 743 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gautschy, A., & Saio, H. 1993, MNRAS, 262, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Gill, A. 1982, AtmosphereOcean Dynamics (New York: Academic Press) [Google Scholar]
 Goldstein, J., & Townsend, R. H. D. 2020, ApJ, 899, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Goupil, M.J., & Buchler, J. R. 1994, A&A, 291, 481 [NASA ADS] [Google Scholar]
 Goupil, M. J., Dziembowski, W. A., & Fontaine, G. 1998, Baltic Astron., 7, 21 [NASA ADS] [Google Scholar]
 Goupil, M., Belkacem, K., Neiner, C., Lignières, F., & Green, J. J. 2013, Studying Stellar Rotation and Convection: Theoretical Background and Seismic Diagnostics, 865 [CrossRef] [Google Scholar]
 Guckenheimer, J., & Holmes, P. 1983, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (New York, NY: Springer) [CrossRef] [Google Scholar]
 Guo, Z. 2020, ApJ, 896, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, Z. 2021, Front. Astron. Space Sci., 8, 67 [NASA ADS] [Google Scholar]
 Hahn, W. 1967, Stability of Motion., Die Grundlehren der mathematischen Wissenschaften, in Einzeldarstellungen mit besonderer Berücksichtigung der Anwendungsgebiete: 138 (Berlin Heidelberg: Springer) [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Henneco, J., Van Reeth, T., Prat, V., et al. 2021, A&A, 648, A97 [EDP Sciences] [Google Scholar]
 Hogg, D. W. 2008, arXiv eprints [arXiv:0807.4820] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Koch, D. G., Borucki, W. J., Basri, G., et al. 2010, ApJ, 713, L79 [Google Scholar]
 Kubicek, M., & Marek, M. 1983, Computational Methods in Bifurcation Theory and Dissipative Structures (Berlin, Heidelberg: Springer) [CrossRef] [Google Scholar]
 Kumar, P., & Goodman, J. 1996, ApJ, 466, 946 [NASA ADS] [CrossRef] [Google Scholar]
 Lai, D., & Wu, Y. 2006, Phys. Rev. D, 74, 024007 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U. 2001, ApJ, 557, 311 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U. 2012, MNRAS, 420, 2387 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U. 2022, MNRAS, 513, 2522 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [Google Scholar]
 LonguetHiggins, M. S. 1968, Philos. Trans. R. Soc. London Ser. A, 262, 511 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293 [Google Scholar]
 Mathis, S. 2013, in Transport Processes in Stellar Interiors, eds. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, 865, 23 [Google Scholar]
 Mathis, S., & Prat, V. 2019, A&A, 631, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56 [Google Scholar]
 Michielsen, M., Pedersen, M. G., Augustson, K. C., Mathis, S., & Aerts, C. 2019, A&A, 628, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michielsen, M., Aerts, C., & Bowman, D. M. 2021, A&A, 650, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miglio, A., Montalbán, J., Noels, A., & Eggenberger, P. 2008, MNRAS, 386, 1487 [Google Scholar]
 Moravveji, E., Aerts, C., Pápics, P. I., Triana, S. A., & Vandoren, B. 2015, A&A, 580, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moravveji, E., Townsend, R. H. D., Aerts, C., & Mathis, S. 2016, ApJ, 823, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Morsink, S. M. 2002, ApJ, 571, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Moskalik, P. 1985, Acta Astron., 35, 229 [NASA ADS] [Google Scholar]
 Mourabit, M., & Weinberg, N. N. 2023, ApJ, 950, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Nayfeh, A. H. 1973, Perturbation Methods (New York: Wiley) [Google Scholar]
 Nayfeh, A. H. 1981, Introduction to Perturbation Techniques (New York: Wiley) [Google Scholar]
 Nayfeh, A. H., & Mook, D. T. 1979, Nonlinear Oscillations (Wiley) [Google Scholar]
 Newville, M., Otten, R., Nelson, A., et al. 2020, https://doi.org/10.5281/zenodo.3814709 [Google Scholar]
 Nieva, M. F., & Przybilla, N. 2012, A&A, 539, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 O’Leary, R. M., & Burkart, J. 2014, MNRAS, 440, 3036 [CrossRef] [Google Scholar]
 Osaki, Y. 1971, PASJ, 23, 485 [NASA ADS] [Google Scholar]
 Pamyatnykh, A. A. 1999, Acta Astron., 49, 119 [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Pedersen, M. G. 2022, ApJ, 930, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Pedersen, M. G., Aerts, C., Pápics, P. I., & Rogers, T. M. 2018, A&A, 614, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nat. Astron., 5, 715 [NASA ADS] [CrossRef] [Google Scholar]
 Prat, V., Mathis, S., Buysschaert, B., et al. 2019, A&A, 627, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Przybilla, N., Nieva, M. F., Irrgang, A., & Butler, K. 2013, EAS Publ. Ser., 63, 13 [CrossRef] [EDP Sciences] [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [Google Scholar]
 Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 848, L1 [CrossRef] [Google Scholar]
 Schenk, A. K., Arras, P., Flanagan, É. É., Teukolsky, S. A., & Wasserman, I. 2001, Phys. Rev. D, 65, 024001 [CrossRef] [Google Scholar]
 Schutz, B. F. 1979, ApJ, 232, 874 [NASA ADS] [CrossRef] [Google Scholar]
 Schutz, B. F. 1980a, MNRAS, 190, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Schutz, B. F. 1980b, MNRAS, 190, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Seaton, M. J. 2005, MNRAS, 362, L1 [Google Scholar]
 Seydel, R. 2009, Practical Bifurcation and Stability Analysis (New York: SpringerVerlag) [Google Scholar]
 Szewczuk, W., & DaszyńskaDaszkiewicz, J. 2017, MNRAS, 469, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Szewczuk, W., & DaszyńskaDaszkiewicz, J. 2018, MNRAS, 478, 2243 [Google Scholar]
 Szewczuk, W., Walczak, P., & DaszyńskaDaszkiewicz, J. 2021, MNRAS, 503, 5894 [NASA ADS] [CrossRef] [Google Scholar]
 Szewczuk, W., Walczak, P., DaszyńskaDaszkiewicz, J., & Moździerski, D. 2022, MNRAS, 511, 1529 [NASA ADS] [CrossRef] [Google Scholar]
 Takeuti, M., & Buchler, J. R. 1993, Ap&SS, 210, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, R. H. D. 2003, MNRAS, 340, 1020 [Google Scholar]
 Townsend, R. H. D. 2005, MNRAS, 360, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, R. H. D., & Teitler, S. A. 2013, MNRAS, 435, 3406 [Google Scholar]
 Townsend, R. H. D., Goldstein, J., & Zweibel, E. G. 2018, MNRAS, 475, 879 [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: University of Tokyo Press) [Google Scholar]
 Ushomirsky, G., & Bildsten, L. 1998, ApJ, 497, L101 [NASA ADS] [CrossRef] [Google Scholar]
 Van Beeck, J. 2023, Ph.D. Thesis, KU Leuven, Belgium https://fys.kuleuven.be/ster/pub/phdthesisjordanvanbeeck/phdthesisjordanvanbeeck [Google Scholar]
 Van Beeck, J., Bowman, D. M., Pedersen, M. G., et al. 2021, A&A, 655, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Hoolst, T. 1994a, A&A, 292, 471 [NASA ADS] [Google Scholar]
 Van Hoolst, T. 1994b, A&A, 286, 879 [NASA ADS] [Google Scholar]
 Van Hoolst, T. 1995, A&A, 295, 371 [Google Scholar]
 Van Hoolst, T. 1996, A&A, 308, 66 [NASA ADS] [Google Scholar]
 Van Hoolst, T., & Smeyers, P. 1993, A&A, 279, 417 [NASA ADS] [Google Scholar]
 Van Hoolst, T., Dziembowski, W. A., & Kawaler, S. D. 1998, MNRAS, 297, 536 [NASA ADS] [CrossRef] [Google Scholar]
 Van Reeth, T., De Cat, P., Van Beeck, J., et al. 2022, A&A, 662, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vandakurov, Y. V. 1981, Sov. Astron. Lett., 7, 128 [NASA ADS] [Google Scholar]
 Vick, M., Lai, D., & Anderson, K. R. 2019, MNRAS, 484, 5645 [NASA ADS] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Waelkens, C. 1991, A&A, 246, 453 [NASA ADS] [Google Scholar]
 Walczak, P., DaszyńskaDaszkiewicz, J., Pigulski, A., et al. 2019, MNRAS, 485, 3544 [Google Scholar]
 Weinberg, N. N. 2016, ApJ, 819, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, N. N., & Arras, P. 2019, ApJ, 873, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, N. N., & Quataert, E. 2008, MNRAS, 387, L64 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, N. N., Arras, P., & Burkart, J. 2013, ApJ, 769, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, N. N., Arras, P., & Pramanik, D. 2021, ApJ, 918, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, T., & Li, Y. 2019, ApJ, 881, 86 [Google Scholar]
 Wu, Y., & Goldreich, P. 2001, ApJ, 546, 469 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, T., Li, Y., Deng, Z.M., et al. 2020, ApJ, 899, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, H., Weinberg, N. N., & Fuller, J. 2020, MNRAS, 496, 5482 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, H., Fuller, J., & Burdge, K. B. 2021a, MNRAS, 501, 1836 [Google Scholar]
 Yu, H., Weinberg, N. N., & Arras, P. 2021b, ApJ, 917, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, H., Weinberg, N. N., & Arras, P. 2022, ApJ, 928, 140 [CrossRef] [Google Scholar]
 Zanazzi, J. J., & Wu, Y. 2021, AJ, 161, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Zong, W., Charpinet, S., & Vauclair, G. 2016a, A&A, 594, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zong, W., Charpinet, S., Vauclair, G., Giammichele, N., & Van Grootel, V. 2016b, A&A, 585, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Phasespace representation of the oscillation model within the TAR
Equation (1) can be written as a dynamical system in phase space that satisfies (following S01)
$$\begin{array}{c}\hfill \dot{\mathit{\zeta}}=\mathsf{T}\mathbf{\xb7}\mathit{\zeta}+\mathsf{F}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(A.1)
in which the nonHermitian operator T is defined as
$$\begin{array}{c}\hfill \mathsf{T}=\left[\begin{array}{cc}\phantom{\rule{0.166667em}{0ex}}\frac{1}{2}\mathit{B}& \mathbf{1}\phantom{\rule{0.166667em}{0ex}}\\ \phantom{\rule{0.166667em}{0ex}}\mathit{C}+\frac{1}{4}{\mathit{B}}^{2}& \frac{1}{2}\mathit{B}\phantom{\rule{0.166667em}{0ex}}\\ \end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(A.2)
where ζ and F are defined as
$$\begin{array}{c}\hfill \mathit{\zeta}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)\equiv \left[\begin{array}{c}\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)\\ \mathit{\pi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\mathsf{F}=\left[\begin{array}{c}0\\ {\mathit{a}}_{\mathrm{ext}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(A.3)
with x the position vector, t the time coordinate, and
$$\begin{array}{c}\hfill \mathit{\pi}=\dot{\mathit{\xi}}+\frac{1}{2}\mathit{B}\left(\mathit{\xi}\right)\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(A.4)
By making a time dependence Ansatz similar to Eq. (2),
$$\begin{array}{c}\hfill \mathit{\zeta}(\mathit{x},t)=\mathit{\zeta}\left(\mathit{x}\right)\phantom{\rule{0.166667em}{0ex}}{e}^{i\omega t}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(A.5)
the equations of motion for free linear oscillations yield the eigenvalue problem
$$\begin{array}{c}\hfill [\mathsf{T}+i\omega ]\mathbf{\xb7}\phantom{\rule{0.166667em}{0ex}}\mathit{\zeta}\left(\mathit{x}\right)=0\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(A.6)
Because T is nonHermitian, one should distinguish between its right and left eigenvectors. The right eigenvectors ζ_{φ} of operator T are the solutions of
$$\begin{array}{c}\hfill [\mathsf{T}+i{\omega}_{\phi}]\mathbf{\xb7}\phantom{\rule{0.166667em}{0ex}}{\mathit{\zeta}}_{\phi}\left(\mathit{x}\right)=0\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(A.7)
These right eigenvectors defined by Eq. (A.7) are the eigenvectors of nonconjugated eigenmodes satisfying Eq. (1). We refer the reader to S01 for details on the left eigenvector problem.
Appendix B: Proof of the orthogonality condition
We give an explicit proof of Eq. (4) that employs the rotatingframe symplectic product W(ξ_{φ},ξ_{β}) , which was defined by Friedman & Schutz (1978a) as
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta})=& \phantom{\rule{0.166667em}{0ex}}\langle {\mathit{\xi}}_{\phi},\phantom{\rule{0.166667em}{0ex}}\dot{{\mathit{\xi}}_{\beta}}+\frac{1}{2}\mathit{B}({\mathit{\xi}}_{\beta})\rangle \langle \dot{{\mathit{\xi}}_{\phi}}+\frac{1}{2}\mathit{B}({\mathit{\xi}}_{\phi}),{\mathit{\xi}}_{\beta}\rangle \phantom{\rule{0.166667em}{0ex}},\hfill \end{array}\end{array}$$(B.1)
where the definition of the inner product in Eq. (5) was used. Because of the time dependence Ansatz (2), Eq. (B.1) becomes
$$\begin{array}{c}\hfill iW({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta})=({\omega}_{\beta}+{\omega}_{\phi}^{\ast})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}\rangle +\langle {\mathit{\xi}}_{\phi},i\mathit{B}({\mathit{\xi}}_{\beta})\rangle \phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.2)
after multiplication with i, and using the antiHermiticity of B. The symplectic product W is a conserved quantity (e.g., Friedman & Schutz 1978b), hence, we have that
$$\begin{array}{c}\hfill \frac{\partial W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta})}{\partial t}=0=i({\omega}_{\phi}^{\ast}{\omega}_{\beta})W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta})\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.3)
because of Ansatz (2) and the assumed timeindependence of Ω. For ${\omega}_{\beta}\ne {\omega}_{\phi}^{\ast}$, we obtain from Eq. (B.3) that iW(ξ_{φ},ξ_{β}) = 0, which proves that in that case Eq. (4) holds. For modes with real eigenfrequencies Ω_{β} ≠ Ω_{φ}, Eqs. (B.2) and (B.3) imply that
$$\begin{array}{c}\hfill ({\mathrm{\Omega}}_{\beta}+{\mathrm{\Omega}}_{\phi})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}\rangle +\langle {\mathit{\xi}}_{\phi},i\mathit{B}({\mathit{\xi}}_{\beta})\rangle =0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.4)
which is equal to Eq. (6) for φ ≠ β. In the case of complex conjugate degeneracy (i.e., ${\omega}_{\beta}={\omega}_{\phi}^{\ast}$), or, in the case of degeneracy for realvalued eigenfrequencies (i.e., Ω_{β} = Ω_{φ}), Eqs. (B.3) and (B.4) yield no orthogonality conditions. We do not consider degenerate modes in this work, and instead refer the reader to S01 for information on the degenerate eigenvalue problem.
For the mode ξ_{φ}, we obtain from Eq. (B.1) that
$$\begin{array}{c}\hfill W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\phi})=i({\omega}_{\phi}+{\omega}_{\phi}^{\ast})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\phi}\rangle +\langle {\mathit{\xi}}_{\phi},\mathit{B}({\mathit{\xi}}_{\phi})\rangle \phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.5)
Because of the conserved nature of W(ξ_{φ},ξ_{φ}) and the assumed timeindependence of Ω, we write the time derivative of the symplectic product (B.5) as
$$\begin{array}{c}\hfill \frac{\partial W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\phi})}{\partial t}=i({\omega}_{\phi}^{\ast}{\omega}_{\phi})W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\phi})=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.6)
indicating that ${\omega}_{\phi}^{\ast}={\omega}_{\phi}$ (i.e., ω_{φ} is real) or iW(ξ_{φ},ξ_{φ}) = 0. For real eigenfrequencies Ω_{φ}, Eq. (B.6) is always fulfilled. This proves Eq. (6) for φ = β, because Eq. (B.5) multiplied by i is equal to the realvalued constant b_{φ} defined in Eq. (7) for a real eigenfrequency Ω_{φ}. For complex nondegenerate eigenfrequencies ω_{φ}, Eq. (B.6) is fulfilled only if b_{φ} = 0, because ${\omega}_{\phi}\ne {\omega}_{\phi}^{\ast}$.
An entirely equivalent proof can be constructed using left contractions of Eq. (3) and its complex conjugate. Lee (2022) provided an example of such a proof, relying on a different time dependence Ansatz. Their derived Eq. (11) is similar to Eq. (6), and is valid for nondegenerate realvalued eigenfrequencies. The lack of an orthogonality relation for the nonadiabatic complexvalued eigenfunctions of oscillations in rotating stars (with associated complexvalued eigenfrequencies) that does not involve Jordan chain modes (see Appendix A of Schenk et al. 2001) motivates the use of adiabatic realvalued eigenfunctions (with associated realvalued eigenfrequencies) in the computations for the mode coupling coefficient (22).
Appendix C: Expressions for the explicit terms of the threemode coupling coefficients
We base ourselves on explicit expressions derived by L12 for the terms that make up the threemode coupling coefficients κ_{ABC} defined in their Eq. (26), and correct some of these expressions before using them to compute the coupling coefficient ${\kappa}_{\phi}^{\beta \gamma}$ defined in Eq. (22). Additionally, we show how we compute the adiabatic derivative of the adiabatic exponent Γ_{1}, ${\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}$, which appears in one of the terms of ${\kappa}_{\phi}^{\beta \gamma}$.
C.1. Information on the explicit terms of ${\kappa}_{\phi}^{\beta \gamma}$
Equations (14), (15) and (16) contain covariant derivatives of the displacements ξ and the gravitational potential Φ, which are used to compute the nonlinear coupling term in the Cowling approximation. These covariant derivatives are computed in spherical coordinates using Eqs. (B6) to (B16) in appendix B of L12, in which we replace the Hough functions defined in L12, $\stackrel{\sim}{\mathrm{\Theta}}$, ${\stackrel{\sim}{\mathrm{\Theta}}}^{\theta}$ and ${\stackrel{\sim}{\mathrm{\Theta}}}^{\varphi}$, by their equivalents H_{r}(θ) e^{iϕφ}, H_{θ}(θ) e^{iϕφ} and H_{ϕ}(θ) e^{iϕφ} defined for a mode φ in this work. We further note that equation (B12) of L12 should read
$$\begin{array}{c}\hfill {\mathrm{\nabla}}_{r}\phantom{\rule{0.166667em}{0ex}}{\xi}^{\varphi}=\frac{1}{rsin\theta}\frac{\partial {\xi}^{\varphi}}{\partial r}=\frac{i}{r}\frac{\partial}{\partial r}\left(r\phantom{\rule{0.166667em}{0ex}}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}}\right)\frac{{H}_{\varphi}}{sin\theta}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(C.1)
in which z_{2} denotes the dimensionless horizontal displacement vector of a mode in the Cowling approximation, $\overline{\omega}$ is the dimensionless mode frequency of that mode, and c_{1} is a dimensionless ratio; all defined in L12.
Following S01 and L12, the terms in the coupling coefficient integral can be split up in four different integrals ${\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(u)}$ (u ∈ ⟦4⟧), similar to the definitions of ${\kappa}_{\mathit{ABC}}^{(u)}$ in Eqs. (B1) to (B5) of L12, or their equivalents ${\eta}_{1}^{(u)}\equiv {\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(u)}\phantom{\rule{0.166667em}{0ex}}/\phantom{\rule{0.166667em}{0ex}}{\u03f5}_{\phi}$. Equations (B2) and (B5) of L12 describe how to compute ${\kappa}_{\mathit{ABC}}^{(1)}$ and ${\kappa}_{\mathit{ABC}}^{(4)}$, and are employed in this work to compute the coupling integrals ${\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(1)}$ and ${\eta}_{1}^{(1)}$, and ${\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(4)}$ and ${\eta}_{1}^{(4)}$, respectively. We note that in Eq. (B3) of L12 for ${\kappa}_{\mathit{ABC}}^{(2)}$, the derivative $\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}$ should be replaced by its adiabatic equivalent, ${\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}$, before using that equation to compute the coupling integrals ${\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(2)}$ and ${\eta}_{1}^{(2)}$. We compute this adiabatic derivative as described in Appendix C.2. Each of the integrands of the coupling integrals ${\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(u)}$ or ${\eta}_{1}^{(u)}$ consists of a factor composed only of quantities derived from the stellar equilibrium structure, multiplied with some expression that either involves covariant derivatives of the eigenfunctions, or covariant derivatives of the gravitational potential Φ and the eigenfunctions themselves (for ${\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(4)}$ and ${\eta}_{1}^{(4)}$). The explicit expressions of these latter factors of the integrands are constructed from the expressions in Eqs. (B17) to (B22) of L12. We note that the last four terms of Eq. (B19) in L12 need to be multiplied with a factor 2 before using them to compute the coupling integrals ${\kappa}_{\phi}^{\beta \gamma \phantom{\rule{0.166667em}{0ex}}(3)}$ and ${\eta}_{1}^{(3)}$, and that the first two terms of their Eq. (B19) can be simplified, because
$$\begin{array}{c}\hfill \begin{array}{cc}& S(\frac{\partial (r\phantom{\rule{0.166667em}{0ex}}{z}_{1})}{\partial r}({z}_{1}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}})\frac{\partial}{\partial r}\left(r\phantom{\rule{0.166667em}{0ex}}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}}\right)\phantom{\rule{0.166667em}{0ex}}:\phantom{\rule{0.166667em}{0ex}}{H}_{r}[{H}_{\theta}{H}_{\theta}{H}_{\varphi}{H}_{\varphi}])\phantom{\rule{0.166667em}{0ex}}+\phantom{\rule{0.166667em}{0ex}}\hfill \\ \hfill & S(\frac{\partial (r\phantom{\rule{0.166667em}{0ex}}{z}_{1})}{\partial r}\phantom{\rule{0.166667em}{0ex}}{z}_{1}\phantom{\rule{0.166667em}{0ex}}\frac{\partial}{\partial r}\left(r\phantom{\rule{0.166667em}{0ex}}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}}\right)\phantom{\rule{0.166667em}{0ex}}:{H}_{r}[\delta {H}_{r}{H}_{\theta}\frac{m\phantom{\rule{0.166667em}{0ex}}{H}_{r}\phantom{\rule{0.166667em}{0ex}}{H}_{\varphi}}{sin\theta}]\phantom{\rule{0.166667em}{0ex}})\phantom{\rule{0.166667em}{0ex}}\hfill \\ \hfill & S(\frac{\partial (r\phantom{\rule{0.166667em}{0ex}}{z}_{1})}{\partial r}\phantom{\rule{0.166667em}{0ex}}{z}_{1}\phantom{\rule{0.166667em}{0ex}}\frac{\partial}{\partial r}\left(r\phantom{\rule{0.166667em}{0ex}}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}}\right)\phantom{\rule{0.166667em}{0ex}}:\phantom{\rule{0.166667em}{0ex}}{H}_{r}[{H}_{\theta}{H}_{\theta}{H}_{\varphi}{H}_{\varphi}])=\hfill \\ \hfill & S(\frac{\partial (r\phantom{\rule{0.166667em}{0ex}}{z}_{1})}{\partial r}\phantom{\rule{0.166667em}{0ex}}{z}_{1}\phantom{\rule{0.166667em}{0ex}}\frac{\partial}{\partial r}\left(r\phantom{\rule{0.166667em}{0ex}}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}}\right)\phantom{\rule{0.166667em}{0ex}}:{H}_{r}[\delta {H}_{r}{H}_{\theta}\frac{m\phantom{\rule{0.166667em}{0ex}}{H}_{r}\phantom{\rule{0.166667em}{0ex}}{H}_{\varphi}}{sin\theta}]\phantom{\rule{0.166667em}{0ex}})\phantom{\rule{0.166667em}{0ex}}\hfill \\ \hfill & S(\frac{\partial (r\phantom{\rule{0.166667em}{0ex}}{z}_{1})}{\partial r}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}}\frac{\partial}{\partial r}\left(r\phantom{\rule{0.166667em}{0ex}}\frac{{z}_{2}}{{c}_{1}\phantom{\rule{0.166667em}{0ex}}{\overline{\omega}}^{2}}\right)\phantom{\rule{0.166667em}{0ex}}:\phantom{\rule{0.166667em}{0ex}}{H}_{r}[{H}_{\theta}{H}_{\theta}{H}_{\varphi}{H}_{\varphi}])\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}\end{array}$$(C.2)
in which we separated the factors containing the azimuthal angle ϕ, and where δH_{r} replaces $\delta \stackrel{\sim}{\mathrm{\Theta}}$, which was defined together with the symmetrical operator S in L12.
Finally, as noted by for example L12, the coupling coefficient integrals over the spherical star can be transformed into products of integrals over the radial part of the integrand with integrals over the angular part of the integrand. The selection rules discussed in Sect. 2.4 originate from the integral over the angular part of the integrand of these coupling integrals (see appendix E for the explicit dependence of the selection rules on the Hough functions).
C.2. A computationfriendly form of ${\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}$
In this section, we show that the adiabatic derivative ${\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}$, which should be used in Eq. (B3) of L12, can be written as a function of nonadiabatic derivatives that are computed in stellar evolution codes such as MESA (Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019). First, let us write a twodimensional Jacobian as
$$\begin{array}{c}\hfill \frac{\partial (u,\phantom{\rule{0.166667em}{0ex}}v)}{\partial (x,\phantom{\rule{0.166667em}{0ex}}y)}\equiv \left\begin{array}{cc}{\left(\frac{\partial u}{\partial x}\right)}_{y}& {\left(\frac{\partial u}{\partial y}\right)}_{x}\\ {\left(\frac{\partial v}{\partial x}\right)}_{y}& {\left(\frac{\partial v}{\partial y}\right)}_{x}\\ \end{array}\right\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(C.3)
Using Eq. (C.3), the expression for ${\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}$ becomes
$$\begin{array}{c}\hfill {\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}=\rho \phantom{\rule{0.166667em}{0ex}}\left\begin{array}{cc}{\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial \rho}\right)}_{S}& {\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial S}\right)}_{\rho}\\ 0& 1\\ \end{array}\right=\rho \phantom{\rule{0.166667em}{0ex}}\frac{\partial ({\mathrm{\Gamma}}_{1},\phantom{\rule{0.166667em}{0ex}}S)}{\partial (\rho ,\phantom{\rule{0.166667em}{0ex}}S)}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(C.4)
If we then multiply Eq. (C.4) with a Jacobian equal to 1, we can write
$$\begin{array}{c}\hfill {\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}=\rho \phantom{\rule{0.166667em}{0ex}}\frac{\partial ({\mathrm{\Gamma}}_{1},\phantom{\rule{0.166667em}{0ex}}S)}{\partial (\rho ,\phantom{\rule{0.166667em}{0ex}}S)}\frac{\partial (\rho ,\phantom{\rule{0.166667em}{0ex}}T)}{\partial (\rho ,\phantom{\rule{0.166667em}{0ex}}T)}=\rho \phantom{\rule{0.166667em}{0ex}}\frac{\frac{\partial ({\mathrm{\Gamma}}_{1},\phantom{\rule{0.166667em}{0ex}}S)}{\partial (\rho ,\phantom{\rule{0.166667em}{0ex}}T)}}{\frac{\partial (\rho ,\phantom{\rule{0.166667em}{0ex}}S)}{\partial (\rho ,\phantom{\rule{0.166667em}{0ex}}T)}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(C.5)
which after using the definition of the Jacobian in Eq. (C.3) becomes
$$\begin{array}{c}\hfill {\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial ln\rho}\right)}_{S}=\rho \phantom{\rule{0.166667em}{0ex}}[{\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial \rho}\right)}_{T}\left\{\frac{{\left(\frac{\partial S}{\partial \rho}\right)}_{T}}{{\left(\frac{\partial S}{\partial T}\right)}_{\rho}}\right\}\phantom{\rule{0.166667em}{0ex}}{\left(\frac{\partial {\mathrm{\Gamma}}_{1}}{\partial T}\right)}_{\rho}]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(C.6)
whose right hand side only contains thermodynamic quantities that are computed in stellar evolution codes such as MESA.
An adiabatic secondorder derivative of adiabatic exponent Γ_{1} with respect to mass density ρ appears in coupling coefficients that describe fourmode interactions (or higherorder derivatives when higherorder coupling terms are considered). Higherorder derivatives of the thermodynamic quantities used in Eq. (C.6) are not readily computed in stellar evolution codes such as MESA, but would be necessary to obtain an analytical expression for the adiabatic higherorder derivatives of Γ_{1}. Such higherorder derivatives can alternatively be estimated from numerical derivatives of Eq. (C.6) when using the output of current stateoftheart stellar evolution codes.
Appendix D: Defining c_{φ}(t) and ${c}_{\phi}^{\ast}(t)$
We prove the validity of the expressions (18) for c_{φ}(t) and ${c}_{\phi}^{\ast}(t)$ by making use of the orthogonality relations among the modes. First, let us derive the complex orthogonality relations containing complex conjugate modes using the approach in Appendix B. For $W({\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta})$, $W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}^{\ast})$ and $W({\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}^{\ast})$) we obtain from the definition of the symplectic product in Eq. (B.1) that
$$\begin{array}{cc}\hfill W({\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta})=& \phantom{\rule{0.166667em}{0ex}}i({\omega}_{\beta}{\omega}_{\phi})\langle {\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}\rangle +\langle {\mathit{\xi}}_{\phi}^{\ast},\mathit{B}({\mathit{\xi}}_{\beta})\rangle \phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(D.1a)
$$\begin{array}{cc}\hfill W({\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}^{\ast})=& \phantom{\rule{0.166667em}{0ex}}i({\omega}_{\phi}^{\ast}{\omega}_{\beta}^{\ast})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}^{\ast}\rangle +\langle {\mathit{\xi}}_{\phi},\mathit{B}({\mathit{\xi}}_{\beta}^{\ast})\rangle \phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(D.1b)
$$\begin{array}{cc}\hfill W({\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}^{\ast})=& \phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}i({\omega}_{\phi}+{\omega}_{\beta}^{\ast})\langle {\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}^{\ast}\rangle +\langle {\mathit{\xi}}_{\phi}^{\ast},\mathit{B}({\mathit{\xi}}_{\beta}^{\ast})\rangle \phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}$$(D.1c)
The orthogonality relation derived from Eq. (D.1a) becomes
$$\begin{array}{c}\hfill ({\omega}_{\beta}{\omega}_{\phi})\langle {\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}\rangle +\langle {\mathit{\xi}}_{\phi}^{\ast},i\mathit{B}({\mathit{\xi}}_{\beta})\rangle =0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.2)
when ω_{φ} ≠ −ω_{β}, which has a realvalued equivalent
$$\begin{array}{c}\hfill ({\mathrm{\Omega}}_{\beta}{\mathrm{\Omega}}_{\phi})\langle {\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}\rangle +\langle {\mathit{\xi}}_{\phi}^{\ast},i\mathit{B}({\mathit{\xi}}_{\beta})\rangle =0\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(D.3)
Orthogonality relation (D.3) is valid if Ω_{φ} ≠ −Ω_{β}. Similarly, we obtain for Eq. (D.1b) that
$$\begin{array}{c}\hfill ({\omega}_{\phi}^{\ast}{\omega}_{\beta}^{\ast})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}^{\ast}\rangle +\langle {\mathit{\xi}}_{\phi},i\mathit{B}({\mathit{\xi}}_{\beta}^{\ast})\rangle =0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.4)
when ${\omega}_{\phi}^{\ast}\ne {\omega}_{\beta}^{\ast}$. The realvalued equivalent of Eq. (D.4) is
$$\begin{array}{c}\hfill ({\mathrm{\Omega}}_{\phi}{\mathrm{\Omega}}_{\beta})\langle {\mathit{\xi}}_{\phi},{\mathit{\xi}}_{\beta}^{\ast}\rangle +\langle {\mathit{\xi}}_{\phi},i\mathit{B}({\mathit{\xi}}_{\beta}^{\ast})\rangle =0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.5)
when Ω_{φ} ≠ −Ω_{β}. Finally, we have that
$$\begin{array}{c}\hfill ({\omega}_{\beta}^{\ast}+{\omega}_{\phi})\langle {\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}^{\ast}\rangle \langle {\mathit{\xi}}_{\phi}^{\ast},i\mathit{B}({\mathit{\xi}}_{\beta}^{\ast})\rangle =0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.6)
for Eq. (D.1c) and ${\omega}_{\phi}\ne {\omega}_{\beta}^{\ast}$. Its realvalued equivalent
$$\begin{array}{c}\hfill ({\mathrm{\Omega}}_{\beta}+{\mathrm{\Omega}}_{\phi})\langle {\mathit{\xi}}_{\phi}^{\ast},{\mathit{\xi}}_{\beta}^{\ast}\rangle \langle {\mathit{\xi}}_{\phi}^{\ast},i\mathit{B}({\mathit{\xi}}_{\beta}^{\ast})\rangle ={\delta}_{\beta}^{\phi}\phantom{\rule{0.166667em}{0ex}}{b}_{\phi}^{\ast}={\delta}_{\beta}^{\phi}\phantom{\rule{0.166667em}{0ex}}{b}_{\phi}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.7)
is the complex conjugate of Eq. (6), and is valid for nondegenerate eigenfrequencies (i.e., Ω_{φ} ≠ Ω_{β}). The second equality in Eq. (D.7) is valid because b_{β} is a real number for the adiabatic eigenfunctions used in this work.
Multiplying the term ${\mathrm{\Omega}}_{\beta}\phantom{\rule{0.166667em}{0ex}}\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)+i\phantom{\rule{0.166667em}{0ex}}\dot{\mathit{\xi}}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)+i\mathit{B}\left(\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)\right)$ with ${\mathit{\xi}}_{\beta}^{\ast}$ and integrating over the stellar mass yields
$$\begin{array}{c}\hfill \langle {\mathit{\xi}}_{\beta},\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\beta}\phantom{\rule{0.166667em}{0ex}}\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)+i\phantom{\rule{0.166667em}{0ex}}\dot{\mathit{\xi}}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)+i\mathit{B}\left(\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)\right)\rangle \equiv \mathcal{P}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(D.8)
By virtue of the phase space mode expansion (17), we write the quantity 𝒫 as
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{P}=& {\displaystyle \sum _{\phi}}\phantom{\rule{0.166667em}{0ex}}({\mathrm{\Omega}}_{\beta}+{\mathrm{\Omega}}_{\phi})\phantom{\rule{0.166667em}{0ex}}{c}_{\phi}\left(t\right)\langle {\mathit{\xi}}_{\beta}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\rangle \hfill \\ \hfill & \phantom{\rule{2em}{0ex}}+({\mathrm{\Omega}}_{\beta}{\mathrm{\Omega}}_{\phi})\phantom{\rule{0.166667em}{0ex}}{c}_{\phi}^{\ast}\left(t\right)\langle {\mathit{\xi}}_{\beta}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right)\rangle \hfill \\ \hfill & \phantom{\rule{2em}{0ex}}+{c}_{\phi}\left(t\right)\langle {\mathit{\xi}}_{\beta}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}i\mathit{B}\left({\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\right)\rangle \hfill \\ \hfill & \phantom{\rule{2em}{0ex}}+{c}_{\phi}^{\ast}\left(t\right)\langle {\mathit{\xi}}_{\beta}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}i\mathit{B}\left({\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right)\right)\rangle \phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}\end{array}$$(D.9)
If we then substitute the orthogonality relation (6) and a relabeled version of orthogonality relation (D.5) in Eq. (D.9), we get for nondegenerate modes φ and β that
$$\begin{array}{c}\hfill \mathcal{P}={c}_{\beta}\left(t\right)\phantom{\rule{0.166667em}{0ex}}{b}_{\beta}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.10)
assuming Ω_{φ} ≠ −Ω_{β}, and using the Kronecker delta in Eq. (6). Equation (D.10) thus proves Eq. (18a) for nondegenerate modes φ and β that have Ω_{β} ≠ −Ω_{φ}, and a nonzero b_{β}.
Similarly, if we multiply ${\mathrm{\Omega}}_{\beta}\phantom{\rule{0.166667em}{0ex}}\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)i\phantom{\rule{0.166667em}{0ex}}\dot{\mathit{\xi}}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)i\mathit{B}\left(\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)\right)$ with ξ_{β} and integrate over the stellar mass, we obtain
$$\begin{array}{c}\hfill \langle {\mathit{\xi}}_{\beta}^{\ast},\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\beta}\phantom{\rule{0.166667em}{0ex}}\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)i\phantom{\rule{0.166667em}{0ex}}\dot{\mathit{\xi}}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)i\mathit{B}\left(\mathit{\xi}(\mathit{x},\phantom{\rule{0.166667em}{0ex}}t)\right)\rangle \equiv \mathcal{Q}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.11)
which, by virtue of the phase space mode expansion (17) can be written as
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{Q}=& \phantom{\rule{0.166667em}{0ex}}{\displaystyle \sum _{\phi}}\phantom{\rule{0.166667em}{0ex}}({\mathrm{\Omega}}_{\beta}{\mathrm{\Omega}}_{\phi})\phantom{\rule{0.166667em}{0ex}}{c}_{\phi}\left(t\right)\langle {\mathit{\xi}}_{\beta}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\rangle \hfill \\ \hfill & \phantom{\rule{2em}{0ex}}+({\mathrm{\Omega}}_{\beta}+{\mathrm{\Omega}}_{\phi})\phantom{\rule{0.166667em}{0ex}}{c}_{\phi}^{\ast}\left(t\right)\langle {\mathit{\xi}}_{\beta}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right)\rangle \hfill \\ \hfill & \phantom{\rule{2em}{0ex}}{c}_{\phi}\left(t\right)\langle {\mathit{\xi}}_{\beta}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}i\mathit{B}\left({\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\right)\rangle \hfill \\ \hfill & \phantom{\rule{2em}{0ex}}{c}_{\phi}^{\ast}\left(t\right)\langle {\mathit{\xi}}_{\beta}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}i\mathit{B}\left({\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right)\right)\rangle \phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}\end{array}$$(D.12)
The following relations
$$\begin{array}{cc}& \langle {\mathit{\xi}}_{\beta}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\rangle =\langle {\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\beta}\left(\mathit{x}\right)\rangle ,\hfill \end{array}$$(D.13a)
$$\begin{array}{cc}& \langle {\mathit{\xi}}_{\beta}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}i\mathit{B}\left({\mathit{\xi}}_{\phi}\left(\mathit{x}\right)\right)\rangle \begin{array}{cc}& =\langle {\mathit{\xi}}_{\phi}^{\ast}\left(\mathit{x}\right),\phantom{\rule{0.166667em}{0ex}}i\mathit{B}\left({\mathit{\xi}}_{\beta}\left(\mathit{x}\right)\right)\rangle ,\hfill \end{array}\hfill \end{array}$$(D.13b)
are valid because of the definition of the inner product and the antiHermiticity of the B operator. Substituting Eq. (D.13) in orthogonality relation (D.3) yields
$$\begin{array}{c}\hfill ({\mathrm{\Omega}}_{\beta}{\mathrm{\Omega}}_{\phi})\langle {\mathit{\xi}}_{\beta}^{\ast},{\mathit{\xi}}_{\phi}\rangle \langle {\mathit{\xi}}_{\beta}^{\ast},i\mathit{B}({\mathit{\xi}}_{\phi})\rangle =0\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(D.14)
By substituting the orthogonality relations (D.7) and (D.14) in Eq. (D.12), we obtain for Ω_{β} ≠ −Ω_{φ} that
$$\begin{array}{c}\hfill \mathcal{Q}={c}_{\beta}^{\ast}\left(t\right)\phantom{\rule{0.166667em}{0ex}}{b}_{\beta}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(D.15)
where we use the Kronecker delta in Eq. (D.7). This proves Eq. (18b) for nonzero b_{β} and Ω_{β} ≠ −Ω_{φ}.
Appendix E: Deriving the meridional selection rule
We derive the meridional selection rule mentioned in Sect. 2.4 based on the fact that a rotating star is invariant under the map f: (r, θ, ϕ)→(r, π − θ, ϕ), in analogy to what was done in S01 for their coupling coefficient. This invariance implies that modes are either even or odd under the pullback operator f_{*} (e.g., S01)
$$\begin{array}{c}\hfill {f}_{\ast}\mathit{\xi}={\xi}^{\widehat{r}}(r,\pi \theta ,\varphi )\phantom{\rule{0.166667em}{0ex}}{\mathit{e}}_{\widehat{r}}{\xi}^{\widehat{\theta}}(r,\pi \theta ,\varphi )\phantom{\rule{0.166667em}{0ex}}{\mathit{e}}_{\widehat{\theta}}+{\xi}^{\widehat{\varphi}}(r,\pi \theta ,\varphi )\phantom{\rule{0.166667em}{0ex}}{\mathit{e}}_{\widehat{\varphi}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(E.1)
We then have
$$\begin{array}{c}\hfill {f}_{\ast}{\mathit{\xi}}_{\phi}={\mathcal{Z}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(E.2)
with 𝒵_{φ} (called the zparity in S01) equal to 1 if mode φ is an even mode and −1 if it is an odd mode. The radial Hough function H_{r}(θ), the latitudinal Hough function H_{θ}(θ), and the azimuthal Hough function H_{ϕ}(θ) defined in Prat et al. (2019) are realvalued. Based on Eq. (10) we then have that
$$\begin{array}{c}\hfill {f}_{\ast}{\mathit{\xi}}_{\phi}^{\ast}={\mathcal{Z}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\phi}^{\ast}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(E.3)
By applying the pullback operator f_{*} on the coupling coefficient, and taking into account the commutation of f with geometrical operations such as the computation of covariant derivatives (e.g., S01), we can write
$$\begin{array}{c}\hfill {f}_{\ast}{\kappa}_{\phi}^{\beta \gamma}=\kappa ({f}_{\ast}{\mathit{\xi}}_{\phi}^{\ast},\phantom{\rule{0.166667em}{0ex}}{f}_{\ast}{\mathit{\xi}}_{\beta},\phantom{\rule{0.166667em}{0ex}}{f}_{\ast}{\mathit{\xi}}_{\gamma})={\mathcal{Z}}_{\phi}\phantom{\rule{0.166667em}{0ex}}{\mathcal{Z}}_{\beta}\phantom{\rule{0.166667em}{0ex}}{\mathcal{Z}}_{\gamma}\phantom{\rule{0.166667em}{0ex}}\kappa ({\mathit{\xi}}_{\phi}^{\ast},\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\beta},\phantom{\rule{0.166667em}{0ex}}{\mathit{\xi}}_{\gamma})\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(E.4)
Because ${\kappa}_{\phi}^{\beta \gamma}$ is a number and therefore invariant under f_{*}, it follows that
$$\begin{array}{c}\hfill (1{\mathcal{Z}}_{\phi}{\mathcal{Z}}_{\beta}{\mathcal{Z}}_{\gamma}){\kappa}_{\phi}^{\beta \gamma}=0\phantom{\rule{0.166667em}{0ex}},\end{array}$$(E.5)
resulting in the meridional or zparity selection rule (S01)
$$\begin{array}{c}\hfill {i}_{\mathrm{odd}}=\mathrm{mod}(0\phantom{\rule{0.166667em}{0ex}},\phantom{\rule{0.166667em}{0ex}}2)\phantom{\rule{0.166667em}{0ex}},\end{array}$$(E.6)
where i_{odd} denotes the number of odd modes. Hence, for threemode coupling, the zparity selection rule or meridional selection rule implies that two or zero odd modes need to be involved in the coupling, as mentioned in Sect. 2.4, where this is cast into a condition on the sum of the ordering numbers of the interacting modes.
Appendix F: Equivalence of the amplitude equations for Ω_{1d} ≈ Ω_{2d} − Ω_{3d}
In this section, we show the equivalence of the AEs determined in Sect. 2.6 for the sumfrequencies Ω_{1} ≈ Ω_{2} + Ω_{3} with those derived for differencefrequencies Ω_{1d} ≈ Ω_{2d} − Ω_{3d} using the same methods employed in Sect. 2.6. This equivalence indicates that the stationary properties of resonant isolated sumfrequency triads that we derive in the main body of this paper are the same as those properties that one would derive for the (relabeled) difference frequencies for quadratic nonlinear mode couplings.
To prove this equivalence, we first derive the differencefrequency AEs using the same procedure as in Sect. 2.6. Here, we define the resonance condition Ω_{1d} ≈ Ω_{2d} − Ω_{3d} in terms of the differencefrequency detuning parameter δω_{d}:
$$\begin{array}{c}\hfill \mathfrak{J}\phantom{\rule{0.166667em}{0ex}}\delta {\omega}_{d}={\mathrm{\Omega}}_{{1}_{d}}\approx {\mathrm{\Omega}}_{{2}_{d}}{\mathrm{\Omega}}_{{3}_{d}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(F.1)
In Eq. (F.1), the subscript d refers to the fact that we are considering a differencefrequency resonance. If we then set the seculartermgenerating terms in Eq. (33) equal to zero, use the differencefrequency resonance condition (F.1), and introduce the linear growth or linear damping rates (similar to what was done in Sect. 2.6), we obtain the extended complex AEs for the differencefrequency resonance Ω_{1d} ≈ Ω_{2d} − Ω_{3d}:
$$\begin{array}{c}\hfill \frac{\partial {\mathit{a}}_{d}}{\partial {t}_{1}}={\mathit{\gamma}}_{d}\xb0{\mathit{a}}_{d}+2\phantom{\rule{0.166667em}{0ex}}i\phantom{\rule{0.166667em}{0ex}}({\mathit{\eta}}_{\mathit{d}}\xb0{\mathit{a}}_{{M}_{d}}\xb0{\mathbf{\Omega}}_{{M}_{d}}\xb0{\mathit{e}}_{d})\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(F.2)
In the extended complex AEs (F.2), we use the vectors
$$\begin{array}{c}\hfill {\mathit{a}}_{d}=\left[\begin{array}{c}{a}_{{1}_{d}}\\ {a}_{{2}_{d}}\\ {a}_{{3}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{a}}_{{M}_{d}}=\left[\begin{array}{c}{a}_{{2}_{d}}\phantom{\rule{0.166667em}{0ex}}{a}_{{3}_{d}}^{\ast}\\ {a}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}{a}_{{3}_{d}}\\ {a}_{{1}_{d}}^{\ast}\phantom{\rule{0.166667em}{0ex}}{a}_{{2}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{e}}_{d}=\left[\begin{array}{c}exp(i\phantom{\rule{0.166667em}{0ex}}\delta {\omega}_{d}\phantom{\rule{0.166667em}{0ex}}{t}_{1})\\ exp(i\phantom{\rule{0.166667em}{0ex}}\delta {\omega}_{d}\phantom{\rule{0.166667em}{0ex}}{t}_{1})\\ exp(i\phantom{\rule{0.166667em}{0ex}}\delta {\omega}_{d}\phantom{\rule{0.166667em}{0ex}}{t}_{1})\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(F.3)
and
$$\begin{array}{c}\hfill {\mathit{\gamma}}_{d}=\left[\begin{array}{c}{\gamma}_{{1}_{d}}\\ {\gamma}_{{2}_{d}}\\ {\gamma}_{{3}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{\eta}}_{d}=\left[\begin{array}{c}{\eta}_{{1}_{d}}^{\ast}\\ {\eta}_{{1}_{d}}\\ {\eta}_{{1}_{d}}^{\ast}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathbf{\Omega}}_{{M}_{d}}=\left[\begin{array}{c}{\mathrm{\Omega}}_{{1}_{d}}\\ {\mathrm{\Omega}}_{{2}_{d}}\\ {\mathrm{\Omega}}_{{3}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(F.4)
in which we use the (differencefrequency) coupling coefficient η_{1d} that is defined as:
$$\begin{array}{c}\hfill {\eta}_{{2}_{d}}^{{1}_{d}{3}_{d}}={\eta}_{{2}_{d}}^{{3}_{d}{1}_{d}}={\left({\eta}_{{1}_{d}}^{{2}_{d}{\overline{3}}_{d}}\right)}^{\ast}={\left({\eta}_{{1}_{d}}^{{\overline{3}}_{d}{2}_{d}}\right)}^{\ast}={\left({\eta}_{{3}_{d}}^{{2}_{d}{\overline{1}}_{d}}\right)}^{\ast}={\left({\eta}_{{3}_{d}}^{{\overline{1}}_{d}{2}_{d}}\right)}^{\ast}\equiv {\eta}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(F.5)
By introducing real amplitudes A_{φd} and phases ϕ_{φd} as in Sect. 2.6, a_{φd} = A_{φd} exp(i ϕ_{φd}) (with φ ∈ ⟦3⟧), we obtain the (realvalued) AEs for difference frequencies,
$$\begin{array}{cc}\hfill \frac{\partial {\mathit{A}}_{d}}{\partial {t}_{1}}& ={\mathit{\gamma}}_{d}\xb0{\mathit{A}}_{d}+2\phantom{\rule{0.166667em}{0ex}}{\eta}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}sin\left({\mathrm{{\rm Y}}}_{d}\right)\phantom{\rule{0.166667em}{0ex}}({\mathit{A}}_{{N}_{d}}\xb0{\mathbf{\Omega}}_{{M}_{d}})\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(F.6a)
$$\begin{array}{cc}\hfill {\mathit{A}}_{d}\xb0\frac{\partial {\mathit{\varphi}}_{d}}{\partial {t}_{1}}& =2\phantom{\rule{0.166667em}{0ex}}{\eta}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}cos\left({\mathrm{{\rm Y}}}_{d}\right)\phantom{\rule{0.166667em}{0ex}}({\mathit{A}}_{{P}_{d}}\xb0{\mathbf{\Omega}}_{{M}_{d}})\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(F.6b)
where
$$\begin{array}{c}\hfill {\eta}_{{1}_{d}}={\eta}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}{e}^{i\phantom{\rule{0.166667em}{0ex}}{\delta}_{{1}_{d}}}\phantom{\rule{0.166667em}{0ex}},{\mathit{A}}_{d}=\left[\begin{array}{c}{A}_{{1}_{d}}\\ {A}_{{2}_{d}}\\ {A}_{{3}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{\varphi}}_{d}=\left[\begin{array}{c}{\varphi}_{{1}_{d}}\\ {\varphi}_{{2}_{d}}\\ {\varphi}_{{3}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(F.7)
as well as
$$\begin{array}{c}\hfill {\mathit{A}}_{{N}_{d}}=\left[\begin{array}{c}{A}_{{2}_{d}}\phantom{\rule{0.166667em}{0ex}}{A}_{{3}_{d}}\\ {A}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}{A}_{{3}_{d}}\\ {A}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}{A}_{{2}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},{\mathit{A}}_{{P}_{d}}=\left[\begin{array}{c}{A}_{{2}_{d}}\phantom{\rule{0.166667em}{0ex}}{A}_{{3}_{d}}\\ {A}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}{A}_{{3}_{d}}\\ {A}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}}{A}_{{2}_{d}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(F.8)
In Eq. (F.6), we define the generic (frequencydifference) phase coordinate Υ_{d} as
$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}_{d}\equiv \delta {\omega}_{d}\phantom{\rule{0.166667em}{0ex}}{t}_{1}+{\varphi}_{{1}_{d}}{\varphi}_{{2}_{d}}+{\varphi}_{{3}_{d}}{\delta}_{{1}_{d}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(F.9)
which contains the combination phase Φ_{d} = ϕ_{1d} − ϕ_{2d} + ϕ_{3d} for difference frequencies. The timedependence of Υ_{d} is described by
$$\begin{array}{c}\hfill \frac{\partial {\mathrm{{\rm Y}}}_{d}}{\partial {t}_{1}}=\delta {\omega}_{d}+cot\left({\mathrm{{\rm Y}}}_{d}\right)({\gamma}_{{\u229e}_{d}}+\frac{\partial ln{A}_{{1}_{d}}}{\partial {t}_{1}}+\frac{\partial ln{A}_{{2}_{d}}}{\partial {t}_{1}}+\frac{\partial ln{A}_{{3}_{d}}}{\partial {t}_{1}})\phantom{\rule{0.166667em}{0ex}},\end{array}$$(F.10)
because the coupling coefficient η_{1d} does not depend on time. In Eq. (F.10), γ_{⊞d} = γ_{1d} + γ_{2d} + γ_{3d}. The equation for the time dependence of Υ_{d} is of the exact same form as Eq. (42), which describes the time dependence of Υ. This suggests that the phase coordinates Υ and Υ_{d} are connected by some direct relation.
Finally, we prove the equivalence of the AEs derived in Sect. 2.6 and the AEs derived in this Appendix using a relabeling operation. The difference frequency Ω_{1d} ≈ Ω_{2d} − Ω_{3d} can be written as a sum frequency Ω_{2d} ≈ Ω_{1d} + Ω_{3d}. Swapping labels in the definition of the coupling coefficient η_{1d} defined in Eq. (F.5) to describe the coupling among modes in this sumfrequency analogue (dropping the subscript d in the process), leads to the equivalence of its expression with that of the sumfrequency coupling coefficient η_{1} defined in Eq. (22). Performing the same relabeling operation for the frequencydifference phase coordinate Υ_{d} defined in Eq. (F.9), we obtain that it is equal to −Υ. Substituting this relation between Υ_{d} and Υ in the differencefrequency AEs (F.6) converts them into the sumfrequency AEs (39a). The temporal dependence of the individual mode phases is therefore the same as those described by the AEs (39b). Hence, the AEs (39) and (42) describe the temporal evolution of sum and difference frequencies, and the properties of the modes in these resonances are derived in the main body of this article.
Appendix G: Materials available online
The MESA and GYRE inlists used to generate our models, as well as the final model data products themselves can be accessed using the MESA market place, https://cococubed.com/mesa_market/inlists.html , which contains a link to a Zenodo repository. This repository^{4} contains the respective inlists and additional information on the modeling process.
We developed a Python code that computes the stationary mode amplitude ratios using the formalism described in this work. This code generates the figures included in this work. You may download the latest version of the code from the GitHub repository at https://github.com/JVB11/AESolver. The online documentation may be found at https://jvb11.github.io/AESolver/, and contains examples on how to manipulate different parts of the code. The inlists for this Python code that generate the coupling models discussed in this work, as well as the final data products generated by this code, can be accessed using the same Zenodo repository that stores the MESA and GYRE inlists.
All Tables
Radial order ranges {n_{u}} (u ∈ ⟦3⟧) of the (k, m)=(0, 2) parent and (k, m)=(0, 1) daughter modes computed with GYRE and based on the MESA models listed in Table 1, angular rotation rates Ω expressed in percentages of the Roche critical rotation rate Ω_{Roche} and in d^{−1}, model radius R, and the order of magnitude of the bounds for the ranges of ϑ_{2} and ϑ_{3} (denoted by ϑ_{2, 3}), q_{1} and the quality factors Q_{φ} (φ ∈ ⟦3⟧).
Overview of the number of mode triads that satisfy different model validity and stability criteria.
Linear corotatingframe frequencies Ω_{φ}, largest absolute spin parameters s_{φ} (of the triad modes), radial orders n_{φ} and linear driving or linear damping rates γ_{φ} of g modes φ in all identified isolated (g) mode triads with ordering numbers (k_{1}, k_{2}, k_{3}) = (0, 0, 0) and azimuthal orders (m_{1}, m_{2}, m_{3}) = (2, 1, 1) for the models in Table 2 that yield stable and valid stationary solutions.
Dimensionless quantities for the mode triads listed in Table 4: dimensionless AE validity parameter Ψ_{AE}, absolute detuningdamping ratio q and absolute detuningdriving ratio q_{1}, the drivingdamping rate ratios ϑ_{2} and ϑ_{3}, as well as the order of magnitude of the energynormalized coupling coefficient, 𝒪_{η1}.
Maximal absolute nonlinear frequency shifts as a fraction of the linear inertialframe angular frequency $\mid \delta {\mathrm{\Omega}}_{\phi ,\mathrm{m}}^{s}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\Omega}}_{\phi ,\phantom{\rule{thinmathspace}{0ex}}i}\mid $ (in partsperthousand), zeropointcorrected nonlinear combination phases ${\stackrel{\sim}{\mathrm{\Phi}}}_{s,0}^{s}$, and expected stationary surface luminosity fluctuation ratios $\mid {\mathfrak{L}}_{2}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ and $\mid {\mathfrak{L}}_{3}^{s}\mid /\mid {\mathfrak{L}}_{1}^{s}\mid $ for the mode triads listed in Table 4.
All Figures
Fig. 1. Representations of isolated resonant threemode coupling networks for modes φ, β and γ. They represent direct resonances (i), parametric resonances (ii), and driven resonances (iii). A linearly unstable (i.e., excited) mode is pictured as an orange circle, whereas a linearly stable (i.e., damped) mode is pictured as a blue square. 

In the text 
Fig. 2. Stability domains as a function of ϑ_{2} and ϑ_{3}, indicated by the dark, shaded area, and evaluated using Eqs. (61a) and (64) for different values of q_{1}. Specifically, q_{1} is set equal to 0.1 (left top panel), 1 (middle top panel), 10 (right top panel), 10^{2} (left bottom panel), 10^{3} (middle bottom panel) and 10^{4} (right bottom panel). The stability domain determined by the stability conditions defined in Dziembowski (1982) is larger and is also indicated in these panels by the hatched areas. Unhatched white areas of the figure panels indicate unstable domains. Note the different axis scales for the different panels, which indicates the importance of the value of q_{1} in determining the stability of stationary solutions. 

In the text 
Fig. 3. Empirical probability distributions of δω / (γ_{2} + γ_{3}) for the 21 identified (stable and AE valid) isolated mode triads (blue and dotted) and their AE valid and stable counterparts (yellow) that do not fulfill the isolation criterion. We choose the bin width in a similar way to what is done in Hogg (2008). 

In the text 
Fig. 4. Coupling coefficient and squared BruntVäisälä frequency (N^{2}) profiles as a function of the fractional radius for the isolated mode triad with g mode radial orders (n_{1}, n_{2}, n_{3}) = (−22, −15, −53) in the midMS ΔX_{c, 1} model of Table 4. 

In the text 
Fig. 5. Same as Fig. 4, but for the isolated mode triad with g mode radial orders (n_{1}, n_{2}, n_{3}) = (−14, −26, −10) in the ΔX_{c, 2}_{(a)} model of Table 4. 

In the text 
Fig. 6. Observed daughterparent amplitude ratios ${A}_{d}^{\mathrm{\prime}}/{A}_{p}^{\mathrm{\prime}}$ of candidate resonances of SPB stars detected in models of light curves measured by Kepler, compared to their theoretically predicted equivalent observables, the monochromatic stationary daughterparent surface luminosity perturbation ratios 𝔏_{d} / 𝔏_{p}, which are computed with the theoretical oscillation model described in Sect. 2. These quantities are shown as a function of the minimal inertialframe frequency ν_{min} of the three signals in the considered candidate resonance, or, the minimal mode frequency observable ${\mathrm{\Omega}}_{\mathrm{min},\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}^{\mathrm{NL},\phantom{\rule{0.166667em}{0ex}}s}\phantom{\rule{0.166667em}{0ex}}/\phantom{\rule{0.166667em}{0ex}}2\pi $ among the three mode frequency observables of the considered triad, respectively. Left panel: ${A}_{d}^{\mathrm{\prime}}/{A}_{p}^{\mathrm{\prime}}$ as a function of ν_{min}. The light curve models were generated using strategy 3 of V21, which uses the signaltonoise ratio of the detected signals to determine their significance. Colors indicate the pseudoclass defined in V21: highf_{sv} (‘hfsv’), highmodedensity (‘hmd’), or outbursting (‘outb’). Right panel: 𝔏_{d} / 𝔏_{p} as a function of ${\mathrm{\Omega}}_{\mathrm{min},\phantom{\rule{0.166667em}{0ex}}\mathfrak{i}}^{\mathrm{NL},\phantom{\rule{0.166667em}{0ex}}s}/2\pi $. Colors indicate whether all stability criteria are fulfilled (‘Stab.’), or if additionally the AE validity criterion (‘Stab. + AE’) or all validity criteria (‘Stab. + Val.’; i.e., the isolated triads listed in Table 6) are fulfilled. We show alias signals in the frequency range of 0 to 24.4512 d^{−1} as fainter orangered (‘Stab.’) symbols, because we do not account for the splitting of the alias signals (which will affect the amplitude ratios). 

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