© The Authors 2023

1. Introduction

Slowly pulsating B (SPB) stars are non-radial pulsators with spectral types between B9 and B3, effective temperatures ranging from 11 000 K to 22 000 K, and masses from about 3 to 9 M_⊙ (Waelkens 1991). They are main-sequence stars that exhibit high radial order gravity (g-) mode oscillations, which allows the use of asteroseismology to investigate the physical processes taking place in their interiors. In order to achieve this, space-based monitoring lasting for years is necessary to resolve the frequencies of the g modes, which have individual periods of roughly 0.5 to 3 days, with a good enough precision for asteroseismology (Aerts et al. 1999; Mathias et al. 2001; De Cat & Aerts 2002; De Cat et al. 2007). Despite these high demands, the monitoring efforts are worthwhile, as SPB stars have the potential to provide the much needed calibration of the stellar structure and evolution theory for massive stars with convective cores and radiative envelopes. Indeed, the well-known asteroseismic scaling relations used for the low-mass stars based upon their solar-like oscillations cannot be extrapolated to stars with a convective core. This leaves poorly calibrated physical properties of their convective core boundary layers, such as its thermal structure or size (often imposed via a free parameter).

Gravity-mode asteroseismology of main-sequence stars saw its birth thanks to the five-month light curves assembled by the CoRoT space telescope (Degroote et al. 2010; Neiner et al. 2012) and underwent a major boost from the four-year datasets from the Kepler mission (e.g. Aerts 2021, for a review). It has meanwhile been applied to numerous cases, both in SPB stars and less massive γ Doradus stars to probe a variety of physical processes such as near-core rotation rates (e.g. Van Reeth et al. 2016; Li et al. 2020; Takata et al. 2020), chemical mixing in the radiative envelope (e.g. Mombarg et al. 2020, 2022), magnetic fields (e.g. Buysschaert et al. 2018; Prat et al. 2019; Lecoanet et al. 2022), opacities (e.g. Szewczuk & Daszyńska-Daszkiewicz 2018; Walczak et al. 2019), and convective boundary mixing (CBM) whether looking solely at chemical element transport (e.g. Moravveji et al. 2016) or also altering the thermal structure in this CBM region (e.g. Pedersen et al. 2021; Michielsen et al. 2021). These modelling efforts typically involve a parameter study in a multidimensional space, which is generally computationally intense and time-consuming. In order to be able to perform asteroseismic modelling of a sample of SPB stars covering the entire core-hydrogen burning stage and rotational frequencies from almost zero up to the critical rate within a reasonable computation time, Pedersen et al. (2021) developed a statistical approach to get an initial estimate for the parameters and internal properties of 26 sample SPB stars.

Our work concerns a novel methodological framework, which we applied to one of these 26 sample stars, namely the double-lined spectroscopic binary SPB KIC 4930889. In order to achieve this, we relied on the most recent spectroscopic analyses and frequency determinations done since the study by Pedersen et al. (2021). The emphasis of this work is on the CBM processes including the thermal structure and the near-core rotation rate of this binary g-mode pulsator. We first provide an overview of its observed properties in the next section.

2. Gravito-inertial pulsations in the double-lined spectroscopic binary KIC 4930889

Pápics et al. (2017) analysed the Kepler light curve of the system and identified a g-mode period series consisting of 20 pulsations with mode identification (ℓ,m) = (1, 1) and consecutive radial orders, where a positive m-value denotes prograde modes. Additionally they gathered high-resolution spectra of the target using the HERMES spectrograph (Raskin et al. 2011). Their spectral synthesis of the disentangled spectra of the system resulted in the parameters listed in Table 1. The analysis places both components in the SPB instability strip, which raises the question as to which star the prograde dipole mode pattern belongs.

Asteroseismic modelling, assuming that the primary is the pulsator and based on a statistical approximation for both stellar model properties and mode frequencies, was performed by Pedersen et al. (2021). The parameters of their best stellar models are listed in Table 2. While such an approximative statistical approach was chosen to keep simultaneous treatment of 26 pulsators within reasonable computation time, it can never be as detailed as an approach tuned to a particular star, as offered here. Moreover, the results listed in Table 2 were based on the frequency list and spectroscopic parameters derived by Pápics et al. (2017). Johnston et al. (2019) revisited and renormalised the spectra obtained by Pápics et al. (2017), deriving a new spectroscopic solution listed in Table 3. They found their solution to be in better agreement with the evolutionary expectations of their isochrone-cloud modelling. We rely on these spectroscopic parameters to guide the modelling.

To compute the stellar luminosity, we used

$$\begin{array}{c}\hfill log\frac{L}{L_{\odot }}=-0.4(M_{S_{\lambda }}+{\mathrm{BC}}_{S_{\lambda }}-M_{\mathrm{bol},\odot })\end{array}$$(1)

with

$$\begin{array}{c}\hfill M_{S_{\lambda }}=m_{S_{\lambda }}-5log\frac{d}{10\mathrm{pc}}-R_{S_{\lambda }}E(B-V).\end{array}$$(2)

In this expression, m_Sλ represents the apparent magnitude from the Gaia G-band, and d is the Gaia eDR3 distance from Bailer-Jones et al. (2021). The E(B − V) reddening value was obtained with the 3D reddening map (Bayestar19) from Green et al. (2019), and the reddening vector R_Sλ equals 3.002 for the G band (Pedersen et al. 2020). The bolometric correction BC_Sλ was calculated adopting the prescription from Pedersen et al. (2020) (model 3, LTE+non-LTE), and M_bol, ⊙ = 4.74.

However, the system is a close binary, yet it was treated as a single object in Gaia DR3. To derive the luminosity for both components separately, we therefore had to adjust their apparent magnitudes based on their respective light factors (here denoted as f_i). Given that the apparent magnitude of the system depends on the total observed flux divided by the zero-point (for simplicity denoted as F)

$$\begin{array}{c}\hfill m_{S_{\lambda },\mathrm{system}}=-2.5log(F),\end{array}$$(3)

we can calculate the apparent magnitude of the individual components by adjusting the total apparent magnitude based on their respective light factors

$$\begin{array}{cc}\hfill m_{S_{\lambda },{\mathrm{component}}_i}=-2.5log(f_i·F)& =-2.5log(F)-2.5log(f_i)\hfill \end{array}$$(4)

$$\begin{array}{cc}& =m_{S_{\lambda },\mathrm{system}}-2.5log(f_i).\hfill \end{array}$$(5)

To derive the light factors for the spectroscopic solution from Johnston et al. (2019), we followed Eq. (3) from Tkachenko (2015). We calculated the ratio of continuum intensities for atmospheric models with parameters given in Table 3 weighed by the Gaia eDR3 G passband transmissivity, and used the radii ratio of ℛ = R_secondary/R_primary = 0.76 from Johnston et al. (2019). This lead us to light factors of 0.67  ±  0.01 and 0.33  ±  0.01 for the primary and secondary, respectively. Using these corrections to the apparent magnitude, we derived a luminosity $log\frac{L}{L_{\odot }}=2.61\pm 0.04$ for the primary, and $log\frac{L}{L_{\odot }}=2.23\pm 0.08$ for the secondary.

Both components of the binary fall within the SPB instability strip, and they both contribute a significant amount to the total observed flux of the system. To aid the determination of the pulsator, we divided the Kepler light curve in 20 segments according to orbital phase, using the methodology from Van Reeth et al. (2023), in search of pulsation amplitude and phase modulations as a function of the binary orbital phase. In close binaries, these can be caused by either tidal perturbations (e.g. Reyniers & Smeyers 2003a,b; Samadi Ghadim et al. 2018; Bowman et al. 2019; Steindl et al. 2021) or tidal tilting (e.g. Springer & Shaviv 2013; Handler et al. 2020; Kurtz et al. 2020; Fuller et al. 2020). In intermediate to wide binaries, these can be caused by the light travel time effect (e.g. Shibahashi & Kurtz 2012; Murphy et al. 2014). Only one dominant mode turned out to be part of a multiplet separated by the orbital frequency in the Fourier transform of the star’s light curve. This mode is f₄₃ from Table B.1 (1.1949135 ± 0.0000019 d⁻¹), and the multiplet consists of only two modes. Moreover, none of the detected pulsation phase and amplitude modulations are consistent with each other, except when the modulation is zero within the uncertainty. This is contrary to what is expected for frequency modulations caused by binarity (Shibahashi & Kurtz 2012). The tidal perturbation of pulsations can affect individual modes more strongly than others, depending on their geometry. However, for g-modes that form a period-spacing pattern (i.e. have the same geometry) we expect the orbital multiplet structure to be very similar for all modes, and most easily detectable for the dominant g-modes (Van Reeth et al. 2023). This was also not the case. From this we conclude that neither tidal perturbations or tilting nor the light travel time effect can be detected for this target.

Additionally, we attempted to determine the pulsator via line profile variations in the 26 spectra observed by Pápics et al. (2017). However, due to the relevant lines being either very weak or broad and shallow, the S/N of the spectra proved insufficient to draw any conclusions from this approach. We therefore make no assumption regarding which component hosts the pulsations, in contrast to the previous study by Pedersen et al. (2021) who assumed the primary to be the pulsator. We subsequently show that this is the least likely of the two scenarios.

The analysis of the Kepler light curve and the subsequent frequency extraction was revisited by Van Beeck et al. (2021). We used the frequency list obtained by their ‘strategy 3’ for the prewhitening procedure, as this method explains the highest fraction of variance for this light curve. This full list of frequencies is provided in Table B.1. The mode period pattern that we identified from this frequency list consists of 22 prograde dipole modes of consecutive radial order, which are listed in Table B.2 and shown in Fig. 1. This provides us with the asteroseismic input for our modelling, Y^obs composed of the individual mode frequencies $Y_i^{\text{Obs}}$ with i = 1, …, 22. Two additional modes occur in our pattern, while they were not present in the one based on the 20 modes found by Pápics et al. (2017). These are the two modes with the longest periods in Fig. 1.

Fig. 1. Prograde dipole g-mode pattern of KIC 4930889. The top panel shows the amplitude spectrum in grey and the frequencies extracted by Van Beeck et al. (2021) in blue. The frequencies selected to be part of the prograde dipole mode pattern are indicated by dashed red lines. The bottom panel shows the period-spacing pattern (ΔPn ≡ Pn + 1 − Pn) of the selected prograde dipole modes.

We additionally found a high amount of modes at lower frequencies that form two similar period series, listed in Table B.3 and shown in Fig. 2. They have an upward tilt, typical for retrograde modes, and were also detected by Pápics et al. (2017). We also found one isolated peak in the amplitude spectrum at a period of 18.297 ± 0.014 d (f₁₆₅ in Table B.1 with frequency 0.05465 ± 0.00004 d⁻¹), which is in perfect agreement with the orbital period found by Pápics et al. (2017). We searched the full list of detected frequencies in Table B.1 for frequencies that correspond to multiples of the orbital frequency within a 2σ uncertainty interval. Frequencies f₁₆₀ (0.10934 ± 0.00005 d⁻¹), f₁₅₅ (0.16388 ± 0.00006 d⁻¹), f₁₅₀ (0.21865 ± 0.00008 d⁻¹), and f₁₄₄ (0.27325 ± 0.00009 d⁻¹) coincide with two, three, four, and five times the orbital frequency, respectively. f₁₆₀ and f₁₄₄ correspond to f₁₁ and f₄ from the first additional pattern in Table B.3, and f₁₅₀ corresponds to f₆ from the second additional pattern Table B.3. Frequencies f₁₁₅ (0.81955 ± 0.00013 d⁻¹) and f₇₀ (1.09295  ±  0.00005 d⁻¹) coincide with 15 and 20 times the orbital frequency, respectively, with f₇₀ corresponding to f₁₃ from the main dipole mode period series in Table B.2. The latter of these two may not be a tidally excited oscillation, since coinciding with a relatively high multiple of the orbital frequency could be coincidental, and since it falls in line with the dipole mode period series.

Fig. 2. Same as Fig. 1, but for the possible additional series of the low frequency peaks.

Since some of the signals in the secondary patterns are low multiples of the orbital frequency, those signals might be caused by proximity effects instead of actual pulsations, and those patterns as a whole are in any case likely influenced by the binary orbit. Additionally, we could not unambiguously determine the degree of the involved modes in these patterns from our single-star approach to the asteroseismology. We therefore do not include them in our modelling at this stage of the work. These two extra patterns offer potential for future more in-depth modelling based on close binary evolution models. Inclusion of these patterns requires developing a dedicated method to include tidal interactions, which is beyond the scope of the current study. Here, we restrict ourselves to asteroseismology based on single-star evolution models, as discussed in Sect. 4.

3. Orbital modelling of proximity effects

Figure 3 shows the phase-folded light curve after prewhitening for all detected significant frequencies apart from the multiples of the orbital frequency. The residual intrinsic variability in the light curve that is left after the stopping criterion employed by Van Beeck et al. (2021) is dominant over the one caused by the orbital motion. We employ PHOEBE version 2.4.11 (Conroy et al. 2021) to construct models for the orbital harmonics. We fix the mass ratio, eccentricity, argument of periapsis, and projected semi-major axis to the values listed in Table 1. Furthermore we fix the orbital period to 18.297 d, since this is the orbital frequency retrieved from the light curve itself, it falls within the uncertainty of the value in Table 1, and the phase-folded light curve has less scatter than if 18.296 d from Table 1 is used in the phase fold. The effective temperatures are fixed to the values from Table 3 since they only have a minor impact on the simulated light curve compared to the other parameters. The surface gravity of the stars are left as free parameters, with the values from Table 3 as initial guesses. The inclination is a free parameter as well, as no prior information is available on its value.

Fig. 3. Phase-folded residual light curve. Prewhitened for all frequencies listed in Table B.1, except those corresponding to the first five orbital harmonics. The signal from those first five orbital harmonics is shown in red.

The Nelder-Mead algorithm (Nelder & Mead 1965) is employed to optimise this initial setup. Afterwards we compute a small parameter study around the retrieved solution, of which the projected goodness of fit can be seen in Fig. 4. As discussed before, the leftover intrinsic variability in the light curve is dominant, we therefore employ a mask during this parameter study to only model the bump in the light curve around orbital phase 0.4 since this is the clearest signal that is present in the harmonics. As can be seen from Fig. 4, the best fitting surface gravity of the primary star deviates from the value of 3.55 ± 0.24 from Johnston et al. (2019). Their value of 4.38 ± 0.10 for the surface gravity of the secondary agrees very well with the best fitting values we retrieve in Fig. 4. The distribution for the χ² values of the inclination, as shown in Fig. 4, is however much flatter than those of the surface gravities. Figure 5 demonstrates that inclination is indeed not very well constrained. The figure shows models with inclination angles of 60° and 74°, which both reproduce the modelled signal comparably well. The absence of a clear minimum in the χ² distribution of the models, combined with the wide range of inclinations that produce visually similarly good models, leaves us unable to confidently use this as a constraint on our asteroseismic modelling.

Fig. 4. χ2 values of the PHOEBE models. The values are projected along the log(g) of the primary (left), the log(g) of the secondary (center), and along the inclination of the binary (right).

Fig. 5. Phase-folded PHOEBE model for the five orbital harmonics from Van Beeck et al. (2021). The dots in blue are included in the mask and modelled in the parameter study, whereas the dots in grey are excluded by the mask.

4. Computation of theoretical mode frequencies

4.1. Stellar equilibrium models

Following a similar setup as Michielsen et al. (2021), we compute two grids of single star models as input for the pulsation computations. The input physics for these two grids is the same, apart from the adopted temperature gradient in the core CBM region. The first, called the radiative grid, adopts the radiative temperature gradient in that transition zone. The other is termed the Péclet grid and adopts a temperature gradient based on the Péclet number in that transition zone, following the same prescription as in Eq. (5) of Michielsen et al. (2021). This prescription includes a convective penetration zone extending the convective core, which entails that (at least a part of) the CBM region is fully adiabatic. The mixing coefficient in the CBM region is governed by two parameters, α_CBM and f_CBM which dictate the step-like and exponentially decaying parts of the region respectively. The diffusive mixing in the radiative envelope is implemented to increase going further outwards due to internal gravity waves as deduced by Rogers & McElwaine (2017), following Pedersen et al. (2018). The level of this mixing in the radiative envelope at its inner boundary with the CBM region is set by the parameter D_env. The parameter ranges for the two grids with different temperature gradients are identical and listed in Table 4. We note that the upper bound of the central hydrogen fraction is the initial fraction at the zero-age main sequence, which can vary depending on the initial metallicity of that model. Figure 6 illustrates the stellar structure of a model with the maximum amount of mixing included in our grid, and compares it with the structure of one of the models with considerably less mixing. We can clearly see that strong mode trapping occurs in the CBM region when the amount of mixing is on the lower end. An increased amount of mixing in the envelope causes the chemical gradient to be less steep, entailing a much less pronounced peak of the Brunt-Väisälä frequency, and a greatly reduced (or even absent) mode trapping.

Fig. 6. Radial profiles of a 4 M⊙star from the Péclet grid with a central hydrogen content Xc = 0.5. The top panel shows the temperature gradients and the mean molecular weight per gas particle (μ). The middle panel shows the Brunt-Väisälä frequency (N), as well as the shape of the mixing profiles, divided in convective core (grey), near-core mixing (blue), and diffusive mixing in the outer radiative envelope (green). The bottom panel shows the mode inertia of two g modes with different radial orders. All dashed lines and transparent colours correspond to a model with the maximum amount of mixing included in our grid (αCBM = 0.3,  fCBM = 0.03,  log(Denv) = 4), whereas the solid lines and darker colours correspond to a model with a considerably lower amount of mixing (αCBM  =  0.1,  fCBM  =  0.01,  log(Denv)  =  1).

The two grids of stellar evolution models are computed using the stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) version r15140. The models employ an Eddington grey atmosphere as atmospheric boundary condition and make use of the OP opacity tables (Seaton 2005). They contain the standard chemical mixture of OB stars in the solar neighbourhood deduced by Nieva & Przybilla (2012) and Przybilla et al. (2013). We determine the initial helium fraction by adopting an enrichment law Y_ini = Y_p + (ΔY/ΔZ)Z_ini. We set the primordial helium abundance Y_p = 0.2465, as determined by Aver et al. (2013). Since there is currently no consensus on the value of $\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}$ (e.g. Verma et al. 2019, and references therein), we require that the galactic enrichment ratio, $\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}$, is able to reproduce the mass fractions of the adopted chemical mixture (X = 0.71, Y = 0.276, Z = 0.014) derived by Nieva & Przybilla (2012). This leads us to adopt ΔY/ΔZ = 2.1. After Y_ini is determined according to this enrichment law, X_ini is set following X_ini = 1 − Y_ini − Z_ini.

We adopt the mixing length theory as developed by Cox & Giuli (1968) with a mixing length parameter α_mlt = 2.0, and use the Ledoux criterion for convection without allowing for semi-convection. This is warranted since this form of slow mixing is absent in the presence of CBM (e.g. Kaiser et al. 2020), which is included in the vast majority of our models. The exact location where the transition from core to near-core mixing is made, is determined by the f₀ parameter in MESA. We fix f₀ = 0.005, except for setting f₀ = 0 in the models where both α_CBM and f_CBM are equal to zero as there is no CBM region for this case. A link to the detailed MESA setup is provided in Appendix A.

4.2. Pulsation computations

The pulsation mode properties of the MESA equilibrium models are computed employing the stellar oscillation code GYRE (Townsend & Teitler 2013; Townsend et al. 2018), version 6.0.1. Since non-adiabatic effects mainly become important in the outer stellar envelope, the adiabatic approximations are sufficient for our modelling work due to the mode inertias of the g modes being dominant near the stellar core. For computational reasons and given that it does not affect the mode frequencies at the level of measurement errors, we only perform non-adiabatic computations for some of our best models after the forward modelling is finished, in order to evaluate their mode excitation. We compute the dipole g modes for all our equilibrium models for an initial guess of the rotation frequency, assuming rigid rotation and relying on the traditional approximation of rotation (TAR; e.g. Eckart 1960; Bildsten et al. 1996; Lee & Saio 1997), following its implementation in GYRE (as described in Townsend 2020, Sect. 4).

The stellar rotation frequency required to optimally reproduce the observed stellar pulsations differs for the varying equilibrium models. We therefore start from the same initial guess for each equilibrium model and rescale the g-mode frequencies for each model separately, following the TAR and assuming rigid rotation, to reproduce the observed pulsations as closely as possible. This optimisation is performed using the Levenberg-Marquardt method implemented in LMFIT (Newville et al. 2020). To reduce the chances of the optimisation method returning a local minimum, we start the optimisation procedure from two separate initial values for the rotation. The first one being the initial guess, ω_{initial, 1} = ω_guess, used to calculate the GYRE model. The second initial value is taken by adjusting the first one by twice the difference between the initial value and its solution, so ω_{initial, 2} = ω_guess − 2 ⋅ (ω_guess − ω_{optimised, 1}) = 2 ⋅ ω_{optimised, 1} − ω_guess. This way the global minimum of the initial value problem is approached both from a higher and a lower initial value. In the case where these solutions do not converge, we take the best of the two returned solutions since it indicates that the other one returns a local minimum.

Figure 7 illustrates the rescaling of the period-spacing pattern due to a change in the rotation rate. It also shows the relative differences between the periods of the rescaled modes and the periods obtained by repeating the GYRE computation using that same optimised rotation rate. We find that the rescaled mode periods agree well with the periods computed by GYRE for the new rotation rate. The differences are of order 10⁻³% in the asymptotic mode frequency regime where the observed pulsations occur, and even the largest differences at low radial orders are still relatively small (< 0.05%). Rescaling the g-mode frequencies to the optimised rotation frequency and selecting a set of the theoretical frequencies to match the observations yields for each equilibrium model a list of theoretically predicted dipole mode frequencies, Y^Theo composed of ${\mathbf{Y}}_i^{\text{Theo}}$, where i stands for the radial order. The GYRE inlist to compute the initial frequency lists is provided through the link in Appendix A.

Fig. 7. Rescaling period spacing patterns to an optimised rotation rate. Top panel shows the period-spacing pattern as calculated by GYRE with the initial guess for the rotation in blue, and the pattern rescaled to the optimised rotation frequency in orange. The inset figures are zoomed in on the region with the observed pulsations, with the modes selected to match the observations circled in red. The bottom panel shows the relative difference between the mode periods calculated by GYRE given the optimised rotation rate, and the mode periods from the rescaled pattern. The grey region denotes the observational uncertainties.

5. Modelling approach

We utilise the same asteroseismic modelling procedure as Michielsen et al. (2021). A brief overview is provided here for convenience without going too much into the details.

5.1. General mathematical framework

We employ the Mahalanobis distance as a merit function for the maximum likelihood estimation in the asteroseismic modelling (see Aerts et al. 2018, for its application to asteroseismic modelling),

$$\begin{array}{c}\hfill {\text{MD}}_j={({\mathbf{Y}}_j^{\text{theo}}-{\mathbf{Y}}^{\text{obs}})}^T{(V+\mathrm{\Sigma })}^{-1}({\mathbf{Y}}_j^{\text{theo}}-{\mathbf{Y}}^{\text{obs}}),\end{array}$$(6)

with Y^obs the vector of observations and ${\mathbf{Y}}_j^{\text{theo}}$ the corresponding vector of predicted values in gridpoint j. Σ is the variance matrix due to the measurement errors of Y^obs and V is the variance-covariance matrix of Y^theo capturing the theoretical uncertainties in the mode frequency predictions caused by the limited knowledge of the physical ingredients in the input physics of the equilibrium models, taking as well the correlations among the free parameters used to describe these ingredients into account.

The modelling involves both statistical models that are non-nested, comparing models within one grid of equilibrium models, and statistical models that are nested, comparing equilibrium models across different grids where none, one, or both of the CBM parameters are fixed at zero. This allows for a comparison between different numbers of free parameters, including an evaluation whether the increase in goodness of fit outweighs the entailed punishment by the selection criterion for having an increased number of free parameters.

We use the Akaike Information Criterion corrected for small sample size (AICc; Claeskens & Hjort 2008, Chap. 2) since it rewards fit quality but penalises complexity. It is defined as

$$\begin{array}{c}\hfill \text{AICc}=-2ln\mathcal{L}+\frac{2kN}{N-k-1},\end{array}$$(7)

with N and k the number of observables and free parameters, respectively, and ℒ the likelihood of a stellar model. In our framework for this star, N = 22 when fitting periods, or 21 when fitting period spacings (that is the differences in period between two pulsation modes of consecutive radial order, ΔP_n ≡ P_n + 1 − P_n). The number of free parameters k is 4, 5, or 6 depending on whether two, one, or zero of the CBM parameters (α_CBM, f_CBM) are fixed in the nested grids. In case k = 6, the list of parameters consists of (M_ini, Z_ini, α_CBM, f_CBM, log(D_env),X_c). Rewriting the AICc for the likelihood function of the Mahalanobis Distance yields

$$\begin{array}{c}\hfill \text{AICc}=ln(|V+\mathrm{\Sigma }|)+kln(2\pi )+\text{MD}+\frac{2kN}{N-k-1}.\end{array}$$(8)

The performance of two nested models can be compared through their difference in AICc values ΔAICc = AICc_A − AICc_B. Model B is preferred over model A if ΔAICc > 2, with a (very) strong preference if ΔAICc > 6 (10).

We determine the uncertainty region of the best solution by employing Bayes’ theorem, stating that the probability of a parameter θ^m occurring in the interval |${\theta }_a^m$, ${\theta }_b^m$| is given by

$$\begin{array}{cc}\hfill P({\theta }_a^m<{\theta }^m<{\theta }_b^m|\mathbf{D})=& \frac{{\sum }_i^{q}P(\mathbf{D}|{\mathit{\theta }}_i)P({\mathit{\theta }}_i)}{{\sum }_j^{Q}P(\mathbf{D}|{\mathit{\theta }}_j)P({\mathit{\theta }}_j)}\hfill \\ \hfill =& \frac{{\sum }_i^{q}P(\mathbf{D}|{\mathit{\theta }}_i){\prod }_l^{k}P({\theta }_i^l)}{{\sum }_j^{Q}P(\mathbf{D}|{\mathit{\theta }}_j){\prod }_l^{k}P({\theta }_j^l)}.\hfill \end{array}$$(9)

Index j is summed over all Q equilibrium models in the grid that are consistent within 3σ of the spectroscopic log g, T_eff, and stellar luminosity, that are also consistent with the observed metallicity and constraints from the binarity of the system. These binary constraints are explained in more detail in Sect. 5.2 Index i is summed over the q models with the highest likelihood so that $P({\theta }_a^m<{\theta }^m<{\theta }_b^m|D)=0.95$.

We consider three approaches to match theoretical mode periods to the observed ones, and analyse the results of the method that performs best for each grid. In the first two we begin matching mode periods starting from the theoretical period that is closest to the either the mode with the highest observed amplitude or the highest-frequency detected in the observed pattern. The third option is to match each observed mode period to its best matching theoretical counterpart, and adopt the longest sequence of consecutive modes that we get in this way. The rest of the pattern is then build consecutively in radial order starting from this sequence. These three options of pattern construction will henceforth be referred to as highest amplitude, highest frequency, and longest sequence.

Apart from just the mode periods, we also consider the period-spacing values as a set of observables to be used in our modelling procedure. The condition numbers of the variance-covariance matrices V + Σ are used to determine the best of these sets of observables. The condition number κ is defined as the ratio of its maximum to minimum eigenvalue,

$$\begin{array}{c}\hfill \kappa (V+\mathrm{\Sigma })=\frac{|{\lambda }_{\max }(V+\mathrm{\Sigma })|}{|{\lambda }_{\min }(V+\mathrm{\Sigma })|}.\end{array}$$(10)

This gives an indication of how well- or ill-conditioned the matrix is with respect to the inversion to be computed, with lower values being better conditioned.

5.2. Isochrone clouds

The methodology from Michielsen et al. (2021) as summarised in Sect. 5.1 considers the system’s asteroseismic and spectroscopic data from a single-star perspective. KIC 4930889 is however a double-lined spectroscopic eccentric (e = 0.32) binary. Hence, we can utilise the information obtained from the binarity of the system to put additional constraints on the models in an attempt to lift some of the degeneracies that are present. We employ the use of isochrone clouds (Johnston et al. 2019), which is in this application the collection of isochrones of a given age but for all combinations of α_CBM, f_CBM, log D_env present in our grid. Constructing an isochrone cloud coupled to a model of a certain grid, we enforce all models in that cloud to be of an age that differs less than one gridstep in age from this model, have the same initial metallicity, and have a mass that is compatible within the error margin of the mass ratio of the system (listed in Table 1). We computed some additional evolutionary tracks for masses above and below the grid range listed in Table 4 to allow all masses that we could expect for the companion from the observed mass ratios to be present in the isochrone clouds.

Figure 8 illustrates how the constraints of the isochrone clouds are applied for the case where we assume the secondary star to be the pulsator. It shows the isochrone clouds of the primary star for three different models of the secondary, which are arbitrarily chosen for the purpose of this visual representation. These three models are among the best, that is having the lowest AICc values, and they would have been included in the error ellipses if the system were modelled from a single-star perspective without any constraints from binarity. The light-grey tracks in the background show all models with the same metallicity and a mass within the observed mass ratio, and only the models that have an age difference smaller than one gridstep are shown in colour. Two of the isochrone clouds fall partially within the 1σ or 3σ errors of the companion star and are thus accepted as viable solutions. The third isochrone cloud falls completely outside of the 3σ spectroscopic error region, and is hence not accepted as a solution. Although Fig. 8 only showcases the constraints on log g and T_eff, the stellar luminosity is also used in these isochrone-cloud constraints.

Fig. 8. Isochrone clouds of the primary star matching three different models of the secondary. The black lines show the 1σ and 3σ spectroscopic Teff and log g error boxes of the primary star.

6. Modelling results

The condition numbers of the variance-covariance matrices computed via Eq. (10) are of order κ(A)∼10³ to 10⁴ when considering mode periods, but are significantly smaller, down to κ(A)∼10¹, when considering period spacings. We therefore primarily consider the period spacings as the set of observables to fit, but still list the results from using the periods as observables as well. Although the individual parameters of the best models may differ between these two sets of observables, they are in most cases quite similar if not the same, and always fall within the other error ellipse. Our conclusions are therefore independent from the chosen observable.

We model the observed pulsation pattern twice. Once with the spectroscopic constraints of the primary, for which the corner plots of the radiative and Péclet grid are shown in Figs. 9 and 10, and once with the spectroscopic constraints of the secondary, with Figs. 11 and 12 showing the corner plots for the radiative and Péclet grid respectively. The models included in the 2σ error ellipse of the MD according to Eq. (9) are shown in colour, while the models in grey scale fall outside of this error ellipse. Additionally we make a comparison between the AICc values of the best models of the full grid and each partial grid with fewer free parameters, for both prescriptions of the temperature gradient in the CBM region. We hereafter refer to the grids with six free parameters as the radiative and Péclet grid, but specify when talking about grids with fewer free CBM parameters. Among the nested grids with five free parameters, we have α_CBM = 0 but varying f_CBM, henceforth denoted as exponential radiative or exponential Péclet grid, and f_CBM = 0 but varying α_CBM, henceforth denoted as step radiative or step Péclet grid. The nested grid with four free parameters, having both α_CBM and f_CBM set to zero, is referred to as the grid without CBM. Comparing all these grids with various numbers of free parameters not only enables us to investigate which temperature gradient is preferred in the CBM region, but also to examine if the increased fit quality outweighs the penalties for higher model complexity.

Fig. 9. Corner plot for the radiative grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the primary star. The 50% best models are shown, colour-coded according to the log of their merit function value (at right). The models in colour fall within the 2σ error ellipse of the MD constructed using Eq. (9), whilst the models in grey fall outside of this error ellipse. The figures on the diagonal show binned parameter distributions of the models in the error ellipse, and the panel at the top right shows a Hertzsprung–Russell diagram with the 1 and 3σTeff and log L error boxes.

Fig. 10. Corner plot (as in Fig. 9) for the Péclet grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the primary star.

Fig. 11. Corner plot (as in Fig. 9) for the radiative grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the secondary star.

Fig. 12. Corner plot (as in Fig. 9) for the Péclet grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the secondary star.

The model parameters and AICc values of these best models are listed in Tables 5 and 6 when enforcing the constraints on the luminosity and spectroscopic T_eff and log g of the primary and secondary, respectively. We cannot distinguish the preferred temperature gradient whilst modelling the primary star, since there is no preference between the exponential radiative, exponential Péclet, or the grid without any CBM, given that ΔAICc < 2 between their best models. We do however find a preference of these three grids over the grids with a step-like mixing profile, or the ones with a combined step and exponential mixing. The period spacings of the best models from these indistinguishable grids are shown in Fig. 13.

Fig. 13. Period-spacing patterns of the observations, and of the best models of the preferred grids that are not distinguishable from one another. These are the models in bold in Tables 5 and 6 that use ΔP as observables. The formal errors on the observations are smaller than the symbol sizes. The largest of the observational errors is enlarged ten times and shown for comparison. The vertical bars in the bottom left corner of the top panel show the maximum considered uncertainty for the theoretical predictions approximated by the variance-covariance matrix of that particular grid. The middle and bottom panels show the relative difference in period spacing and period, respectively, between the observation and the model. The narrow grey areas indicate the formal 1σ observational uncertainty from Table B.2.

Modelling the secondary star also yields no possibility to distinguish which temperature gradient is preferred. It is in this case not possible to differentiate between both grids with an exponentially decaying mixing in the CBM region, the grid without any CBM, and the radiative grid with six free parameters. The period spacings of the best models from these indistinguishable grids are shown in Fig. 13.

We clearly see from both Fig. 13 that the variance of the theoretical predictions is much larger than the uncertainties on the observations. Our solutions are therefore dominated by the theoretical uncertainties, rather than the observational ones. The larger variance of the grid without CBM is one of the reasons why there is no selection capacity between this grid and those with exponential CBM. A reduction of the variance of the theoretical predictions would lead to a much stronger preference for the presence of CBM over the absence of CBM according to the AICc. This would in particular be the case should we ignore the (co)variance due to limits in the theoretical predictions (V = 0), that is when reducing the merit function from a Mahalanobis distance to a χ².

6.1. Near-core rotation rate and convective core mass

From the 2σ MD error ellipses on the best models, we find the near-core rotation rate of the star to be well constrained. The values are listed alongside the best models of each nested grid in Tables 5 and 6. If we take the grids that are indistinguishable according to the AICc, we find all of them to be consistent with the result of our best model grid; ${\mathrm{\Omega }}_{\text{rot}}=0.{73}_{-0.06}^{+0.02}{\mathrm{d}}^{-1}$. This is about 13.4 times the orbital frequency of this eccentric binary, and about 37% of the best model’s Roche critical rotation rate. Even when considering the nested grids that are not preferred, we can see that their near-core rotation rates are also consistent and agree very well with one another.

Defining the convective core mass as the one determined by the Ledoux criterion without including the CBM region (as is visualised by the grey area in Fig. 6), we constrain it to $M_{\mathrm{cc}}=0.{67}_{-0.20}^{+0.17}$ M_⊙. This value is consistent across all grids we considered in the modelling, keeping in mind its uncertainties.

6.2. Mode excitation

We are left with multiple different solutions that cannot be distinguished from each other based on modelling the period-spacing values and using the spectroscopic, astrometric, and isochrone-cloud restraints. Therefore we look at which of these models performs best at reproducing the mode excitation of our observed period-spacing pattern.

Figure 14 shows the normalised growth rates, η, of the modes (Stellingwerf 1978). These indicate an excited or damped mode for a positive or negative value of η, respectively. For the primary, the models from the exponential Péclet and the exponential radiative grid have ten modes excited out of the 22 observed modes in our pulsation pattern. The model from the grid without CBM has sixteen of the observed modes as excited. This higher amount of excited modes is an effect of its more evolved nature compared to the other best models, rather than a direct effect of the absence of CBM (e.g. Fig. 1. of Pápics et al. 2017, which shows an increasing number excited modes during the first part of the main-sequence evolution). All of these models show some excited modes at shorter periods that were not observed in our pattern. The model for the secondary star from the exponential Péclet grid shows twelve out of the 22 modes excited. The models from the exponential radiative grid and the grid without CBM show seven excited modes, while the model from the radiative grid shows no excited modes at all.

Fig. 14. Instability parameter η. The parameter is shown as a function of mode period for the best models of the primary (top panel) and secondary (bottom panel). The period spacing patterns of these models are shown in Fig. 13, coloured circles indicate excited modes, while empty circles indicate the non-excited ones. The vertical lines show the observed mode periods and their amplitude.

Accurately reproducing the excitation of high radial order g modes in B-type pulsators often requires opacity enhancements which were not considered in this study (e.g. Moravveji 2016; Daszyńska-Daszkiewicz et al. 2017; Walczak et al. 2019; Szewczuk et al. 2022). In our results, we also see that the increased number of excited modes in general corresponds to an increased metallicity, and hence elevated opacity of the iron- and nickel-group chemical elements. We therefore confirm that the standard OP opacity tables are insufficient to accurately reproduce the observed mode excitations, rather than using this result to constrain our solutions.

7. Conclusions

In this work we investigated the gravito-inertial modes in the double-lined B-type binary KIC 4930889. We explored which of the components hosts the pulsations, what the preferred temperature gradient and mixing profiles are in the CBM region of that star, and constrained its near-core rotation rate. We employed asteroseismic, spectroscopic and astrometric information, and constraints obtained from the binarity through isochrone clouds.

The quality of our best asteroseismic solutions are better for the spectroscopic and astrometric constraints of the secondary star, although not statistically significantly better than when using those of the primary star. The difference in luminosity between both stars, with 67% and 33% of the light contribution from the primary and secondary respectively, is not large enough to assign the pulsation signal to one or the other. Furthermore, we do not find a preference for the temperature gradient based on the 22 mode periods or 21 period spacings. However, we are able to constrain the near-core rotation rate of the pulsating component to ${\mathrm{\Omega }}_{\text{rot}}=0.{73}_{-0.06}^{+0.02}{\mathrm{d}}^{-1}$. We also obtain better solutions for models with some type of exponentially decaying CBM over those with a step-like mixing profile in the CBM region.

We find that the model from the exponential Péclet grid performs better at explaining the mode excitation than the models with a radiative temperature gradient. The Péclet model shows twelve out of the 22 observed modes excited, whereas the radiative models have at most seven of the observed modes excited. The larger number of excited modes is due to the higher metallicity of the stellar equilibrium model selected from this grid; that is, the different temperature gradient only indirectly influences the predicted mode excitation. The model that best explains the excited modes is however the one without CBM present that assumes the spectroscopic and astrometric constraints for the primary star. This absence of CBM also only indirectly influences the mode excitation, since this model is more evolved than the others along the first half of its main sequence, entailing a higher number of theoretically predicted excited modes.

Comparing our detailed follow-up treatment to the earlier performed statistical modelling approach of Pedersen et al. (2021) and Pedersen (2022b), we note a few key differences in the modelling setup. Whilst Pedersen et al. (2021) investigated two different prescriptions for the CBM region and four for the envelope mixing, our study considers only one case of envelope mixing but seven different ones for the CBM region. We opted for such an approach since the g modes that we are considering have a much higher probing power in the CBM region than in the stellar envelope, as can be seen from their mode inertia in Fig. 6. Pedersen (2022b) was able to distinguish between different shapes of the envelope mixing, but found no preference between their considered CBM prescriptions. This would equal a comparison between our exponential radiative and step Péclet grids, where we do find a preference for the exponential grids over the step grids.

As far as the retrieved model parameters are concerned, we list both the results from Pedersen (2022b) and from this work (using the spectroscopic and astrometric constraints of the primary and secondary) in Table 7. All parameters are in agreement when considering our error estimation and considering the spectroscopic and astrometric constraints for the primary star. In particular, the stellar rotation rate aligns well. However, our best asteroseismic model was found when we consider the constraints for the secondary star. For this case the mass and central hydrogen content are no longer compatible within the projected error ellipses. Our best point estimators do deliver a younger, less massive star with more CBM and less envelope mixing. We note that the uncertainties on the two sets of results differ substantially, where our uncertainties on the parameters encompass the results from Pedersen (2022b), but not vice versa. These different results are influenced slightly by the different spectroscopic constraints and set of prograde dipole modes that are employed, but stem dominantly from the modelling approach, where we used a more detailed treatment of the asteroseismic modelling as compared to the approximative statistical modelling approach of Pedersen et al. (2021, 2022b,a), who used statistical approximations for the pulsations rather than detailed GYRE computations. We find their uncertainties obtained from the approximative statistical modelling to be underestimated for this star.

Although the error ellipses of our solutions contain less than 3% of our initial models, the projections of the six dimensional error ellipse in one dimension results in uncertainties on each individual parameter that range over most of the initial model grid. The vast majority of the parameter combinations that are included in these one-dimensional projections are however not part of the actual higher dimensional error ellipse. An indication of this can be seen in the two-dimensional projections of the error ellipse in Figs. 9–12. In contrast, the results obtained by Michielsen et al. (2021) for KIC 7760680 yielded much smaller error ellipses so that the uncertainties remained small even when they were projected on one dimension. This difference is due to the number of modes in the observed prograde dipole mode pattern, which amounted to 36 modes for KIC 7760680 and to 22 modes for KIC 4930889, where a larger number of observed modes entails a better probing power of the stellar interior. With this in mind, constraints from spectroscopic data, Gaia astrometric data, and from the binarity of the system are valuable to complement asteroseismic information. These complementary constraints become all the more beneficial for stars with lower asteroseismic probing power due to fewer observed modes.

Similar to Michielsen et al. (2021), so regardless of the amount of pulsations in our observed period pattern, we find that the uncertainties on the theoretically predicted pulsation patterns are much larger than the uncertainties on the observed patterns. The uncertainties in our modelling are therefore dominated by the theoretical model uncertainties, rather than the observational ones. Additionally, the theoretical variance is largest in the grid without CBM, causing the lack of selection capacity between these models and the ones with CBM. A reduction of the variances would lead to a stronger preference of model grids with CBM over the ones without it. Hence, future work should prioritise improving stellar evolutionary models by both refining and expanding the physical processes that are included in them. KIC 4930889 is a good target to evaluate tidal effects in close binary evolution models, given our detection of multiples of the orbital frequency in its secondary period spacing patterns.

Acknowledgments

The authors thank Jordan van Beeck for having provided the frequencies he found for KIC 4930889 from his 2021 paper in electronic format, and Sarah Gebruers and Alex Kemp for the insightful scientific discussions. The authors are grateful to the MESA and GYRE developers teams for their efforts and for releasing their software publicly; this study would not have been possible without their codes. The research leading to these results has received funding from the Research Foundation Flanders (FWO) by means of a PhD scholarship to MM under project No. 11F7120N, a postdoctoral fellowship to TVR with grant agreement No. 12ZB620N, and to AT through grant agreement No. G089422N, from the KU Leuven Research Council (grant C16/18/005: PARADISE to PI Aerts), as well as from the BELgian federal Science Policy Office (BELSPO) through PRODEX grant PLATO.

References

Appendix A: MESA and GYRE Inlists

The example MESA and GYRE inlists used for this work are available from the MESA inlists section of the MESA marketplace: https://cococubed.com/mesa_market/inlists.html.

Appendix B: Frequency list

We provide lists of frequencies.

All Tables

All Figures

	Fig. 1. Prograde dipole g-mode pattern of KIC 4930889. The top panel shows the amplitude spectrum in grey and the frequencies extracted by Van Beeck et al. (2021) in blue. The frequencies selected to be part of the prograde dipole mode pattern are indicated by dashed red lines. The bottom panel shows the period-spacing pattern (ΔPn ≡ Pn + 1 − Pn) of the selected prograde dipole modes.
In the text

	Fig. 2. Same as Fig. 1, but for the possible additional series of the low frequency peaks.
In the text

	Fig. 3. Phase-folded residual light curve. Prewhitened for all frequencies listed in Table B.1, except those corresponding to the first five orbital harmonics. The signal from those first five orbital harmonics is shown in red.
In the text

	Fig. 4. χ2 values of the PHOEBE models. The values are projected along the log(g) of the primary (left), the log(g) of the secondary (center), and along the inclination of the binary (right).
In the text

	Fig. 5. Phase-folded PHOEBE model for the five orbital harmonics from Van Beeck et al. (2021). The dots in blue are included in the mask and modelled in the parameter study, whereas the dots in grey are excluded by the mask.
In the text

Fig. 6. Radial profiles of a 4 M⊙star from the Péclet grid with a central hydrogen content Xc = 0.5. The top panel shows the temperature gradients and the mean molecular weight per gas particle (μ). The middle panel shows the Brunt-Väisälä frequency (N), as well as the shape of the mixing profiles, divided in convective core (grey), near-core mixing (blue), and diffusive mixing in the outer radiative envelope (green). The bottom panel shows the mode inertia of two g modes with different radial orders. All dashed lines and transparent colours correspond to a model with the maximum amount of mixing included in our grid (αCBM = 0.3,  fCBM = 0.03,  log(Denv) = 4), whereas the solid lines and darker colours correspond to a model with a considerably lower amount of mixing (αCBM  =  0.1,  fCBM  =  0.01,  log(Denv)  =  1).

In the text

Fig. 7. Rescaling period spacing patterns to an optimised rotation rate. Top panel shows the period-spacing pattern as calculated by GYRE with the initial guess for the rotation in blue, and the pattern rescaled to the optimised rotation frequency in orange. The inset figures are zoomed in on the region with the observed pulsations, with the modes selected to match the observations circled in red. The bottom panel shows the relative difference between the mode periods calculated by GYRE given the optimised rotation rate, and the mode periods from the rescaled pattern. The grey region denotes the observational uncertainties.

In the text

	Fig. 8. Isochrone clouds of the primary star matching three different models of the secondary. The black lines show the 1σ and 3σ spectroscopic Teff and log g error boxes of the primary star.
In the text

Fig. 9. Corner plot for the radiative grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the primary star. The 50% best models are shown, colour-coded according to the log of their merit function value (at right). The models in colour fall within the 2σ error ellipse of the MD constructed using Eq. (9), whilst the models in grey fall outside of this error ellipse. The figures on the diagonal show binned parameter distributions of the models in the error ellipse, and the panel at the top right shows a Hertzsprung–Russell diagram with the 1 and 3σTeff and log L error boxes.

In the text

	Fig. 10. Corner plot (as in Fig. 9) for the Péclet grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the primary star.
In the text

	Fig. 11. Corner plot (as in Fig. 9) for the radiative grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the secondary star.
In the text

	Fig. 12. Corner plot (as in Fig. 9) for the Péclet grid. Made using period spacings in a Mahalanobis distance merit function and spectroscopic and luminosity constraints from the secondary star.
In the text

Fig. 13. Period-spacing patterns of the observations, and of the best models of the preferred grids that are not distinguishable from one another. These are the models in bold in Tables 5 and 6 that use ΔP as observables. The formal errors on the observations are smaller than the symbol sizes. The largest of the observational errors is enlarged ten times and shown for comparison. The vertical bars in the bottom left corner of the top panel show the maximum considered uncertainty for the theoretical predictions approximated by the variance-covariance matrix of that particular grid. The middle and bottom panels show the relative difference in period spacing and period, respectively, between the observation and the model. The narrow grey areas indicate the formal 1σ observational uncertainty from Table B.2.

In the text

	Fig. 14. Instability parameter η. The parameter is shown as a function of mode period for the best models of the primary (top panel) and secondary (bottom panel). The period spacing patterns of these models are shown in Fig. 13, coloured circles indicate excited modes, while empty circles indicate the non-excited ones. The vertical lines show the observed mode periods and their amplitude.
In the text