4pt 4pt 11institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan
200D, 3001 Leuven, Belgium
e-mail<EMAIL_ADDRESS>22institutetext: Kavli Institute for
Theoretical Physics, University of California, Santa Barbara, CA 93106, USA
33institutetext: Reference Systems and Planetology, Royal Observatory of
Belgium, Brussels, Belgium 44institutetext: Dept. of Astrophysics, IMAPP,
Radboud University Nijmegen, 6500 GL, Nijmegen, The Netherlands
55institutetext: Max Planck Institute for Astronomy, Koenigstuhl 17, 69117
Heidelberg, Germany

# Detection of nonlinear resonances among gravity modes of slowly pulsating B
stars: results from five iterative prewhitening strategies

J. Van Beeck Detection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening strategies D. M.
Bowman Detection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening strategies M. G.
Pedersen Detection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening strategies T. Van
Reeth Detection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening strategies T. Van
Hoolst Detection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening strategies C.
Aerts Detection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening
strategiesDetection of nonlinear resonances among gravity modes of slowly
pulsating B stars: results from five iterative prewhitening strategies

(Received June 17, 2021 / Accepted August 5, 2021)

###### Abstract

Context. Slowly pulsating B (SPB) stars are main-sequence multi-periodic
oscillators that display non-radial gravity modes. For a fraction of these
pulsators, 4-year photometric light curves obtained with the Kepler space
telescope reveal period spacing patterns from which their internal rotation
and mixing can be inferred. In this inference, any direct resonant mode
coupling has usually been ignored so far.

Aims. We re-analysed the light curves of a sample of 38 known Kepler SPB
stars. For 26 of those, the internal structure, including rotation and mixing,
was recently inferred from their dipole prograde oscillation modes. Our aim is
to detect direct nonlinear resonant mode coupling among the largest-amplitude
gravity modes.

Methods. We extract up to 200 periodic signals per star with five different
iterative prewhitening strategies based on linear and nonlinear regression
applied to the light curves. We then identify candidate coupled gravity modes
by verifying whether they fulfil resonant phase relations.

Results. For 32 of 38 SPB stars we find at least 1 candidate resonance that is
detected in both the linear and the best nonlinear regression model fit to the
light curve and involves at least one of the two largest-amplitude modes.

Conclusions. The majority of the Kepler SPB stars reveal direct nonlinear
resonances based on the largest-amplitude modes. These stars are thus prime
targets for nonlinear asteroseismic modelling of intermediate-mass dwarfs to
assess the importance of mode couplings in probing their internal physics.

###### Key Words.:

asteroseismology – stars: oscillations (including pulsations) – stars: early-
type – stars: variables: general – stars: rotation – methods: data analysis

##  1 Introduction

Three decades ago Waelkens (1991) first coined the group of “slowly pulsating
B” (SPB) stars to describe mid-B variable stars that display low-frequency
oscillations with physical properties similar to 53 Persei variables (see
e.g., Waelkens & Rufener 1985, and references therein). The SPB stars are mid-
to-late B main-sequence (MS) stars with masses ranging from $\sim 3$ to $\sim
9$ M☉ (De Cat & Aerts 2002; Aerts et al. 2010; Pedersen et al. 2021)
displaying all levels of rotation, from very slow up to the critical rotation
rate (Pedersen et al. 2021). A distinct feature is their multi-periodic
photometric and spectroscopic variability, with periods typically ranging from
$\sim 0.5$ to $\sim 5.0$ d and amplitudes typically less than $\sim 10$ mmag
(e.g., Aerts et al. 2010; Pápics et al. 2017; Pedersen et al. 2021). Their
variability is attributed to low-degree, high-radial order gravity ($g$)
modes, which are excited by the heat engine (i.e. $\kappa$) mechanism driven
by an opacity enhancement from the iron-group elements, also called the Z
Bump, at a temperature of $\sim 200\leavevmode\nobreak\ 000$ K (Gautschy &
Saio 1993; Dziembowski et al. 1993; Pamyatnykh 1999).

The oscillations in SPB stars are fairly well-characterised from ground-based
photometric and spectroscopic observations (e.g., Aerts et al. 1999; De Cat et
al. 2000; Mathias et al. 2001; De Cat & Aerts 2002; De Cat et al. 2005, 2007),
but their asteroseismic potential only became clear when CoRoT (Auvergne et
al. 2009) and Kepler (Borucki et al. 2010) space photometry became available
(see e.g., Degroote et al. 2010; Balona et al. 2011, 2015; Pápics et al. 2012,
2014, 2015, 2017; Szewczuk & Daszyńska-Daszkiewicz 2018; Pedersen et al. 2021;
Szewczuk et al. 2021, for recent examples of light curves). Because of their
wealth in number of excited and identified oscillation modes, they have been
the target of various asteroseismic modelling efforts, which focused on
inferring the internal mixing and rotation (e.g., Moravveji et al. 2015, 2016;
Triana et al. 2015; Pápics et al. 2017; Pedersen et al. 2021; Michielsen et
al. 2021) and on mode excitation (e.g., Szewczuk & Daszyńska-Daszkiewicz 2018;
Szewczuk et al. 2021). There have also been efforts to characterise the
magnetic field of a small number of SPB stars (e.g. HD 43317) by analysing the
frequency shifts of its $g$-modes (Buysschaert et al. 2018; Prat et al. 2019,
2020). The influence of rotation on the $g$ modes in SPB stars is more
important than the magnetic shifts and imply that almost all modes are in the
gravito-inertial regime (Aerts et al. 2017, 2019). The effects of rotation,
therefore, cannot be incorporated from a perturbative treatment of the
Coriolis acceleration, which led many authors to adopt the so-called
traditional approximation of rotation (TAR; see e.g., Longuet-Higgins 1968;
Lee & Saio 1997; Townsend 2003; Degroote et al. 2009b; Mathis 2013) for
pulsation computations (e.g., Szewczuk & Daszyńska-Daszkiewicz 2018; Pedersen
et al. 2021; Szewczuk et al. 2021; Michielsen et al. 2021).

Rapid rotation also deforms a star, causing it to become oblate. The rotation
rates of SPB stars cover the cases from slight to major deformation. However,
the detected and identified modes have their dominant mode energy deep inside
the star, close to the convective core, where spherical symmetry still applies
well. This justifies a treatment of g modes in the TAR (Mathis & Prat 2019;
Henneco et al. 2021; Dhouib et al. 2021). An important consequence of the
rapid rotation of some SPB stars is the occurrence of mode coupling between
overstable convective modes based on frozen-in core convection and envelope
$g$ modes excited by the $\kappa$ mechanism (Lee 2021).

SPB stars have so far been modelled from linear pulsation theory. Period
spacing patterns constructed from the frequencies of identified $g$ modes of
consecutive radial order and the same angular degree and azimuthal order, are
the most important tool used in linear asteroseismic modelling. Such patterns
can display oscillatory features attributable to chemical gradients inside the
star (e.g., Miglio et al. 2008). Moreover, they display a slope due to the
rotation rate near the convective core (e.g., Bouabid et al. 2013; Pápics et
al. 2017). Some SPB stars, however, can show structure in their g-mode period
spacing patterns that is hard to explain without incorporating additional
physical effects in the oscillation model, such as rotational deformation
(e.g., Dhouib et al. 2021) or (internal) magnetism (see e.g., Van Beeck et al.
2020), or resonances of envelope $g$ modes with inertial modes in the core
(e.g., Ouazzani et al. 2020; Saio et al. 2021).

The input to asteroseismic modelling based on linear oscillation theory are
the independent frequencies extracted from an observed time series (e.g., a
light curve), where frequency extraction is usually done using an iterative
prewhitening procedure (e.g., Degroote et al. 2009b). However, non-sinusoidal
light curves give rise to harmonics and combination frequencies (dependent
frequencies) in the amplitude spectrum (Pápics 2012; Kurtz et al. 2015). These
combination frequencies are therefore also included in the extracted
frequencies. The phases and amplitudes of the independent frequencies, as well
as the frequencies, amplitudes, and phases of the combination frequencies are
subsequently neglected during linear asteroseismic modelling. On the other
hand, frequency perturbations within g-mode period spacing patterns are
expected if nonlinear mode coupling is accounted for (e.g., Buchler & Goupil
1984; Van Hoolst 1996). Hence, this may be an alternative explanation for the
deviations in the period spacing patterns of SPB stars. Such frequency shifts
may be small if the $g$ modes are only moderately nonlinear, as is the case
for most ‘well-observed’ nonradial pulsators such as white dwarfs or $\delta$
Sc stars (Buchler et al. 1997), yet this remains to be verified for SPB stars.

As detailed in Buchler et al. (1997), resonances among (excited) non-radial
oscillation modes also play an important role in mode selection and
interaction, and can lead to additional constraints on mode identification.
Resonant interactions among modes can further hamper pattern identification
due to so-called frequency or phase locking. To complicate pattern
identification the exact resonance relation, $\sum_{i}n_{i}\,\omega_{i}=0$,
need not be satisfied for linear mode frequencies. It suffices that linear
mode frequencies are in near resonance, such that the exact resonance relation
between angular frequencies $\omega_{j}$ is satisfied aside from a small
angular frequency difference $\delta\omega$. Nonlinear frequency shifts can
then lead to a nonlinear frequency locking or phase locking, where the
nonlinear locked frequencies satisfy the resonance condition exactly, which
can distort the period spacing patterns relied upon in linear
asteroseismology.

The theoretical framework of resonant mode coupling, the amplitude equation
(AE) formalism, is described in numerous publications (e.g., Dziembowski 1982;
Buchler & Goupil 1984; Dziembowski & Krolikowska 1985; Buchler 1993; Buchler
et al. 1997; Van Hoolst & Smeyers 1993; Van Hoolst 1994a, b, 1995). It has the
inherent validity assumption that the (linear) growth rates, $\gamma$, of the
modes are small compared to their angular frequencies, which is generally
assumed to be satisfied for observed MS oscillators (see e.g., Barceló Forteza
et al. 2015; Bowman et al. 2016). The AEs constitute a set of coupled first-
order complex nonlinear ordinary differential equations and predict the
temporal behaviour of the mode amplitudes and phases. Three main regimes of
mode interaction are distinguished (as explained in Buchler et al. 1997):

  1. 1.

If two or more modes are close to being in resonance, such that
$D_{g}\equiv\delta\omega/\gamma$ is of order $1$ or smaller (where
$\delta\omega$ is the angular frequency mismatch or angular frequency
detuning), and a stable fixed point (FP) solution exists for the AEs, for
which the mode amplitudes and phases are constant over time, the frequencies
are locked. These locked frequencies can be substantially different from their
non-resonant counterparts.

  2. 2.

With an increasing angular frequency mismatch, $\delta\omega$, the FP solution
becomes unstable or disappears, and a bifurcation to another solution takes
place, with amplitude (and phase) modulation. The larger the value of $D_{g}$,
the smaller the amplitude of the modulation. This is called the ‘intermediate’
regime (Goupil et al. 1998).

  3. 3.

At large $\delta\omega$ one encounters the solution where frequencies are
approximately equal to their (unlocked) non-resonant values. The nonlinear
nonresonant frequencies have a small frequency shift compared to their linear
counterparts: this is the mildly nonlinear regime in which linear
asteroseismology operates, also called the ‘nonresonant’ regime (Goupil et al.
1998).

In addition, a narrow hysteresis (transitory) regime exists in between the
frequency-locked and ‘intermediate’ regime, in which amplitudes are modulated,
but frequencies (or equivalently, phases) remain constant (Buchler et al.
1997).

The study of the relationships among amplitudes and phases, which provide the
observables to be matched to the theoretical predictions from the AEs, is a
standard approach for the analysis of light curves of RR Lyr stars and
Cepheids (as pioneered by Simon & Lee 1981). It is also a well-established
technique for unravelling nonlinearities in the light curves of oscillating
white dwarf stars (e.g., Wu 2001; Zong et al. 2016b), and hot B subdwarf stars
(e.g., Zong et al. 2016a). Yet, it has not received much attention for
oscillating MS stars (Lee 2012; Bowman & Kurtz 2014; Kurtz et al. 2015; Bowman
et al. 2016, 2021). In this work, we provide the identification of candidate
direct resonances among the extracted oscillation frequencies in the light
curves of 38 SPB stars, and test the impact of different iterative
prewhitening strategies on the results. In Sect. 2 we summarise the
characteristics of our sample, the strategies we employ to analyse the light
curves, and discuss how we search for candidate direct resonances. We discuss
the light curve analysis strategy selection process in Sect. 3, and the
results of the candidate resonance search in Sect. 4. A summary of our
conclusions and future prospects can be found in Sect. 5.

##  2 Light curves and their analysis

###  2.1 The Pedersen (2020) SPB sample

Our sample consists of 38 stars originally identified as SPB stars based on
_Kepler_ photometry (Balona et al. 2011, 2015; McNamara et al. 2012; Pápics et
al. 2013, 2014, 2015, 2017; Zhang et al. 2018; Szewczuk & Daszyńska-
Daszkiewicz 2018). We use the light curves extracted from the 30-min cadence
target pixel data by Pedersen (2020), who used customised pixel masks (see
Pápics et al. 2017 for further details). After light curve extraction,
Pedersen (2020) manually removed outliers, and separately corrected each
quarter for instrumental trends by fitting a low order polynomial to the data.
Pedersen (2020) then converted the light curves to have flux units of ppm, and
normalised them to have an average flux of zero.

Contrary to previous analyses of Kepler SPB stars by Szewczuk & Daszyńska-
Daszkiewicz (2018), Pedersen (2020), Szewczuk et al. (2021) and Pedersen et
al. (2021), we use five different iterative prewhitening strategies, each of
which produce a final regression model fit for the light curve. Our approach
is motivated by the fact that iterative prewhitening techniques vary in the
literature, while it is not assessed how these various approaches influence
frequency extraction results. Our comparative capacity assessment of the
different iterative prewhitening procedures is an important aspect of the hunt
for nonlinear resonances in SPB stars.

###  2.2 SPB light curve analysis

Table 1: Five iterative prewhitening strategies employed in this work.

Strategy | Hinting | Stop criterion | Final optimisation aa$a$Type of least-squares optimisation of the entire model at the end of each iterative prewhitening step.
---|---|---|---
$1$ | S/N | S/N $<4.0$ | nonlinear
$2$ | A |  $p_{A}$ bb$b$The $p$-value for a Z-test based on extracted amplitude. $>0.05$ | nonlinear
$3$ | A |  $p_{\rm LRT}$ cc$c$The $p$-value for a likelihood ratio test (LRT) based on Bayesian information criterion (BIC) values for the nested light curve regression models. $>0.05$ | nonlinear
$4$ | S/N |  $p_{\rm LRT}$ cc$c$The $p$-value for a likelihood ratio test (LRT) based on Bayesian information criterion (BIC) values for the nested light curve regression models. $>0.05$ | nonlinear
$5$ dd$d$The frequency extracted during the iterative prewhitening step is optimised by selecting the largest-amplitude signal in an over-sampled LS periodogram in a 1 d-1 region around the frequency hint. | A |  $p_{\rm LRT}$ cc$c$The $p$-value for a likelihood ratio test (LRT) based on Bayesian information criterion (BIC) values for the nested light curve regression models. $>0.05$ | linear
111 Parameter hinting is either performed by selecting the highest-amplitude
peak (A) in the residual LS periodogram, or the peak with the highest
amplitude S/N value using a 1 d-1 window (S/N). The selected peak is not
closer than 2.5 / $T$ to earlier extracted frequencies, where $T$ is the total
time span of the data (Loumos & Deeming 1978). Linear interpolation between
the over-sampled frequency grid points of the Lomb-Scargle periodogram (Lomb
1976; Scargle 1982) is employed to determine the S/N value associated with the
optimised frequency.

We employ five prewhitening strategies to derive frequencies, amplitudes and
phases from the light curves, which emulate various approaches taken in the
literature. The first timestamp of each light curve is subtracted from all
timestamps of that light curve to provide a consistent zero point in time for
phase calculation. Strategy 1 uses a nonlinear least-squares fit of the entire
multi-parameter model in the time domain at each prewhitening step (Bowman et
al. 2016). It requires each signal to have an amplitude S/N $\geq$ 4 (Breger
et al. 1993) computed in a window of 1 d-1 around the target frequency, and
fits a model comprising a sum of sinusoids, $F(t_{i})$, to the light curve:

$F(t_{i})=\beta_{0}+\sum_{j=1}^{n_{f}}A_{j}\sin\left[2\pi\nu_{j}t_{i}+\phi_{j}\right]\,,$
(1)

for each $t_{i}$, $i=1,\ldots,n_{t}$, with $n_{t}$ the total number of
timestamps of the light curve. In this expression, $n_{f}$ is the number of
fitted frequencies, $\beta_{0}$ is the y-intercept, $A_{j}$ denotes the
amplitude, $\nu_{j}\left(\equiv\omega_{j}/2\pi\right)$ is the temporal
frequency, and $\phi_{j}$ is the phase of fitted sinusoid $j$. The other four
prewhitening strategies differ from this first one in the way they determine
the initial parameters for the optimisation (i.e. parameter hinting; see the
notes of Table 1), stop criterion, and whether a linear or nonlinear (multi-
parameter) least-squares optimisation is performed at the end of each
iterative prewhitening step. These differences are listed in Table 1.

As is commonly done, we calculate the uncertainties for the extracted
parameters, under the assumption that they are uncorrelated, using

$\hat{\sigma}_{\nu}=D\,\dfrac{\sqrt{6}\,\sigma_{n_{t}}}{\pi\,\sqrt{n_{t}}\,\hat{A}\,{T}},\hskip
8.5359pt\hat{\sigma}_{A}={D}\,\sqrt{\dfrac{2}{n_{t}}}\sigma_{n_{t}},\hskip
8.5359pt\hat{\sigma}_{\phi}={D}\,\sqrt{\dfrac{2}{n_{t}}}\dfrac{\sigma_{n_{t}}}{\hat{A}}\,,$
(2)

where $\sigma_{n_{t}}$ is the standard deviation of the residual signal (after
the regression model has been subtracted from the light curve), $\hat{A}$ is
the optimised amplitude, $T$ is the total time span of the data, and where $D$
is a correction factor for the correlated nature of the light curve data,
which can be estimated as the square root of the average number of consecutive
data points of the same sign in the residual light curve (e.g., Schwarzenberg-
Czerny 1991; Montgomery & O’Donoghue 1999; Schwarzenberg-Czerny 2003; Degroote
et al. 2009b). During each prewhitening step we compute the Lomb-Scargle (LS)
periodogram (Lomb 1976; Scargle 1982) of the (residual) light curve, with
frequency grid step equal to $0.1\mathfrak{R}_{\nu}$ (where
$\mathfrak{R}_{\nu}\equiv\tfrac{1}{T}$ is the Rayleigh limit). The start and
end values of this grid are equal to $1.5\,\mathfrak{R}_{\nu}$ and the Nyquist
frequency $\nu_{\rm Nyq}$ ($\nu_{\rm Nyq}=1/(2\,\tilde{\Delta t})$, where
$\tilde{\Delta t}$ is the median cadence of the time series), respectively.
The LS periodogram, together with a specific parameter hinting rule listed in
Table 1, are used to provide the preliminary amplitude and frequency hint,
taking into account the frequency resolution (Loumos & Deeming 1978),
hereafter referred to as the LD78 resolution. The preliminary phase hint is
subsequently obtained from a linear least-squares fit to the residual light
curve, where frequencies and amplitudes in Eqn. (1) are fixed at either their
frequency and amplitude hint values (the frequency being extracted), or the
values they had in the previous prewhitening step (all other frequencies). The
phase hinting optimisation employs the Trust-region reflective least-squares
regression method implemented in the Python package lmfit (Newville et al.
2020), whose exact configuration is specified in Appendix A. For strategies 1
to 4 the optimised parameter hints are then fed to the nonlinear least-squares
optimiser at the end of this prewhitening step to optimise the frequencies,
amplitudes and phases of the entire model. For strategy 5, however, these
parameter hints are only used in a linear least-squares optimisation of the
amplitudes and phases, where the frequencies are fixed at the values
determined from an over-sampled LS periodogram in a 1 d-1 region around the
preliminary frequency hint (signal extracted during this iterative
prewhitening step), and the values they had in the previous iteration step
(other signals). The Levenberg-Marquardt method implemented in lmfit (Newville
et al. 2020) is used for this optimisation, whose exact configuration is
specified in Appendix A.

Each iterative prewhitening process continues until 200 frequencies have been
extracted in this way, or until the specified stop criterion as in Table 1 is
reached. After the frequency extraction, we verify that the amplitude of each
extracted frequency is significant at the 95% level using a Z-test, whilst
taking the LD78 resolution into account. If at least one insignificant
extracted frequency is detected, the last extracted insignificant frequency is
removed. After that frequency is removed, the entire model is re-optimised,
and a new significance check is performed. This iterative process, which we
refer to as the frequency filtering step, continues until no further
insignificant frequencies occur in the regression model.

The regression models for each of the prewhitening strategies that have
undergone frequency filtering are then compared with one another. To quantify
the quality of the regression model, we rely on the residual sum of squares,
weighted by the total variance in the observed light curve and by a factor
including the total number of free parameters of the regression model,
$n_{p}$, versus $n_{t}$. To end up with a maximum likelihood estimator, we
subsequently use the scaled fraction of variance, $f_{\rm sv}$, which we
define as

$\displaystyle f_{\rm sv}\equiv\
1-\left(\dfrac{\sum_{i}^{n_{t}}\left(y_{i}-F(t_{i})\right)^{2}}{\sum_{i}^{n_{t}}\left(y_{i}-\bar{y}\right)^{2}}\right)\left\\{\dfrac{{n_{t}}-1}{{n_{t}}-n_{p}}\right\\}\,,$
(3)

where $y_{i}$ denotes the time series signal at timestamp $t_{i}$, and
$\bar{y}=\sum_{i}^{n_{t}}y_{i}/n_{t}$ is the mean signal of the time series.
The higher the quality of the regression model, the closer $f_{\rm sv}$
approaches unity while we punish for a higher number of free parameters.

###  2.3 Candidate resonant Kepler SPB oscillations

In our identification of direct resonant mode interactions we focus on FP
solutions of the AEs, for which frequency locking persists (Buchler et al.
1997). Any such solutions have parent-daughter frequency combinations that
approximately satisfy the (direct) resonance relation. Our aim is to extract
information on the least complicated mode interactions, hence, we use strict
candidate resonance identification criteria. The first criterion requires that
a frequency combination of a candidate resonance satisfies the following
relation:

$\exists\nu_{\rm daughter}:\hskip 5.69046pt\delta\nu\,\equiv\nu_{\rm
daughter}-\sum_{i}n_{i}\,\nu_{\rm parent,\,i}\approx 0\,\text{d}^{-1}\,,$ (4)

where $\nu_{\rm daughter}$ is the daughter frequency, $\nu_{\rm parent,\,i}$
is its $i$-th parent frequency, $n_{i}\in\mathbb{Z}$, and where $\delta\nu$ is
the frequency detuning, for which $|\delta\nu|<\nu_{\rm daughter}$ holds. This
follows the parent-daughter terminology introduced in Degroote et al. (2009b).
To further characterise and verify the direct candidate resonance, and
following Buchler et al. (1997), Vuille (2000) and Vuille & Brassard (2000),
we compute relative phases and amplitudes. The second candidate resonance
identification criterion then consists of verifying whether the corresponding
relative phases $\Phi$ of the combinations that satisfy Eqn. (4) also satisfy
the following (resonant) relative phase relation (see e.g., Fig. $13$ in Van
Hoolst 1995, in which $\Gamma^{\,0}$ and $\sigma$ denote the FP relative phase
and frequency detuning, respectively, as well as Fig. $1$ in Buchler et al.
1997):

$\exists k\in\left[-4,\,-3,\,\ldots,\,4\right]:\hskip
4.26773pt\Phi\equiv\phi_{\rm daughter}\,-\sum_{i}n_{i}\,\phi_{\rm
parent,\,i}\approx k\cdot\dfrac{\pi}{2}\,,$ (5)

where $\phi_{\rm daughter}$ is the extracted phase of the daughter frequency,
and $\phi_{\rm parent,\,i}$ is the extracted phase of its $i$-th parent
frequency. If for a given candidate combination of parent(s) and daughter
frequencies Eqns (4) and (5) are satisfied, this pair is classified as a
candidate resonance. We define the combination order $o$ of this candidate
resonance as: $o=\sum_{i}\left|n_{i}\right|$. Following Vuille & Brassard
(2000) we compute the relative amplitude of the candidate resonance, $A_{r}$,
in the following way:

$A_{r}=\dfrac{A_{\rm daughter}}{\displaystyle\prod_{i}\left(A_{\rm
parent,\,i}\right)^{\left|n_{i}\right|}}\,.$ (6)

In practice, we limit ourselves to identifying up to second-order combinations
of parent frequencies ($o=2$). For each potential candidate direct resonance,
we then verify whether Eqn. (4) is satisfied to within 1 Rayleigh limit
$\mathfrak{R}_{\nu}$ at the 1-$\sigma$ level:

$\exists\nu_{\rm daughter}:\hskip 5.69046pt\left|\nu_{\rm
daughter}-\nu\right|-\sigma_{\nu_{\rm
daughter}}-\sigma_{\nu}\leq\mathfrak{R}_{\nu}\,,$ (7)

where $\nu\equiv\sum_{i}n_{i}\,\nu_{\rm parent,\,i}$ is a combination
frequency composed of other extracted frequencies, with propagated uncertainty
$\sigma_{\nu}$. The comparison with $\mathfrak{R}_{\nu}$ hence measures
whether $\nu$ is resolved with respect to $\nu_{\rm daughter}$. We use this
practical approximation for assigning candidate resonance identifiers, to
avoid use of the theoretical criterion by Buchler et al. (1997),
$|\delta\omega|/\gamma\lesssim O(1)$, setting the approximate limits of
parameter space for which an FP solution with frequency locking exists, as it
relies on model-dependent nonadiabatic oscillation calculations for $\gamma$,
while we want to establish a data-driven approach. If Eqn. (7) holds, we
further verify whether the relative phase $\Phi$ of the candidate combination
satisfies Eqn. (5) to within 1 $\sigma$:

$\exists k\in\left[-4,\,-3,\,\ldots,\,4\right]:\hskip
5.69046pt\Phi-\sigma_{\Phi}\leq k\cdot\dfrac{\pi}{2}\leq\Phi+\sigma_{\Phi}\,.$
(8)

##  3 Results: light curve regression model capacity

In this section the resulting regression models obtained for each of the five
prewhitening strategies applied to all 38 stars in our sample are compared.
The metrics used for this comparison are the $f_{\rm sv}$ value defined in
Eqn. (3), and the (residual) light curves and their LS periodogram. The
$f_{\rm sv}$ metric measures how well the five regression models describe the
light curve of any given star, where we recall that Eqn. (1) assumes that each
oscillation mode has constant amplitude, frequency, and phase. Under this
assumption and using the criteria in Table 1, we define the best regression
model as the one with the highest $f_{\rm sv}$. Any violation of the assumed
constant amplitudes, phases, or frequencies results in lower values of $f_{\rm
sv}$ along with the occurrence of more harmonics and combination frequencies.

Based on the $f_{\rm sv}$ metric, we divide the stars in our sample into two
pseudoclasses: (i) stars for which the light curve regression models attain
$f_{\rm sv}$ values above $0.9$ after frequency filtering for all five
iterative prewhitening strategies; and (ii) stars for which at least one of
these models attains a value of $f_{\rm sv}<0.9$. For brevity, we instead
refer to strategies ‘having’ $f_{\rm sv}$ values above $0.9$ in the remainder
of the text. This classification accentuates the apparent bimodality in the
distribution of the $f_{\rm sv}$ metric values that can be seen in Fig. 1.
Half of the 38 SPB stars in our sample have $f_{\rm sv}\geq 0.9$ for all five
prewhitening strategies and thus belong to pseudoclass 1, hereafter referred
to as the ‘high-$f_{\rm sv}$’ stars. We discuss their light curve regression
model selection process in subsection 3.1. We additionally subdivide the stars
belonging the second pseudoclass (the ‘low-$f_{\rm sv}$’ stars) into two
additional pseudoclasses based on whether outburst-like features are detected
in the light curves, and thus end up with three pseudoclasses that describe
the whole sample. The light curve regression model selection process for these
two additional pseudoclasses, the high mode density starands and the
outbursting stars, are discussed in subsection 3.2, both guided by particular
example stars. The percentages of explained scaled variance $f_{\rm sv}$
averaged over all stars in our sample, and averaged over the stars belonging
to each of these pseudoclasses are shown in Fig. 2. According to this ‘sample’
metric, high-$f_{\rm sv}$ stars are explained well by the regression models
belonging to the five prewhitening strategies. However, significantly
differing values of $f_{\rm sv}$ can be noted for low-$f_{\rm sv}$ stars. For
the sample of $38$ SPB stars as a whole, the overall highest-$f_{\rm sv}$
regression model is the one delivered by prewhitening strategy $3$, as will be
discussed further below. The best regression model and other properties for
all individual stars in our sample are discussed in Appendix B.

Figure 1: Comparison of the scaled fraction of variance $f_{\rm sv}$ values of
the resulting light curve regression models obtained for the different
prewhitening strategies applied to the stars in the sample, with colouring as
a function of $f_{\rm sv}$. Each cell contains the $f_{\rm sv}$ value as
computed using Eqn. (3), as well as the total number of frequencies extracted
using this approach, $n_{\rm e}$. The number $n_{\rm e}$ and the $f_{\rm sv}$
value are computed after the frequency filtering step has removed any
unresolved frequencies with respect to the LD78 criterion and any signals that
have insignificant amplitudes (at $95$% confidence), as discussed in section
2.2.

Figure 2: Total percentage of explained scaled variance $f_{\rm sv}$ averaged
over all 38 stars in the sample (‘all stars’) for each of the five
prewhitening strategies. This averaged percentage is also displayed for each
of the three pseudoclasses introduced in Section 3.

###  3.1 High-$f_{\rm sv}$ stars

Figure 3: Top: Full light curve of the high-$f_{\rm sv}$ star KIC005941844
(left), and a small part of the light curve indicating its multi-periodic
constant-frequency nature (right). Middle: Comparison of the different stop
criteria used for the five prewhitening strategies throughout the iterative
prewhitening procedure, applied to the light curve of KIC005941844. The BIC
values are relevant for strategies 3, 4, and 5, whereas the S/N value and the
estimated amplitude $A$ expressed in estimated standard deviations
$\hat{\sigma}_{A}$ are relevant for strategies 1 and 2, respectively. Bottom:
Lomb-Scargle periodogram (in orange) obtained after subtracting the light
curve model obtained by prewhitening strategy 3 for KIC005941844, from the
original light curve. The extracted frequencies for this model are displayed
as grey lines with height equal to the extracted amplitude. Note the different
y-axis scale for the inset compared to the main panel.

We discuss the light curve regression model selection process for high-$f_{\rm
sv}$ stars guided by the example of KIC005941844, classified as an SPB star by
McNamara et al. (2012) and analysed by Pedersen (2020) and Pedersen et al.
(2021). Its light curve is well explained by sinusoids of constant amplitude
and frequency, as demonstrated by its high $f_{\rm sv}$ values and the light
curve in the top panel of Fig. 3. The light curves of high-$f_{\rm sv}$ stars
are, in general, (nearly) sinusoidal, because a large degree of variance is
explained by the model generated by Eqn. (1). The most prominent data gaps in
the light curve of KIC005941844 are caused by the loss of Kepler’s CCD module
3, one year into the mission.

The middle panel of Fig. 3 shows a comparison of the values of the metrics
important for the stop criteria of all five strategies, throughout the
iterative prewhitening process. The values of these metrics as a function of
the prewhitening step number indicate the progress of that process before
frequency filtering is carried out, and are henceforth referred to as the
‘iteration progress curves’. If the gradient of the Bayesian information
criterion (BIC) iteration progress curve (strategy 3, 4, or 5) converges to a
value of zero, the regression model is approaching the optimal model according
to the likelihood ratio test (LRT). The LRT compares the nested light curve
regression models of this and the previous step based on the likelihoods
inferred from the BIC values. A significant upward trend in the BIC iteration
progress curve increases the p-value of the LRT to a value greater than the
stop criterion, with the significance of that trend depending on the
difference in BIC values and number of optimisation parameters of the nested
regression models. The S/N iteration progress curve (strategy 1) represents
the signal-to-noise ratio of the frequency extracted during each iterative
prewhitening step. If it falls below the pre-set limit $S/N=4.0$ (Breger et
al. 1993), the stop criterion is triggered and iterative prewhitening stops.
Finally, the $A/\hat{\sigma}_{A}$ iteration progress curve represents the
amplitude of the frequency extracted during each iterative prewhitening step,
in terms of its standard deviation. Iterative prewhitening for strategy 2 is
stopped if the amplitude of the frequency extracted is consistent with $0.0$
ppm at $95$% confidence level (i.e. $A/\hat{\sigma}_{A}\leq 1.96$).

The iteration progress curves provide a metric that dictates whether a light
curve regression model is nearing its optimal parameters. This needs to be
investigated on a star-to-star basis, because it intrinsically depends on the
oscillation mode density in the (residual) LS periodograms and amplitudes
attained by any of the periodic signals in the light curves. For KIC005941844,
the iteration progress curves indicate that for all five iterative
prewhitening strategies, the optimal or near-optimal model has been obtained,
before any additional post-processing frequency filtering. The frequency
filtering step removes unresolved frequencies and frequencies that do not have
significant amplitudes (at a $95$% confidence level), so that the final light
curve regression models are obtained. The numbers of frequencies retained in
the regression models after this frequency filtering step (i.e. the numbers of
final resolved significant extracted frequencies) are indicated in Fig. 1 by
the $n_{e}$ parameter for each strategy applied to the light curve of each
star.

The differences among the $f_{\rm sv}$ values of the regression models for all
five strategies are small for the high-$f_{\rm sv}$ stars, as can be noticed
when comparing these values listed for KIC005941844 in Fig. 1. Yet, the
numbers of extracted frequencies, $n_{\rm e}$, differ strongly. Good
correspondence with respect to the number of extracted frequencies for
strategies 2 and 3 occurs for almost all high-$f_{\rm sv}$ stars. For
strategies 1 and 4, the parameter hinting selects the signal in the residual
LS periodogram with the highest S/N value. Because we analyse $g$-mode
oscillators, the parameter hinting for these strategies renders them prone to
select a new frequency in a mode-scarce region of the periodogram, which is
illustrated in Fig. 4, with the mode-scarce regions typically found at higher
frequencies. Signals in mode-scarce regions of a typical LS periodogram of a
$g$-mode oscillator are of lower amplitude than those in mode-dense regions.
These signals therefore contribute less to $f_{\rm sv}$, explaining why for
many high-$f_{\rm sv}$ stars, including KIC005941844, the stop criterion is
triggered earlier for strategy 4 than for the other strategies. Conversely,
the parameter hinting for strategies 2, 3, and 5 selects the largest-amplitude
signal, and thus targets the mode-dense regions of such a typical periodogram.
This hence reiterates the important role that mode density plays in frequency
extraction, even for Kepler light curves, which offer the most optimal
frequency resolution for space photometry available today. Moreover, because
different strategies focus their extraction efforts in different frequency
regimes of the periodogram, the physical information extracted from these
differing light curve regression models is different and can be complementary.

Figure 4: Distribution of the amplitudes and frequencies of (independent)
signals not identified as combination frequencies for all 38 SPB stars in our
sample, for each of the prewhitening strategies defined in Table 1. The colour
indicates if these signals are involved in a candidate resonance (‘red’) or
not (‘blue’), and its gradient indicates the number of detected frequencies:
colours are brighter if more signals are detected within a specific amplitude-
frequency bin. The ‘red and blue distributions’ are plotted independently, so
that overlap between these distributions is the cause for the appearance of
darker, maroon-like coloured bins belonging to the ‘red distribution’. The
maximal frequency displayed in each panel is $8$ d-1, which is the highest
frequency of any of the detected (robust two-signal) candidate resonances (see
Section 4, and especially, Fig. 10). We display independent signals with an
amplitude $>10$ ppm because very few independent signals have amplitudes
smaller than this limit.

The effect of a high mode density in the periodogram of KIC005941844 is
clearly visible in the bottom panel of Fig. 3: dozens of small-amplitude
signals persist in the residual LS periodogram, which are unresolved with
respect to the LD78 criterion. Such unresolved signals can cause modulation of
amplitudes and frequencies because they can interfere with other modes (e.g.,
Bowman et al. 2016) or be residual artefacts from imperfect prewhitening under
the assumption that all oscillations have constant frequencies, amplitudes and
phases. In this case the amplitudes of the extracted frequencies are large
compared to the amplitudes of these unresolved signals, hence $f_{\rm sv}$
values are not strongly affected by such effects and remain high.

A sometimes neglected point in frequency analysis is that some of the free
parameters in the fit to the light curve covary with other free parameters in
the fit. Particularly strong covariances exist between the extracted
frequencies and their corresponding phases. As an example, neglecting such
frequency-phase covariances leads to differences in the frequency values of
the two largest-amplitude extracted frequencies of KIC005941844 at
$2$-$\sigma$ level when comparing their amplitudes for the same mode
frequencies extracted with strategy 5. Covariances among regression parameters
not only influence the estimation of relative phases that are crucial to
identifying candidate resonances, but can also exert influence on
asteroseismic modelling relying on the frequency values (Bowman & Michielsen,
submitted).

From a purely mathematical point of view, parsimony dictates that any nested
regression model that explains the same scaled fraction of variance as a more
complex model (i.e. from the same model strategy but containing more free
parameters) should be preferred over that more complex model, because of bias-
variance trade-off (Claeskens & Hjort 2008). For this reason, one can use
information criteria for model selection when the regression models are
nested, as per each model strategy separately. However, one cannot use such
criteria to distinguish between the capacity of the different regression
models based on different stopping criteria and optimisation methods.
Moreover, the aim of any asteroseismic data analysis is to deduce as many
significant and resolved oscillation frequencies as possible (i.e. we wish to
maximise $n_{\rm e}$ as listed in Fig. 1). Since $n_{\rm e}$ is a discrete
unknown that we want to maximise, within the adopted conditions on frequency
resolution and significance of the modes’ amplitudes, one should not use
classical statistical model selection criteria to choose among the strategies
1 to 5. The case of KIC005941844 illustrates this: its five $f_{\rm sv}$
values differ by less than $\sim 6\cdot 10^{-3}$ but the different
prewhitening strategies lead to vastly different $n_{\rm e}$ values. A
classical model selection procedure via a BIC or LRT applied to the models
based on the different strategies by just using their $f_{\rm sv}$ alone would
fail, as it would favour the model obtained by strategy 4 for KIC005941844.
However, it is clear that the other four regression models with the vastly
larger numbers of extracted modes capture more astrophysical information on,
for example, the number of excited modes, potential mode interactions or
rotational multiplets. Hence, the aim of finding a good regression model (high
$f_{\rm sv}$) along with a maximum number of resolved oscillation frequencies
of significant amplitude (high $n_{\rm e}$) favour strategies 3, 4 or 5 for
this star. These three regression models preferentially select signals at low
frequencies, as shown in Fig. 4. The interplay between the physical and
(purely) mathematical considerations in selecting the optimal iterative
prewhitening model is challenging. For the SPB stars with high $f_{\rm sv}$,
Fig. 2 shows all five prewhitening strategies to be almost equivalent by just
looking at the summed $f_{\rm sv}$ for the 19 stars. However, adding the quest
to maximise the number of resolved modes gives preference to strategies 2 or
3. In any case, KIC005941844’s large-amplitude oscillations and small-
amplitude signals unresolved with respect to the LD78 resolution, as well as
its high $f_{\rm sv}$ values, make it a prime target for nonlinear
asteroseismic analysis.

###  3.2 Low-$f_{\rm sv}$ stars

Figure 5: Box-and-whisker plot of the rotation rates of 26 SPB stars obtained
from averaging over 8 theories of core-boundary and envelope mixing by
Pedersen et al. (2021). These are grouped into the two classes defined in this
work: low and high $f_{\rm sv}$, as indicated schematically on the x-axis.
Whiskers are extended past the quartiles by maximally $1.5$ times the
interquartile range. A single outlier for the low $f_{\rm sv}$ class is
indicated with a rhombus mark. Individual rotation rates and estimated
standard deviations are indicated in grey. The horizontal offset for these
rotation rates was imposed for reasons of visibility without any additional
meaning.

High mode densities are observed in the LS periodograms of all low-$f_{\rm
sv}$ stars, because dozens of potential signals are left unresolved with
respect to the LD78 resolution. Six of the nineteen low-$f_{\rm sv}$ stars
furthermore display outburst-like features: marked departures from the
baseline brightness of the original light curve, which are in phase with
departures from the baseline in the residual light curve. We make a specific
distinction between those stars that display these features in their light
curves, and stars that do not. The former are referred to as ‘outbursting
stars’, whereas the latter are called ‘high mode density stars’. The
discussion of the light curve regression model selection for low-$f_{\rm sv}$
stars is therefore guided by examples from both types of low-$f_{\rm sv}$
stars and the overall performance of the regression models in Fig. 2 is split
up accordingly.

Oscillation mode density is expected to increase with increasing rotation rate
due to rotational frequency shifts. No statistically significant difference is
however found between the rotation rates of 26 of the high- and low-$f_{\rm
sv}$ stars, derived by Pedersen et al. (2021), as is illustrated by the box
plot in Fig. 5. Weighted averages, on the other hand, indicate a difference,
as illustrated in Tables 2 and 3. The rotation rate for the high-$f_{\rm sv}$
star we discussed in the previous subsection is located around the centre of
the respective class rotation rate distributions, whereas the rotation rates
of the the low-$f_{\rm sv}$ stars we discuss in this subsection are at the
lower boundaries of their distributions (see Tables 2 and 3 for the explicit
values). The rotation rate of the high mode density star example,
KIC008255796, is in fact an outlier in the rotation rate box-and-whisker plot
for the low-$f_{\rm sv}$ class, next to the slowly-rotating SPB KIC008459899,
which is suspected to be a double lined spectroscopic binary due to slightly
larger O-C residuals of the fit to its observed spectrum (Lehmann et al.
2011). The latter star is the only (suspected) binary within the pseudoclass
of low-$f_{\rm sv}$ stars. Whether and how its binary nature is correlated
with its low rotation rate and low $f_{\rm sv}$ is unknown. With the exception
of this binary and KIC008255796, the low $f_{\rm sv}$ values of the
low-$f_{\rm sv}$ stars thus seem to be connected to rapid rotation.

Table 2: Rotation rates obtained from averaging over eight theories of core-
boundary and envelope mixing by Pedersen et al. (2021) for the 19 high-$f_{\rm
sv}$ stars, if applicable. The prewhitening strategy that delivers the
highest-$f_{\rm sv}$ light curve regression model among the models generated
by all five strategies, is noted in the ‘Strategy’ column. The (approximate)
amplitude of the largest-amplitude unresolved signal in the LS periodogram of
the light curve subtracted by the highest-$f_{\rm sv}$ regression model is
noted in the $A_{\rm unres/res}$ column, where it is symbolically divided by
the amplitude of the largest-amplitude resolved signal in the LS periodogram
of the original light curve.

KIC |  $\Omega_{\rm rot}/\Omega_{\rm crit}$ (%) | Strategy |  $A_{\rm unres/res}$ (ppm)
---|---|---|---
$003459297$ | $55\pm 6$ | $2$ | $\sim 135/4600$
$003839930$ | / | $2$ | $\sim 330/10000$
$003865742$ | $88\pm 10$ | $2$ | $\sim 160/8750$
$\mathbf{004930889}$ | $60\pm 12$ | $3$ | $\sim 310/8000$
$004936089$ | $16\pm 7$ | $3$ | $\sim 80/5700$
$005084439$ | / | $2$ | $\sim 175/3650$
$005309849$ | $46\pm 14$ | $3$ | $\sim 65/3400$
$005941844$ | $39\pm 4$ | $3$ | $\sim 300/14000$
$\mathbf{006352430}$ | $50\pm 11$ | $2$ | $\sim 350/7400$
$006780397$ | $70\pm 24$ | $3$ | $\sim 26/440$
$007760680$ | $25\pm 5$ | $5$ | $\sim 390/9800$
$008324482$ | / | $2$ | $\sim 175/2000$
$008714886$ | $60\pm 18$ | $2$ | $\sim 155/2150$
$008766405$ | $77\pm 10$ | $3$ | $\sim 110/1750$
$009020774$ | $51\pm 10$ | $2$ | $\sim 75/7700$
$009227988$ | / | $3$ | $\sim 170/2200$
$010220209$ | / | $2$ | $\sim 44/1300$
$010526294$ | $4\pm 3$ | $3$ | $\sim 2000/14500$
$010658302$ | / | $3$ | $\sim 620/3000$
Average |  $29$ aa$a$Average weighted rotation rate, where the inverse variances ($1/\sigma^{2}_{\Omega_{\rm rot}/\Omega_{\rm crit}}$) are used as weights. | N.A. | $300/5807$
222 Known binary stars (Pápics et al. 2013, 2017) have their KIC identifiers
marked in boldface.

Table 3: Same as Table 2, but for the 19 low-$f_{\rm sv}$ stars in our sample.

KIC |  $\Omega_{\rm rot}/\Omega_{\rm crit}$ (%) | Strategy |  $A_{\rm unres/res}$ (ppm)
---|---|---|---
$003240411$ | $67\pm 5$ | $3$ | $\sim 500/570$
$003756031$ | / | $3$ | $\sim 130/690$
$004939281$**$*$Outbursting star. | $78\pm 17$ | $3$ | $\sim 350/1850$
$006462033$ | $85\pm 6$ | $2$ | $\sim 60/600$
$007630417$ | $66\pm 16$ | $2$ | $\sim 160/540$
$008057661$ | $59\pm 16$ | $3$ | $\sim 50/310$
$008087269$ | / | $3$ | $\sim 175/1150$
$008255796$ | $3\pm 5$ | $3$ | $\sim 280/970$
$008381949$ | $76\pm 6$ | $3$ | $\sim 135/650$
$\mathbf{008459899}$ | $8\pm 4$ | $2$ | $\sim 145/640$
$009715425$**$*$Outbursting star. | $83\pm 9$ | $1$ | $\sim 2400/4700$
$009964614$ | / | $2$ | $\sim 130/520$
$010285114$**$*$Outbursting star. | / | $3$ | $\sim 380/730$
$010536147$ | $67\pm 8$ | $3$ | $\sim 145/2075$
$011152422$**$*$Outbursting star. | / | $3$ | $\sim 135/540$
$011360704$**$*$Outbursting star. | $96\pm 4$ | $3$ | $\sim 1225/2075$
$011454304$ | / | $3$ | $\sim 350/535$
$011971405$**$*$Outbursting star. | $72\pm 8$ | $2$ | $\sim 950/3500$
$012258330$ | $73\pm 9$ | $5$ | $\sim 155/500$
Average | $57$ | N.A. | $441/1218$
\- Outb. aa$a$Average sample properties for the outbursting stars. | $90$ | N.A. | $907/2232$
\- H.m.d. bb$b$Average sample properties for the high mode density stars. | $44$ | N.A. | $186/750$
333 A suspected binary star (Lehmann et al. 2011) has its KIC identifier
marked in bold. The weighted averages for the rotation rates are computed in
the same way as in Table 2. The numbers in the ‘Strategy’ column denote the
strategy for which the highest $f_{\rm sv}$ value is obtained.

####  3.2.1 High mode density stars

Figure 6: Same as Fig. 3, but for high mode density star KIC008255796.

Figure 7: Same as Fig. 3, but for outbursting star KIC011971405.

An example of a high mode density star, a low-$f_{\rm sv}$ star that does not
display outburst-like features in its light curve, is KIC008255796. It was
first identified as a misclassified B star by Zhang et al. (2018). Its light
curve is shown in the top panel of Fig. 6 and displays periodic variations on
longer time scales than KIC005941844, the high-$f_{\rm sv}$ star example. The
characteristic properties of the light curves of high mode density stars such
as KIC008255796 are not dissimilar to those of high-$f_{\rm sv}$ stars. The
lower $f_{\rm sv}$ values attained by high mode density stars can be
attributed to the fact that on average maximal excursions from the baseline
brightness are smaller, so that small-amplitude unresolved signals have a
larger impact on the $f_{\rm sv}$ value. The maximal excursion from the light
curve baseline for KIC005941844 is $\sim 35$ ppt, whereas this value is only
$\sim 6$ ppt for KIC008255796, as can be deduced from the top panels of Fig. 3
and Fig. 6.

The iteration progress curves of high mode density stars are also not
dissimilar to those of high-$f_{\rm sv}$ stars, as is illustrated by comparing
the middle panels of Figs 3 and 6, obtained for KIC005941844 and KIC008255796,
respectively. The strategy 2 and 3 iteration progress curves of KIC008255796
indicate that more frequencies need to be extracted to reach the ‘optimal
model’ for that strategy. This is not unexpected, because these strategies
focus their extraction efforts on high mode density regions of the
periodogram, which are strongly present in the LS periodograms of high mode
density stars. The other strategies have triggered their stop criteria, and
those in which hinting is S/N-based are especially strict, as expected.

The eponymous high mode density within the periodogram is well-illustrated in
the bottom panel of Fig. 6 for KIC008255796. Higher mode densities cause more
signals to be unresolved with respect to the LD78 criterion, and therefore
increase the probability that a larger-amplitude signal is unresolved.
Combined with the on average smaller amplitude signals observed in high mode
density stars compared to high-$f_{\rm sv}$ stars, which are a corollary of
the observed differences in maximal excursions from the light curve baseline,
a larger influence of unresolved signals on $f_{\rm sv}$ is expected. Indeed,
as can be derived from Tables 2 and 3, the average unresolved/resolved
amplitude for high mode density stars is $\sim 186/750$ ppm, whereas that
value is $\sim 300/5807$ ppm for high-$f_{\rm sv}$ stars.

The similarities noted in the light curve and Fourier properties of high mode
density low-$f_{\rm sv}$ stars and high-$f_{\rm sv}$ stars render their
classification somewhat artificial. This is why we refer to the high- and
low-$f_{\rm sv}$ stars as pseudoclasses. It is however expected that amplitude
modulation (and frequency modulation) more frequently occurs for high mode
density stars, as they on average exhibit more unresolved signals. Whether
that implies that the observed differences in $f_{\rm sv}$ between the two
pseudoclasses are not of a purely instrumental origin (i.e. not being able to
resolve high mode densities due to limited time span of the photometry),
remains to be validated. What is clear, however, is that the amplitude
modulation observed in some high mode density stars could originate from
beating of close frequency pairs (Bowman et al. 2016), or nonlinear mode
interactions in the ‘intermediate’ regime (e.g., Buchler et al. 1997; Goupil
et al. 1998). Hence, some high mode density stars also are prime targets for
in-depth nonlinear asteroseismic analysis.

####  3.2.2 Outbursting stars

An example of a low-$f_{\rm sv}$ star that displays the eponymous outburst-
like features in its light curve, an outbursting star, is KIC011971405. It is
considered to be a prototypical oscillating Be star by Kurtz et al. (2015),
and was first identified by McNamara et al. (2012) as an SPB oscillator.
Pápics et al. (2017) spectroscopically confirmed it as a Be star, and studied
its variability and outbursts in detail. The days-long outburst-like features
in its light curve are readily visible in the upper panel of Fig. 7. Similar
features can be observed in the light curves of the other outbursting stars,
with differing lengths and amplitudes. However, as illustrated in the same
panel of Fig. 7, the parts of the light curve of KIC011971405 away from these
outburst-like features are very similar to that of a non-outbursting SPB
oscillator. This similarity is also noted for the other outbursting stars, and
fits within the angular momentum transport model put forward by Neiner et al.
(2020) that is used to describe Be star outbursts. Within this model gravito-
inertial modes stochastically excited by core convection dominate the LS
periodogram at the time of outbursts, whereas $g$ modes excited by the
$\kappa$ mechanism prevail for the entire time span of the light curve. Hence,
amplitude and frequency modulation is expected for all outbursting stars,
which gives rise to high mode densities in the periodogram. For KIC011971405
modulation has been detected and characterised by Pápics et al. (2017).
Interestingly, while for five of six outbursting stars these high mode
densities are distinctly grouped, KIC011152422 displays a high mode density in
its LS periodogram that spans from $0$ to $1$ d-1.

The S/N stop criteria trigger an early stop to the prewhitening (i.e. before
200 frequencies are extracted) for almost all outbursting stars, an example of
which can be observed for KIC01197405 in the middle panel of Fig. 7. High mode
densities are the main cause for these early triggers. The exception is
KIC010285114, whose outburst amplitude, in comparison to the typical
excursions from the baseline in its light curve, is the smallest for all
outbursting stars. The stop criteria of strategies 4 and 5 trigger early
abortion of the prewhitening for all outbursting stars, whereas the iteration
progress curves for strategies 2 and 3 indicate more frequencies can be
extracted.

The high mode densities in the LS periodograms and the outburst-like features
in the light curves are also the prime factors to which the low $f_{\rm sv}$
values for outbursting stars can be attributed. Because of these high mode
densities we do not expect any significant changes in $f_{\rm sv}$ after
incorporating additional frequencies extracted by strategies 2 and 3 for
KIC01197405, and the other outbursting stars. The observed amplitude
modulation and outbursts in the outbursting stars in our sample could be
linked to many $g$ modes coming into phase with each other, some of which by
nonlinear interaction (Kurtz et al. 2015). If this is the sole origin of the
outbursts, it corroborates the idea that oscillating Be stars are complex
analogues of SPB stars (Aerts 2006), and supports the idea of the
artificiality of the high- and low-$f_{\rm sv}$ pseudoclasses that was
discussed in the previous subsection. The Kurtz et al. (2015) theory relies
upon a light curve model similar to the one utilised in this work, consisting
of (interacting) sinusoids of constant amplitude and frequency. Kurtz et al.
(2015) show that the bulk of the variability of KIC011971405 can be modelled
in such a manner. Whether such a light curve model is sufficient to explain
the bulk of the observed modulation for all of the outbursting B stars
observed by Kepler remains to be verified, but this is outside of the scope of
this work. As seen from Figs 1 and 2, the choice in prewhitening strategy
clearly exerts a stronger influence for outbursting stars than for the other
SPB stars.

##  4 Results: candidate resonance characteristics

We now discuss the results of our search for candidate nonlinear $g$-mode
resonances within the light curve models obtained for the prewhitening
strategies applied to all stars in our sample. The degree to which nonlinear
interactions are important increases with pulsation amplitude. We therefore
limit our search for candidate resonances to the combination frequencies that
include as a parent frequency at least one of the two largest-amplitude
oscillation signals detected in the light curves, or both. We hereafter refer
to such candidate resonances as ‘two-signal’ candidate resonances. We also
verify whether the number of two-signal candidate resonances is correlated to
the binary or outbursting nature of the SPB stars, or their $f_{\rm sv}$
pseudo-classification. More candidate resonances that do not involve these two
largest-amplitude signals may be detectable. We discuss these additional
candidate resonances for each star in the sample in Appendix B, in which we
also indicate the priority for nonlinear asteroseismic follow-up.

Table 4: Number of two-signal candidate resonances for each of the five
prewhitening strategies applied to all stars in our sample, and the number of
two-signal candidate resonances defined as ‘Robust’, i.e. fulfilling Eqn. (9).

KIC | Prewhitening strategy
---|---
1 | 2 | 3 | 4 | 5 | Robust
$\mathbf{003240411}$ |  $2$/$3$ |  $2$/$17$ |  $\mathbf{1}$/$\mathbf{16}$ |  $0$/$0$ |  $0$/$0$ |  $0$/$0$
$003459297$ |  $3$/$3$ |  $\mathbf{3}$/$\mathbf{11}$ |  $4$/$10$ |  $3$/$3$ |  $2$/$2$ |  $0$/$0$
$\mathbf{003756031}$ |  $6$/$7$ |  $5$/$20$ |  $\mathbf{5}$/$\mathbf{20}$ |  $0$/$0$ |  $2$/$13$ |  $1$/$6$
$003839930$ |  $5$/$5$ |  $\mathbf{4}$/$\mathbf{16}$ |  $4$/$16$ |  $3$/$3$ |  $1$/$1$ |  $1$/$1$
$003865742$ |  $4$/$7$ |  $\mathbf{1}$/$\mathbf{14}$ |  $1$/$13$ |  $4$/$4$ |  $1$/$16$ |  $0$/$6$
$004930889$ |  $4$/$6$ |  $2$/$9$ |  $\mathbf{2}$/$\mathbf{9}$ |  $1$/$1$ |  $1$/$6$ |  $0$/$4$
$004936089$ |  $4$/$7$ |  $2$/$7$ |  $\mathbf{2}$/$\mathbf{7}$ |  $1$/$1$ |  $4$/$6$ |  $2$/$2$
$\mathit{004939281}$ |  $1$/$2$ |  $2$/$12$ |  $\mathbf{2}$/$\mathbf{12}$ |  $0$/$0$ |  $0$/$4$ |  $0$/$0$
$005084439$ |  $8$/$14$ |  $\mathbf{9}$/$\mathbf{27}$ |  $9$/$27$ |  $10$/$20$ |  $7$/$9$ |  $6$/$6$
$005309849$ |  $10$/$10$ |  $4$/$7$ |  $\mathbf{4}$/$\mathbf{7}$ |  $7$/$7$ |  $3$/$6$ |  $3$/$5$
$005941844$ |  $12$/$14$ |  $9$/$18$ |  $\mathbf{9}$/$\mathbf{18}$ |  $5$/$5$ |  $5$/$15$ |  $5$/$12$
$006352430$ |  $1$/$1$ |  $\mathbf{6}$/$\mathbf{16}$ |  $6$/$16$ |  $0$/$0$ |  $5$/$15$ |  $5$/$12$
$\mathbf{006462033}$ |  $2$/$2$ |  $\mathbf{2}$/$\mathbf{8}$ |  $2$/$8$ |  $0$/$0$ |  $1$/$1$ |  $0$/$0$
$006780397$ |  $6$/$7$ |  $2$/$12$ |  $\mathbf{2}$/$\mathbf{12}$ |  $0$/$0$ |  $0$/$10$ |  $0$/$8$
$\mathbf{007630417}$ |  $2$/$3$ |  $\mathbf{2}$/$\mathbf{4}$ |  $2$/$4$ |  $3$/$3$ |  $3$/$5$ |  $2$/$4$
$007760680$ |  $4$/$5$ |  $3$/$3$ |  $3$/$3$ |  $0$/$0$ |  $\mathbf{5}$/$\mathbf{8}$ |  $3$/$4$
$\mathbf{008057661}$ |  $2$/$2$ |  $0$/$6$ |  $\mathbf{0}$/$\mathbf{6}$ |  $2$/$2$ |  $0$/$6$ |  $0$/$4$
$\mathbf{008087269}$ |  $1$/$3$ |  $2$/$9$ |  $\mathbf{4}$/$\mathbf{9}$ |  $2$/$2$ |  $2$/$7$ |  $2$/$4$
$\mathbf{008255796}$ |  $0$/$3$ |  $0$/$17$ |  $\mathbf{0}$/$\mathbf{17}$ |  $0$/$0$ |  $0$/$14$ |  $0$/$10$
$008324482$ |  $6$/$6$ |  $\mathbf{4}$/$\mathbf{13}$ |  $4$/$13$ |  $2$/$2$ |  $3$/$4$ |  $2$/$3$
$\mathbf{008381949}$ |  $1$/$1$ |  $2$/$16$ |  $\mathbf{2}$/$\mathbf{16}$ |  $0$/$0$ |  $0$/$2$ |  $0$/$2$
$\mathbf{008459899}$ |  $1$/$4$ |  $\mathbf{3}$/$\mathbf{29}$ |  $3$/$29$ |  $0$/$0$ |  $1$/$3$ |  $1$/$2$
$008714886$ |  $6$/$6$ |  $\mathbf{6}$/$\mathbf{17}$ |  $6$/$17$ |  $1$/$1$ |  $3$/$11$ |  $3$/$11$
$008766405$ |  $5$/$5$ |  $7$/$16$ |  $\mathbf{7}$/$\mathbf{16}$ |  $5$/$5$ |  $7$/$16$ |  $6$/$12$
$009020774$ |  $4$/$5$ |  $\mathbf{4}$/$\mathbf{11}$ |  $4$/$11$ |  $4$/$5$ |  $1$/$2$ |  $0$/$1$
$009227988$ |  $3$/$4$ |  $1$/$13$ |  $\mathbf{1}$/$\mathbf{13}$ |  $3$/$4$ |  $2$/$13$ |  $0$/$2$
$\mathit{009715425}$ |  $\mathbf{3}$/$\mathbf{5}$ |  $1$/$4$ |  $1$/$4$ |  $1$/$1$ |  $2$/$7$ |  $1$/$1$
$\mathbf{009964614}$ |  $1$/$3$ |  $\mathbf{0}$/$\mathbf{13}$ |  $0$/$13$ |  $0$/$0$ |  $0$/$10$ |  $0$/$2$
$010220209$ |  $2$/$5$ |  $\mathbf{2}$/$\mathbf{8}$ |  $2$/$8$ |  $1$/$9$ |  $0$/$0$ |  $0$/$0$
$\mathit{010285114}$ |  $6$/$6$ |  $4$/$10$ |  $\mathbf{4}$/$\mathbf{10}$ |  $0$/$0$ |  $0$/$1$ |  $1$/$1$
$010526294$ |  $7$/$8$ |  $7$/$8$ |  $\mathbf{7}$/$\mathbf{8}$ |  $4$/$4$ |  $9$/$10$ |  $7$/$8$
$\mathbf{010536147}$ |  $0$/$1$ |  $1$/$7$ |  $\mathbf{1}$/$\mathbf{7}$ |  $0$/$0$ |  $1$/$11$ |  $1$/$7$
$010658302$ |  $4$/$4$ |  $3$/$9$ |  $\mathbf{4}$/$\mathbf{17}$ |  $1$/$1$ |  $3$/$3$ |  $1$/$1$
$\mathit{011152422}$ |  $3$/$3$ |  $2$/$15$ |  $\mathbf{2}$/$\mathbf{16}$ |  $3$/$6$ |  $2$/$13$ |  $2$/$4$
$\mathit{011360704}$ |  $6$/$7$ |  $\mathbf{3}$/$\mathbf{19}$ |  $3$/$19$ |  $2$/$2$ |  $1$/$11$ |  $1$/$8$
$\mathbf{011454304}$ |  $3$/$5$ |  $3$/$11$ |  $\mathbf{3}$/$\mathbf{11}$ |  $0$/$0$ |  $2$/$2$ |  $1$/$1$
$\mathit{011971405}$ |  $5$/$9$ |  $\mathbf{3}$/$\mathbf{16}$ |  $4$/$17$ |  $2$/$2$ |  $1$/$2$ |  $0$/$0$
$\mathbf{012258330}$ |  $5$/$5$ |  $4$/$10$ |  $4$/$11$ |  $3$/$3$ |  $\mathbf{5}$/$\mathbf{8}$ |  $3$/$6$
444 Result format: S/N $\geq 5$ restriction imposed / no S/N restriction
imposed. Boldface for the overall highest-$f_{\rm sv}$ strategy.
KIC number style: outbursting stars - italics, high mode density stars -
boldface, high-$f_{\rm sv}$ stars - regular.

Figure 8: Left panel: box-and-whisker plot of the number of identified two-
signal candidate resonances, grouped into the two pseudoclasses defined in
this work: low and high $f_{\rm sv}$. The colour depicts the prewhitening
strategy employed to generate the light curve models. Right panel: box-and-
whisker plot of the number of identified robust two-signal candidate
resonances (defined by Eqn. (9)), grouped by pseudoclass. Whiskers and
outliers for both panels are indicated in a similar way to how they are
indicated in Fig. 5.

The strategy discussed in subsection 2.3 is employed to identify two-signal
candidate resonances in the frequency-locked regime of mode interaction for
the models generated by all five strategies for each star. To do so, we verify
whether the highest-$f_{\rm sv}$ strategy that performs nonlinear least-
squares optimisation of the frequencies (i.e. one of strategies $1$-$4$)
contains unique matches in frequency to within 1 Rayleigh limit at the
1-$\sigma$ level, with the frequencies in the linear regression model
generated by prewhitening strategy 5:

$\exists!\nu_{f_{\rm sv}}:\hskip 5.69046pt\left|\nu_{\rm lin}-\nu_{f_{\rm
sv}}\right|-\sigma_{\nu_{\rm lin}}-\sigma_{\nu_{f_{\rm
sv}}}\leq\mathfrak{R}_{\nu}\,.$ (9)

In Eqn. (9) $\nu_{f_{\rm sv}}$ denotes the frequency of the model that attains
the highest-$f_{\rm sv}$ among all strategies that employ nonlinear least-
squares optimisation, and $\nu_{\rm lin}$ denotes the frequency obtained for
prewhitening strategy $5$. We refer to candidate resonances for which Eqn. (9)
is fulfilled as ‘robust’ candidate resonances and discuss how their occurrence
covaries with (pseudo)classification in subsection 4.1. We also discuss how
this is reflected in the characteristic properties of these resonances in
subsection 4.2.

Baran et al. (2015) and Zong et al. (2016b) note that the Breger et al. (1993)
S/N detection threshold of $4.0$ is low as a detection threshold for signals
in K2 data and Kepler short-cadence data, respectively. On the other hand,
Degroote et al. (2009b) and Degroote et al. (2009a) conclude that not all
signals failing the Breger et al. (1993) threshold correspond to noise peaks
(i.e. S/N-based methods are not clearcut in their optimal cut-off value as
significance criterion). For the detection of resonances, we therefore make an
additional distinction between frequencies with S/N $\geq 5$ and S/N $<5$,
when analysing candidate resonances in the following subsections.

###  4.1 Numbers of identified candidate resonances

The numbers of two-signal candidate resonances, with an indication of how many
of those candidate resonances have S/N values $\geq 5$, are displayed in Table
4 for all five strategies applied to each star in our sample. The number of
robust two-signal candidate resonances and its S/N $\geq 5$ counterpart are
also shown. It is apparent that for every star in our sample there is at least
one prewhitening strategy that delivers a final light curve model in which
two-signal candidate resonances are identified. For 32 out of 38 stars ($\sim
84$%) we find robust resonance identifications in the definition in Eqn. (9).
Imposing the additional S/N $\geq 5$ restriction lowers this number to 24 out
of 38 stars ($\sim 63$%). More candidate resonances are identified if we
loosen the restriction on involvement of the two largest-amplitude candidate
resonances. Candidate resonances are thus omnipresent in the light curves of
the SPB stars we consider in this work, and SPB stars in general are therefore
intrinsic nonlinear nonradial oscillators.

The covariance matrices of the light curve models obtained from applying the
prewhitening strategies to all stars in the sample vary from strategy to
strategy and from star to star. The figures of Appendix A display
visualisations of the covariance matrices obtained for the highest-$f_{\rm
sv}$ prewhitening strategy (i.e. the best regression model), in a block-
normalised form, which is discussed in detail in Appendix A.4. An example of
such a visualisation of a covariance matrix is given in Fig. 9, which displays
the covariance matrix of the best regression model of the light curve of the
high-$f_{\rm sv}$ star KIC005941844. Large numbers of weak inter-parameter
covariations and the strong frequency-phase covariations mentioned in Sect.
3.1 can be seen, similar to what is observed in all covariance matrices
displayed in Appendix A.4. These covariations differ from star to star.
Several covariance and correlation matrices display a structured pattern
within their different blocks (i.e. they are structured) that indicates the
strongest inter-parameter covariations and correlations.

Figure 9: Covariance matrix of the best regression model of the light curve of
KIC005914844 (obtained from prewhitening strategy 3) in the block-normalised
form described in Appendix A.4 and Fig. A.1. The colour map shows the
normalised (co)variance values ($S_{N}$)${}_{ij,\,{\rm norm}}$ and the change
in darkness between the upper and lower triangular part of the matrix is only
for visualisation purposes.

The covariance matrix propagates into the relative frequencies and phases
defined in Eqns (4) and (5), which are used to identify a candidate resonance.
The occurrence of robust candidate resonances as indicated in Table 4 thus
connects to the covariance matrix structure among differing light curve
models.

In asteroseismic forward modelling one compares observed mode frequencies with
frequencies obtained from models of stellar evolution. As Davies et al. (2016)
note for asteroseismic modelling based on solar-like oscillations, it is
essential to account for regression parameter covariances and correlations in
this process. Regression parameter covariances and correlations can also
impact asteroseismic forward modelling based on gravity modes, as discussed in
Sect. 3.1. These covariances and correlations also impact nonlinear
asteroseismic modelling, which utilises extracted amplitudes, phases, and
derived quantities such as the frequency detuning $\delta\nu$ and the relative
phase $\Phi$, next to the extracted frequencies. The strong frequency-phase
covariances detected in all covariance matrices displayed in Appendix A.4
therefore render it necessary to include such correlations and covariances in
future nonlinear asteroseismic modelling of the SPB stars in our sample.

Figure 10: Highest-$f_{\rm sv}$ strategy characteristic properties of the
identified robust candidate resonances, which involve the two largest-
amplitude signals, as a function of each other, with $1$-$\sigma$
uncertainties indicated: amplitude $A$, frequency $\omega/2\pi$, frequency
detuning $\delta\nu$ (see Eqn. (4)), relative amplitude $A_{r}$ (see Eqn.
(6)), and relative phase $\Phi$ (see Eqn. (5)). Colours indicate the different
pseudoclasses to which the stars that exhibit the candidate resonances belong.
The transparent symbols indicate robust candidate resonances having $4<$ S/N
$<5$, whereas the full symbols indicate robust candidate resonances for which
S/N $\geq 5$. Grey lines indicate the zero-point for the frequency detuning,
and the faint grey lines in the last row of panels indicate the relative
phases satisfying Eqn. (5).

No differences in the number of identified two-signal candidate resonances are
found when comparing the models of stars belonging to the low- and
high-$f_{\rm sv}$ pseudoclasses, as is shown in the left panel of Fig. 8. This
conclusion is not changed when only considering two-signal candidate
resonances with S/N $\geq 5$. The box-and-whisker plots in that panel
furthermore show that no significant dependence on prewhitening strategy can
be identified. Any additional uncertainty in parameter values due to
covariance among parameters is neglected in the construction of the box-and-
whisker plots. However, including such uncertainties in their construction
would only support the conclusion that no significant dependencies are
observed.

One might suspect that the differing (rich) objective function landscape for
which the global minimum is sought during optimisation is a reason for a
difference between stars belonging to the two pseudoclasses, high and low
$f_{\rm sv}$. Amplitude modulation is more important in shaping this landscape
for low-$f_{\rm sv}$ stars than their high-$f_{\rm sv}$ counterparts, although
this varies from star to star and from signal to signal. However, no
difference in number of identified two-signal robust candidate resonances
(i.e. those that fulfil the condition in Eqn. (9)) can be noted when comparing
the two pseudoclasses, as is shown on the right panel of Fig. 8. This lack of
difference is also noted for robust two-signal candidate resonances with S/N
$\geq 5$. Besides a few exceptions, amplitude modulation of low-$f_{\rm sv}$
star signals thus seems to have only a minor influence on the identification
of signals relevant for (two-signal) robust candidate resonances.

We find no evidence for a distinction in number of identified candidate
resonances when comparing the numbers strategy per strategy for outbursting
and non-outbursting stars, even if we only compare candidate resonances for
which S/N $\geq 5$. A class imbalance exists when making the latter comparison
(8 outbursting stars vs. 30 non-outbursting stars). However, the presence of
robust candidate resonances in stars of all (pseudo)classes supports the lack
of distinction between outbursting and non-outbursting stars. It also lends
support to the interpretation of Be star variability in which the times of
outbursts signify several modes coming into phase (e.g., Kurtz et al. 2015).

Finally, we also check for correlation of the number of candidate resonances
with binarity. Two stars in our sample are confirmed binaries, KIC006352430
(Pápics et al. 2013) and KIC004930889 (Pápics et al. 2017), whereas Lehmann et
al. (2011) suspect KIC008459899 to be a binary as well. All other stars in our
sample do not have confirmed or suspected binary companions. We find no
difference in the number of identified two-signal candidate resonances between
(suspected) binary and single SPB stars, and only including candidate
resonances for which S/N $\geq 5$ does not make a difference.

###  4.2 Characteristic properties of identified candidate resonances

Amplitudes, frequencies, frequency detunings $\delta\nu$, relative amplitudes
$A_{r}$, and relative phases $\Phi$ define candidate resonances. These
characteristic properties are for example discussed by Degroote et al. (2009b)
for a $\beta$ Cep (CoRoT) star, HD 180642, whereas Vuille (2000) and Vuille &
Brassard (2000) discuss them for an oscillating DA white dwarf, G29-38. We
present how these properties covary in Fig. 10 for the identified robust two-
signal candidate resonances of all stars in our sample. The amplitudes and
frequencies of signals involved in these resonances are not significantly
different from those not involved in them, as shown in Fig. 4.

Most robust (two-signal) candidate resonances originate from high $f_{\rm sv}$
stars, as can be seen in Table 4. The panels in the left column of Fig. 10
show that most robust two-signal candidate resonances exhibit frequencies
$\lesssim 2.5$ d-1. Several of those low-frequency ($\omega/2\pi\lesssim 2.5$
d-1) candidate resonances have $4<$ S/N $<5$. The bulk of the high-frequency
robust two-signal candidate resonances ($\omega/2\pi\gtrsim 2.5$ d-1) depicted
in the top panel of Fig. 10 exhibit amplitudes $<100$ ppm, whereas the bulk of
the low-frequency robust two-signal candidate resonances exhibit amplitudes
varying from $\sim 30$ to $\sim 300$ ppm. Only a few of the detected robust
two-signal candidate resonances have amplitudes $\geq 200$ ppm, predominantly
at low frequencies. Similarly, only a few robust two-signal candidate
resonances are detected with amplitudes of the order of $10$ ppm. Few of the
frequency detunings of the robust candidate resonances are significantly
different from $0$ d-1, as can be deduced from the panels in Fig. 10 that
display the frequency detunings ($\delta\nu$). This is expected because of the
way we select robust candidate resonances (see Eqn. (9)). The few outliers we
observe likely either belong to the frequency-locked or to the ‘intermediate’
regime, although definitive conclusions can only be made if the linear growth
rate $\gamma$ of the mode is known, which requires nonlinear asteroseismic
models. No distinct trends are present in any of the panels of Fig. 10 that
display the relative phases ($\Phi$) and the relative amplitudes ($A_{r}$).

Variances for the characteristic properties are larger than average for robust
two-signal candidate resonances that have extreme amplitudes. By comparing the
individual properties of the two-signal candidate resonances for each star in
our sample, we find that the largest-amplitude signals shown in Fig. 10 are
the robust two-signal candidate resonances of KIC009715425 and KIC010526294.
The larger uncertainties in their properties are attributable to the high mode
density for the outbursting star KIC009715425, and the unresolved rotationally
split triplet signals in the light curve of KIC010526294 (Pápics et al. 2014).
The covariance matrices of the best light curve regression models of these
stars are strongly structured, as shown in Figs A.6 and A.7, indicating
significant covariance among some of the regression parameters. For the small-
amplitude signals shown in Fig. 10, the increased variance of the
characteristic properties is explained by the fact that the signals are closer
to the noise level, as is indicated by half of these signals having $4<$ S/N
$<5$.

The approach taken in this section to characterise candidate resonances
contains the conservative approximation of cataloguing only the candidate
resonances influenced by the two largest-amplitude signals. It is thus limited
in capturing the characteristic properties of candidate resonances. This
raises questions that need to be solved:

  1. 1.

Is the conservative approach that only considers robust candidate resonances
defined by Eqn. (9) satisfactory for initial target selection for nonlinear
asteroseismic studies?

  2. 2.

Assuming the conservative approach is satisfactory, how do the number and
characteristic properties of candidate resonances change if more, lower-
amplitude signals are allowed to be involved in the candidate resonances?

  3. 3.

Is it expected from nonlinear stellar oscillation theory that no significant
trends exist among the characteristic candidate resonance properties?

The second question is dealt with in Appendix B, whereas the other questions
require nonlinear asteroseismic models in order to be answered, which we plan
to take up in follow-up studies.

##  5 Conclusions and prospects

We analysed the Kepler light curves of an ensemble of $38$ SPB stars by making
use of five different iterative prewhitening strategies. Amplitude-based
prewhitening strategies tend to extract larger-amplitude signals from the
lower-frequency, mode-dense regions of the LS periodograms in $g$-mode
pulsators, whereas S/N-based prewhitening strategies tend to extract more
smaller-amplitude signals at higher-frequency, mode-scarce regions, as
illustrated in Fig. 4.

The scaled fraction of variance $f_{\rm sv}$, defined in Eqn. (3), is one of
the metrics used to compare the frequency analysis results, in order to select
the optimal model for the light curve. Based on that metric, as well as the
(residual) LS periodogram and light curve, we divide the stars in our sample
in three different pseudoclasses: the high $f_{\rm sv}$, the high mode
density, and the outbursting pseudoclass. The first pseudoclass consists of
$19$ stars for which the models of all five prewhitening strategies attain
$f_{\rm sv}$ values $>0.9$. The two other pseudoclasses are discerned by
detection of outburst-like features. We refer to them as pseudoclasses because
we found that some of the characterising properties of these pseudoclasses are
attributable to the limited frequency resolution that does not allow one to
resolve several of the signals observed in the Kepler light curves. This
limited frequency resolution can be attributed to the high mode densities that
are detected in the LS periodograms, which also leads to modulation of some of
the detected signals. The resonant excitation of envelope $g$ modes with
inertial modes propagating in the core, which requires rapid rotation to be
efficient (Lee 2021), may also be a reason why SPB stars exhibit such mode-
dense LS periodograms.

An important point in the frequency analysis of SPB stars is the covariance
among the fitting parameters. We demonstrate the particularly strong
covariances between extracted frequencies and phases in the figures in
appendix A. These covariances are sometimes neglected in the uncertainty
budget of the light curve models, so that regression parameter uncertainties
are underestimated. It is necessary to include such regression parameter
covariations in future asteroseismic modelling of the SPB stars in our sample.

The role of resonances among $g$ modes was described eloquently by Buchler et
al. (1997). Few to no space-based photometric data were available at the time
to match the mainly theoretical discussions of its importance. Lee (2012)
provided account of the importance of rotation on nonlinear mode interactions
in the rapidly rotating B stars by studying $g$ modes excited in stellar
models near the zero-age main sequence. Applications of this theory were
limited by the lack of available photometric data. Because we now have the
full Kepler data set at our disposal, the research field of nonlinear
asteroseismology applied to main sequence $g$-mode pulsators can be
reinvigorated. Although intermediate-mass dwarfs have primarily been the
target of linear asteroseismic studies for widely differing purposes, such as
their sizing, weighing, ageing (Serenelli et al. 2021), as well as the
characterisation of their internal rotation and mixing (e.g., Pedersen et al.
2021; Michielsen et al. 2021), no observation-based modelling studies have
focused on nonlinear mode interactions.

Given the observed deficit of angular momentum transport during the main
sequence in current stellar evolution models (Aerts et al. 2019), the
importance of nonlinear mode interactions for angular momentum transport
deserves follow-up studies. In this work we provide an observational basis for
this. Based on the light curve models constructed during frequency analysis,
we identified second-order candidate direct resonant Kepler SPB oscillations.
We guided ourselves in this identification process by the frequency detuning
and relative phase relations defined in Eqns (4) and (5), inspired by Van
Hoolst (1995) and Buchler et al. (1997). Future nonlinear asteroseismic models
describing frequency-locked signals based on the amplitude equations, will
make use of such characteristic properties of candidate resonances (as well as
the relative amplitude defined in Eqn. (6)), or derivates thereof.

It is more likely that a resonant $g$ mode interaction occurs if one of the
parent frequencies has a large amplitude. As a first approximation we look at
the number of identified candidate resonances that involve the two largest-
amplitude signals in the light curve, which we refer to as ‘two-signal’
candidate resonances, and the characteristic properties of such resonances. We
find that 32 of the 38 SPB stars turn out to be intrinsic nonlinear $g$-mode
oscillators, reiterating the need for quantification of the nonlinear $g$ mode
interactions.

One should take into account that the practical implementation of the
criterion indicating that the frequency detuning relation is fulfilled is only
approximate, because the growth rates needed to verify the Buchler et al.
(1997) criterion cannot be derived directly from observations. Moreover, as
Goupil et al. (1998) explain, the boundaries of the Buchler et al. (1997)
regimes of resonant mode interaction might not be that distinct for
observational data, and they are undetermined for SPB stars A step towards the
determination of the Buchler et al. (1997) regimes consists of verifying
whether the candidate resonances detected in this work display the oscillatory
behaviour that is deemed characteristic of the ‘intermediate regime’ (Goupil
et al. 1998). This is a next logical step in characterising $g$ mode
resonances in SPB stars.

Our current observational analyses constitute the ground work to perform
forward asteroseismic modelling of the nonlinear mode properties, with the aim
to investigate whether it leads to different or additional results on internal
rotation and mixing compared to the current state-of-the-art based on the
linear asteroseismic approach in stellar modelling of SPB stars (Pedersen et
al. 2021).

###### Acknowledgements.

The research leading to these results received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (grant agreement No. 670519: MAMSIE with PI CA.), and from
the KU Leuven Research Council (grant C16/18/005: PARADISE, with PI CA.). DMB
is grateful for a senior post-doctoral fellowship from the Research Foundation
Flanders (FWO) with grant agreement No. 1286521N. TVR gratefully acknowledges
support from the Research Foundation Flanders (FWO) under grant agreement nr.
12ZB620N. This research was supported in part by the National Science
Foundation under Grant No. NSF PHY-1748958. A part of the computational
resources and services used in this work were provided by the VSC (Flemish
Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and
the Flemish Government – department EWI. JVB thanks J. De Ridder for his
advice on visualising correlations among light curve regression model
parameters, and thanks W. Zong for his useful comments on the manuscript.
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 (Wes McKinney 2010), Matplotlib (Hunter 2007), and Seaborn
(Waskom & the seaborn development team 2020). We thank their authors for
making these great software packages open source.

## References

  * Aerts (2006) Aerts, C. 2006, in ESA Special Publication, Vol. 624, Proceedings of SOHO 18/GONG 2006/HELAS I, Beyond the spherical Sun, ed. K. Fletcher & M. Thompson, 131
  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology
  * Aerts et al. (1999) Aerts, C., De Cat, P., Peeters, E., et al. 1999, A&A, 343, 872
  * Aerts et al. (2019) Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35
  * Aerts et al. (2017) Aerts, C., Van Reeth, T., & Tkachenko, A. 2017, ApJ, 847, L7
  * Auvergne et al. (2009) Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411
  * Balona et al. (2015) Balona, L. A., Baran, A. S., Daszyńska-Daszkiewicz, J., & De Cat, P. 2015, MNRAS, 451, 1445
  * Balona et al. (2011) Balona, L. A., Pigulski, A., De Cat, P., et al. 2011, MNRAS, 413, 2403
  * Baran et al. (2015) Baran, A. S., Koen, C., & Pokrzywka, B. 2015, MNRAS, 448, L16
  * Barceló Forteza et al. (2015) Barceló Forteza, S., Michel, E., Roca Cortés, T., & García, R. A. 2015, A&A, 579, A133
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Bouabid et al. (2013) Bouabid, M. P., Dupret, M. A., Salmon, S., et al. 2013, MNRAS, 429, 2500
  * Bowman et al. (2021) Bowman, D. M., Hermans, J., Daszyńska-Daszkiewicz, J., et al. 2021, MNRAS, 504, 4039
  * Bowman & Kurtz (2014) Bowman, D. M. & Kurtz, D. W. 2014, MNRAS, 444, 1909
  * Bowman et al. (2016) Bowman, D. M., Kurtz, D. W., Breger, M., Murphy, S. J., & Holdsworth, D. L. 2016, MNRAS, 460, 1970
  * Breger et al. (1993) Breger, M., Stich, J., Garrido, R., et al. 1993, A&A, 271, 482
  * Buchler (1993) Buchler, J. R. 1993, Ap&SS, 210, 9
  * Buchler & Goupil (1984) Buchler, J. R. & Goupil, M. J. 1984, ApJ, 279, 394
  * Buchler et al. (1997) Buchler, J. R., Goupil, M. J., & Hansen, C. J. 1997, A&A, 321, 159
  * Buysschaert et al. (2018) Buysschaert, B., Aerts, C., Bowman, D. M., et al. 2018, A&A, 616, A148
  * Claeskens & Hjort (2008) Claeskens, G. & Hjort, N. L. 2008, Model Selection and Model Averaging, Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press)
  * Davies et al. (2016) Davies, G. R., Silva Aguirre, V., Bedding, T. R., et al. 2016, MNRAS, 456, 2183
  * De Cat & Aerts (2002) De Cat, P. & Aerts, C. 2002, A&A, 393, 965
  * De Cat et al. (2000) De Cat, P., Aerts, C., De Ridder, J., et al. 2000, A&A, 355, 1015
  * De Cat et al. (2007) De Cat, P., Briquet, M., Aerts, C., et al. 2007, A&A, 463, 243
  * De Cat et al. (2005) De Cat, P., Briquet, M., Daszyńska-Daszkiewicz, J., et al. 2005, A&A, 432, 1013
  * Degroote et al. (2010) Degroote, P., Aerts, C., Baglin, A., et al. 2010, Nature, 464, 259
  * Degroote et al. (2009a) Degroote, P., Aerts, C., Ollivier, M., et al. 2009a, A&A, 506, 471
  * Degroote et al. (2009b) Degroote, P., Briquet, M., Catala, C., et al. 2009b, A&A, 506, 111
  * Dhouib et al. (2021) Dhouib, H., Prat, V., Van Reeth, T., & Mathis, S. 2021, A&A, in press, arXiv:2104.09302
  * Dziembowski (1982) Dziembowski, W. 1982, Acta Astron., 32, 147
  * Dziembowski & Krolikowska (1985) Dziembowski, W. & Krolikowska, M. 1985, Acta Astron., 35, 5
  * Dziembowski et al. (1993) Dziembowski, W. A., Moskalik, P., & Pamyatnykh, A. A. 1993, MNRAS, 265, 588
  * Gautschy & Saio (1993) Gautschy, A. & Saio, H. 1993, MNRAS, 262, 213
  * Goupil et al. (1998) Goupil, M. J., Dziembowski, W. A., & Fontaine, G. 1998, Baltic Astronomy, 7, 21
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362
  * Henneco et al. (2021) Henneco, J., Van Reeth, T., Prat, V., et al. 2021, A&A, 648, A97
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
  * Kurtz et al. (2015) Kurtz, D. W., Shibahashi, H., Murphy, S. J., Bedding, T. R., & Bowman, D. M. 2015, MNRAS, 450, 3015
  * Lee (2012) Lee, U. 2012, MNRAS, 420, 2387
  * Lee (2021) Lee, U. 2021, MNRAS, 505, 1495
  * Lee & Saio (1997) Lee, U. & Saio, H. 1997, ApJ, 491, 839
  * Lehmann et al. (2011) Lehmann, H., Tkachenko, A., Semaan, T., et al. 2011, A&A, 526, A124
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447
  * Longuet-Higgins (1968) Longuet-Higgins, M. S. 1968, Philosophical Transactions of the Royal Society of London Series A, 262, 511
  * Loumos & Deeming (1978) Loumos, G. L. & Deeming, T. J. 1978, Ap&SS, 56, 285
  * Mathias et al. (2001) Mathias, P., Aerts, C., Briquet, M., et al. 2001, A&A, 379, 905
  * Mathis (2013) Mathis, S. 2013, Transport Processes in Stellar Interiors, ed. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, Vol. 865, 23
  * Mathis & Prat (2019) Mathis, S. & Prat, V. 2019, A&A, 631, A26
  * McNamara et al. (2012) McNamara, B. J., Jackiewicz, J., & McKeever, J. 2012, AJ, 143, 101
  * Michielsen et al. (2021) Michielsen, M., Aerts, C., & Bowman, D. M. 2021, A&A, 650, A175
  * Miglio et al. (2008) Miglio, A., Montalbán, J., Noels, A., & Eggenberger, P. 2008, MNRAS, 386, 1487
  * Montgomery & O’Donoghue (1999) Montgomery, M. H. & O’Donoghue, D. 1999, Delta Scuti Star Newsletter, 13, 28
  * Moravveji et al. (2015) Moravveji, E., Aerts, C., Pápics, P. I., Triana, S. A., & Vandoren, B. 2015, A&A, 580, A27
  * Moravveji et al. (2016) Moravveji, E., Townsend, R. H. D., Aerts, C., & Mathis, S. 2016, ApJ, 823, 130
  * Neiner et al. (2020) Neiner, C., Lee, U., Mathis, S., et al. 2020, A&A, 644, A9
  * Newville et al. (2020) Newville, M., Otten, R., Nelson, A., et al. 2020, lmfit/lmfit-py 1.0.1
  * Ouazzani et al. (2020) Ouazzani, R. M., Lignières, F., Dupret, M. A., et al. 2020, A&A, 640, A49
  * Pamyatnykh (1999) Pamyatnykh, A. A. 1999, Acta Astron., 49, 119
  * Pápics (2012) Pápics, P. I. 2012, Astronomische Nachrichten, 333, 1053
  * Pápics et al. (2012) Pápics, P. I., Briquet, M., Baglin, A., et al. 2012, A&A, 542, A55
  * Pápics et al. (2014) Pápics, P. I., Moravveji, E., Aerts, C., et al. 2014, A&A, 570, A8
  * Pápics et al. (2013) Pápics, P. I., Tkachenko, A., Aerts, C., et al. 2013, A&A, 553, A127
  * Pápics et al. (2015) Pápics, P. I., Tkachenko, A., Aerts, C., et al. 2015, ApJ, 803, L25
  * Pápics et al. (2017) Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, A&A, 598, A74
  * Pedersen (2020) Pedersen, M. G. 2020, PhD thesis, KU Leuven
  * Pedersen et al. (2021) Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nature Astronomy, 5, 715
  * Prat et al. (2019) Prat, V., Mathis, S., Buysschaert, B., et al. 2019, A&A, 627, A64
  * Prat et al. (2020) Prat, V., Mathis, S., Neiner, C., et al. 2020, A&A, 636, A100
  * Saio et al. (2021) Saio, H., Takata, M., Lee, U., Li, G., & Van Reeth, T. 2021, MNRAS, 502, 5856
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  * Schwarzenberg-Czerny (1991) Schwarzenberg-Czerny, A. 1991, MNRAS, 253, 198
  * Schwarzenberg-Czerny (2003) Schwarzenberg-Czerny, A. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 292, Interplay of Periodic, Cyclic and Stochastic Variability in Selected Areas of the H-R Diagram, ed. C. Sterken, 383
  * Serenelli et al. (2021) Serenelli, A., Weiss, A., Aerts, C., et al. 2021, A&A Rev., 29, 4
  * Simon & Lee (1981) Simon, N. R. & Lee, A. S. 1981, ApJ, 248, 291
  * Szewczuk & Daszyńska-Daszkiewicz (2018) Szewczuk, W. & Daszyńska-Daszkiewicz, J. 2018, MNRAS, 478, 2243
  * Szewczuk et al. (2021) Szewczuk, W., Walczak, P., & Daszyńska-Daszkiewicz, J. 2021, MNRAS, 503, 5894
  * Townsend (2003) Townsend, R. H. D. 2003, MNRAS, 340, 1020
  * Triana et al. (2015) Triana, S. A., Moravveji, E., Pápics, P. I., et al. 2015, ApJ, 810, 16
  * Van Beeck et al. (2020) Van Beeck, J., Prat, V., Van Reeth, T., et al. 2020, A&A, 638, A149
  * Van Hoolst (1994a) Van Hoolst, T. 1994a, A&A, 292, 471
  * Van Hoolst (1994b) Van Hoolst, T. 1994b, A&A, 286, 879
  * Van Hoolst (1995) Van Hoolst, T. 1995, A&A, 295, 371
  * Van Hoolst (1996) Van Hoolst, T. 1996, A&A, 308, 66
  * Van Hoolst & Smeyers (1993) Van Hoolst, T. & Smeyers, P. 1993, A&A, 279, 417
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Vuille (2000) Vuille, F. 2000, MNRAS, 313, 179
  * Vuille & Brassard (2000) Vuille, F. & Brassard, P. 2000, MNRAS, 313, 185
  * Waelkens (1991) Waelkens, C. 1991, A&A, 246, 453
  * Waelkens & Rufener (1985) Waelkens, C. & Rufener, F. 1985, A&A, 152, 6
  * Waskom & the seaborn development team (2020) Waskom, M. & the seaborn development team. 2020, mwaskom/seaborn
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61
  * Wu (2001) Wu, Y. 2001, MNRAS, 323, 248
  * Zhang et al. (2018) Zhang, C., Liu, C., Wu, Y., et al. 2018, ApJ, 854, 168
  * Zong et al. (2016a) Zong, W., Charpinet, S., & Vauclair, G. 2016a, A&A, 594, A46
  * Zong et al. (2016b) Zong, W., Charpinet, S., Vauclair, G., Giammichele, N., & Van Grootel, V. 2016b, A&A, 585, A22