© The Authors 2024

1. Introduction

Slowly pulsating B (SPB) stars are mid-to-late-B variable stars on the main sequence that display a variety of low-frequency 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 multi-periodic 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 gravito-inertial 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ńska-Daszkiewicz 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 non-linear response of the stellar medium to the pulsation wave (see e.g., Bowman et al. 2016), which is referred to as non-linear distortion by Degroote et al. (2009). In this work, however, we focus on combination frequencies that are explained by non-linear coupling among oscillation modes (e.g., Buchler & Goupil 1984 and Van Hoolst 1994b). The amplitudes of heat-driven non-radial oscillations in SPB stars cannot be explained in the linear approximation. Observed amplitudes are therefore currently not used in asteroseismic inference. Non-linear 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 resource-costly numerical integration of the non-linear hydrodynamical equations that govern the oscillation dynamics, we consider weakly non-linear effects of the lowest order. We therefore consider isolated weak non-linear 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 non-linear 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 non-radial 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 non-rotating 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 non-linear three-mode 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 non-linear tides in multiple-star 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 star-exoplanet systems (Essick & Weinberg 2016; Vick et al. 2019; Yu et al. 2021b, 2022), non-linear 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), non-linear 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 so-called 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., Longuet-Higgins 1968; Lee & Saio 1997; Townsend 2003; Mathis 2013). For high-radial 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 Brunt-Väisälä frequency N in those deep near-core 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 low-frequency (Poincaré) modes¹. 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 non-linear 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 zero-age 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 non-linear three-mode 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 non-linear 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 co-rotating reference frame as (e.g., Frieman & Rotenberg 1960, Lynden-Bell & Ostriker 1967 or S01)

$$\begin{array}{c}\hfill \ddot{\mathit{\xi }}+\mathit{B}(\dot{\mathit{\xi }})+\mathit{C}(\mathit{\xi })={\mathit{a}}_{\mathrm{ext}},\end{array}$$(1)

where ξ denotes the Lagrangian displacement, the superscripted dot indicates a partial time derivative, $\mathit{B}(\dot{\mathit{\xi }})\equiv 2\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 anti-Hermitian and Hermitian, respectively (Lynden-Bell & 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})e^{-i\omega t},\end{array}$$(2)

we can rewrite Eq. (1) for free oscillations as

$$\begin{array}{c}\hfill -{\omega }^2\mathit{\xi }-i\omega \mathit{B}(\mathit{\xi })+\mathit{C}(\mathit{\xi })=0,\end{array}$$(3)

where ω is the (real- or complex-valued) angular frequency in the co-rotating 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 time-independent, 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.\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 }{\mathit{\xi }}^{\ast }·{\mathit{\xi }}^{\prime }\mathrm{d}V.\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 complex-valued ω_φ described by Eq. (3), the linear heat-driven growth (damping) rate γ_φ is defined as γ_φ ≡ Im[ω_φ]. Because of Ansatz (2), linear growth (damping) occurs when γ_φ is positive (negative).

We describe the coupling among non-degenerate 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 }{b}_{\phi },\end{array}$$(6)

where ${\delta }_{\beta }^{\phi }$ denotes a Kronecker delta, Re[ω_φ] = Ω_φ (and similar for mode β), and the real-valued 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{\mathrm{\Omega }}_{\phi }\langle {\mathit{\xi }}_{\phi },{\mathit{\xi }}_{\phi }\rangle .\end{array}$$(7)

The constant b_φ defined in Eq. (7) is related to the rotating-frame mode energy ϵ_φ at unit complex amplitude, if Ω_φ ≠ 0 (S01):

$$\begin{array}{c}\hfill {ϵ}_{\phi }={\mathrm{\Omega }}_{\phi }{b}_{\phi }.\end{array}$$(8)

Because we describe uniformly rotating non-magnetic 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 ,\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, Lynden-Bell & 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)H_r(\theta ),{\xi }_h^{\phi }(r)H_{\theta }(\theta ),i{\xi }_h^{\phi }(r)H_{\varphi }(\theta )]e^{i(m_{\phi }{\varphi }_{\phi }-{\mathrm{\Omega }}_{\phi }t)},\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 real-valued 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\pi {\displaystyle {\int }_{-1}^1}{|H_r(\theta )|}^2\mathrm{d}\mu =1,\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. Non-linear oscillations within the TAR

Non-linear mode interactions add acceleration (i.e., forcing) terms to the governing equation of motion, so that the higher-order equation of motion becomes the following quadratic eigenfunction problem (see e.g., Frieman & Rotenberg 1960)

$$\begin{array}{c}\hfill -{\omega }^2\mathit{\xi }-i\omega \mathit{B}\left(\mathit{\xi }\right)+\mathit{C}\left(\mathit{\xi }\right)=\mathit{a}[\mathit{\xi }],\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 non-linear coupling term a[ξ], which can be split up into contributions from pressure gradients (a_P[ξ]) and gravity (a_G[ξ]).

Because we limit ourselves to non-linear three-mode 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⁽²⁾[ξ, ξ′], 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,n}^{(2)}[\mathit{\xi },{\mathit{\xi }}^{\prime }]=& -\frac{1}{2}{\xi }^k{\xi }^l{\mathrm{\nabla }}_k{\mathrm{\nabla }}_l{\mathrm{\nabla }}_n\mathrm{\Phi },\hfill \end{array}$$(14a)

$$\begin{array}{cc}\hfill {\mathit{a}}_{P,n}^{(2)}[\mathit{\xi },{\mathit{\xi }}^{\prime }]=& -\frac{1}{\rho }{\mathrm{\nabla }}_j[p({\mathrm{\Gamma }}_1-1){\mathrm{\Theta }}_n^j+p{\mathrm{\Xi }}_n^j+\mathrm{\Psi }{\delta }_n^j],\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 Γ₁ = (∂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=& \left({\mathrm{\nabla }}_i{\xi }^j\right)\left({\mathrm{\nabla }}_k{\xi }^k\right)=\left({\mathrm{\nabla }}_i{\xi }^j\right)(\mathrm{\nabla }·\mathit{\xi }),\hfill \end{array}$$(15a)

$$\begin{array}{cc}\hfill {\mathrm{\Xi }}_i^j=& \left({\mathrm{\nabla }}_k{\xi }^j\right)\left({\mathrm{\nabla }}_i{\xi }^k\right),\hfill \end{array}$$(15b)

and of the quantities

$$\begin{array}{cc}& \mathrm{\Psi }=\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)\},\hfill \end{array}$$(16a)

$$\begin{array}{cc}& \mathrm{\Theta }={\delta }_j^i{\mathrm{\Theta }}_i^j={(\mathrm{\nabla }·\mathit{\xi })}^2,\hfill \end{array}$$(16b)

$$\begin{array}{cc}& \mathrm{\Xi }={\delta }_j^i{\mathrm{\Xi }}_i^j=\left({\mathrm{\nabla }}_k{\xi }^i\right)\left({\mathrm{\nabla }}_i{\xi }^k\right).\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. Coupled-mode 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 coupled-mode 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{\mathrm{\Omega }}_{\phi }{\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{\mathrm{\Omega }}_{\phi }{\mathit{\xi }}_{\phi }^{\ast }\left(\mathit{x}\right)\end{array}\right],\end{array}$$(17)

under the assumption that mode φ is not degenerate and has a real-valued 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 real-valuedness of the eigenfrequencies further justifies the use of eigenmode orthogonality relation (6). In principle, the expansion (17) should also include modes with non-ordinary 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 r-modes in S01, we assume that such modes with non-ordinary 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 non-ordinary 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)=& \frac{1}{b_{\phi }}\langle {\mathit{\xi }}_{\phi },{\mathrm{\Omega }}_{\phi }\mathit{\xi }(t)+i\dot{\mathit{\xi }}(t)+i\mathit{B}(\mathit{\xi }(t))\rangle ,\hfill \end{array}$$(18a)

$$\begin{array}{cc}\hfill c_{\phi }^{\ast }(t)=& \frac{1}{b_{\phi }}\langle {\mathit{\xi }}_{\phi }^{\ast },{\mathrm{\Omega }}_{\phi }\mathit{\xi }(t)-i\dot{\mathit{\xi }}(t)-i\mathit{B}(\mathit{\xi }(t))\rangle ,\hfill \end{array}$$(18b)

for non-zero b_φ and Ω_φ ≠ −Ω_β, where the subscript β refers to a different non-degenerate mode in expansion (17). The rotating-frame mode energy E_φ for mode φ with rotating-frame frequency Ω_φ and (complex) amplitude c_φ is (see Appendix K in S01)

$$\begin{array}{c}\hfill E_{\phi }={ϵ}_{\phi }|c_{\phi }{|}^2={\mathrm{\Omega }}_{\phi }{b}_{\phi }{|c_{\phi }|}^2.\end{array}$$(19)

We disregard non-linear 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 inter-g mode coupling. For faster-rotating 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 inter-g mode couplings. As discussed by Buchler et al. (1995), stable (constant), periodic or irregular time-dependent behavior of the oscillation amplitudes (and frequencies) may be observed. In this work, we focus on the stable (time-independent) 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².

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 coupled-mode equations in the quadratic approximation (see Appendix A of S01 for the technical details),

$$\begin{array}{c}\hfill {\dot{c}}_{\phi }+i{\mathrm{\Omega }}_{\phi }{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 }),\end{array}$$(20)

in which we employ the Einstein summation convention. If we combine Eqs. (8) and (20) and assume non-zero Ω_φ, we get

$$\begin{array}{c}\hfill \frac{{\dot{c}}_{\phi }}{i}+{\mathrm{\Omega }}_{\phi }{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 }).\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 energy-scaled 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 }={ϵ}_{\phi }{\eta }_{\phi }^{\beta \gamma }=\langle {\mathit{\xi }}_{\phi },{\mathit{a}}^{(2)}[{\mathit{\xi }}_{\beta },{\mathit{\xi }}_{\gamma }]\rangle ,\end{array}$$(22)

which deviates from the definition in S01: ξ_φ is used instead of its complex conjugate. The (real-valued) total displacement is thus expanded in terms of non-conjugated 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⁽²⁾[ξ_β,  ξ_γ]. We normalize the eigenfunctions ξ_φ, ξ_β and ξ_γ such that ϵ_φ = ϵ_β = ϵ_γ = GM²/R, following L12. Additional information on the explicit expressions used in Eq. (22) can be found in Appendix C.

The coupled-mode Eqs. (20) and (21) thus describe the temporal dynamics of the amplitudes of all (coupled) modes in a coupling network. To study specific three-mode non-linear interactions among g modes in SPB stars, we derive amplitude equations (AEs) in Sect. 2.6, based on the coupled-mode 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 non-zero, 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\pi }}exp\left(i(-H_{\phi }{m}_{\phi }+H_{\beta }{m}_{\beta }+H_{\gamma }{m}_{\gamma })\varphi \right)d\varphi .\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 }{m}_{\phi }+H_{\beta }{m}_{\beta }+H_{\gamma }{m}_{\gamma }=0,\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(|k_{\phi }|+|k_{\beta }|+|k_{\gamma }|,2\right)=0,\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 three-mode 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.

In addition to these three isolated resonant coupling networks, non-isolated 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 multiple-mode or multi-mode 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 great-granddaughter modes have also been constructed. The reader is referred to Kumar & Goodman (1996) and Weinberg et al. (2021) for additional information on such networks.

Higher-order coupling terms, such as those that appeared in the cubic four-wave coupling networks developed to describe four interacting modes in non-rotating stars (by for example Van Hoolst & Smeyers 1993 and Van Hoolst 1994a) do not adhere to the selection rules in Sect. 2.4. The four-wave selection rules can however be determined from similar arguments. For example, the four-wave 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 multi-mode or higher-order coupling terms are necessary to explain the (stationary) amplitudes and phases of the three-mode 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 three-mode sum-frequency coupling networks that produce stable stationary solutions, if those solutions effectively describe observed (time-independent) amplitudes, phases and frequencies. We limit ourselves in the remainder of this work to describing such coupling networks.

2.6. Amplitude equations for Ω₁ ≃ Ω₂ + Ω₃

Under the assumption that the linear driving or linear damping rate γ_φ of a mode φ fulfills |γ_φ/Ω_φ|≪1 for the modes that are weakly non-linearly coupled, we can derive the AEs. In this section we use the coupled-mode Eqs. (20) and (21) valid for quadratically non-linear oscillations in rotating stars to derive the specific AEs for the sum-frequency resonant interaction Ω₁ ≃ Ω₂  +  Ω₃ 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 difference-frequency 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 non-linear coupling at the quadratic level. We introduce time variables

$$\begin{array}{c}\hfill t_0=t,t_1=\mathfrak{J}t,\dots ,\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}\frac{\partial }{\partial t_1}+O({\mathfrak{J}}^2).\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}{c}_{\phi 1}(t_0,t_1,\dots )+{\mathfrak{J}}^2{c}_{\phi 2}(t_0,t_1,\dots )+O({\mathfrak{J}}^3).\end{array}$$(28)

Perturbation series (28) needs to be uniform: each of the higher-order terms should be a small correction to lower-order terms.

Substituting expansions (27) and (28) into the coupled-mode equations (20), subsequently dividing by the non-zero 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{\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 & +O({\mathfrak{J}}^2).\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 𝔍⁰ is

$$\begin{array}{c}\hfill \frac{\partial c_{\phi 1}}{\partial t_0}+i{\mathrm{\Omega }}_{\phi }{c}_{\phi 1}\equiv L_{\phi }(c_{\phi 1})=0,\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{\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{\mathrm{\Omega }}_{\phi }{t}_0},\end{array}$$(31)

where a_φ is a complex amplitude factor that varies only on time scales slower than the time scale t₀ 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 & +{\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 \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}{e}^{-i{\mathrm{\Omega }}_{\phi }{t}_0}\hfill \\ \hfill & +\frac{i}{b_{\phi }}(a_{\beta }{a}_{\gamma }{\kappa }_{\phi }^{\beta \gamma }{e}^{-i({\mathrm{\Omega }}_{\beta }+{\mathrm{\Omega }}_{\gamma })t_0}+a_{\beta }^{\ast }{a}_{\gamma }{\kappa }_{\phi }^{\overline{\beta }\gamma }{e}^{-i(-{\mathrm{\Omega }}_{\beta }+{\mathrm{\Omega }}_{\gamma })t_0}\hfill \\ \hfill & +a_{\beta }{a}_{\gamma }^{\ast }{\kappa }_{\phi }^{\beta \overline{\gamma }}{e}^{-i({\mathrm{\Omega }}_{\beta }-{\mathrm{\Omega }}_{\gamma })t_0}+a_{\beta }^{\ast }{a}_{\gamma }^{\ast }{\kappa }_{\phi }^{\overline{\beta }\overline{\gamma }}{e}^{i({\mathrm{\Omega }}_{\beta }+{\mathrm{\Omega }}_{\gamma })t_0}).\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 secular-term-generating terms have a multiplying factor exp(−iΩ_φt₀) 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 Ω₁ ≃ Ω₂ + Ω₃, defined in terms of the detuning parameters δω and ΔΩ^l,

$$\begin{array}{c}\hfill \mathfrak{J}\delta \omega ={\mathrm{\Omega }}_1-{\mathrm{\Omega }}_2-{\mathrm{\Omega }}_3\equiv \mathrm{\Delta }{\mathrm{\Omega }}^l,\end{array}$$(34)

in Eq. (33). The condition that ensures that the sum of these secular-term-generating terms vanishes leads to the (extended complex) AEs for the sum-frequency resonance Ω₁ ≃ Ω₂ + Ω₃:

$$\begin{array}{c}\hfill \frac{\partial \mathit{a}}{\partial t_1}=\mathit{\gamma }°\mathit{a}+2i(\mathit{\eta }°{\mathit{a}}_M°{\mathbf{\Omega }}_M°\mathit{e}),\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 Hadamard-Schur (i.e., element-wise, 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],{\mathit{a}}_M=\left[\begin{array}{c}{a}_2{a}_3\\ a_1{a}_3^{\ast }\\ a_1{a}_2^{\ast }\end{array}\right],\mathit{e}=\left[\begin{array}{c}exp(i\delta \omega t_1)\\ exp(-i\delta \omega t_1)\\ exp(-i\delta \omega t_1)\end{array}\right],\end{array}$$(36)

and

$$\begin{array}{c}\hfill \mathit{\gamma }=\left[\begin{array}{c}{\gamma }_1\\ {\gamma }_2\\ {\gamma }_3\end{array}\right],\mathit{\eta }=\left[\begin{array}{c}{\eta }_1\\ {\eta }_1^{\ast }\\ {\eta }_1^{\ast }\end{array}\right],{\mathbf{\Omega }}_M=\left[\begin{array}{c}{\mathrm{\Omega }}_1\\ {\mathrm{\Omega }}_2\\ {\mathrm{\Omega }}_3\end{array}\right],\end{array}$$(37)

in Eq. (35). In Eq. (37), we define the isolated three-mode coupling coefficient η₁:

$$\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,\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 non-linear energy transfer is expected, and one enters the regime of intermediate to strong non-linear 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∈ℕ₀ | 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 }°\mathit{A}+2|{\eta }_1|sin\left(\mathrm{{\rm Y}}\right)\left({\mathit{A}}_N^{°}{\mathbf{\Omega }}_M\right),\hfill \end{array}$$(39a)

$$\begin{array}{cc}\hfill \mathit{A}°\frac{\partial \mathit{\varphi }}{\partial t_1}& =2|{\eta }_1|cos\left(\mathrm{{\rm Y}}\right)({\mathit{A}}_P°{\mathbf{\Omega }}_M),\hfill \end{array}$$(39b)

where η₁ = |η₁| e^−i δ₁, and

$$\begin{array}{c}\hfill \mathit{A}=\left[\begin{array}{c}{A}_1\\ A_2\\ A_3\end{array}\right],\mathit{\varphi }=\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ {\varphi }_3\end{array}\right],{\mathit{A}}_N=\left[\begin{array}{c}{A}_2{A}_3\\ -A_1{A}_3\\ -A_1{A}_2\end{array}\right],{\mathit{A}}_P=\left[\begin{array}{c}{A}_2{A}_3\\ A_1{A}_3\\ A_1{A}_2\end{array}\right].\end{array}$$(40)

In the AEs (39) we also introduce the generic phase coordinate

$$\begin{array}{c}\hfill \mathrm{{\rm Y}}\equiv -\delta \omega t_1+{\varphi }_1-{\varphi }_2-{\varphi }_3+{\delta }_1,\end{array}$$(41)

which contains the combination phase Φ = ϕ₁ − ϕ₂ − ϕ₃. 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 }_{⊞}+\frac{\partial ln{A}_1}{\partial t_1}+\frac{\partial ln{A}_2}{\partial t_1}+\frac{\partial ln{A}_3}{\partial t_1}),\end{array}$$(42)

for non-zero resonant mode amplitudes, because of Eq. (39a) and the definition γ_⊞ ≡ γ₁ + γ₂ + γ₃. Equations (39a) and (42) thus form an autonomous four-dimensional system equivalent to the six-dimensional system in Eqs. (39a) and (39b), in which the individual mode phases ϕ₁, ϕ₂ and ϕ₃ do not explicitly appear.

Non-linear interactions lead to non-linear 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{t}_0+2{\displaystyle \int ({\mathit{A}}_P\oslash \mathit{A}°{\mathbf{\Omega }}_M)}|{\eta }_1|cos\left(\mathrm{{\rm Y}}\right)d{t}_1)\right\},\hfill \end{array}\end{array}$$(43)

in which the symbol ‘⊘’ denotes a Hadamard-Schur (element-wise) 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],{\mathit{c}}_1=\left[\begin{array}{c}{c}_{11}\\ c_{21}\\ c_{31}\\ \end{array}\right].\end{array}$$(44)

The third term in the exponential factor of Eq. (43) defines the quadratic non-linear frequency shift multiplied with the elapsed time. If the linear angular mode frequency in the co-rotating frame – the second term in that factor – is added to the non-linear frequency shift, we obtain a shifted frequency, hereafter referred to as the non-linear frequency, in the same reference frame. The expression for the (second-order) non-linear correction to the complex amplitude, c_φ2, can be found in Appendix B.1 of Van Beeck (2023). The harmonic resonance equivalents of the non-linear frequency shifts and second-order non-linear 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 three-mode couplings identified by V21 are time-independent, with the modes in a non-linear frequency (and phase) lock (i.e., they are synchronized). We therefore derive the time-independent 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 }°{\mathit{A}}^s=-2|{\eta }_1|sin\left({\mathrm{{\rm Y}}}^s\right)({\mathit{A}}_N^s°{\mathbf{\Omega }}_M),\end{array}$$(45)

in which we define the stationary phase coordinate Υ^s as

$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}^s=-\delta \omega t_1+{\varphi }_1^s-{\varphi }_2^s-{\varphi }_3^s+{\delta }_1,\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){\gamma }_{⊞},\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)}^{°2}=\frac{\mathit{Q}}{4|{\eta }_1{|}^2}[1+q^2],\end{array}$$(48)

where the superscript ‘°n’ indicates a Hadamard-Schur (element-wise) nth power. In Eq. (48), we use the detuning-damping ratio q and quality factor vector Q, which we define as

$$\begin{array}{c}\hfill \mathit{Q}=\left[\begin{array}{c}{\displaystyle {\gamma }_2{\gamma }_3/{\mathrm{\Omega }}_2{\mathrm{\Omega }}_3}\\ {\displaystyle -{\gamma }_1{\gamma }_3/{\mathrm{\Omega }}_1{\mathrm{\Omega }}_3}\\ {\displaystyle -{\gamma }_1{\gamma }_2/{\mathrm{\Omega }}_1{\mathrm{\Omega }}_2}\\ \end{array}\right]=\left[\begin{array}{c}{\displaystyle 1/Q_2{Q}_3}\\ {\displaystyle -1/Q_1{Q}_3}\\ {\displaystyle -1/Q_1{Q}_2}\end{array}\right],q\equiv \frac{\delta \omega }{{\gamma }_{⊞}}=-cot{\mathrm{{\rm Y}}}^s,\end{array}$$(49)

with quality factor Q_φ ≡ Ω_φ/γ_φ for φ ∈ ⟦3⟧. The resonant stationary amplitudes are real if

$$\begin{array}{c}\hfill {\mathbf{\Omega }}_N{>}_{°}0\vee {\mathbf{\Omega }}_N{<}_{°}0\Rightarrow {\mathit{A}}^s{\in }_{°}\mathbb{R},\end{array}$$(50)

for a driven resonance scenario, or if

$$\begin{array}{c}\hfill {\mathbf{\Omega }}_P{>}_{°}0\vee {\mathbf{\Omega }}_P{<}_{°}0\Rightarrow {\mathit{A}}^s{\in }_{°}\mathbb{R},\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 element-wise 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 non-linear interaction described by coupled-mode 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 non-linear coupling coefficients |η₁|, the efficiency for non-linear 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 non-linear effect leading to amplitude saturation. Quality factor Q_φ expresses the ratio of the mode φ’s e-folding 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²/R) to damp or grow by a factor of e due to linear heat-driven 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 non-linearly excite (i.e., a larger parent mode energy is required to have a measurable non-linear effect).

Increasing values of |q| lead to increasing values of the stationary amplitudes. Two factors affect the value of the detuning-damping ratio q: the detuning δω and γ_⊞ (≡γ₁ + γ₂ + γ₃). For larger detunings the stationary amplitudes increase, because the increased detuning reduces the efficiency of the non-linear mode coupling (similar to decreasing |η₁|), requiring a larger amplitude for a non-linear effect. For smaller values of |γ_⊞|, that is, for the case where the parent and daughter mode linear heat-driven 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 non-linear 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},{\left(\frac{A_3^s}{A_1^s}\right)}^2=-\frac{Q_3}{Q_1},{\left(\frac{A_3^s}{A_2^s}\right)}^2=\frac{Q_3}{Q_2}.\end{array}$$(52)

The energy ratios (52) do not depend on the non-linear 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 non-linear energy transfer from the linearly excited mode, because it loses energy more slowly.

Non-linear 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,\end{array}$$(53)

with the non-linear 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}q\left[\begin{array}{c}-{\gamma }_1\\ {\gamma }_2\\ {\gamma }_3\\ \end{array}\right]=\mathrm{\Delta }{\mathrm{\Omega }}^l\left[\begin{array}{c}-{\gamma }_1/{\gamma }_{⊞}\\ {\gamma }_2/{\gamma }_{⊞}\\ {\gamma }_3/{\gamma }_{⊞}\\ \end{array}\right],\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 non-linear stationary frequencies ${\mathrm{\Omega }}_{\phi }^{nl,s}$ are frequency-locked, 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,\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 difference-frequency 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 Hartman-Grobman 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 non-linear coupling coefficients because the hyperbolicity condition considers the geometry of the non-linear 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+\delta A_{\phi })exp\left[i({\varphi }_{\phi }^s+\delta {\varphi }_{\phi })\right],\forall \phi \in [[3]].\end{array}$$(56)

Substituting the perturbed complex amplitude factor (56) into the complex-valued AEs (35), subsequently linearizing and separating the real-valued and complex-valued 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}\mathsf{Z}.\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{\gamma }_1\\ -{\gamma }_2& {\gamma }_2& -{\gamma }_2& q{\gamma }_2\\ -{\gamma }_3& -{\gamma }_3& {\gamma }_3& q{\gamma }_3\\ q{\gamma }_{{⊞}_1}& q{\gamma }_{{⊞}_2}& q{\gamma }_{{⊞}_3}& {\gamma }_{⊞}\\ \end{array}\right),\end{array}$$(58)

and Z is the perturbation tensor defined as

$$\begin{array}{c}\hfill \mathsf{Z}={(\delta A_1/A_1^s,\delta A_2/A_2^s,\delta A_3/A_3^s,\delta \mathrm{{\rm Y}})}^{\mathrm{T}},\end{array}$$(59)

for which

$$\begin{array}{c}\hfill \delta \mathrm{{\rm Y}}=\delta {\varphi }_1-\delta {\varphi }_2-\delta {\varphi }_3,{\gamma }_{{⊞}_{\phi }}\equiv {\gamma }_{⊞}-2{\gamma }_{\phi },\forall \phi \in [[3]].\end{array}$$(60)

We infer an analytical stability domain of the stationary solutions of the AEs based on the Routh-Hurwitz 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=& w_1=-2{\gamma }_{⊞}>0\iff {\gamma }_{⊞}<0,\hfill \end{array}$$(61a)

$$\begin{array}{cc}\hfill H_2=& w_1{w}_2-w_3>0,\hfill \end{array}$$(61b)

$$\begin{array}{cc}\hfill H_3=& w_3{H}_2-w_1^2{w}_4>0,\hfill \end{array}$$(61c)

$$\begin{array}{cc}\hfill H_4=& w_4{H}_3>0,\hfill \end{array}$$(61d)

in which the coefficients w_u  (u ∈ ⟦4⟧₀, with ⟦n⟧₀ ≡ {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}_{4-u}{\lambda }^u=w_4+w_3\lambda +w_2{\lambda }^2+w_1{\lambda }^3+w_0{\lambda }^4=0,\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=& 1,\hfill \end{array}$$(63a)

$$\begin{array}{cc}\hfill w_1=& -2{\gamma }_{⊞},\hfill \end{array}$$(63b)

$$\begin{array}{cc}\hfill w_2=& {\gamma }_{⊞}^2[1+q^2]-4q^2[{\gamma }_1{\gamma }_2+{\gamma }_1{\gamma }_3+{\gamma }_2{\gamma }_3],\hfill \end{array}$$(63c)

$$\begin{array}{cc}\hfill w_3=& 4{\gamma }_1{\gamma }_2{\gamma }_3[1+3q^2],\hfill \end{array}$$(63d)

$$\begin{array}{cc}\hfill w_4=& -4{\gamma }_{⊞}{\gamma }_1{\gamma }_2{\gamma }_3[1+q^2].\hfill \end{array}$$(63e)

These are the same coefficients as the ones that were determined by Dziembowski (1982) in his formalism for quadratic non-linear mode coupling of three distinct modes in non-rotating stars. Because H₃ >  0 due to the stability condition (61c), stability condition (61d) becomes w₄ >  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 third-order couplings that can limit the amplitude. The direct resonant three-mode 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 time-variable 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 three-mode 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).

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 three-mode 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 }_{⊞}}[{\mathfrak{d}}_4+{\mathfrak{d}}_2{q}^2+{\mathfrak{d}}_0{q}^4]>0\iff {\mathfrak{d}}_4+{\mathfrak{d}}_2{q}^2+{\mathfrak{d}}_0{q}^4>0.\end{array}$$(64)

The factor $-{\gamma }_1^3/{\gamma }_{⊞}$ 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=& -18{\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)],\hfill \end{array}\hfill \end{array}$$(65a)

$$\begin{array}{cc}& \begin{array}{cc}\hfill {\mathfrak{d}}_2=& -12{\vartheta }_2{\vartheta }_3-(1-{\vartheta }_2-{\vartheta }_3)\hfill \\ \hfill & [2+{({\vartheta }_2-{\vartheta }_3)}^2+{({\vartheta }_2+{\vartheta }_3)}^2],\hfill \end{array}\hfill \end{array}$$(65b)

$$\begin{array}{cc}& {\mathfrak{d}}_4=-2{\vartheta }_2{\vartheta }_3+{(1-{\vartheta }_2-{\vartheta }_3)}^3,\hfill \end{array}$$(65c)

where we define ϑ_2, 3 ≡ −γ_2, 3 / γ₁, 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 ϑ₂ and ϑ₃, in accordance with the symmetry of the coupling coefficient. Coefficients 𝔡₂ and 𝔡₄ are the same as in Eq. (6.14) of Dziembowski (1982). Although the coefficient 𝔡₀ 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 ϑ₂ and ϑ₃ only, whereas q can be expanded as a function of ϑ₂, ϑ₃, and the ratio of the linear frequency detuning to the parent’s linear driving rate q₁ ≡ δω/γ₁. We therefore can explore the stability domain of the (hyperbolic) stationary solutions by varying only three dimensionless ratios of linear variables: q₁ = δω/γ₁, ϑ₂ and ϑ₃. The dimensionless ratio q₁ can be written as

$$\begin{array}{c}\hfill q_1\equiv Q_1\left(\frac{\delta \omega }{{\mathrm{\Omega }}_1}\right)=\frac{Q_1\delta \omega }{2\pi }{P}_1,\end{array}$$(66)

where the cyclic co-rotating frame period P₁ = 2 π / Ω₁. It is therefore a period-weighted combined measure of the driving time scale of the linearly excited parent mode (Q₁) and the efficiency with which non-linear energy transfer occurs (δω). Hence, when comparing values of q₁ for triads with linearly excited modes of similar period, one might expect that the triad with larger |q₁| has larger stationary mode amplitudes because of decreased efficiency of non-linear energy transfer (larger δω) and/or smaller stationary mode energy ratios $(A_{2,3}^s/A_1^s{)}^2$ (larger Q₁, when values of Q_2,  3 are comparable).

Figure 2 displays the domains of stability of stationary solutions in the three-dimensional phase space of the parameters q₁, ϑ₂ and ϑ₃, 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 ϑ₂ + ϑ₃ >  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 − ϑ₂ − ϑ₃) in that stability condition is always negative due to Eq. (61a) and can be interpreted as an effective total-damping-to-linear-driving ratio γ_⊞/γ₁. 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 𝔡₀. Condition (64) can also be expressed in terms of q₁, ϑ₂ and ϑ₃. 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 ϑ₂ and ϑ₃ 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 ϑ₂ − ϑ₃ plane moves towards larger values of ϑ₂ and ϑ₃ for larger values of q₁, as shown in the different panels of Fig. 2. Conversely, the stability domain decreases considerably in size for smaller values of q₁. 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 non-linear 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 higher-order non-linear mode coupling terms, such as the cubic self-coupling terms that were described in Van Hoolst (1996). Alternatively, limit cycles characterized by non-stationary 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₁/2) (for k ∈ {2 , 3}, following Dziembowski 1982) the explicit time-dependence disappears. Further assuming that the complex amplitude factor a₁ 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\frac{\delta \omega }{2})S_2+2i{\eta }^{\ast }{a}_1{S}_3^{\ast }{\mathrm{\Omega }}_2,\hfill \end{array}$$(67a)

$$\begin{array}{cc}\hfill \frac{\partial S_3^{\ast }}{\partial t_1}& =({\gamma }_3-i\frac{\delta \omega }{2})S_3^{\ast }-2i\eta a_1^{\ast }{S}_2{\mathrm{\Omega }}_3.\hfill \end{array}$$(67b)

Under the assumption that S_k ∼ exp(σ t₁) (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\delta \omega )}^2+16{|{\eta }_1|}^2{A}_1^2{\mathrm{\Omega }}_2{\mathrm{\Omega }}_3}}{2},\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₂ and A₃. 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₁ 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)}.\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|{\eta }_1|\sqrt{Q_2{Q}_3}}\sqrt{1+{\left(\frac{\delta \omega }{{\gamma }_2+{\gamma }_3}\right)}^2}.\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₁ ≥ A_t, and $A_1^s$ is always larger than A_t. In the limit of very small (non-zero) detuning δω, the threshold amplitude is solely dependent on the coupling coefficient and the quality factors. In that case, A_t increases with decreasing |η₁| 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 non-linear energy transfer, thus requiring a larger parent mode energy for a visible non-linear effect. A larger detuning increases the threshold amplitude because of the less efficient energy transfer.

3. Theoretically predicted observables

An important observable in linear g-mode asteroseismic modeling is a g-mode 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⁻⁷ d⁻¹, based on long-term ground-based 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: model-generated 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 disc-averaged 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}}{H}_r(i_s;k),\end{array}$$(71)

within the TAR, where

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{\Delta }{L}_{\phi }}{L}\right)}_{\mathrm{mode}}=t_k\frac{\mathrm{\Delta }{L}_{\phi ,R}(R)}{L(R)}-e_k\frac{{\xi }_{\phi ,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 limb-darkening 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}\mu H_r(\mu ;k)h(\mu )\mathrm{d}\mu ,\hfill \end{array}$$(73a)

$$\begin{array}{cc}\hfill e_k& ={\displaystyle {\int }_0^1}[2{\mu }^2\frac{\mathrm{d}{H}_r(\mu ;k)}{\mathrm{d}\mu }-(\mu -{\mu }^3)\frac{{\mathrm{d}}^2{H}_r(\mu ;k)}{\mathrm{d}{\mu }^2}]h(\mu )\mathrm{d}\mu ,\hfill \end{array}$$(73b)

where μ ≡ cosθ and h(μ) is a limb-darkening 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 limb-darkening 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 limb-darkening laws by changing the limb-darkening 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 disc-averaged 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,\phi }=\left|{\left(\frac{\mathrm{\Delta }{L}_{\phi }}{L}\right)}_{\mathrm{mode}}\right|H_r(i_s;k).\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}{A}_t,|{\mathfrak{L}}_{\phi }^s|=o_{A,\phi }{A}_{\phi }^s.\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 limb-darkening 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 daughter-parent 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 limb-darkening 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,2}}{o_{A,1}}\right)\sqrt{\frac{|Q_2|}{|Q_1|}},\frac{|{\mathfrak{L}}_3^s|}{|{\mathfrak{L}}_1^s|}=\left(\frac{o_{A,3}}{o_{A,1}}\right)\sqrt{\frac{|Q_3|}{|Q_1|}}.\end{array}$$(76)

The daughter-parent surface luminosity fluctuation ratios (76) are the most robust amplitude-based theoretically predicted observables that can be used in resonant non-linear asteroseismic modeling. To derive the expressions for these ratios, we assume a parametric resonant mode triad (for a three-mode sum-frequency 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 non-adiabatic parameters. The equivalent expression for the daughter-parent 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 multi-color 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 daughter-parent ratios considered here. Such integrations will not be considered in this work but will be considered in follow-up application papers.

3.2. Frequencies and phases: frequency detuning and combination phase for Ω₁ ≈ Ω₂ + Ω₃

The inherent assumption made in linear asteroseismic inference is that any non-linear 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 non-linear 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,\mathfrak{i}}-{\stackrel{\sim }{\mathrm{\Omega }}}_{2,\mathfrak{i}}-{\stackrel{\sim }{\mathrm{\Omega }}}_{3,\mathfrak{i}},\end{array}$$(77)

where the observed frequencies of modes φ are defined as ${\stackrel{\sim }{\mathrm{\Omega }}}_{\phi ,\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 co-rotating 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,\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\pi {\mathfrak{R}}_{\nu },\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 \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 ,\mathfrak{i}}^{\mathrm{NL},\mathrm{s}}\equiv {\mathrm{\Omega }}_{\phi }^s+m_{\phi }\mathrm{\Omega }+\delta {\mathrm{\Omega }}_{\phi }^s={\mathrm{\Omega }}_{\phi ,\mathfrak{i}}^s+\delta {\mathrm{\Omega }}_{\phi }^s,\end{array}$$(80)

which includes the quadratic stationary non-linear 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]-\left[\begin{array}{c}{\varphi }_{L1}\\ {\varphi }_{L2}\\ {\varphi }_{L3}\\ \end{array}\right],\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,\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 }_{L1}-{\varphi }_{L2}-{\varphi }_{L3}$. The coupling coefficient η₁ is real-valued 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 three-g-mode 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 B-type 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}}}=\frac{X_{\ast }}{Y_{\ast }}$, with X_* and Y_* equal to the Galactic standard values for B-type 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 proton-proton 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).

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 core-boundary mixing (CBM) processes, we employ an approach similar to the diffusive exponential overshooting model with efficiency parameters f_CBM and f₀ 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² s⁻¹. 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 mid-MS and near-ZAMS 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_⊙ mid-MS 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 non-linear quadratic mode coupling coefficients (22), we need the adiabatic derivative of Γ₁ 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 non-adiabatic 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 non-adiabatic linear eigenfrequencies computed by GYRE.

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 non-adiabatic GYRE computations – which use trial frequencies based on the adiabatic GYRE computation results – do not show evidence of missing non-adiabatic mode radial orders (see e.g., Goldstein & Townsend 2020), so that there is a one-to-one mapping between the adiabatic and non-adiabatic 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 Ω₁ ≈ Ω₂ + Ω₃ (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 γ₁ >  0, and restrict the computations of potential daughter modes to those (k, m)=(0, 1) modes with linear damping rates γ₂ <  0 and γ₃ <  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_⊙ near-ZAMS and 8 M_⊙ mid-MS and near-ZAMS models with a rotation rate of 20% of Ω_Roche, we obtain no linearly excited potential parent modes.

We compute mode coupling coefficients η₁ for mode triads whose stationary AE solutions satisfy the (linear) stability conditions derived in Sect. 2.8. To guarantee non-linear 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 non-linear mode coupling computations, because the computed stationary Jacobian matrix contains the coupling coefficient |η₁| (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 non-resonant 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 multi-mode fixed point is reached with the same modes for resonant and non-resonant AEs (Buchler et al. 1997). If that were not the case, the non-resonant 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 non-linear couplings considered in this work, stable multi-mode 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 co-rotating-frame 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,{\mathrm{\Omega }}_2,{\mathrm{\Omega }}_3)|},\end{array}$$(83)

which compares |ΔΩ^l| with the minimum absolute value of the co-rotating-frame angular mode frequencies of the triad. We require that mode triads considered for modeling have Ψ_AE ≤ 0.1, so that |ΔΩ^l| ≪ |min(Ω₁, Ω₂, Ω₃)|. Many mode triads satisfy this model validity criterion, as is shown in Table 3.

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 second-lowest |𝔏_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 second-lowest |𝔏_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 multi-mode 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.

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 |δω / (γ₂+γ₃)| to be small. Figure 3 shows the empirical probability distribution of |δω / (γ₂+γ₃)| for the 21 identified (stable, AE-valid) isolated mode triads, and their AE-valid and stable counterparts that do not fulfill the isolation criterion. While it is certainly true that |δω / (γ₂+γ₃)| is small for the isolated mode triads, there are non-isolated mode triads that have similarly low values of |δω / (γ₂+γ₃)|, which re-iterates 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_⊙ near-ZAMS 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 η₁(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 η₁ ≡ η₁(R)). The profile η₁(r) then indicates the zones in the stellar interior for which contributions to the coupling coefficient are significant for these mode triads.

In the near-core regions the squared Brunt-Väisälä frequency (N²) 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² which changes the eigenfunctions: their modal inertia becomes more confined to the near-core 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 TAR-valid near-core mixing zone with model age. The width of that maximum depends on the evolutionary stage: more evolved SPB models have wider N² maxima, and some modes become trapped in the near-core 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 near-core regions. Both these effects likely explain why the contributions to coupling coefficients are more concentrated in the near-core regions for more evolved models.

Examples of such radial coupling coefficient profiles are given in Figs. 4 and 5, which display the N² profile along with the coupling coefficient profile η₁(r). The profile in Fig. 4 describes the coupling contributions for the mode triad characterized by g mode radial orders (n₁, n₂, n₃) = (−22, −15, −53) in the mid-MS ΔX_c,  1 model in Table 4. The near-ZAMS example displayed in Fig. 5 describes the coupling contributions for the mode triad characterized by g mode radial orders (n₁, n₂, n₃) = (−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 mid-MS model are typically more confined to the near-core region in comparison to those contributions in the near-ZAMS model.

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 |η₁|, 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 sub-inertial (i.e., |s_φ|> 1), requiring the non-perturbative (TAR) description of the Coriolis force.

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 non-sectoral modes (due to geostrophic balance; see e.g., Gill 1982; Townsend 2005; Daszyńska-Daszkiewicz et al. 2015; Szewczuk & Daszyńska-Daszkiewicz 2017).

The co-rotating-frame mode frequencies also (on average) decrease for higher-mass 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 co-rotating frame is instrumental for determining linear, heat-driven 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³ – 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ńska-Daszkiewicz 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 non-linear 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 slower-rotating 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₁| for the identified isolated mode triads is [10⁴, 10⁰] (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 Ω₁ ≈ Ω₂ + Ω₃. 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 proof-of-concept study and do not study the dependence of non-linear parameters – such as the coupling coefficient – on stellar parameters in detail.

5. Impact on asteroseismic modeling

In this section, we discuss how the non-linear theoretically predicted observables defined in Sect. 3 can aid or complement current frequency-based asteroseismic modeling of SPB stars. The various non-linear observables of interest, or derived quantities thereof, are listed in Table 6.

5.1. Oscillation frequency and combination phase observables

We find that the non-linear 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{\mathrm{\Omega }}_{\phi ,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⁻³ to 10⁻⁴, 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⁻¹, this means that there are several frequencies whose (non-linear) frequency shifts are of order 10⁻⁵ d⁻¹ or smaller, rendering them of similar order as the typical errors for the frequencies derived from 4-year Kepler light curves (e.g., Bowman & Michielsen 2021). The largest shifts are obtained for the modes in the faster-rotating models, and are of order 10⁻² or 10⁻³ d⁻¹, which is approximately two orders of magnitude larger than typical uncertainties. In practice, the resonant non-linear frequency shifts are hard to distill from frequency lists deduced from prewhitening analyses similar to what was done in V21. The non-resonant frequency shifts are even smaller, justifying the approximation of using linear frequencies in asteroseismic modeling.

We show the zero-point-corrected non-linear 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, ϕ₀, 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 stationary-state (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_⊙ near-ZAMS models ΔX_c,  2(a) and ΔX_c,  2(b), as well as the 6 M_⊙ near-TAMS 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 ϑ₂ or ϑ₃ and proportional to the frequency ratio Ω₂ / Ω₁ or Ω₃ / Ω₁. We find that the frequency ratios for these specific isolated mode triads that have larger-amplitude daughter modes are not particularly different from those in other identified isolated mode triads, as can be rationalized from the co-rotating mode frequency values in Table 4. The observed difference in amplitude ratio is thus is mostly due to different values of ϑ₂. The smaller ϑ₂, 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 ϑ₂ in Table 5, this effect seems to become important when ϑ₂ ≲ 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 ϑ₂ at which this effect becomes important by performing a linear regression between the ϑ₂ 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 ϑ₂ and $\mid {\mathfrak{L}}_2^s\mid /\mid {\mathfrak{L}}_1^s\mid$ values obtained from the two mode triads with parent-daughter amplitude ratios greater than unity but approaching the unity limit the closest, and the values of these quantities obtained from those mode triads with parent-daughter amplitude ratios smaller than unity that approach the unity limit the closest. The linear regression yields a linear slope of −0.08 ± 0.04 (ϑ₂ 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 non-AE valid triads have both daughter-parent 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₁, the necessary but not sufficient condition for both stationary daughter-parent rotating-frame 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₁ are of the same order of magnitude) and likely obtain daughter-parent 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 daughter-parent 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).

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 signal-to-noise 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 pseudo-classes to the different members of the ensemble of analyzed SPB stars to denote the difficulties encountered during the analysis process: high-f_sv stars had light curves that were described adequately by the harmonic light curve models; high mode density stars had many close-spaced frequencies in their Fourier transforms or Lomb-Scargle periodograms; and outbursting stars were found to have distinct variability ‘outbursts’ in their light curves, resulting in many close-spaced and high-amplitude frequency groups (see V21). To select the relevant observed signals for our comparison, we look for sum-frequency resonance combinations of three signals that are part of the light curve models and fulfill the resonance criterion |ν₁ − ν₂ − ν₃|≤ℜ_ν + σ_{ν,  prop} (with σ_{ν,  prop} the propagated uncertainty of the observed frequency difference ν₁ − ν₂ − ν₃ and ν₁, ν₂ and ν₃ the observed frequencies). These combinations also need to contain (i) at least one of the two highest-amplitude signals in the light curve model, and (ii) lower-amplitude components (i.e., not one of the two highest-amplitude signals) that are not present in other considered combinations. The low-over-highest 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 non-linear 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 daughter-parent 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 zero-point-corrected non-linear 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 daughter-parent amplitude ratios with the relevant computed daughter-parent (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 }}_{\min ,\mathfrak{i}}^{\mathrm{NL},s}$ (used as the x-axis variable in the right panel of Fig. 6) is smaller than the non-linear 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⁻¹ (using the notation in Table 2; i.e., up to the notional Kepler long-cadence 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 (orange-red) symbols.

The ensemble of stationary surface luminosity fluctuation ratios of computed triads seems to be able to explain the low-amplitude-ratio part of the ensemble of observed daughter-parent 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 daughter-parent 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ńska-Daszkiewicz et al. (2017) used in modeling the linear stability of the hybrid main-sequence early B-type 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 co-rotating 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 daughter-parent amplitude ratios shown in Fig. 6. For comparison with the theoretically predicted daughter-parent 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 multi-mode coupling or higher-order 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 higher-order couplings might be needed to explain the high amplitude ratios of high-amplitude SPB pulsators.

6. Conclusions

We derive a theoretical oscillation modeling framework that describes non-linear three-mode coupling among gravito-inertial (g) modes of Slowly Pulsating B (SPB) stars within the Traditional Approximation for Rotation (TAR; e.g., Longuet-Higgins 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 non-linear energy exchange between the modes, by computing three-mode (energy-scaled) coupling coefficients |η₁| based on a phase space mode decomposition inserted in the relevant coupled equations of motion. To describe three-mode resonant couplings, we derive amplitude equations (AEs) from the coupled equations of motion for a sum-frequency resonance Ω₁ ≃ Ω₂ + Ω₃, 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 three-mode resonant couplings among g modes in rapidly rotating g-mode 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₁, k₂, k₃) = (0, 0, 0) and azimuthal orders (m₁, m₂, m₃) = (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 three-mode coupling scenarios in such stars. To ensure stability of the three-mode 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 multi-mode 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 near-core zones of the SPB models, where N² peaks and modes can become (partially) trapped (Miglio et al. 2008; Michielsen et al. 2021). For more evolved models, this N² maximum is wider and the κ mechanism typically excites higher radial order modes, which are increasingly confined to the near-core zones. Conversely, for less evolved models, significant contributions to |η₁| can be found outside of the near-core zones.

Linear heat-driven vibrational instability defines which modes are available for the non-linear (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 non-linear 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 detuning-damping 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 non-linear frequency shifts (54) for isolated triads are small enough to neglect when compared with the typical errors for the frequencies derived from 4-year 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 daughter-parent amplitude ratios are greater than unity, because of the typical values for the quality factors and detuning-damping 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 non-linear 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 non-linear 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 higher-order and r − g mode coupling can saturate some of the unexplained observed candidate couplings in SPB stars. Higher-order coupling will definitely already include self-saturation 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 non-linear 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 non-linear saturation is crucially dependent on the linear excitation properties, in-depth investigations of the linear vibrational instability of g modes in SPB stars using the latest improvements in opacity computations are crucial for future non-linear excitation studies (e.g., Daszyńska-Daszkiewicz 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.

---

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: 4D-STAR). 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. PHY-1748958. 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

Appendix A: Phase-space 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{·}\mathit{\zeta }+\mathsf{F},\end{array}$$(A.1)

in which the non-Hermitian operator T is defined as

$$\begin{array}{c}\hfill \mathsf{T}=\left[\begin{array}{cc}-\frac{1}{2}\mathit{B}& \mathbf{1}\\ -\mathit{C}+\frac{1}{4}{\mathit{B}}^2& -\frac{1}{2}\mathit{B}\\ \end{array}\right],\end{array}$$(A.2)

where ζ and F are defined as

$$\begin{array}{c}\hfill \mathit{\zeta }(\mathit{x},t)\equiv \left[\begin{array}{c}\mathit{\xi }(\mathit{x},t)\\ \mathit{\pi }(\mathit{x},t)\end{array}\right],\mathsf{F}=\left[\begin{array}{c}0\\ {\mathit{a}}_{\mathrm{ext}}\end{array}\right],\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).\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)e^{-i\omega t},\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{·}\mathit{\zeta }\left(\mathit{x}\right)=0.\end{array}$$(A.6)

Because T is non-Hermitian, 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{·}{\mathit{\zeta }}_{\phi }\left(\mathit{x}\right)=0.\end{array}$$(A.7)

These right eigenvectors defined by Eq. (A.7) are the eigenvectors of non-conjugated 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 rotating-frame 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 })=& \langle {\mathit{\xi }}_{\phi },\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 ,\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 ,\end{array}$$(B.2)

after multiplication with i, and using the anti-Hermiticity 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 }),\end{array}$$(B.3)

because of Ansatz (2) and the assumed time-independence 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,\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 real-valued 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 .\end{array}$$(B.5)

Because of the conserved nature of W(ξ_φ,ξ_φ) and the assumed time-independence 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,\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 real-valued constant b_φ defined in Eq. (7) for a real eigenfrequency Ω_φ. For complex non-degenerate 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 non-degenerate real-valued eigenfrequencies. The lack of an orthogonality relation for the non-adiabatic complex-valued eigenfunctions of oscillations in rotating stars (with associated complex-valued eigenfrequencies) that does not involve Jordan chain modes (see Appendix A of Schenk et al. 2001) motivates the use of adiabatic real-valued eigenfunctions (with associated real-valued eigenfrequencies) in the computations for the mode coupling coefficient (22).

Appendix C: Expressions for the explicit terms of the three-mode coupling coefficients

We base ourselves on explicit expressions derived by L12 for the terms that make up the three-mode 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 Γ₁, ${\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 non-linear 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{\xi }^{\varphi }=\frac{1}{rsin\theta }\frac{\partial {\xi }^{\varphi }}{\partial r}=\frac{i}{r}\frac{\partial }{\partial r}\left(r\frac{z_2}{c_1{\overline{\omega }}^2}\right)\frac{H_{\varphi }}{sin\theta },\end{array}$$(C.1)

in which z₂ 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₁ 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 (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 (u)}/{ϵ}_{\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 (1)}$ and ${\eta }_1^{(1)}$, and ${\kappa }_{\phi }^{\beta \gamma (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 (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 (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 (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 (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{z}_1)}{\partial r}(z_1-\frac{z_2}{c_1{\overline{\omega }}^2})\frac{\partial }{\partial r}\left(r\frac{z_2}{c_1{\overline{\omega }}^2}\right):H_r[H_{\theta }{H}_{\theta }-H_{\varphi }{H}_{\varphi }])+\hfill \\ \hfill & S(\frac{\partial (r{z}_1)}{\partial r}{z}_1\frac{\partial }{\partial r}\left(r\frac{z_2}{c_1{\overline{\omega }}^2}\right):H_r[\delta H_r{H}_{\theta }-\frac{m{H}_r{H}_{\varphi }}{sin\theta }])-\hfill \\ \hfill & S(\frac{\partial (r{z}_1)}{\partial r}{z}_1\frac{\partial }{\partial r}\left(r\frac{z_2}{c_1{\overline{\omega }}^2}\right):H_r[H_{\theta }{H}_{\theta }-H_{\varphi }{H}_{\varphi }])=\hfill \\ \hfill & S(\frac{\partial (r{z}_1)}{\partial r}{z}_1\frac{\partial }{\partial r}\left(r\frac{z_2}{c_1{\overline{\omega }}^2}\right):H_r[\delta H_r{H}_{\theta }-\frac{m{H}_r{H}_{\varphi }}{sin\theta }])-\hfill \\ \hfill & S(\frac{\partial (r{z}_1)}{\partial r}\frac{z_2}{c_1{\overline{\omega }}^2}\frac{\partial }{\partial r}\left(r\frac{z_2}{c_1{\overline{\omega }}^2}\right):H_r[H_{\theta }{H}_{\theta }-H_{\varphi }{H}_{\varphi }]),\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 computation-friendly 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 non-adiabatic 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 two-dimensional Jacobian as

$$\begin{array}{c}\hfill \frac{\partial (u,v)}{\partial (x,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|.\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 \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 \frac{\partial ({\mathrm{\Gamma }}_1,S)}{\partial (\rho ,S)}.\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 \frac{\partial ({\mathrm{\Gamma }}_1,S)}{\partial (\rho ,S)}\frac{\partial (\rho ,T)}{\partial (\rho ,T)}=\rho \frac{\frac{\partial ({\mathrm{\Gamma }}_1,S)}{\partial (\rho ,T)}}{\frac{\partial (\rho ,S)}{\partial (\rho ,T)}},\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 [{\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\}{\left(\frac{\partial {\mathrm{\Gamma }}_1}{\partial T}\right)}_{\rho }],\end{array}$$(C.6)

whose right hand side only contains thermodynamic quantities that are computed in stellar evolution codes such as MESA.

An adiabatic second-order derivative of adiabatic exponent Γ₁ with respect to mass density ρ appears in coupling coefficients that describe four-mode interactions (or higher-order derivatives when higher-order coupling terms are considered). Higher-order 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 higher-order derivatives of Γ₁. Such higher-order derivatives can alternatively be estimated from numerical derivatives of Eq. (C.6) when using the output of current state-of-the-art 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 })=& -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 ,\hfill \end{array}$$(D.1a)

$$\begin{array}{cc}\hfill W({\mathit{\xi }}_{\phi },{\mathit{\xi }}_{\beta }^{\ast })=& -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 ,\hfill \end{array}$$(D.1b)

$$\begin{array}{cc}\hfill W({\mathit{\xi }}_{\phi }^{\ast },{\mathit{\xi }}_{\beta }^{\ast })=& 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 .\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,\end{array}$$(D.2)

when ω_φ ≠ −ω_β, which has a real-valued 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.\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,\end{array}$$(D.4)

when ${\omega }_{\phi }^{\ast }\ne -{\omega }_{\beta }^{\ast }$. The real-valued 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,\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,\end{array}$$(D.6)

for Eq. (D.1c) and ${\omega }_{\phi }\ne {\omega }_{\beta }^{\ast }$. Its real-valued 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 }{b}_{\phi }^{\ast }={\delta }_{\beta }^{\phi }{b}_{\phi },\end{array}$$(D.7)

is the complex conjugate of Eq. (6), and is valid for non-degenerate 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 }\mathit{\xi }(\mathit{x},t)+i\dot{\mathit{\xi }}(\mathit{x},t)+i\mathit{B}\left(\mathit{\xi }(\mathit{x},t)\right)$ with ${\mathit{\xi }}_{\beta }^{\ast }$ and integrating over the stellar mass yields

$$\begin{array}{c}\hfill \langle {\mathit{\xi }}_{\beta },{\mathrm{\Omega }}_{\beta }\mathit{\xi }(\mathit{x},t)+i\dot{\mathit{\xi }}(\mathit{x},t)+i\mathit{B}\left(\mathit{\xi }(\mathit{x},t)\right)\rangle \equiv \mathcal{P}.\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 }}({\mathrm{\Omega }}_{\beta }+{\mathrm{\Omega }}_{\phi })c_{\phi }\left(t\right)\langle {\mathit{\xi }}_{\beta }\left(\mathit{x}\right),{\mathit{\xi }}_{\phi }\left(\mathit{x}\right)\rangle \hfill \\ \hfill & +({\mathrm{\Omega }}_{\beta }-{\mathrm{\Omega }}_{\phi })c_{\phi }^{\ast }\left(t\right)\langle {\mathit{\xi }}_{\beta }\left(\mathit{x}\right),{\mathit{\xi }}_{\phi }^{\ast }\left(\mathit{x}\right)\rangle \hfill \\ \hfill & +c_{\phi }\left(t\right)\langle {\mathit{\xi }}_{\beta }\left(\mathit{x}\right),i\mathit{B}\left({\mathit{\xi }}_{\phi }\left(\mathit{x}\right)\right)\rangle \hfill \\ \hfill & +c_{\phi }^{\ast }\left(t\right)\langle {\mathit{\xi }}_{\beta }\left(\mathit{x}\right),i\mathit{B}\left({\mathit{\xi }}_{\phi }^{\ast }\left(\mathit{x}\right)\right)\rangle .\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 non-degenerate modes φ and β that

$$\begin{array}{c}\hfill \mathcal{P}=c_{\beta }\left(t\right)b_{\beta },\end{array}$$(D.10)

assuming Ω_φ ≠ −Ω_β, and using the Kronecker delta in Eq. (6). Equation (D.10) thus proves Eq. (18a) for non-degenerate modes φ and β that have Ω_β ≠ −Ω_φ, and a non-zero b_β.

Similarly, if we multiply ${\mathrm{\Omega }}_{\beta }\mathit{\xi }(\mathit{x},t)-i\dot{\mathit{\xi }}(\mathit{x},t)-i\mathit{B}\left(\mathit{\xi }(\mathit{x},t)\right)$ with ξ_β and integrate over the stellar mass, we obtain

$$\begin{array}{c}\hfill \langle {\mathit{\xi }}_{\beta }^{\ast },{\mathrm{\Omega }}_{\beta }\mathit{\xi }(\mathit{x},t)-i\dot{\mathit{\xi }}(\mathit{x},t)-i\mathit{B}\left(\mathit{\xi }(\mathit{x},t)\right)\rangle \equiv \mathcal{Q},\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}=& {\displaystyle \sum _{\phi }}({\mathrm{\Omega }}_{\beta }-{\mathrm{\Omega }}_{\phi })c_{\phi }\left(t\right)\langle {\mathit{\xi }}_{\beta }^{\ast }\left(\mathit{x}\right),{\mathit{\xi }}_{\phi }\left(\mathit{x}\right)\rangle \hfill \\ \hfill & +({\mathrm{\Omega }}_{\beta }+{\mathrm{\Omega }}_{\phi })c_{\phi }^{\ast }\left(t\right)\langle {\mathit{\xi }}_{\beta }^{\ast }\left(\mathit{x}\right),{\mathit{\xi }}_{\phi }^{\ast }\left(\mathit{x}\right)\rangle \hfill \\ \hfill & -c_{\phi }\left(t\right)\langle {\mathit{\xi }}_{\beta }^{\ast }\left(\mathit{x}\right),i\mathit{B}\left({\mathit{\xi }}_{\phi }\left(\mathit{x}\right)\right)\rangle \hfill \\ \hfill & -c_{\phi }^{\ast }\left(t\right)\langle {\mathit{\xi }}_{\beta }^{\ast }\left(\mathit{x}\right),i\mathit{B}\left({\mathit{\xi }}_{\phi }^{\ast }\left(\mathit{x}\right)\right)\rangle .\hfill \end{array}\end{array}$$(D.12)

The following relations

$$\begin{array}{cc}& \langle {\mathit{\xi }}_{\beta }^{\ast }\left(\mathit{x}\right),{\mathit{\xi }}_{\phi }\left(\mathit{x}\right)\rangle =\langle {\mathit{\xi }}_{\phi }^{\ast }\left(\mathit{x}\right),{\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),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),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 anti-Hermiticity 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.\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)b_{\beta },\end{array}$$(D.15)

where we use the Kronecker delta in Eq. (D.7). This proves Eq. (18b) for non-zero 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 ){\mathit{e}}_{\widehat{r}}-{\xi }^{\widehat{\theta }}(r,\pi -\theta ,\varphi ){\mathit{e}}_{\widehat{\theta }}+{\xi }^{\widehat{\varphi }}(r,\pi -\theta ,\varphi ){\mathit{e}}_{\widehat{\varphi }}.\end{array}$$(E.1)

We then have

$$\begin{array}{c}\hfill f_{\ast }{\mathit{\xi }}_{\phi }={\mathcal{Z}}_{\phi }{\mathit{\xi }}_{\phi },\end{array}$$(E.2)

with 𝒵_φ (called the z-parity 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 real-valued. Based on Eq. (10) we then have that

$$\begin{array}{c}\hfill f_{\ast }{\mathit{\xi }}_{\phi }^{\ast }={\mathcal{Z}}_{\phi }{\mathit{\xi }}_{\phi }^{\ast }.\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 },f_{\ast }{\mathit{\xi }}_{\beta },f_{\ast }{\mathit{\xi }}_{\gamma })={\mathcal{Z}}_{\phi }{\mathcal{Z}}_{\beta }{\mathcal{Z}}_{\gamma }\kappa ({\mathit{\xi }}_{\phi }^{\ast },{\mathit{\xi }}_{\beta },{\mathit{\xi }}_{\gamma }).\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,\end{array}$$(E.5)

resulting in the meridional or z-parity selection rule (S01)

$$\begin{array}{c}\hfill i_{\mathrm{odd}}=\mathrm{mod}(0,2),\end{array}$$(E.6)

where i_odd denotes the number of odd modes. Hence, for three-mode coupling, the z-parity 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 sum-frequencies Ω₁ ≈ Ω₂ + Ω₃ with those derived for difference-frequencies Ω_1d ≈ Ω_2d − Ω_3d using the same methods employed in Sect. 2.6. This equivalence indicates that the stationary properties of resonant isolated sum-frequency 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 non-linear mode couplings.

To prove this equivalence, we first derive the difference-frequency AEs using the same procedure as in Sect. 2.6. Here, we define the resonance condition Ω_1d ≈ Ω_2d − Ω_3d in terms of the difference-frequency detuning parameter δω_d:

$$\begin{array}{c}\hfill \mathfrak{J}\delta {\omega }_d={\mathrm{\Omega }}_{1_d}\approx {\mathrm{\Omega }}_{2_d}-{\mathrm{\Omega }}_{3_d}.\end{array}$$(F.1)

In Eq. (F.1), the subscript d refers to the fact that we are considering a difference-frequency resonance. If we then set the secular-term-generating terms in Eq. (33) equal to zero, use the difference-frequency 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 difference-frequency resonance Ω_1d ≈ Ω_2d − Ω_3d:

$$\begin{array}{c}\hfill \frac{\partial {\mathit{a}}_d}{\partial t_1}={\mathit{\gamma }}_d°{\mathit{a}}_d+2i({\mathit{\eta }}_{\mathit{d}}°{\mathit{a}}_{M_d}°{\mathbf{\Omega }}_{M_d}°{\mathit{e}}_d).\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],{\mathit{a}}_{M_d}=\left[\begin{array}{c}{a}_{2_d}{a}_{3_d}^{\ast }\\ a_{1_d}{a}_{3_d}\\ a_{1_d}^{\ast }{a}_{2_d}\end{array}\right],{\mathit{e}}_d=\left[\begin{array}{c}exp(i\delta {\omega }_d{t}_1)\\ exp(-i\delta {\omega }_d{t}_1)\\ exp(i\delta {\omega }_d{t}_1)\end{array}\right],\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],{\mathit{\eta }}_d=\left[\begin{array}{c}{\eta }_{1_d}^{\ast }\\ {\eta }_{1_d}\\ {\eta }_{1_d}^{\ast }\end{array}\right],{\mathbf{\Omega }}_{M_d}=\left[\begin{array}{c}{\mathrm{\Omega }}_{1_d}\\ {\mathrm{\Omega }}_{2_d}\\ {\mathrm{\Omega }}_{3_d}\end{array}\right],\end{array}$$(F.4)

in which we use the (difference-frequency) 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}.\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 (real-valued) AEs for difference frequencies,

$$\begin{array}{cc}\hfill \frac{\partial {\mathit{A}}_d}{\partial t_1}& ={\mathit{\gamma }}_d°{\mathit{A}}_d+2|{\eta }_{1_d}|sin\left({\mathrm{{\rm Y}}}_d\right)({\mathit{A}}_{N_d}°{\mathbf{\Omega }}_{M_d}),\hfill \end{array}$$(F.6a)

$$\begin{array}{cc}\hfill {\mathit{A}}_d°\frac{\partial {\mathit{\varphi }}_d}{\partial t_1}& =2|{\eta }_{1_d}|cos\left({\mathrm{{\rm Y}}}_d\right)({\mathit{A}}_{P_d}°{\mathbf{\Omega }}_{M_d}),\hfill \end{array}$$(F.6b)

where

$$\begin{array}{c}\hfill {\eta }_{1_d}=|{\eta }_{1_d}|e^{-i{\delta }_{1_d}},{\mathit{A}}_d=\left[\begin{array}{c}{A}_{1_d}\\ A_{2_d}\\ A_{3_d}\end{array}\right],{\mathit{\varphi }}_d=\left[\begin{array}{c}{\varphi }_{1_d}\\ {\varphi }_{2_d}\\ {\varphi }_{3_d}\end{array}\right],\end{array}$$(F.7)

as well as

$$\begin{array}{c}\hfill {\mathit{A}}_{N_d}=\left[\begin{array}{c}{A}_{2_d}{A}_{3_d}\\ -A_{1_d}{A}_{3_d}\\ A_{1_d}{A}_{2_d}\end{array}\right],{\mathit{A}}_{P_d}=\left[\begin{array}{c}{A}_{2_d}{A}_{3_d}\\ A_{1_d}{A}_{3_d}\\ A_{1_d}{A}_{2_d}\end{array}\right].\end{array}$$(F.8)

In Eq. (F.6), we define the generic (frequency-difference) phase coordinate Υ_d as

$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}_d\equiv -\delta {\omega }_d{t}_1+{\varphi }_{1_d}-{\varphi }_{2_d}+{\varphi }_{3_d}-{\delta }_{1_d},\end{array}$$(F.9)

which contains the combination phase Φ_d = ϕ_1d − ϕ_2d + ϕ_3d for difference frequencies. The time-dependence 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 }_{{⊞}_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}),\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 sum-frequency analogue (dropping the subscript d in the process), leads to the equivalence of its expression with that of the sum-frequency coupling coefficient η₁ defined in Eq. (22). Performing the same relabeling operation for the frequency-difference phase coordinate Υ_d defined in Eq. (F.9), we obtain that it is equal to −Υ. Substituting this relation between Υ_d and Υ in the difference-frequency AEs (F.6) converts them into the sum-frequency 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⁴ 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

All Figures