© ESO 2021

1. Introduction

Stellar evolutionary models are the realisation of our understanding of stellar structure and evolution, and they are the foundation for more complex models and simulations of stellar populations, clusters, and galaxies (Salaris & Cassisi 2005). As such, any inaccuracies in the mathematical description or calibration of the various input physics that are included in these models inevitably propagate to the results of studies that rely on them as input. In addition, the one-dimensional (1D) nature of stellar evolutionary models necessitates the inclusion of parametric simplifications of multi-dimensional physical processes, such as the 1D approximation of turbulent convection through the incorporation of the mixing length theory (Böhm-Vitense 1958).

Stellar evolutionary models are therefore limited by the framework of these simplifications, complicating the task of improving these models. However, it has become increasingly clear in recent decades that one of the dominant sources of uncertainty regulating the evolutionary pathways of stars is the amount and prescriptions of the various mixing processes occurring in stellar interiors due to their ability to modify the amount of fuel available for nuclear fusion (Maeder 2009; Brott et al. 2011; Langer 2012; Meynet et al. 2013). The influence of these mechanisms is particularly significant for the evolution of stars that are born with or that develop a convective core on the main sequence (M ≳ 1.15 M_⊙): Uncertainties in the calibration and implementation of these mixing processes propagate to later evolutionary stages, most ubiquitously to the red-giant branch for intermediate-mass stars with 1.15 M_⊙ ≲ M ≲ 2.5 M_⊙ (Constantino et al. 2015; Silva Aguirre et al. 2020) and to the pre-supernova stage for the massive stars (Martins et al. 2015; Bowman 2020).

Internal chemical mixing is a result of independent or coupled physical mechanisms operating in different regions of the stellar interior (see Salaris & Cassisi 2017 for a review) and can be broadly divided into two classes for the M ≳ 1.15 M_⊙ regime: core-boundary mixing (CBM) and envelope mixing. CBM is an collective term that includes all prescriptions for the transport of material from the convective core into the radiative envelope (or vice versa), such as convective entrainment (Meakin & Arnett 2007; Cristini et al. 2019), convective penetration and overshooting (Viallet et al. 2015). Of these various prescriptions, it is convective penetration and overshooting that are most commonly implemented in most 1D stellar evolutionary codes (Viallet et al. 2015), and both are referred to in the literature as ‘overshooting’. Both of these prescriptions are based on the traversing of convective fluid parcels beyond the limits set by the Schwarzchild or Ledoux criterion due to their inertia. In convective penetration (also known as ‘step overshooting’; Zahn 1991), the region in which these fluid particles traverse is subjected to the adiabatic temperature gradient, resulting in an extended region around the convective core being instantaneously mixed and modifying the thermal structure of the star. In convective overshooting (also known as ‘diffusive exponential overshooting’; Freytag et al. 1996; Herwig 2000), the extended region is subjected to the radiative temperature gradient, and hence only the chemical structure of star is changed. In both types of CBM regions, different functional forms of temperature gradient have been implemented, but distinguishing between them observationally is non-trivial (Michielsen et al. 2019, 2021).

The net effect of overshooting is an increase in the convective core mass, and hence for the star to appear more luminous (i.e. brighter). This effect is therefore degenerate with the stellar mass, and classical observables would not be able to distinguish between both types of overshooting. However, the analyses of stars that show gravity-mode (g-mode) pulsations (e.g., Pedersen et al. 2018; Michielsen et al. 2019) have demonstrated that the type and amount of CBM can be constrained by modelling them asteroseismically. This has been demonstrated by a number of observational studies of g-mode pulsators (e.g., Moravveji et al. 2015, 2016; Buysschaert et al. 2018; Walczak et al. 2019; Wu & Li 2019; Fedurco et al. 2020).

Envelope mixing, similar to CBM, is an umbrella term that includes a number of rotational and pulsational mechanisms (see Aerts 2021 for a detailed discussion), including meridional circulation, hydrodynamical instabilities, magnetorotational instabilities and internal gravity waves. Envelope mixing in 1D stellar models is typically implemented by means of diffusive mixing coefficients in the transport equations for the regions outside of the convective core and CBM zones. There are notable exceptions, such as the advective implementation in the Geneva (Ekström et al. 2012) and CESTAM (Marques et al. 2013) stellar evolutionary codes. Due to this implementation, the envelope mixing is artifically decoupled from the CBM. It was noted by Moravveji et al. (2015, 2016) that the inclusion of envelope mixing significantly improved the fit of the g-mode pulsational frequencies and therefore they confirm that such mixing is an important component that should be included in 1D stellar models.

Due to the sensitivity of g modes to stellar interior mixing profiles, the analysis of g-mode pulsating stars allows for the calibration of these interior mixing profiles. For pulsating F-type stars, these modes are low-frequency pulsations that are excited by the convective flux-blocking mechanism (Guzik et al. 2000; Dupret et al. 2005). At high radial order, the g modes that have the same spherical harmonic degree ℓ and azimuthal order m (see Aerts et al. 2010 for a complete description) have the characteristic of being equally spaced in period in the non-rotating approximation. This gives rise to the so-called asymptotic approximation of the period spacing (Tassoul 1980):

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{\ell }=\frac{{\mathrm{\Pi }}_0}{\sqrt{\ell (\ell +1)}};\end{array}$$(1)

where,

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_0=2{\pi }^2{\left({{\displaystyle \int }}_{r_1}^{r_2}N\frac{\mathrm{d}r}{r}\right)}^{-1}·\end{array}$$(2)

In these equations, r is the distance from the stellar centre, N is the Brunt–Väisälä frequency and r₁ and r₂ are the radial boundaries of the g-mode propagation cavity in the star. Deviations from this equidistant spacing of g modes are a result of either (1) mode trapping due to the near-core chemical gradient (Miglio et al. 2008), or (2) near-core rotation, which introduces a slope into the pattern due to the pulsational periods being affected by the Coriolis force (Bouabid et al. 2013). The morphology of the observed g-mode period-spacing patterns can therefore be linked to mixing and rotational characteristics of the g mode propagation cavity (e.g., Van Reeth et al. 2015a,b, 2016; Ouazzani et al. 2017, 2020; Li et al. 2018, 2019a, 2020a).

The potential for g modes in unravelling the interior mixing profiles in stars is unfortunately hampered by the highly correlated nature of the various free parameters used in the evolutionary models. Aerts et al. (2018) proposed a forward modelling scheme for the purposes of taking these correlations, as well as uncertainties due to the imperfections in the input physics of the equilibrium models, into account when performing asteroseismic modelling. In addition, restricting the ranges of the multidimensional parameter space would significantly alleviate these degeneracies in stellar modelling, resulting in more precise and accurate solutions.

The largest number and most stringent of constraints to limit the parameter space can be obtained when studying double-lined eclipsing binary systems. One can extract the surface metallicities, effective temperatures, masses, radii, and surface gravities and rotational frequencies of the individual components. In addition to all of these parameters, binarity also enforces a co-evolutionary scenario, demanding equal ages of the individual components. These systems have been used to calibrate the amount of CBM included in stellar evolutionary models, testing either the convective penetration prescription (e.g., Ribas et al. 2000; Torres et al. 2014; Tkachenko et al. 2014a; Stancliffe et al. 2015; Claret & Torres 2016), the convective overshooting prescription (e.g., Tkachenko et al. 2014b; Stancliffe et al. 2015; Claret & Torres 2017, 2018, 2019; Constantino & Baraffe 2018; Johnston et al. 2019a,b) or both (e.g., Claret & Torres 2017; Tkachenko et al. 2020). A secondary advantage of studying such systems is the ability to compare three different sets of parameters: the dynamical and spectroscopic parameters, those derived from the evolutionary models, and those derived from the pulsational models.

There have been a total of 34 eclipsing binary systems with g-mode period-spacing patterns discovered so far in the Kepler data (Li et al. 2020b), and only three have been modelled (Zhang et al. 2018, 2020; Guo & Li 2019), none of which from the viewpoint of assessing internal mixing. As such, it is the goal of our study to add to the dearth of literature on the evolutionary and asteroseismic modelling of g-mode pulsators in eclipsing binary systems. In this paper, we present the evolutionary and asterosesismic modelling of the pulsating eclipsing binary KIC 9850387 that was identified as the most promising candidate in terms of g-mode asteroseismic potential in Sekaran et al. (2020), hereafter Paper I. It is a short period (P_orb = 2.74 d) eclipsing binary in a near-circular orbit with an intermediate-mass primary star and a solar-like secondary star. The primary star displayed remarkable ℓ = 1 and ℓ = 2 period-spacing patterns as detailed in Paper I. Section 2 details the isochrone-cloud construction and fitting methodology adopted in this work, Sect. 3 details our asteroseismic modelling results, and we present a discussion of our results and conclusions in Sect. 5.

2. Isochrone-cloud construction and fitting

In Paper I, a full observational analysis of KIC 9850387 was performed, resulting in the extraction of the orbital and fundamental parameters of the system as well as of the individual components (listed in Table 1). The ≲1% errors on many of these parameters are not unusual for detached eclipsing systems: The individual analyses of several tens of such systems (listed in DEBCat, Southworth 2015) boast similar precisions. Using these parameters, we can calculate evolutionary models that best match our observations. To achieve this goal, we constructed isochrone-clouds (hereafter isoclouds) as detailed in Johnston et al. (2019b). Each isocloud is comprised of isochrones generated from equivalent evolutionary phase (EEP) tracks (see Dotter 2016 for a full description) generated from single-star main-sequence evolutionary tracks computed with the stellar evolutionary code MESA (version 10348; Paxton et al. 2011, 2018). These evolutionary tracks are generated using the same input physics as Johnston et al. (2019b), fixing the initial helium abundance Y_ini at 0.274 as per the cosmic B-star standard (Nieva & Przybilla 2012), fixing α_MLT at 1.8 (Joyce & Chaboyer 2018), and using an initial metallicity (Z_ini) of 0.010 to match the spectroscopic metallicity derived for the primary star.

We adopted the Ledoux criterion to position the convective boundary and diffusive exponential overshooting as the CBM prescription according to Freytag et al. (1996) and Herwig (2000). Restricting to one type of CBM is justified because asteroseismic analyses of a sample of A/F-type γ Dor pulsators have shown the results using this CBM to be indistinguishable from convective penetration (Mombarg et al. 2019). This prescription assumes an exponential decrease of the mixing efficiency with distance from the convective core as convective fluid parcels propagate further away from the core into the radiative zone. The amount of diffusive mixing in the CBM region [D_CBM(r)] is calculated using the following equation:

$$\begin{array}{c}\hfill D_{\mathrm{CBM}}(r)=D_0\exp \left[\frac{-2(r-r_0)}{f_{\mathrm{ov}}{H}_{\mathrm{p}}}\right]·\end{array}$$(3)

Here, the reference radius r₀ = r_cc − f₀H_p, where r_cc is the radius of the convective core according to the Ledoux criterion, f₀ = 0.002 and H_p is the local pressure scale height. D₀ is the amount of diffusive mixing within the inner edge of the convective-core boundary at the reference radius r₀ and f_ov is the length scale over which the diffusive exponential overshooting applies.

The prescriptions of Johnston et al. (2019b) include ranges for f_ov and the amount of diffusive envelope mixing D_mix, taken from the asteroseismic calibrations of Aerts (2015), Bowman (2020), and Moravveji et al. (2015, 2016) respectively. More specifically, the diffusive envelope mixing adopted in Johnston et al. (2019b) takes the form of D_mix ∝ ρ^−1/2, which is an approximation of the mixing due to internal gravity waves taken from 2D hydrodynamical simulations (Rogers & McElwaine 2017) and implemented in MESA by Pedersen et al. (2018). However, these parameters have been calibrated using g-mode pulsators, and are therefore only appropriate for stars that develop a substantial convective core over the main sequence (M_ini ≳ 1.15 M_⊙ as noted by Dotter 2016). The low dynamical mass of the secondary ($M_{\mathrm{s}}=1.{0617}_{-0.005}^{+0.003}{M}_{\odot }$) necessitates the calculation of evolutionary tracks of stars below this limit. Choi et al. (2016) computed a grid of solar-calibrated stellar models and isochrones but only include rotational mixing in their tracks with M_ini ≥ 1.15 M_⊙ to reflect the slow rotational frequencies observed in low-mass stars. Considering that rotational mixing is the only form of envelope mixing present in their tracks, it can be concluded that they effectively set the diffusive envelope mixing at zero for their evolutionary tracks with M_ini < 1.15 M_⊙. In addition, a significant portion of the envelope of low-mass stars is convective, implying that mixing in the radiative region outside of the core would have a minimal impact on their evolution. We therefore fix both f_ov and D_mix at 0 for our evolutionary tracks for M_ini < 1.15 M_⊙. The ranges of M_ini, f_ov, and log D_mix used to create our evolutionary tracks are listed in Table 2. It should be noted that our use of an upper bound of 4.0 for log D_mix is atypically high for intermediate-mass g-mode pulsators (e.g., Mombarg et al. 2019 had used an upper bound of 1.0 for log D_mix). Such high values allow for the testing of the existence of tidal mixing mechanisms in the stellar interior.

An isochrone is traditionally constructed from evolutionary model tracks with the same input physics and the same free-parameter (e.g., f_ov and D_mix) values. However, due to our choice of fixing f_ov and log D_mix at 0 for M_ini < 1.15 M_⊙, we effectively have to fuse the isochrones created from the M_ini ≤ 1.10 M_⊙ and M_ini ≥ 1.15 M_⊙ evolutionary tracks together: A single isochrone therefore has an f_ov and log D_mix value of 0 for 0.80 M_⊙ ≤ M_ini ≤ 1.10 M_⊙ and a fixed combination of f_ov and log D_mix as per the ranges listed in Table 2 for 1.15 M_⊙ ≤ M_ini ≤ 2.00 M_⊙, with a total of 1000 datapoints with 0.80 ≤ M_ini ≤ 2.00. Similarly, each isocloud, which is a combination of all isochrones at an age τ for all combinations of f_ov and log D_mix in our grid, exhibits the same behaviour. Each of our isoclouds comprise a total of 72 000 datapoints: 8 f_ov values × 9 log D_mix values × 1000 datapoints per isochrone.

In order to obtain the best-fitting isocloud to the dynamical parameters of KIC 9850387, we perturbed the primary and secondary T_eff and log g within their 1σ observational errors in a Monte-Carlo framework for 1000 iterations, and retained the best-fitting isocloud in each iteration. Due to the fact that an isocloud is effectively a set of curves, we calculated the reduced χ² (${\chi }_{\mathrm{red}}^2$) value for each of the 72 000 datapoints in an isocloud with respect to the primary and secondary T_eff and log g (i.e. 144 000 total ${\chi }_{\mathrm{red}}^2$ values), and took the sum of the 50 smallest ${\chi }_{\mathrm{red}}^2$ values for both the primary and secondary star (i.e. the sum of 100 ${\chi }_{\mathrm{red}}^2$ values across both components) as our goodness-of-fit metric. In addition, to account for the fact that the stellar models are in much better agreement with dynamical parameters for low-mass stars (e.g., Chaplin et al. 2014; Higl & Weiss 2017), we chose to bias the isocloud fitting towards the secondary component. This was done by decreasing the errors on the secondary parameters by a factor of five, effectively weighting the secondary parameters five times more heavily than the primary ones. We then sample this distribution of the 1000 best-fitting isoclouds and determined the evolutionary parameters (τ, M, f_ov, log D_mix, the central hydrogen fraction X_c, and the inferred convective-core mass M_cc) based on the intersection of the dynamical T_eff¹ and log g constraints with the isocloud parameters corresponding to the 95% highest posterior density (HPD) interval² of the Monte-Carlo age distribution, corresponding to ${\tau }_{\mathrm{MC}}=1.3_{-0.2}^{+1.5}$ Gyr. These parameters are listed in Table 3, and the log T_eff − log g (or Kiel) diagram displaying our isocloud fitting results is shown in Fig. 1. The large asymmetry in the errors of τ_MC is a consequence of older isoclouds covering a larger fraction of the T_eff − log g parameter space at higher masses, as the varying amounts of core-boundary and envelope mixing included in the models cause the corresponding evolutionary tracks to diverge more strongly with increasing age. These isoclouds therefore start to overlap more greatly at older ages, allowing for a greater range in the upper age bound. It should be noted that the isocloud diverges into three distinct sub-clouds due to the logarithmic scaling of the log D_mix values in our grid.

Fig. 1. Isocloud fit to the components of KIC 9850387 on the log Teff − log g diagram. The red regions correspond to the best-fitting isocloud, and the pink regions represent the isoclouds in the 95% HPD interval of the fit. The black regions on the isocloud represent the log Teff − log g values corresponding to the dynamical masses of the individual components. The grey curves are MESA evolutionary tracks with fov = 0.005 and log Dmix = 0.0 (solid), fov = 0.04 and log Dmix = 0.0 (dashed), and fov = 0.04 and log Dmix = 4.0 (dotted) with their corresponding masses (in units of M⊙) indicated at the ZAMS of each track. The inset plots are magnified regions around the position of the primary (top) and secondary (bottom) component.

As noted in Paper I, there is a slight discrepancy between the evolutionary and dynamical masses of the primary star. However, our isoclouds were calculated from evolutionary tracks using a fixed Z_ini = 0.010, which is close to the best-fitting spectroscopic surface metallicity but its associated errors are not taken into account. Varying the initial metallicity within the spectroscopic errors (i.e. 0.008 ≲ Z_ini ≲ 0.012) resolves this discrepancy, but introducing Z_ini as another free parameter leads to additional degeneracies (e.g., between Z_ini and M, and between Z_ini and f_ov). These degeneracies would have further propagated to our subsequent asteroseismic analyses, and as such we chose not to include Z_ini as a free parameter for our study.

It can be seen in Table 3 that both f_ov and D_mix are unconstrained in our fit, although log D_mix > 3 is incompatible with the position of the star in the Kiel diagram (cf. Fig. 1). This has been noted in other studies as well (e.g., Johnston et al. 2019a,b; Mombarg et al. 2019; Tkachenko et al. 2020). However, as noted by those very same studies, M_cc is well constrained, and our value of M_cc/M_* = 0.11 ± 0.02 is compatible with the range (0.075 ≲ M_cc/M_* ≲ 0.125) reported for a sample of γ Dor stars by Mombarg et al. (2019).

3. Asteroseismic modelling

Our asteroseismic modelling is based on two different sets of asteroseismic observables: the asymptotic period spacing Π₀ (described in Eqs. (1) and (2) and the individual period-spacing values (ΔP_ℓ) of the ℓ = 1 and ℓ = 2 modes. As noted by Zhang et al. (2020) and Li et al. (2020a), the primary star of KIC 9850387 is an extremely slow rotator, with Li et al. (2020a) reporting a near-core rotational frequency of f_rot, core = 0.0053 d⁻¹. We are therefore able to estimate the observational values of Π₀ from ${\mathrm{\Pi }}_{0,\mathrm{obs}}\simeq \overline{\mathrm{\Delta }{P}_{\ell }}\ast \sqrt{\ell (\ell +1)}$ in the approximation of a non-rotating star. The determination of Π₀ is more complicated for fast rotators (e.g., Van Reeth et al. 2016, 2018; Li et al. 2018, 2019a,b, 2020a,b) as this requires an estimation of the slope of the period spacing pattern, which is absent for KIC 9850387. We are therefore able to perform asteroseismic modelling using the quantity Π₀ directly, without estimation of the rotation frequency, through Eqs. (1) and (2), by comparing our Π_0, obs estimates with values predicted from the non-rotating stellar evolutionary models. One of the advantages of isocloud fitting (see Sect. 2) is that it provides age constraints for our stellar evolutionary models. We therefore restricted our evolutionary models to the ages that were within the 95% HPD interval of the Monte Carlo isocloud-fitting age distribution (${\tau }_{\mathrm{MC}}=1.3_{-0.2}^{+1.5}$ Gyr).

In order to obtain theoretical ΔP_ℓ, values, one would have to calculate pulsational models. These models require stellar structural models (i.e. a model describing the temperature, density, chemical and mixing profiles from the centre to the surface of the star) as an input. Using our isocloud-fitting age constraint, we computed a grid of stellar structural models with masses between 1.2 M_⊙ and 2.0 M_⊙ in steps of 0.02 in X_c as an input to compute pulsational mode predictions in the adiabatic framework using the stellar pulsation code GYRE (revision 5.2; Townsend & Teitler 2013; Townsend et al. 2018). We then use the theoretical pulsational frequencies of the zonal ℓ = 1 and ℓ = 2 modes to construct theoretical period-spacing patterns and confront these to the observed ones, which is a more detailed analysis compared to just using Π₀.

For both approaches, our fitting methodology involves the perturbation of the Π_0, obs and observational ΔP_ℓ values of the period-spacing patterns (ΔP_ℓ, obs) within their respective 1σ observational errors in a Monte-Carlo framework for 10 000 iterations, and then retaining the parameters of the best-fitting model in each iteration. To determine the best-fitting model within this framework, we used two different merit functions: (1) the ${\chi }_{\mathrm{red}}^2$ and (2) the Mahalanobis Distance (MD). As discussed in detail in Aerts et al. (2018), the MD is a maximum-likelihood point estimator that takes into account uncertainties in the theoretical asteroseismic predictions and that treats correlations in the free parameters of our stellar model grid appropriately by incorporating its covariance matrix (V) into the distance calculation:

$$\begin{array}{c}\hfill \mathrm{MD}={({\mathit{Y}}_{\mathrm{theo}}-{\mathit{Y}}_{\mathrm{obs}})}^{\mathrm{T}}{(\mathbf{V}+\mathbf{\Lambda })}^{-1}({\mathit{Y}}_{\mathrm{theo}}-{\mathit{Y}}_{\mathrm{obs}}).\end{array}$$(4)

In this equation, Y_theo and Y_obs are the vectors representing the theoretical asteroseismic predictions and asteroseismic observables that are being compared, and the matrix Λ is a diagonal matrix with squared observational errors as each of the diagonal elements. The free parameters to be estimated, determining the theoretical asteroseismic predictions, are M, X_c, f_ov, and log D_mix. The MD takes uncertainties in the theoretical asteroseismic predictions due to the imperfect input physics of the stellar structural models in the grid into account. It does so via the matrix V, where we assume that its components are well described by the variance induced by the whole range of the four free parameters defining the theoretical asteroseismic predictions.

3.1. Π₀ modelling results

We performed our Π₀ modelling by restricting the size of our grid based on varying types of observational constraints. Due to the tight constraints imposed by the dynamical parameters, a grid setup created by restricting the parameter space based on the dynamical M, T_eff, and log g would necessarily be very small, comprising only a few tens of models in total. In order to investigate the effects of different levels of parameter space restriction on our results, we chose to construct a third grid setup of intermediate detail and size between the full evolutionary grid and the small asteroseismic grid based on the dynamical constraints. To that end, we obtained a second set of T_eff and log g constraints by performing atmospheric parameter determination for the primary star of KIC 9850387 (as detailed in Sect. 3.3 of Paper I) but ignoring any binary or eclipse information (i.e. treating the system as if it was a single star and performing a classical atmospheric analysis).

Using these ‘pseudo-single-star’ spectroscopic and dynamical parameters, we created three grid setups: (1) The full grid, as the name suggests, is the grid resulting from the 95% HPD interval of the Monte-Carlo age distribution of our isocloud fit (as described in Sect. 2) without any other constraints; (2) The spectroscopic grid is a subset of the full grid based on the 3σ intervals of the ‘pseudo-single-star’ spectroscopic T_eff and log g values; and (3) the dynamical grid is a subset of the full grid based on the 3σ intervals of the dynamical M, T_eff, and log g values. The models in the spectroscopic and dynamical grids have 6965 K ≤ T_eff ≤ 7520 K and 3.09 ≤ log g ≤ 4.13 (spectroscopic grid), and 1.63 M_⊙ ≤ M_ini ≤ 1.70 M_⊙, 7081 K ≤ T_eff ≤ 7590 K, and 3.9826 ≤ log g ≤ 4.0021 (dynamical grid) respectively. It should be noted that we use 3σ intervals to construct our grids to ensure that the grids are sufficiently large and to introduce sufficient variance into the grids, which as mentioned significantly impacts our MD calculations.

For each grid setup, we fitted the Π_0,  obs estimates for each of the dipole (ℓ = 1) and quadrupole (ℓ = 2) modes separately, as well as combined into a single vector with two components, resulting in three different solutions. More specifically, for each of our grid setups (1, 2, and 3), we used three different sets of observables and predictions in independent fits: (a) Y = (Π_{0,  ℓ = 1}), (b) Y = (Π_{0,  ℓ = 2}), and (c) Y = (Π_{0,  ℓ = 1}, Π_{0,  ℓ = 2}). The results for each parameter setup using either the ${\chi }_{\mathrm{red}}^2$ or the MD merit functions are listed in Table 4 (full grid), Table 5 (spectroscopic grid) and Table 6 (dynamical grid). Comparisons of the ${\chi }_{\mathrm{red}}^2$ and the MD results in terms of the best-fitting Π₀ values and positions on Kiel diagrams are shown in Figs. 3 and 4 respectively. The errors are based on 68% HPD intervals of the Monte-Carlo parameter distributions (e.g., see Fig. 2). This allows for proper comparison with the evolutionary modelling and dynamical results, which are based on 68% confidence intervals.

Fig. 2. Example of the distribution of M values in the posterior distribution of models selected through our Monte-Carlo fitting framework. The vertical dashed line represents the median value and the vertical solid lines represent the upper and lower bounds of 68% HPD of the M values.

Fig. 3. Best-fitting Π0 values of our Π0-based asteroseismic modelling, based on three different grid setups (full, spectroscopic, and dynamical). Left panel: results of using a ${\chi }_{\mathrm{red}}^2$ merit function and right panel: results of using an MD merit function. The vertical dashed black line represents the zero-point of the difference between the observational and theoretical Π0 (i.e. Π0, obs − Π0, theo = 0). The error bars are based on 68% HPD intervals of the Monte-Carlo parameter distributions.

Fig. 4. Positions of the best-fitting models from Π0-based asteroseismic modelling on Kiel diagrams. The main plots show the results of using the whole grid in the fit, while the inset plots show the results of applying the ‘pseudo-single-star’ spectroscopic Teff and log g (solid box) and dynamical M, Teff, and log g (dotted box) constraints. The observed positions of the star according to the spectroscopic and dynamical parameters are represented by black ‘X’ symbols and stars respectively. The 1σ and 3σ spectroscopic and dynamical error bars are represented with straight end-caps and ball end-caps respectively. The Yobs = (Π0,  ℓ = 1), Yobs = (Π0,  ℓ = 1), and Yobs = (Π0,  ℓ = 1, Π0,  ℓ = 2) solutions for the different grid subsets are represented by the same colours and symbols as in Fig. 3. Left panel: results of using a ${\chi }_{\mathrm{red}}^2$ merit function and right panel: results of using an MD merit function. The error bars on the asteroseismic solutions are based on 68% HPD intervals of the Monte Carlo parameter distributions. The solid and dashed grey curves are the same as in Fig. 1.

It should be noted that in some cases, the models that were selected in every single Monte-Carlo iteration had the same values for the free parameters, due to the limited resolution of the grid and its subsets. In these cases, we adopt upper- and lower-bound errors on the estimated parameters (M, X_c, f_ov, and log D_mix) that are of a single step size in the grid (i.e. 0.5 M_⊙ in M, 0.005 in f_ov and 0.5 in log D_mix) and propagate these errors to the inferred parameters (R, T_eff, log g, M_cc, and τ). These errors (hereafter single-grid-step errors) should be interpreted as a conservative upper limit on the precision of the extracted parameters based on the resolution of our grid. Parameters for which these single-grid-step errors are assumed are indicated by * next to the values in the tables.

It can be seen that there is generally little difference in the results regardless of whether a ${\chi }_{\mathrm{red}}^2$ or the MD merit function is used, regardless of the grid setup or Y configurations. This was an expected result, as the χ² and MD merit functions converge for vectors of unit or near-unit length (due to the term (V + Λ)⁻¹ used in the computation of the MD, see Eq. (4)). The small minima of the respective merit functions shown in Tables 4–6 are due to the near-perfect agreement (sub-second differences in the best cases) of the theoretical predictions with the perturbed observables in the Monte-Carlo framework, given that we fitted only one observable by the means of four free parameters. In addition, there is little difference in the results for the three different Y configurations. This is also unsurprising, as it was already noted in Paper I that the calculated Π₀ values for each of these modes agree within 2σ. The fits using Y_obs = (Π_{0,  ℓ = 1}, Π_{0,  ℓ = 2}) are in agreement with those derived from evolutionary modelling (cf. Table 3). Overall, our results of purely asteroseismic modelling with Y = (Π₀) provide weaker constraints on the external properties (M, R, T_eff, and log g), superior constraints on the interior properties (M_cc, f_ov, and D_mix) and similar constraints on the evolutionary stage (X_c and τ) when compared to evolutionary modelling (cf. Tables 3 and 4). This showcases the superiority of g-mode asteroseismic modelling over evolutionary modelling when it comes to constraining interior properties.

The fits using the spectroscopic and dynamical grids (Tables 5 and 6) resulted in superior constraints on most parameters than using either purely asteroseismic or evolutionary modelling when using both the ${\chi }_{\mathrm{red}}^2$ and MD merit functions. However, the best-fitting Π₀ values and errors are almost identical regardless of the grid subset used in the fit (Fig. 3), demonstrating the value of asteroseismic observational constraints on parameter estimation due to its degenerate nature. It can also be seen that the dynamical grid results cluster at the edge of the grid (see inset plot demarcated by a dotted box in Fig. 4), indicating that the minima of the respective merit functions are outside of the range of the dynamical grid. When compared to the full grid results, the spectroscopic and dynamical grid results also trended towards lower f_ov and D_mix values.

As mentioned previously, the period-spacing patterns of KIC 9850387 are relatively flat (i.e. do not display significant dips), implying either a high amount of mixing at the bottom of the radiative envelope or a young evolutionary stage (i.e. high X_c values). In general, the envelope mixing levels in intermediate-mass g-mode pulsating stars are hard to infer, as investigated and argued by Mombarg et al. (2019). Their sample study utilised principle component analysis to deduce an expectation value of D_mix = 1 cm² s⁻¹ for the 37 pulsators in their study, which is much lower than the wide range of values covered by B-type pulsators with much bigger convective cores (Aerts 2021, Table 1). Here, our results imply an evolved star with flat patterns revealing an above-average level of mixing for F-type g-mode pulsators (D_mix ∼ 3 cm² s⁻¹). In fact, given the flat period-spacing patterns, this value should be seen as a lower limit since higher envelope mixing will not change the already flat spacing patterns and does not change the global parameters of the stars (Mombarg et al. 2019).

The minima of the merit functions listed in Tables 4–6 reveal that using just one observable to estimate the four free parameters (M, X_c, f_ov, and log D_mix) is not very constraining as there are too many degrees of freedom. It is therefore of interest to fit the individual ΔP_ℓ values of the observed dipole and quadrupole modes.

3.2. ΔP_ℓ modelling results

We performed our ΔP_ℓ modelling in an identical way to our Π₀ modelling, based on the three grid subsets (full, spectroscopic and dynamical), using three parameter setups (Y = (ΔP_{ℓ = 1,  i}), Y = (ΔP_{ℓ = 2,  j}), and Y = (ΔP_{ℓ = 1,  i}, ΔP_{ℓ = 2,  j}), where i and j are the indices of the ΔP_ℓ values of the corresponding period-spacing patterns), and computed the same two merit functions (${\chi }_{\mathrm{red}}^2$ and MD). These results are listed in Table 7 (full grid), Table 8 (spectroscopic grid) and Table 9 (dynamical grid). Comparisons of the ${\chi }_{\mathrm{red}}^2$ and the MD results in terms of the best-fitting ΔP_ℓ values, period-spacing patterns and positions on Kiel diagrams are shown in Figs. 5–7 respectively.

Fig. 5. Same as Fig. 3, but representing the results of ΔPℓ-based asteroseismic modelling.

Fig. 6. Best-fitting ℓ = 1 and ℓ = 2 period-spacing patterns using our ΔPℓ-based asteroseismic modelling. Top panel: results of using a ${\chi }_{\mathrm{red}}^2$ merit function and bottom panel: results of using an MD merit function.

Fig. 7. Same as Fig. 4, but representing the results of ΔPℓ-based asteroseismic modelling.

During our best solution and precision estimation, we encountered a number of solutions whose precision could not be assessed, because the same grid point was selected in every Monte-Carlo iteration (as described in Sect. 3). These so-called ‘single-grid-point’ solutions occur more frequently for the ${\chi }_{\mathrm{red}}^2$ (7/9 solutions over the three grids) than for the MD framework (5/9 solutions over the three grids). This is due to the MD incorporating the uncertainties in the theoretical predictions, offering a broader range of acceptable solutions to match the perturbations of the Y_obs vector. In these cases, we once again assigned single-grid-point errors to the estimated parameters and propagated them to the inferred parameters. Similar to our purely asteroseismic Π₀-based results, ΔP_ℓ-based modelling resulted in weaker constraints on the external properties (M, R, T_eff, and log g) and superior constraints on the interior properties (M_cc, f_ov, and D_mix) when compared to evolutionary modelling.

The effect of the choice of merit function is demonstrated by Fig. 5. The MD solutions offer a larger variety of appropriate solutions according the various MD minima, because this metric allows for theoretical uncertainty in the modelling. Since the ${\chi }_{\mathrm{red}}^2$ is a specific solution of the regression based on the assumption that there are no theoretical imperfections, while ignoring correlations among the observables and parameters, its best solutions always occur in the MD solution space as well, but often at much higher MD values. This is graphically illustrated in Fig. 5, where the broad coverage of the best allowed period spacing values selected by the MD is due to the allowance of uncertainty in the theoretical models via the variance covered by the entire grid. For example, the MD value of the model from the dynamical grid (Table 9) with the lowest ${\chi }_{\mathrm{red}}^2$ value (${\chi }_{\mathrm{red}}^2=597$) is 2674, which is larger than the lowest MD value of 2378. However, the overall best ${\chi }_{\mathrm{red}}^2$ and MD solutions from the dynamical grid, demanding compliance with both the ΔP_ℓ values of the dipole and quadrupole modes, binary and the spectroscopy, are very similar. The MD methodology is sensitive to the model-independent constraints imposed upon the problem set, as it allows for the theory of the used grid to be imperfect at the level captured in the variance-covariance matrix induced by the free parameter ranges. Nevertheless, both metrics end up with an almost equal solution when demanding compliance with all observational constraints at the 3σ level. The overall best solution combining all of the asteroseismic and dynamical constraints is indicated by the MD solution in boldface in Table 9.

Similar to our Π₀-based results, our dynamical grid results imply above-average levels of envelope mixing, and in the Y = (ΔP_{ℓ = 1,  i}) and Y = (ΔP_{ℓ = 1,  i}, ΔP_{ℓ = 2,  j}) cases we find high values of D_mix ≃ 25 cm² s⁻¹ compared to those derived for single γ Dor pulsators by Mombarg et al. (2019). While this may indicate extra tidal mixing, casting doubt on the use of eclipsing binaries as testbeds of stellar structural and asteroseismic theory of single stars, several slowly-rotating single stars also display near-flat period-spacing patterns (Li et al. 2019a) as was found for KIC 9850387. Hence we posit that the estimated level of envelope mixing could still be due to intrinsic non-tidal element transport mechanisms.

4. Investigating high-frequency modes

As noted in Paper I, a handful of independent high-frequency modes were observed in the frequency spectrum of KIC 9850387. It was also stated that no frequency splitting or characteristic spacing was observed. Having obtained theoretical models that explain our observed ΔP_ℓ values, we test that claim by computing theoretical frequencies of the ℓ = 0 to ℓ = 4 modes in the 12.8 d⁻¹ to 13.4 d⁻¹ frequency range for the stellar structural models within the overall best MD solution indicated in bold in Table 9. It was found that only two or three frequencies were obtained per structural model within the considered range, and the modes associated with these frequencies have varying ℓ values, with modes of certain ℓ values not occurring in the observed range. For example, the pulsational model with the closest match between the theoretical and observational frequencies, with parameters M = 1.70 M_⊙, X_c = 0.22, f_ov = 0.01, and log D_mix = 1.5, consists an ℓ = 2 mixed mode³ (i.e. n_p = 1, n_g = 3) at a frequency of 13.22 d⁻¹ and a low-order ℓ = 4 mode (i.e. n_g = 5) at a frequency of 12.93 d⁻¹ (see Fig. 8, where we show all modes occurring in the 12.0 d⁻¹ to 14.0 d⁻¹ frequency range).

Fig. 8. Pulsational model with the closest frequency match between the observed and theoretical modes in the high-frequency regime of KIC 9850387. The vertical dotted green lines represent the orbital harmonics.

It can be seen that the theoretical modes are offset from the observed modes. The same phenomenon was also observed for the F-type binary pulsator KIC 10080943 (Schmid & Aerts 2016). Following these authors, we also posit that this offset is due to a so-called surface effect (see Ball 2017 for a recent review) that is commonplace in solar-like stars due to uncertainties in the physical descriptions of the near-surface convection and non-adiabatic effects in the outer envelope. Indeed, while intermediate-mass stars such as the primary component of KIC 9850387 have thin convective envelopes, it has been shown that time-dependent convection theory is required to model the oscillations properly (see Dupret et al. 2005). The MESA and GYRE codes do not include a treatment of time-dependent convection. Hence, frequency shifts between the theoretical predictions and the observed modes are expected given that the best models were fitted to the g modes. Their mode energy is determined by the physics in the deep adiabatic interior of the star, where the approximation of time-independent convection is appropriate.

The pulsational model with the closest match between the theoretical and observational frequencies can only potentially explain two out of the four observational frequencies. In order to explain more of the observational behaviour, we hypothesise that the pair of observed high-amplitude modes at 13.21 d⁻¹ and 13.25 d⁻¹ displayed in Fig. 8 is part of a rotationally split multiplet. However, the closest theoretical frequency is the ℓ = 2 mixed mode at 13.22 d⁻¹, and rotational splitting of such a mode would result in a quintuplet. Therefore, it is more likely that the doublet is a rotationally-split prograde-retrograde (m = 1 and m = −1) doublet with a missing zonal (m = 0) mode. This phenomenon is a result of cancellation effects (Aerts et al. 2010) due to the near 90° inclination of the system, and has been observed previously in hybrid p- and g-mode pulsators in the mass range of KIC 9850387, most notably in Kurtz et al. (2014) where multiple well-resolved complete multiplets were reported.

As such, we further restricted our pulsational models to those that included ℓ = 1 modes within the 12.8 d⁻¹ to 13.4 d⁻¹ frequency range. We then determined the model with the smallest difference in frequency between the ℓ = 1 mode and the midpoint of the observed frequency doublet (13.23 d⁻¹). This best model has parameters M = 1.65 M_⊙, X_c = 0.24, f_ov = 0.005, and log D_mix = 1.5, and three modes were obtained within this range: The fundamental radial mode (i.e. n_p = 1, ℓ = 0) at a frequency of 12.86 d⁻¹, a low-order octupole g mode (i.e. n_g = 4, ℓ = 3) at a frequency of 12.96 d⁻¹ and the g₁ mode (i.e. n_g = 1, ℓ = 1) at a frequency of 13.18 d⁻¹. These theoretical mode predictions, along with the observed modes, are displayed in Fig. 9. It should be noted that modes with ℓ > 2 are generally not observed in space-based photometry due to geometric cancellation effects (Aerts et al. 2010), although there have been claims to the contrary (see e.g., Baran & Østensen 2013). As such, we hypothesise that the observed mode at 12.92 d⁻¹ is the fundamental radial mode. The fundamental radial mode is offset from the nearest observed mode by ∼0.05 d⁻¹ while the dipole g₁ mode is offset from the nearest observed mode by ∼0.03 d⁻¹. The values of these frequency offsets are smaller than those observed for KIC 10080943 (Schmid & Aerts 2016), who were unable to identify the degree of the p modes for KIC 10080943.

Fig. 9. Pulsational model with the closest frequency match between the theoretical ℓ = 1 mode and the midpoint of the high-amplitude doublet at 13.21 d−1 and 13.25 d−1. The vertical dotted green lines represent the orbital harmonics.

We then use this model providing the optimal fit to both the g modes in the deep stellar interior and the observed p modes to calculate the rotational frequency implied by the splitting of the ℓ = 1 mode. Under the assumption of slow rotation, a non-radial mode is split due to the Coriolis force into 2ℓ+1 multiplet components with frequencies f_nlm given by the following equation (see e.g., Aerts et al. 2010):

$$\begin{array}{c}\hfill f_{\mathrm{nlm}}=f_{\mathrm{nl}}+m(1-C_{\mathrm{nl}})f_{\mathrm{rot},\mathrm{puls}}.\end{array}$$(5)

In this equation, f_nl is the unperturbed frequency of the zonal mode, f_rot, puls is the rotational frequency about the pulsational axis and C_nl is the Ledoux constant (Ledoux 1951). We computed it for the best model and obtained a value of C_{n, ℓ = 1} = 0.02. Figure 10 is a display of the rotational kernels of the highest radial order g mode of the best-fitting ℓ = 1 period-spacing pattern, and those of the fundamental radial and dipole g₁ mode. It can be inferred from the figure that the high-order g modes have the highest probing power in the near-core region, and the radial and g₁ modes have the highest probing power in the stellar envelope.

Fig. 10. Rotational kernels of the high- and low-order modes of the best-fitting pulsational model as a function of the fractional radius R/R⋆.

Using Eq. (5) and f_nl = 0.5 ⋅ (13.21 + 13.25) = 13.23 d⁻¹ (the midpoint of the doublet), we obtained a value of f_rot,puls = 0.01689 ± 0.00004 d⁻¹. This value is 3.2 times that of f_rot,core = 0.0053 ± 0.0015 d⁻¹ reported by Li et al. (2020b), and a factor of 7.2 times lower than the surface rotational frequency $f_{\mathrm{rot},\mathrm{surf}}=0.{122}_{-0.008}^{+0.008}$ d⁻¹ reported in Paper I. In that paper, the surface rotation frequency was derived from the projected rotational velocity (v_rot sin i_orb) determined through spectral line modelling, relying on the orbital inclination (i_orb) and radius of the primary star (R_p) determined through eclipse modelling, assuming aligned pulsational and orbital axes. This would imply a sharp decrease in the rotational frequency from the surface to the envelope, and a more gradual decrease from the envelope to the near-core region (as displayed in Fig. 11). However, the surface rotational frequency deduced from v_rot sin i_orb is necessarily an overestimation of the true surface rotational rate as v_rot sin i_orb was used as a proxy for the total velocity broadening required to fit the spectral lines. The contribution from asymmetric line-profile variations due to the oscillation was ignored (as was mentioned in Paper I). Assuming that the envelope rotates rigidly (i.e. f_rot, surf = f_rot, puls), we obtain v_rot sin i_orb = 1.823 ± 0.005 km s⁻¹. Keeping all other atmospheric parameters identical, it was found that the synthetic spectrum generated using a v_rot sin i_orb = 2 km s⁻¹ and including a pulsational velocity broadening in the form of macroturbulence following Aerts et al. (2009, 2014), requires v_macro = 15 km s⁻¹. Such a profile provides a qualitatively similar fit to the observed profile than the best-fitting synthetic spectrum with v_rot sin i_orb = 13 km s⁻¹ and v_macro = 0 km s⁻¹ found in Paper I, as shown in Fig. 12. Therefore, we conclude that the asteroseismic and spectroscopic data are consistent with a rigidly rotating envelope rotating ∼3 times faster than the deep interior near the convective core. The near-core region and envelope of the primary star rotate sub-synchronously with respect to the binary orbit by factors ∼69 and ∼22, respectively.

Fig. 11. Core (taken from Li et al. 2020b, based on the slope of the g-mode period-spacing pattern), envelope (from the rotational splitting of the g1 mode) and surface (from spectroscopic line broadening ignoring pulsational velocity broadening) rotational frequencies (frot) as a function of the fractional radius (R/R⋆) of the primary component of KIC 9850387. The core and envelope rotational frequencies are positioned at the maxima of the rotational kernels of the dipole ng = 39 and ng = 1 modes represented in Fig. 10. The errors on the core and envelope rotational frequencies are smaller than the symbol size. The dashed black line represents the rotational profile from the core to the surface that is compatible with macroturbulent line broadening due pulsations.

Fig. 12. Synthetic spectral fits to a Fe I line of the disentangled primary component spectrum with different compositions of rotational broadening. The solid and dotted black lines represent synthetic spectra with identical atmospheric parameter inputs but with v sin i = 13 km s−1 and vmacro = 0 km s−1 and v sin i = 2 km s−1 and vmacro = 15 km s−1 respectively.

In the above, we have assumed aligned pulsational and orbital axes. Allowing for misalignment (i.e. f_rot, puls = f_rot, surfsin(90 ° −i_offset), where i_offset is the angle between the pulsational and orbital axes), we obtain i_offset = 82 ± 5°. This would imply that the axes are near orthogonal to each other. While this could be an example of a so-called tidally-trapped or single-sided pulsation (Handler et al. 2020; Kurtz et al. 2020), the lack of correlation between these high-frequency modes and the orbital frequency makes this unlikely. Moreover, faster envelope than core rotation has been detected in several other pulsating close binaries (e.g., KIC 819776, Sowicka et al. 2017, and the inner binary of the triple system HD 201433, Kallinger et al. 2017).

One of the potential mechanisms for this effect is the ‘inverse’ tidal mechanism (Fuller 2021), where tidal interaction with unstable pulsation modes can transfer energy and angular momentum in a manner that forces the star away from synchronicity. The differential rotation of the star, coupled with the asynchrononicity of the surface and envelope rotation with respect to the orbit (f_orb = 0.364 d⁻¹) noted in Paper I, seems to reinforce this argument. However, the non-detection of tidally excited or perturbed pulsations (reported in Paper I), and the fact that close binaries with synchronous surface rotation but asynchronous core rotation exist (e.g., KIC 819776, Sowicka et al. 2017) cast doubt on this possibility, though does not rule it out. The second scenario in which a faster envelope than core rotation can develop is by angular momentum transport by internal gravity waves (Rogers 2015), which reinforces our use of mixing profiles calibrated by the theoretical simulations of internal gravity waves in our evolutionary models (see Sect. 2).

5. Discussion and conclusions

Gravity-mode period-spacing series found in the primary of the eclipsing binary KIC 9850387 allowed for high precision estimation of the stellar parameters of this intermediate-mass F-type pulsator. We coupled the period spacing patterns of identified dipole and quadrupole modes with an evolutionary and asteroseismic modelling-based analysis by comparing the observationally and theoretically derived parameters of this star. To achieve this goal, we computed a grid of evolutionary models with 0.80 M_⊙ ≤ M_ini ≤ 2.00 M_⊙ with a range of f_ov and D_mix values typical for this type of pulsator (Mombarg et al. 2019), aside from allowing for a broader range of envelope mixing levels to investigate the potential influence of tidal mixing mechanisms in the envelope. We then performed isocloud fitting in a Monte-Carlo framework similar to Johnston et al. (2019a), and obtained the evolutionary parameters based on the intersection of the dynamical T_eff and log g constraints of each component with the isocloud parameters corresponding to the 95% HPD interval of the Monte Carlo age distribution (${\tau }_{\mathrm{MC}}=1.3_{-0.2}^{+1.5}$ Gyr).

We exploited the slow-rotating nature of the primary star of KIC 9850387 by estimating Π₀ directly from the mean period-spacing $\overline{\mathrm{\Delta }{P}_{\ell }}$ of the ℓ = 1 and ℓ = 2 patterns, and then using it to confront the theoretical Π₀ values extracted from our grids of evolutionary models. Additionally, we modelled the individual ΔP_ℓ values of the observed dipole and quadrupole period-spacing patterns by constructing theoretical period-spacing patterns based on non-rotating stellar structural models. We used two different merit functions (${\chi }_{\mathrm{red}}^2$ and MD) and tested different setups based on the imposition of ‘pseudo-single-star’ spectroscopic and dynamical constraints. It was found that our asteroseismic modelling provided stronger constraints on the interior properties (M_cc, f_ov, and D_mix) of the primary than the evolutionary modelling, demonstrating the probing power of g-modes. We also found an overall agreement between the asteroseismic and evolutionary modelling results.

Our results reinforce the claim of main-sequence binary evolutionary stage made in Paper I, contradicting the conclusion of Zhang et al. (2020) that the system comprises two pre-main-sequence components. Our best-fitting models allowed for strong constraints on the parameters describing the interior mixing profile of the star, comprising a low amount of exponentially-decaying core overshooting (f_ov = 0.006) and a high amount of envelope mixing (D_mix = 25 cm² s⁻¹) for this type of pulsator (Mombarg et al. 2019). These findings led to precise constraints on the evolutionary stage of the primary (X_c,  p = 0.24 ± 0.02) and evolutionary modelling allowed for constraints on the secondary (X_c,  s = 0.57 ± 0.05), corresponding to an age of 1.2 ± 0.1 Gyr (see parameters in bold in Table 9). We also obtained constraints on the mass of the convective core, M_cc = 0.16 ± 0.01, which is within range of expectation values reported for single F-type g-mode pulsators by Mombarg et al. (2019).

We exploited the slow-rotating nature of the primary star of KIC 9850387 by estimating Π₀ directly from the mean period-spacing $\overline{\mathrm{\Delta }{P}_{\ell }}$ of the ℓ = 1 and ℓ = 2 patterns, and then using it to confront the theoretical Π₀ values extracted from our grids of evolutionary models. Additionally, we modelled the individual period spacing values of the consecutive zonal dipole and quadrupole modes via pulsation computations based on non-rotating equilibrium models. We used two different merit functions (${\chi }_{\mathrm{red}}^2$ and MD) and tested different setups based on the imposition of ‘pseudo-single-star’ spectroscopic and dynamical constraints. It was found that our asteroseismic modelling provided stronger constraints on the interior properties (M_cc, f_ov, and D_mix) of the primary than the evolutionary modelling, demonstrating the probing power of g modes. We also found an overall agreement between the asteroseismic and evolutionary modelling results.

We found little difference in the modelling results regardless of whether a ${\chi }_{\mathrm{red}}^2$ or the MD merit function was used and regardless of the grid setup or Y configurations for our Π₀-based modelling. However, the differences are significant for our ΔP_ℓ-based modelling, particularly when the full grid of models is used in the fit. The application ‘pseudo-single-star’ spectroscopic reduces these discrepancies, and the application of dynamical constraints eliminates them altogether. Due to the degenerate nature of the estimated stellar parameters, the best-fitting period spacing values and errors in the ${\chi }_{\mathrm{red}}^2$ framework are almost identical regardless of the grid subset used in the fit, while those in the MD framework vary significantly as expected from the construction of this merit function.

After identifying the overall best model combining the asteroseismic information of the dipole and quadrupole modes with the dynamical constraints, we investigated the few high-frequency modes that were identified in Paper I by calculating the theoretical frequencies of the ℓ = 0 to ℓ = 4 modes in the 12.8 d⁻¹ to 13.4 d⁻¹ frequency range for our best model. The fundamental radial mode and the g₁ mode were found within this frequency range, and frequency offsets of 0.05 d⁻¹ and 0.03 d⁻¹ between these modes and the nearest observed mode were obtained. We posited that these frequency shifts were due to the surface effect (e.g., Ball 2017), following the explanation of similar behaviour in the binary pulsator KIC 10080943 by Schmid & Aerts (2016).

We investigated the hypothesis that the observed high-amplitude frequency peaks at 13.21 d⁻¹ and 13.25 d⁻¹ near the theoretical dipole g₁ mode are a part of a rotationally-split prograde-retrograde doublet with a missing zonal mode, as was found in several p- and g-mode hybrid pulsators (e.g., Kurtz et al. 2014). This corresponds to an envelope rotational of 0.01689 ± 0.00004 d⁻¹ that is thrice as high as the core rotational frequency 0.0053 ± 0.0015 d⁻¹ (Li et al. 2020b). Within the limitations of the data, the surface rotation is compatible with the envelope rotation. Similar behaviour has been observed for other close binaries. For single pulsators of the same mass, the faster envelope than core rotation has been explained in terms of angular momentum transport by internal gravity waves triggered by the convective core (Rogers 2015). Such a mechanism may also be active within the primary star of this binary. However, we reinforce that these conclusions are based on the envelope rotational frequency derived from the splitting of a single frequency doublet with posited identification as a dipole mode, and as such is subject to uncertainty.

Overall, we find that asteroseismic theory and observations are only barely compliant, reinforcing the need for homogeneous analyses of samples of pulsating eclipsing binaries that aim at calibrating interior mixing profiles. Such studies would allow for the investigation of the sources of discrepancy between the various parameters, and address weaknesses in the descriptions of angular momentum transport and interior mixing mechanisms.

Acknowledgments

The research leading to these results has received funding from the Fonds Wetenschappelijk Onderzoek – Vlaanderen (FWO) under the grant agreements G0H5416N (ERC Opvangproject) and G0A2917N (BlackGEM), and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 670519: MAMSIE) and from the KU Leuven Research Council (grant C16/18/005: PARADISE). 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. The authors are much indebted to the MESA and GYRE development teams for making their user-friendly software publicly available. The authors would also like to thank the Leuven MAMSIE team for useful discussions. Last but not least, we would like to acknowledge R. H. D. Townsend for his comments on an earlier version of the manuscript as part of the PhD jury of SS, which led to a more thorough investigation and subsequent refinement of the internal differential rotation results for this star.

References

Appendix A: Correlation plots of the MD values in the period-spacing modelling

Fig. A.1. Correlations between model parameters for our asteroseismic modelling with Y = (ΔPℓ = 1,  i). The points in each subplot are colour-coded according to the MD values of the whole grid. The white circle and error bars represent the maximum-likelihood estimates and half of the 95% HPD of the Monte Carlo parameter distributions.

Fig. A.2. Correlations between model parameters for our asteroseismic modelling with Y = (ΔPℓ = 2,  j). The points in each subplot are colour-coded according to the MD values of the whole grid. The white circle and error bars represent the maximum-likelihood estimates and half of the 95% HPD of the Monte Carlo parameter distributions.

Fig. A.3. Correlations between model parameters for our modelling with Y = (ΔPℓ = 1,  i, ΔPℓ = 2,  j). The points in each subplot are colour-coded according to the MD values of the whole grid. The white circle and error bars represent the maximum-likelihood estimates and half of the 95% HPD of the Monte-Carlo parameter distributions.

Appendix B: Covariance matrices used in the MD calculations in the period-spacing modelling

Fig. B.1. Covariance matrices (V + Λ, see Eq. (4)) of the individual and combined mode-fitting setups, representing the 19 ΔPℓ = 1 values, the 10 ΔPℓ = 2 values of the ℓ = 2 mode, and the 29 total ΔPℓ = 1 and ΔPℓ = 2 values.

Appendix C: Displacement and rotational-kernel plots for the ℓ = 1 and ℓ = 2 modes of KIC 9850387

Fig. C.1. Left and middle panels: radial (ξr(R)) and horizontal (ξh(R)) components of the Lagrangian displacement vectors. Right panel: rotational kernels and the Brunt–Väisälä frequency (N) in the near-core region in which the g modes have the highest probing power. These displacement vectors and kernels are of the lowest and highest radial order modes of the theoretical period-spacing pattern of our overall best model (M = 1.65 M⊙, Xc = 0.24, fov = 0.005, and log Dmix = 1.5). A full derivation of these quantities is provided in Unno et al. (1989) and Aerts et al. (2010).

Appendix D: Propagation diagrams of the ℓ = 1 and ℓ = 2 modes of KIC 9850387

Fig. D.1. Propagation diagrams showing the mode cavities, frequencies and nodes of the theoretical period-spacing pattern of our overall best model (M = 1.65 M⊙, Xc = 0.24, fov = 0.005, and log Dmix = 1.5). The thick solid, dashed and dotted curves represent the Brunt–Väisälä frequency (N), and the ℓ = 1 and ℓ = 2 Lamb frequencies (Sℓ = 1 and Sℓ = 2) respectively. The red region represents the ℓ = 1 g mode cavity that gets extended by the pink region for the ℓ = 2 g-modes. The blue region represents the ℓ = 2 p-mode cavity that gets extended by the light blue region for the ℓ = 1 p modes. The horizontal black and grey dotted lines represent the 20 ℓ = 1 and 11 ℓ = 2 theoretical frequencies, with the nodes of radial Lagrangian displacement vector for each frequency represented by black circles and grey triangles respectively. Right panel: magnification of the frequency region in which the theoretical modes propagate.

All Tables

All Figures

Fig. 1. Isocloud fit to the components of KIC 9850387 on the log Teff − log g diagram. The red regions correspond to the best-fitting isocloud, and the pink regions represent the isoclouds in the 95% HPD interval of the fit. The black regions on the isocloud represent the log Teff − log g values corresponding to the dynamical masses of the individual components. The grey curves are MESA evolutionary tracks with fov = 0.005 and log Dmix = 0.0 (solid), fov = 0.04 and log Dmix = 0.0 (dashed), and fov = 0.04 and log Dmix = 4.0 (dotted) with their corresponding masses (in units of M⊙) indicated at the ZAMS of each track. The inset plots are magnified regions around the position of the primary (top) and secondary (bottom) component.

In the text

	Fig. 2. Example of the distribution of M values in the posterior distribution of models selected through our Monte-Carlo fitting framework. The vertical dashed line represents the median value and the vertical solid lines represent the upper and lower bounds of 68% HPD of the M values.
In the text

Fig. 3. Best-fitting Π0 values of our Π0-based asteroseismic modelling, based on three different grid setups (full, spectroscopic, and dynamical). Left panel: results of using a ${\chi }_{\mathrm{red}}^2$ merit function and right panel: results of using an MD merit function. The vertical dashed black line represents the zero-point of the difference between the observational and theoretical Π0 (i.e. Π0, obs − Π0, theo = 0). The error bars are based on 68% HPD intervals of the Monte-Carlo parameter distributions.

In the text

Fig. 4. Positions of the best-fitting models from Π0-based asteroseismic modelling on Kiel diagrams. The main plots show the results of using the whole grid in the fit, while the inset plots show the results of applying the ‘pseudo-single-star’ spectroscopic Teff and log g (solid box) and dynamical M, Teff, and log g (dotted box) constraints. The observed positions of the star according to the spectroscopic and dynamical parameters are represented by black ‘X’ symbols and stars respectively. The 1σ and 3σ spectroscopic and dynamical error bars are represented with straight end-caps and ball end-caps respectively. The Yobs = (Π0,  ℓ = 1), Yobs = (Π0,  ℓ = 1), and Yobs = (Π0,  ℓ = 1, Π0,  ℓ = 2) solutions for the different grid subsets are represented by the same colours and symbols as in Fig. 3. Left panel: results of using a ${\chi }_{\mathrm{red}}^2$ merit function and right panel: results of using an MD merit function. The error bars on the asteroseismic solutions are based on 68% HPD intervals of the Monte Carlo parameter distributions. The solid and dashed grey curves are the same as in Fig. 1.

In the text

	Fig. 5. Same as Fig. 3, but representing the results of ΔPℓ-based asteroseismic modelling.
In the text

	Fig. 6. Best-fitting ℓ = 1 and ℓ = 2 period-spacing patterns using our ΔPℓ-based asteroseismic modelling. Top panel: results of using a ${\chi }_{\mathrm{red}}^2$ merit function and bottom panel: results of using an MD merit function.
In the text

	Fig. 7. Same as Fig. 4, but representing the results of ΔPℓ-based asteroseismic modelling.
In the text

	Fig. 8. Pulsational model with the closest frequency match between the observed and theoretical modes in the high-frequency regime of KIC 9850387. The vertical dotted green lines represent the orbital harmonics.
In the text

	Fig. 9. Pulsational model with the closest frequency match between the theoretical ℓ = 1 mode and the midpoint of the high-amplitude doublet at 13.21 d−1 and 13.25 d−1. The vertical dotted green lines represent the orbital harmonics.
In the text

	Fig. 10. Rotational kernels of the high- and low-order modes of the best-fitting pulsational model as a function of the fractional radius R/R⋆.
In the text

Fig. 11. Core (taken from Li et al. 2020b, based on the slope of the g-mode period-spacing pattern), envelope (from the rotational splitting of the g1 mode) and surface (from spectroscopic line broadening ignoring pulsational velocity broadening) rotational frequencies (frot) as a function of the fractional radius (R/R⋆) of the primary component of KIC 9850387. The core and envelope rotational frequencies are positioned at the maxima of the rotational kernels of the dipole ng = 39 and ng = 1 modes represented in Fig. 10. The errors on the core and envelope rotational frequencies are smaller than the symbol size. The dashed black line represents the rotational profile from the core to the surface that is compatible with macroturbulent line broadening due pulsations.

In the text

	Fig. 12. Synthetic spectral fits to a Fe I line of the disentangled primary component spectrum with different compositions of rotational broadening. The solid and dotted black lines represent synthetic spectra with identical atmospheric parameter inputs but with v sin i = 13 km s−1 and vmacro = 0 km s−1 and v sin i = 2 km s−1 and vmacro = 15 km s−1 respectively.
In the text

	Fig. A.1. Correlations between model parameters for our asteroseismic modelling with Y = (ΔPℓ = 1,  i). The points in each subplot are colour-coded according to the MD values of the whole grid. The white circle and error bars represent the maximum-likelihood estimates and half of the 95% HPD of the Monte Carlo parameter distributions.
In the text

	Fig. A.2. Correlations between model parameters for our asteroseismic modelling with Y = (ΔPℓ = 2,  j). The points in each subplot are colour-coded according to the MD values of the whole grid. The white circle and error bars represent the maximum-likelihood estimates and half of the 95% HPD of the Monte Carlo parameter distributions.
In the text

	Fig. A.3. Correlations between model parameters for our modelling with Y = (ΔPℓ = 1,  i, ΔPℓ = 2,  j). The points in each subplot are colour-coded according to the MD values of the whole grid. The white circle and error bars represent the maximum-likelihood estimates and half of the 95% HPD of the Monte-Carlo parameter distributions.
In the text

	Fig. B.1. Covariance matrices (V + Λ, see Eq. (4)) of the individual and combined mode-fitting setups, representing the 19 ΔPℓ = 1 values, the 10 ΔPℓ = 2 values of the ℓ = 2 mode, and the 29 total ΔPℓ = 1 and ΔPℓ = 2 values.
In the text

Fig. C.1. Left and middle panels: radial (ξr(R)) and horizontal (ξh(R)) components of the Lagrangian displacement vectors. Right panel: rotational kernels and the Brunt–Väisälä frequency (N) in the near-core region in which the g modes have the highest probing power. These displacement vectors and kernels are of the lowest and highest radial order modes of the theoretical period-spacing pattern of our overall best model (M = 1.65 M⊙, Xc = 0.24, fov = 0.005, and log Dmix = 1.5). A full derivation of these quantities is provided in Unno et al. (1989) and Aerts et al. (2010).

In the text

Fig. D.1. Propagation diagrams showing the mode cavities, frequencies and nodes of the theoretical period-spacing pattern of our overall best model (M = 1.65 M⊙, Xc = 0.24, fov = 0.005, and log Dmix = 1.5). The thick solid, dashed and dotted curves represent the Brunt–Väisälä frequency (N), and the ℓ = 1 and ℓ = 2 Lamb frequencies (Sℓ = 1 and Sℓ = 2) respectively. The red region represents the ℓ = 1 g mode cavity that gets extended by the pink region for the ℓ = 2 g-modes. The blue region represents the ℓ = 2 p-mode cavity that gets extended by the light blue region for the ℓ = 1 p modes. The horizontal black and grey dotted lines represent the 20 ℓ = 1 and 11 ℓ = 2 theoretical frequencies, with the nodes of radial Lagrangian displacement vector for each frequency represented by black circles and grey triangles respectively. Right panel: magnification of the frequency region in which the theoretical modes propagate.

In the text