© ESO 2021

1. Introduction

Asteroseismology, the study of stellar interiors by means of detected identified oscillation modes (Aerts et al. 2010), has two major aims. The first goal is to provide masses, radii, and ages of stars relying on current knowledge of stellar structure, providing these three quantities with far better precision than is available from single-epoch data of surface properties (such as spectroscopy or interferometry). The second goal is to improve the input physics of stellar interiors so as to achieve better models for stellar evolution. The need for better models is obvious when it comes to angular momentum transport across stellar evolution (Gehan et al. 2018; Aerts et al. 2019) and mixing that guides the change in the mass fractions of the chemical species (Salaris & Cassisi 2017). Here our aim is to probe the temperature and chemical structure in the deep interior of stars and to improve our models of those quantities, since these are aspects of the stellar structure that propagate throughout and influence the star’s evolution.

Although angular momentum and element transport are two different physical phenomena, they are intimately related, and as such they are treated by similar diffusion coefficients in stellar evolution codes (Maeder 2009). To date, asteroseismology has mainly been used to evaluate and improve the transport of angular momentum (e.g. Fuller et al. 2019; Eggenberger et al. 2019a,b; den Hartogh et al. 2020; Deheuvels et al. 2020) and to assess the level of chemical mixing (e.g. Deheuvels et al. 2016; Buldgen et al. 2019; Bellinger et al. 2019; Angelou et al. 2020), with a specific focus on red giants and on stars of low mass. These studies are all based on acoustic or mixed modes excited by envelope convection. In the current work we rely on gravito-inertial modes in a rotating star of intermediate mass with a well-developed convective core and a radiative envelope. These modes are excited by a heat mechanism (Pamyatnykh 1999) and have frequencies that are similar to the rotation frequency of the star, such that the modes’ spin parameters s = 2Ω/ω > 1, with Ω and ω the rotation and oscillation frequencies. In this case the Coriolis and buoyancy forces act together as restoring forces of the oscillations. These oscillations correspond to gravito-inertial modes (see Aerts et al. 2019, Fig. 5). In particular, we focus on a B-type pulsator. The extent of and mass in the rotating core-boundary layers are critical unknowns for this type of stars in the current theory of stellar structure, yet these properties drive their evolution. For such stars, the asteroseismic modelling cannot rely on a perturbative treatment of the Coriolis force, in contrast to the case of acoustic or mixed modes in slowly rotating stars (see e.g. Ballot et al. 2010, for the limits of the perturbative method).

Despite their interest from a stellar evolution and chemical yield point of view, few B-type stars with gravito-inertial modes have been modelled asteroseismically, and the focus has been on exceptionally slow rotators (Moravveji et al. 2015; Wu et al. 2020). On the other hand, fast rotators offer the more interesting case to probe stellar interiors with the aim of improving stellar evolution theory (Moravveji et al. 2016; Buysschaert et al. 2018; Szewczuk & Daszyńska-Daszkiewicz 2018). One of the reasons for the limited number of modelled cases is that only the Kepler space telescope (Borucki et al. 2010) offered sufficiently long time-series photometry to do so, and massive stars were avoided by it (see Bowman 2020, for a recent review). While time series covering a few months allow us to detect high-order gravito-inertial modes in main-sequence stars and to identify their mode wavenumbers (l, m) from period spacing patterns (Degroote et al. 2010a; Pápics et al. 2012; Kallinger et al. 2017), only the 4 yr Kepler light curves of B-type pulsators lead to sufficient frequency precision to discriminate among various profiles representing convective boundary mixing (hereafter CBM; Pedersen et al. 2018). Another reason for the scarcity of asteroseismically modelled B-type pulsators is that few of those observed by Kepler have delivered suitable modes to perform modelling of their interiors. Finally, the methodology used to perform the modelling in the case of gravito-inertial modes is more complex than the framework based on modes in the perturbative regime of rotation (see Aerts 2021, for a review of the various modelling regimes).

Several prescriptions to describe the temperature gradient, ∇_T(r), and the diffusion coefficient defining the local chemical mixing profile in the convective core boundary layers, D_mix(r), exist. We consider two main options: firstly convective penetration with an adiabatic temperature gradient, ∇_ad, and a constant mixing profile in the core boundary layers, following Zahn (1991), and secondly diffusive overshooting with a radiative temperature gradient, ∇_rad, and an exponentially decaying CBM profile, following Freytag et al. (1996) and Herwig (2000). Here we are concerned with convective boundary layers adjacent to a convective core in a main-sequence star. We take a data-driven approach based on observed gravito-inertial modes, which are excellent probes of the near-core regions deep inside stars, rather than relying on (uncalibrated) numerical simulations of these regions.

In a proof-of-concept study, Michielsen et al. (2019) showed that observed gravito-inertial modes hold the potential to discriminate between a radiative and an adiabatic temperature gradient, along with assessment of the mixing profile in the core boundary layers, provided that the mass, evolutionary stage, internal rotation, and metallicity are known. In this work we improve this initial theoretical study and apply it to one of the rotating B-type pulsators with a long series of prograde dipole gravito-inertial modes of consecutive radial order, KIC 7760680. Our aim is to assess the temperature and chemical gradients of the core boundary layers in this star from observations, but also to investigate the impact of various choices of asteroseismic observables in the modelling process on the inferred outcome. We do not focus on the observational aspects of the detected gravito-inertial modes because the frequency uncertainties are typically two orders of magnitude below those stemming from unknown input physics of stellar models (see Aerts et al. 2018). We thus consider the case of KIC 7760680 by relying on its observational properties derived previously by Pápics et al. (2015), and briefly recall the observational input of the modelling in the next section.

2. The gravito-inertial B-type pulsator KIC 7760680

We have opted to use the slowly pulsating B star KIC 7760680 as a case study for this paper. Our reasons for this choice are twofold. First, this is the only star rotating at a considerable fraction of its critical rotation rate (about 26%, Moravveji et al. 2016) with an extensive identified mode period pattern (36 prograde dipole modes of consecutive radial order). Second, this star has a remarkably low level of inferred mixing in its envelope compared to slower rotating B-type pulsators, while one would expect its rotation to induce strong mixing beyond the convective core. The radial shear-induced mixing in differentially rotating massive stars is about three to ten orders of magnitude higher according to theoretical predictions than the envelope mixing found by Moravveji et al. (2016). We refer to Fig. 6 in Meynet & Maeder (2000), Fig. 3 in Mathis et al. (2004), and Figs. 15 and 16 in Decressin et al. (2009) for various examples. An explanation for this low envelope mixing proposed by Moravveji et al. (2016) is that the star is nearly a rigid-body rotator.

Following up on the best model parameters for the star found by Moravveji et al. (2016) listed in Table 1, we assess whether a more structured CBM profile is meaningful and leads to a better fit to the star’s observed oscillation modes. Moravveji et al. (2016) considered ∇_rad in the CBM region, limited the parameter estimation to the use of a χ² merit function, and performed mode frequency fitting. Given the mathematical challenges connected to it, no precision estimation for the parameters in Table 1 was done in this previous study. In the case of strong correlations among the parameters to estimate, as is the case here, the very small grid step size listed in Table 1 tends to be an underestimation of the true uncertainties because this neglects the parameter correlation structure. We aim to improve upon this omission of error estimation, while searching for a more accurate interior model of the star. This improvement concerns several aspects, namely the use of more optimized asteroseismic observables, a more appropriate merit function taking the correlation structure of the fitting problem into account, and a more structured D_mix(r) profile compared to the initial study by Moravveji et al. (2016). Another reason to choose this star as a test case is that it has an excellent asteroseismic calibration of its near-core rotation frequency. Moravveji et al. (2016) derived a value of 0.4805 d⁻¹ under the assumption of rigid rotation, which corresponds to 26% of its critical Roche rotation frequency. To reduce dimensionality, we thus fixed this value of the rotation frequency.

The Kepler light curve of this pulsator and its Fourier transform have already been discussed in great detail in Pápics et al. (2015), so we do not repeat this information here for brevity. In practice, we use the period values of the 36 prograde dipole modes listed in Table 1 of Pápics et al. (2015) as asteroseismic input for our modelling, Y^obs composed of ${\text{Y}}_i^{\text{Obs}}$ with i = 1, …, 36. Pápics et al. (2015) list the formal errors from a least-squares fit to the light curve, after pre-whitening all the modes with amplitudes above four times the local noise level. Schwarzenberg-Czerny (1991) showed that the formal least-squares errors are statistically equivalent to the errors derived in the Fourier domain if one takes into account the variance in the data, the signal-to-noise ratio, and the total time base of the data set. However, the author noted that a correction factor to the formal errors is required to take into account the correlated nature of the data, following their covariance structure. In line with the method in Schwarzenberg-Czerny (1991), we multiply the errors listed in Table 1 of Pápics et al. (2015) by a correction factor of four, as was done for the modelling performed by Moravveji et al. (2016). This leads to errors on the mode periods lower than ∼10 s, while the time resolution restriction stemming from the light curve that may have affected the choice of frequencies during the pre-whitening process amounts to ∼60 s.

As additional input to guide the modelling, we rely on T_eff  =  11 650  ±  210 K, log g = 3.97 ± 0.08 dex, and [M/H] = 0.14 ± 0.09 dex, again taken from Pápics et al. (2015). We also rely on log (L/L_⊙) = 2.19 ± 0.06 based on Gaia DR2 astrometry deduced by Pedersen et al. (2020).

3. Chemical mixing profiles

In the diffusion approximation, the transport equation describing the rate of change in the mass fraction of each individual chemical element i in the considered chemical mixture, denoted here as X_i, reads

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial X_i(r)}{\partial t}={\mathcal{R}}_i(r)+\frac{1}{\rho (r)r^2}\frac{\partial }{\partial r}\left[\rho (r)r^2{D}_{\mathrm{mix}}(r)\frac{\partial X_i(r)}{\partial r}\right],}\end{array}$$(1)

with ℛ_i(r) the change due to nuclear reactions and D_mix(r) the local mixing coefficient. In this work this mixing coefficient in the convective boundary layers is implemented as extended convective penetration described by Michielsen et al. (2019), with an extra diffusive mixing profile in the radiative envelope due to internal gravity waves as deduced by Rogers & McElwaine (2017), following Pedersen et al. (2018). The mixing profile outside of the regions determined to be convective by the Ledoux criterion can hence be described by

$$\begin{array}{ccc}& D_{\mathrm{CBM}}(r)=D_0\hfill & \hfill \mathrm{for}{r}_0<r<r_{\mathrm{cp}},\end{array}$$(2)

$$\begin{array}{ccc}\hfill & D_{\mathrm{CBM}}(r)=D_{0}exp\left(\frac{-2(r-r_{\mathrm{cp}})}{f_{\mathrm{CBM}}{H}_{\mathrm{p}}}\right)\hfill & \hfill \mathrm{for}{r}_{\mathrm{cp}}<r<r_{\mathrm{CBM}},\end{array}$$(3)

$$\begin{array}{ccc}\hfill & D_{\mathrm{env}}(r)=D_{\mathrm{env}}\frac{\rho (r_{\mathrm{CBM}})}{\rho (r)}\hfill & \hfill \mathrm{for}{r}_{\mathrm{CBM}}<r,\end{array}$$(4)

with r₀ = r_cc − f₀H_p being the radius where the transition between core mixing and core boundary mixing is made, r_cc being the radius of the convective core as obtained from the Ledoux criterion applied for mixing length theory (MLT; Böhm-Vitense 1958), and H_p the local pressure scale height; D₀ is the mixing coefficient at r₀, and is hence influenced by the choice of f₀. The edge of the step-like CBM in the transition region is given by r_cp = r₀ + α_CBMH_p. Beyond this value the mixing profile takes on an exponentially decaying function with free parameter f_CBM until the core boundary transition region reaches the envelope of the star, whose bottom is situated at r_CBM and where the envelope mixing becomes active rather than CBM. This corresponds to the radius where the convective boundary mixing coefficient becomes smaller than D_env, which stands for the mixing at the bottom of the radiative envelope. The latter mixing coefficient takes over as of this position towards the surface of the star. The free parameters α_CBM and f_CBM determine the extent of the step and exponentially decaying profile part of D_mix(r) in the core boundary region, respectively.

Employing this prescription for the mixing profile, two cases are considered regarding the temperature gradient in the core boundary transition region. The first one (seen in Fig. 1a) adopts the radiative temperature gradient, ∇_T = ∇_rad, outside the convective core (r > r₀). The second one (seen in Fig. 1b) assumes that the overshooting material influences the entropy stratification in the transition layers, as proposed by Zahn (1991). Instead of a transition in temperature gradient based on the mass coordinate in the convective boundary mixing region used in Michielsen et al. (2019), we opt for a transition based on the Péclet number, as this is more physically motivated (see e.g. Viallet et al. 2015, for an overview of the different convective regimes). The Péclet number is defined as the ratio of the advective transport rate to the diffusive transport rate, such that the flow can be considered diffusive when Pe ≪ 1 and adiabatic when Pe ≫ 1. Bridging these two regimes was implemented by making a gradual transition based on the Péclet number, in terms of a factor h = (log(Pe)+2)/4. This factor was constructed to create a small transition region between the regimes instead of introducing a discontinuity at the interface. This implementation causes Pe = 1 to correspond to the temperature structure being halfway between ∇_ad and ∇_rad:

$$\begin{array}{c}\hfill \begin{array}{cc}{\mathrm{\nabla }}_{\mathrm{T}}={\mathrm{\nabla }}_{\mathrm{ad}}\hfill & \mathrm{for}\mathrm{Pe}>100,\hfill \\ {\mathrm{\nabla }}_{\mathrm{T}}=h{\mathrm{\nabla }}_{\mathrm{ad}}+(1-h){\mathrm{\nabla }}_{\mathrm{rad}}\hfill & \mathrm{for}0.01<\mathrm{Pe}<100,\hfill \\ {\mathrm{\nabla }}_{\mathrm{T}}={\mathrm{\nabla }}_{\mathrm{rad}}\hfill & \mathrm{for}\mathrm{Pe}<0.01.\hfill \end{array}\end{array}$$(5)

Fig. 1. Radial profiles of a 3 M⊙ star halfway through the main sequence with a central hydrogen content Xc = 0.4. Top panels: temperature gradients. Middle panels: structure (solid) and composition (dashed) components of the Brunt−Väisälä profile, 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). Bottom panels: mode inertia of two g-modes with different radial orders. (a) Near-core region with radiative temperature gradient. (b) Near-core region with temperature gradient according to Péclet number.

One of our goals is to answer whether the 36 detected dipole gravito-inertial modes of KIC 7760680 allow us to distinguish between these two options of temperature gradients in the transition layer between its convective core and radiative envelope, adopting the observational diagnostics of the star from Pápics et al. (2015).

4. Computation of theoretical mode frequencies

4.1. Stellar equilibrium models

Two grids of stellar models were computed as input for pulsation computations. These two grids provide equilibrium models with the same input physics, except for the adopted temperature gradient in the core boundary transition region. One grid adopts the radiative temperature gradient in that transition zone and is termed the radiative grid. For the other grid we adopt the temperature gradient according to Eq. (5) and we call it the Péclet grid. Both grids are identical in terms of the free parameter ranges listed in Table 2. The ranges of the mass and central hydrogen fraction were guided by the earlier study of Moravveji et al. (2016), while the metallicity range relied on the spectroscopic results taken from Pápics et al. (2015).

The grids were computed using the stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) version r12115. The MESA models contain the standard chemical mixture of OB stars in the solar neighbourhood deduced by Nieva & Przybilla (2012) and Przybilla et al. (2013) and make use of the OP opacity tables (Seaton 2005). An Eddington grey atmosphere is used as atmospheric boundary condition. To set the initial mass fractions (X_ini, Y_ini, Z_ini) we rely on the initial metallicity Z_ini as derived from the spectroscopic measurements derived by Pápics et al. (2015), and vary it along its uncertainty range. The initial helium fraction is then determined by adopting an enrichment law Y_ini = Y_p + (ΔY/ΔZ)Z_ini, where the primordial helium abundance, Y_p, and the galactic enrichment ratio, $\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}$, come into play. There is currently no consensus on the value of $\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}$ (see e.g. Verma et al. 2019, and references therein). Therefore, we have taken Y_p = 0.2465, as determined by Aver et al. (2013), and we required that the enrichment ratio was 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 led 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 refer to Nsamba et al. (2020) for an in-depth comparison between this approach of using an enrichment ratio, and including the initial helium abundance as a free variable.

To determine where the transition from core to near-core mixing is made, we fixed f₀ = 0.005, with the exception of setting f₀ = 0 for the models where both α_CBM and f_CBM are equal to zero as there is no CBM region for this case. The mixing length parameter was set to α_mlt = 2.0 in the mixing length theory as developed by Cox & Giuli (1968). We used the Ledoux criterion for convection without allowing for semi-convection, since this form of slow mixing does not occur in the presence of convective boundary mixing (see e.g. Kaiser et al. 2020), which is included in the vast majority of models in our grids. Appendix A provides the link to the detailed MESA setup.

4.2. Pulsation computations

To model the 36 consecutive dipole g modes detected in KIC 7760680, the pulsation mode properties of the MESA equilibrium models were computed employing the stellar oscillation code GYRE (Townsend & Teitler 2013; Townsend et al. 2018), version 5.2. We computed the dipole g-modes for a fixed rotation frequency of 0.4805 d⁻¹, assuming rigid rotation as discussed previously. We relied on the traditional approximation of rotation (TAR, following the implementation by Townsend et al. 2018) and adopted the adiabatic approximation to compute the dipole gravito-inertial mode frequencies. The mode inertias of the considered g-modes are dominant near the core of the star, as illustrated in Fig. 1. Hence the adiabatic approximation is sufficient for our modelling work, whilst non-adiabatic effects mainly become important in the outer stellar envelope.

Figure 1 reveals the difference in temperature gradient between two models from the two MESA grids. Although it is relatively small and leaves the composition term of the Brunt−Väisälä frequency ($N_{\mathrm{composition}}=\frac{g^2\rho }{P}{\mathrm{\nabla }}_{\mathrm{\mu }}$) unchanged, it does change the thermal structure term ($N_{\mathrm{structure}}=\frac{g^2\rho }{P}$ (∇_ad − ∇_T)), which causes a different mode trapping. This is illustrated for the gravito-inertial modes with radial orders 25 and 45.

It is this type of difference in the mode kernels that brings about the potential to assess the temperature gradient in an observed star, provided that it delivers suitable modes to make this assessment. It was shown by Moravveji et al. (2016, see their Fig. 10) that KIC 7760680 is a suitable candidate to offer this opportunity.

The final GYRE inlist is provided through the link in Appendix A and provides us with a list of theoretically predicted prograde dipole mode frequencies, Y^Theo composed of ${\text{Y}}_j^{\text{Theo}}$, where j stands for the radial order, for each of the equilibrium models in the two grids.

5. Modelling approaches

The asteroseismic modelling of a star involves various aspects, some of statistical nature and others connected with astrophysical considerations. The problem comes down to an estimation effort to assess the free parameters covered in the two model grids listed in Table 2. We achieve this by relying on the list of observed gravito-inertial mode frequencies, or any other observables derived from them and denoted as Y^Obs with ${\text{Y}}_i^{\text{Obs}}$ labelled as i = 1, …, 36, where the radial order of the observed modes is unknown a priori. These observables are matched to a list of theoretically predicted analogues Y^Theo, with ${\text{Y}}_j^{\text{Theo}}$ to be computed from a list of model frequencies of radial order j = 1, …, 100. There are various aspects of this matching process that have to be considered. Ultimately it comes down to solving regression problems, where various choices can be made for the observables to match and the merit functions to minimize. We discuss each of these aspects in the following subsections, noting here that we do not use the spectroscopic and astrometric information as observables in the fitting procedure, but use the measured spectroscopic abundances to determine the range of initial metal and helium mass fractions.

In the case of gravito-inertial modes of B-type pulsators with 4 yr Kepler light curves, we are in the situation where non-seismic observables from spectroscopy or astrometry, such as T_eff, log g, or log (L/L_⊙) have relative observational uncertainties that are typically two orders of magnitude or more above those of the measured frequencies (see Table 1 in Aerts et al. 2019). For this reason, these classical observables are not included in the merit functions, but rather serve to eliminate stellar models in the two grids having values of these non-seismic observables outside the observed regime (see e.g. Buysschaert et al. 2018).

5.1. Model selection criterion

The asteroseismic modelling involves maximum likelihood estimation based on the chosen observables and their theoretically predicted counterparts computed from astrophysical models. The modelling may involve nested or non-nested statistical models (see Claeskens & Hjort 2008). We are dealing with both cases. Firstly a comparison between the performance of models within one grid of equilibrium models where one or both of the CBM parameters are fixed at zero (α_CBM, f_CBM, or both) versus a free parameter to estimate. Such a comparison constitutes a nested model comparison allowing the most appropriate statistical model to be selected. We are doing this model selection as a test to evaluate whether the added complexity in the core boundary transition region of the models improves the fit to the data more than the entailed punishment by the selection criterion for having additional free parameters. Secondly a comparison between the performance of the radiative versus Péclet grids while estimating the same free parameters for both grids deals with a model selection based on non-nested models as only the input physics changes, but not the degrees of freedom.

In the case of model selection for nested statistical models, we use the Akaike Information Criterion corrected for small sample size (AICc, Claeskens & Hjort 2008, Chap. 2). This type of model selection rewards fit quality, but penalizes complexity. Moreover, the AICc is an appropriate criterion to evaluate nested models with a limited number of observables yet considerable dimensionality of the free parameter vector, as we have here (36 modes to estimate six free parameters, see Table 2). The AICc is defined as

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

with N and k being the number of observables and free parameters, respectively, and ℒ being the likelihood of a stellar model.

In our framework of asteroseismic modelling, N = 36 or 35 depending on whether we choose to fit mode frequencies, periods, or period spacings. The number of free parameters k is according to the parameters varied in the grids and is 4, 5, or 6 in our case, which means that (M_ini, Z_ini, α_CBM, f_CBM, log(D_env),X_c) are used to assess the general performance of the two grids, and the same list but with only one or neither of (α_CBM, f_CBM) included to test if increased CBM complexity performs better or worse within one grid. In this latter application, we make use of the property that two nested models A and B can be compared in performance through their difference in AICc; if ΔAICc = AICc_A − AICc_B > 2, model B is favoured over model A, with the evidence being (very) strong if ΔAICc > 6 (10) (Kass & Raftery 1955; Claeskens & Hjort 2008). These regimes are derived from the value of 2ln[P(D|A)/P(D|B)], with P(D|X) the probability of observing D given model X. A difference in ΔAICc of 6 or 10 implies that the probability of reproducing the data better with model A than with model B is roughly 95% or 99.3%, respectively.

5.2. Merit functions: Mahalanobis distance versus χ²

Uncalibrated physical ingredients required as input physics of stellar evolution models are usually dealt with via the inclusion of free parameters. Often these ingredients are correlated. Thus, the systematic uncertainties originating from these limitations in the model input physics cause a variance–covariance structure in the theoretically predicted observables. This, as well as non-linear correlations among the observables and hence in the overall correlated nature of the estimation problem to solve, has to be addressed in the regression problem.

Since the beginnings of space asteroseismology, Bayesian inferences have been developed and applied in the modelling of solar-like oscillations of low-mass stars (e.g. Appourchaux et al. 2009; Serenelli et al. 2013; Silva Aguirre et al. 2013, 2015, 2017; Metcalfe et al. 2014). In such applications the Coriolis acceleration is ignored in the prediction of the oscillation frequencies. The asteroseismic observables of rotating stars with a convective core and gravito-inertial modes demand the inclusion of the Coriolis acceleration in a non-perturbative way, implying that the g-mode frequencies are strongly dependent on the rotational frequency. Moreover, the observables are also correlated in a non-linear way, for example the modes are often strongly trapped in the near-core region (Pápics et al. 2017). In addition, in the mass regime under study here, several of the free parameters of the equilibrium models to be estimated are also strongly correlated, following tight and sometimes non-linear relationships, such as between mass and metallicity, convective-core-mass and mass, and convective-core-mass and age (e.g. Maeder 2009).

The non-linear correlations among the free parameters and observables, along with the inclusion of the uncertainties for the theoretical predictions in the estimation problem lend itself to use a more appropriate merit function instead of the traditional χ² (whose results we include in Appendix B for comparison). Here we take the full variance–covariance structure among the observables and the theoretical predictions into account, thereby mitigating the aforementioned shortcomings. This is achieved by building a regression model via the use of the so-called Mahalanobis distance (see Aerts et al. 2018, for its application to asteroseismic modelling),

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

where Y^obs is the vector of observables, ${\mathbf{Y}}_j^{\mathrm{theo}}$ is the corresponding vector of predicted values in gridpoint j of the grids of equilibrium models, V is the variance–covariance matrix of Y^theo due to the unknown physical ingredients in the input physics of the equilibrium models, and Σ is the variance due to the measurement errors of Y^obs. A major advantage of this approach is that it allows the inclusion of theoretical uncertainties for the predictions Y^theo used in the modelling.

The likelihood function corresponding to the Mahalanobis distance to be used in the statistical model selection is assumed to adhere to normality (Aerts et al. 2018, for arguments why this is a good approach). Using the likelihood function of the multivariate normal distribution, we find the likelihood of parameters θ = (θ¹, θ², …θ^k) given the observed data D,

$$\begin{array}{cc}\hfill \mathcal{L}(\mathit{\theta }|\mathbf{D})& =exp(\frac{1}{2}(ln(|V+\mathrm{\Sigma }|)+kln(2\pi )\hfill \\ \hfill & +{({\mathbf{Y}}^{\mathrm{theo}}-{\mathbf{Y}}^{\mathrm{obs}})}^{\mathrm{T}}{(V+\mathrm{\Sigma })}^{-1}({\mathbf{Y}}^{\mathrm{theo}}-{\mathbf{Y}}^{\mathrm{obs}})))\hfill \\ \hfill & =exp(-\frac{1}{2}(ln(|V+\mathrm{\Sigma }|)+kln(2\pi )+\mathrm{MD})),\hfill \end{array}$$(8)

which is proportional to the probability density function

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{\theta }|\mathbf{D})\propto \mathcal{P}(\mathbf{D}|\mathit{\theta }).\end{array}$$(9)

We can write the AICc in this case as

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

To determine the uncertainty region of the best solution, we use the likelihood in combination with 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}$$(11)

Here the index j is summed over all Q equilibrium models in the grid that are consistent within 3σ of the spectroscopic observables, and index i is summed over the q models with the highest likelihood so that $P({\theta }_a^m<{\theta }^m<{\theta }_b^m|\mathbf{D})=0.95$. Sometimes only one equilibrium model in the grid fulfils the requirement to reach this probability of 0.95. In that case we consider the uncertainty region to span a smaller area than the one defined by the grid step. This leaves us with the step sizes of the grid in each of its dimensions as an upper limit of the uncertainty region in that particular case.

5.3. Observables to fit

We investigate which is the optimal set of observables to be used in the stellar modelling procedure. In addition to just the mode periods, we also investigate the period spacing values (i.e. the differences in period between two pulsation modes of consecutive radial order, ΔP_n ≡ P_n + 1 − P_n). The best of these sets of observables is determined from the condition numbers of the variance–covariance matrices V + Σ following the Mahalanobis distance definition.

In general, the condition number κ of a matrix A is defined as the ratio of its maximum to minimum eigenvalue,

$$\begin{array}{c}\hfill \kappa (\mathbf{A})=\frac{|{\lambda }_{\max }(\mathbf{A})|}{|{\lambda }_{\min }(\mathbf{A})|},\end{array}$$(12)

and gives a handle on how well- or ill-conditioned the matrix is with respect to the inversion to be computed. A rule of thumb is that log(κ) gives an estimate of the number of digits of accuracy lost in addition to the loss due to arithmetic methods. Given this, we search for the vector of asteroseismic observables that gives the lowest condition numbers.

5.4. Constructing the theoretical period spacing pattern

When modelling an observed pattern of gravito-inertial modes of consecutive radial order, the theoretical mode periods are typically matched to the observed ones, as inspired by asymptotic relations for gravito-inertial modes (Aerts et al. 2010). However, different choices can be made to determine the starting point to build the mode pattern. Three such choices are compared in modelling performance in this paper. In the first we begin matching mode periods starting from the theoretical period that is closest to the highest-frequency mode detected in the observed pattern. A second option is to build the pattern starting from the mode with the highest observed amplitude. The third method 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. In the case of multiple mode series with the same length, a final pattern selection is made based on the best match between theory and observations. We henceforth refer to these three options of pattern construction as highest frequency, highest amplitude, and longest sequence, respectively.

6. Modelling results

6.1. Selected observables to fit

When calculating the condition number of the variance–covariance matrices according to Eq. (12), we find values on the order of 10⁶ when we solely consider the observed mode periods to fit. The variance–covariance matrices for this case, which express the coverage of the mode periods across the model grids, are shown for the two grids, and for the three ways of constructing the period spacing pattern in Figs. D.1–D.6. Variances in the theoretical mode period predictions across the entire grids reach up to ∼3400 s in this case. This is about 50 times larger than the observational uncertainties for the mode periods due to the time resolution and noise properties added together.

On the other hand, the variance–covariance matrices for the case of theoretically predicted period spacing values shown in Figs. D.7–D.12 cover up to ∼400 s in variance for the theoretical predictions of the spacings, which is about 20 times larger than the formal uncertainties on the observed period spacing values quoted by Pápics et al. (2015) and Moravveji et al. (2016). When we fit the period spacing values in the detected pattern instead of the periods, the condition numbers significantly improve, yielding κ(A) ∼ 10³. We therefore conclude that the period spacings are the optimal set of observables to fit for our case study of KIC 7760680.

The results for both sets of observables are very similar for this star in terms of the correlations and constraints on the optimal solutions. The correlation structures of the radiative grid are shown in Figs. 2 and 3 for the two sets of observables. While the correlation structures in these figures are similar, the condition numbers for the cases of the period spacing values as observables are far lower than those for the periods as observables, and hence constitute the best modelling setup in terms of numerical stability.

Fig. 2. Correlation plot for the radiative grid, using periods in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the longest sequence method. The 50% best models are shown, colour-coded according to the log of their merit function value (at right). The figures on the diagonal show binned parameter distributions, and the larger panel at the top right shows a Hertzsprung–Russell diagram with the 1 and 3σ spectroscopic Teff and astrometric luminosity error boxes.

Fig. 3. Correlation plot (see Fig. 2) for the radiative grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.

The better suitability of the period spacing values to fit as observables compared to the mode periods themselves also has advantages from the input physics point of view. The mode periods are subject to systematic uncertainty connected to various choices for the chemical mixture and opacities in B-type pulsators, as illustrated in Moravveji et al. (2015, their Fig. 3). This systematic uncertainty decreases drastically when considering mode period or mode frequency differences. Similarly, dipole gravito-inertial modes of γ Doradus pulsators computed from equilibrium models with atomic diffusion are globally shifted compared to those without taking atomic diffusion and radiative levitation into account (Mombarg et al. 2020). The mode period shift due to atomic diffusion is illustrated graphically in Fig. 5 of Aerts (2021).

We note that Moravveji et al. (2016) fitted the mode frequencies in their modelling of KIC 7760680. Just as for the case of considering mode periods, their procedure makes the modelling result more subject to uncertainties in the input physics of the equilibrium models than our new approach in this work. This is particularly so because these authors did not include the covariance among the free parameters to estimate in their modelling. The results in our paper improve the asteroseismic modelling of KIC 7760680 on these two fronts, also providing precision estimation that was not attempted in the first modelling efforts for this g-mode pulsator.

6.2. Theoretical mode period pattern construction

The Mahalanobis distance (MD) yields the same best model regardless of the method used for the theoretical pattern construction in the radiative grid. The MD and AICc values for this grid change slightly (Table C.1) due to the changed correlation structure, as can be seen in the variance–covariance matrices in Figs. D.7–D.9, but they are relatively small differences. The highest amplitude method is slightly disfavoured to compose the pattern compared to the other two in terms of their MD values.

Although the global morphology of the correlation plots remains rather similar in the Péclet grid (Figs. 4–6), the selected best model does change based on the way the theoretical pattern is constructed, as can be seen in the summary plot in Fig. 9 and in Table C.1. Looking at the differences in AICc values, the model returned using the longest sequence method is preferred over the other two methods, and the distributions of the MD values are in general shifted in favour of it (Fig. 7). This method searches for a closely matching sequence of frequencies to start building the pattern from, as opposed to the other methods that match just one frequency and build the pattern from there onward. It thus comes as no surprise that this method yields better matching patterns overall.

Fig. 4. Correlation plot (see Fig. 2) for the Péclet grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.

Fig. 5. Correlation plot (see Fig. 2) for the Péclet grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the highest amplitude method.

Fig. 6. Correlation plot (see Fig. 2) for the Péclet grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the highest frequency method.

Fig. 7. Distribution of the Mahalanobis distances using ΔP as observables to fit. The inset shows the cumulative distribution function of the 20 best models in the grids.

6.3. Radiative grid versus Péclet grid

When comparing the best models of the radiative and Péclet grid, the AICc cannot be coupled to a hypothesis test to determine how strongly one of these two astrophysical models is preferred because we are not dealing with evaluating nested regression models. We have an equal number of free parameters to fit for both of the grids. However, we can still rank the models with the same free parameters based on their AICc values. We find that the radiative grid has systematically lower AICc values than the Péclet grid, both for the best models (Table C.1) and for the distribution of the MD values as a whole (Fig. 7). Without being able to make firm statistical conclusions at the level of a significant difference, we find that the radiative grid provides a better fit to the observations.

In the next section we also consider nested regression models for both the radiative and Péclet grids. These models point out that, statistically speaking, the least complex grid without CBM is preferred (Table C.4). From this point of view, the treatment of the temperature gradient in the CBM region becomes irrelevant since the presence of such a region altogether is statistically less favourable.

6.4. Nested statistical regression models

To verify if the better fit quality outweighs the penalties for increased model complexity, we compare the AICc values of the best regression models based on the full grid with six free parameters (Table C.1) to those based on partial grids with five free parameters where either α_CBM (Table C.2) or f_CBM (Table C.3) is fixed at 0. Additionally, we compare all of these regression model fits to the results from the partial grid composed of stellar models without any presence of CBM having just four free parameters (Table C.4). Since the discovery of their period spacing patterns (Degroote et al. 2010b), slowly rotating B-type g-mode pulsators have been shown in the literature to have good solutions for a very small amount of CBM (Wu et al. 2020). This was our motivation for these tests on nested statistical models.

On the other hand, mixing in the radiative envelope was found to be a necessary ingredient for these stars (e.g. Moravveji et al. 2015). However, in none of these studies with asteroseismic modelling was a fixed zero CBM considered in a way that takes proper account of a smaller number of degrees of freedom in the regression, nor were covariances due to theoretical uncertainties considered in the fitting problem. In particular, Moravveji et al. (2016) did not consider models without CBM in their modelling of KIC 7760680.

The added complexity in the regression models is more heavily penalized statistically than rewarded by the better quality of the fit. We consistently find ΔAICc > 10 when comparing models with five free parameters to those with six, or when comparing four free parameters to five. This implies a strong statistical preference for the less complex regression models in the case of this pulsator with exceptionally low mixing beyond its convective core. This capacity of the best models with and without CBM is illustrated in Fig. 8, where we have indicated the maximum variance in the theoretical predictions according to the complete radiative and Péclet grids.

Fig. 8. Period spacing patterns of the observations, of the best models of the two grids, and of the nested grid where no CBM is included in the models, as extracted by the longest sequence method. The formal errors on the observations are smaller than the symbol sizes. The vertical bars in the bottom left corner show the maximum considered uncertainty for the theoretical predictions approximated by the variance–covariance matrix of that particular grid. Middle and bottom panels: 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 Pápics et al. (2015), including the correction factor of four.

We point out that a different conclusion would be reached if we were to rely on the χ² merit function, as is usually done in the literature of asteroseismic modelling. Appendix B offers all the results of the modelling of the prograde dipole modes in KIC 7760680 when ignoring any variance–covariance in the theoretical predictions (i.e. pretending that the theory adopted in the stellar models is fully correct) as was done by Moravveji et al. (2016). Under this (incorrect) assumption, the quality of the fit increases by adding the complexity of a CBM layer with two free parameters to the models. This is assessed quantitatively by comparing the minima listed in Tables C.5–C.8 with those in Tables C.1–C.4. It is visually shown in Fig. 8 if one ignores the theoretical uncertainties whose maxima are indicated by the coloured lines in the bottom left of the upper panel. A firm conclusion irrespective of the choice of merit function, even in the case where any theoretical uncertainties are ignored, is that the radiative grid delivers the better regression solution.

We conclude that if we allow for theoretical uncertainties in the regression problem at the level of the variance–covariance covered by the model grids, the prograde dipole modes of KIC 7760680 do not provide statistical evidence of the occurrence of CBM in this star, which has an exceptionally low mixing beyond the convective core. For the models without CBM, an increase in the uncertainties on the inferred model parameters occurs compared to the situation where more complex models with CBM are considered. This can be seen by the increased uncertainty of the solutions for the case of the decreasing number of free parameters in Fig. 9.

Fig. 9. Parameters of the selected best model with their uncertainty ranges indicated by the black lines. If no uncertainty was found through Bayes’ theorem, its upper limit is indicated by the dashed grey lines. The blue and red symbols indicate the results for the radiative and Péclet grid, respectively, while the labels on the y-axis indicate how the theoretical frequency pattern was constructed. The top row of figures represents the full grids, while the rows underneath represent the nested grids with fewer parameters.

6.5. Bootstrapping observed frequencies

In this paper we have been relying on the observational input from Pápics et al. (2015), but we wondered how the frequency resolution of the data may have impacted the selection of modes during the extraction and identification of period spacing patterns via iterative pre-whitening (see Sect. 2). This is particularly important for modes of rather low amplitude as choices in the method of iterative pre-whitening can impact the resultant period spacing pattern. To assess if and how our asteroseismic modelling might depend on the robustness of an extracted period spacing pattern, we considered the statistical resampling method of bootstrapping where we replaced one of the modes in the observed period spacing pattern by another value, relying on the frequency resolution of the light curve, R = 1/ΔT = 0.00068 d⁻¹.

We perturbed different observed frequencies within the pattern one at a time by increasing or decreasing its value by R. We performed this experiment for a mode at high period (at 1.45261 d), a trapped mode (at 1.24544 d), a non-trapped mode (at 1.117507 d), and for a mode at low period (at 0.90147 d) within the extracted pattern from Pápics et al. (2015). For each of these cases we then repeated the whole analysis for these patterns based on the one perturbed mode frequency instead of the one from Pápics et al. (2015). This gives us an idea of how sensitive the asteroseismic modelling solutions are to an individual frequency being different in the pattern. Such a sensitivity analysis is standard in modern statistics and is of particular relevance in the case of a trapped gravito-inertial mode leading to a dip in the period spacing pattern.

The results of remodelling all these bootstrapped patterns are shown in Fig. 10. The best solutions sometimes deviate from the one found for the unperturbed pattern, but for the cases where this occurs the uncertainty range almost always includes the unperturbed solution. Our modelling results hence reveal to be mostly robust against the replacement of one mode in the pattern according to the frequency resolution. We recall that this resolution is a gross overestimation of the actual frequency uncertainties, but it must be kept in mind in the context of a pre-whitening procedure (Pápics et al. 2014).

Fig. 10. Same as Fig. 9, but the different lines per extraction method are for each of the different perturbed patterns from the bootstrapping; the upper lines are for the unperturbed pattern.

7. Discussion and conclusion

In this work we investigated the relative importance of various aspects when performing forward asteroseismic modelling of gravito-inertial modes. It is important to not only consider which physical prescriptions or free parameters are used in the computational setup, but to also take different methods for the observational input of the regression problem into account. We highlighted the impact of choices in the procedures via the case study of KIC 7760680, which is the B-type star with the largest number of identified prograde dipole gravito-intertial modes and the best assessed near-core rotation rate among all the published main-sequence gravito-inertial mode pulsators to date (Pápics et al. 2015; Moravveji et al. 2016).

We focused our study on the question of whether it is possible to infer the most likely temperature gradient in the near-core boundary layer from identified high-order gravito-inertial modes. This was done by considering two model grids based on two different prescriptions for this temperature gradient, following the potential capacity of dipole g-mode asteroseismology published in the twin papers by Pedersen et al. (2018) and Michielsen et al. (2019). While doing so, we considered two different merit functions to solve the regression problem, following Aerts et al. (2018). One of these is χ² while the other is the Mahalanobis distance taking into account the variance–covariance properties of the theoretical predictions, including the entire correlation structure of the regression problem at hand. Moreover, we assessed the quality of different sets of observables to perform the regression based on the comparison between observed and theoretically predicted pulsation patterns.

In terms of a comparison between the two merit functions, we found that using the Mahalanobis distance not only changes the preferred model parameters as compared to χ², but it constrains the optimal solution better. Irrespective of the merit function, we found a clear preference for the model grid based on the radiative temperature gradient in the small core-boundary layer for KIC 7760680 in models that include this layer, instead of a gradually changing temperature gradient between the adiabatic and radiative gradient based on the Péclet number. We emphasise that this preference for a radiative temperature gradient in the small core-boundary layer has been obtained for just this one rotating B-type pulsator, and it should not be considered a general result for all B stars given the immense diversity in observed period spacing patterns and the large range of CBM and envelope mixing levels for 26 such stars found by Pedersen et al. (2021), as summarized in Aerts (2021, Table 1). The summary of all available g-mode asteroseismology of B stars has so far revealed that KIC 7760680 has exceptionally low envelope mixing among those stars rotating at a considerable fraction of the critical rotation rate. We also found that models without any convective penetration or overshooting are statistically preferred for this rotating star, as shown by our evaluation of nested regression models.

For the regression problem based on the Mahalanobis distance, the use of different sets of quantities to be compared with the respective observables can heavily influence the condition number of the used variance–covariance matrix. In particular, the way the theoretical pulsation quantities are used deserves attention when adopting the Mahalanobis distance as merit function since it impacts the structure of the variance–covariance matrix to be inverted in the regression. The effects for the regression based on a χ² merit function are smaller since in that case we ignore the covariance among the theoretical prediction or even the theoretical uncertainty as a whole. We should find a compromise between the best numerical stability and the use of maximum observable information in the definition of the regression problem. The application of singular value decomposition or principle component analysis as adopted for some of the modelling efforts of low-mass stars (e.g. Angelou et al. 2017) or intermediate-mass pulsators (Mombarg et al. 2019) might be a useful additional or complementary approach to the one we adopted here for g-mode asteroseismology of B stars.

The best model parameters we found for the cases with and without CBM are provided in Tables C.1 and C.4 and are listed in Table 3. Overall, the model parameters obtained by Moravveji et al. (2016) listed in Table 1 are outside of the uncertainties on the model with CBM, but within those of the statistically preferred model without CBM. The well-known mass–metallicity and mass–CBM degeneracies for pulsators in this mass range (Moravveji et al. 2015, their Fig. 5) are the basis of the uncertainty regions of the estimated parameters, whose projection in one dimension are shown graphically in Fig. 9. The solved regression problem based on the Mahalanobis distance merit function properly takes into account these correlations and places the solution based on the χ² as obtained by Moravveji et al. (2016) within the uncertainty regions of the estimated parameters. This is remarkably consistent given that we used a non-flat D_mix(r) profile and performed the fitting from the ΔP values instead of the mode frequencies as used by Moravveji et al. (2016). Our less refined grid step implied a major reduction in the analysis time of the problem and allowed us to treat the primary goal set for this paper, which was to investigate if we can distinguish between a different thermal and chemical structure in the CBM region while having fixed the shape of the mixing profile as diffusive core overshooting in the core-boundary layer, and as being due to wave mixing in the envelope. With this choice of core overshooting, which is often taken in the literature, we found that a radiative temperature gradient in the CBM region gives a better fit to the data than the models with a temperature gradient based on the Péclet number. However, the presence of such a CBM region altogether turned out to be statistically disfavoured by the AICc for this particular star rotating at about a quarter of its critical rate.

Finally, it is noteworthy that the best model with CBM listed in Table 3 falls within the 3σ error box of the spectroscopic T_eff and astrometrically derived log(L/L_⊙) as shown in Fig. 3, while these observables were not used in the asteroseismic modelling process. This a posteriori evaluation of the asteroseismic solution is reassuring as there was no need to reject models based on these non-asteroseismic observables.

Overall, we conclude that choices in the modelling methodology should be made with great care and consideration. The optimal modelling setup might vary from star to star because some pulsators reveal signatures of trapped modes, while others do not (Pápics et al. 2017). A fully automated modelling strategy may therefore be premature for ensemble modelling of gravito-inertial B-type pulsators covering the entire instability strip of such stars (see Pedersen et al. 2020), particularly as our current case study was done for a star with very high-precision observables and an exceptionally low level of mixing in the core-boundary layer and stellar envelope compared to most other gravito-inertial pulsators of this kind among the B-type stars (Pedersen et al. 2021).

Acknowledgments

The authors are grateful to the MESA and GYRE developers teams led by Bill Paxton and Rich Townsend for their efforts and for releasing their software publicly; this study would not have been possible without their codes. The authors thank Geert Molenberghs for his helpful advice on the statistical aspects of the study, Joey Mombarg for the fruitful discussions, and Cole Johnston for his insights in theoretical pulsation pattern construction. We thank the referee for the comments and suggestions towards improving the manuscript. 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 and a senior postdoctoral fellowship to DMB with grant agreement No. 1286521N, from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 670519: MAMSIE to CA), and from the KU Leuven Research Council (grant C16/18/005: PARADISE).

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¹.

Appendix B: Reduced χ² as merit function

B.1. Model selection criterion

The most commonly used merit function to assess the goodness of fit in forward asteroseismic modelling of SPB pulsators is a χ², which is based on a Euclidian distance. Denoting a set of uncorrelated observables as Y^Obs with errors according to a normal distribution with variance ${\sigma }_{{\mathrm{Y}}^{\mathrm{Obs}}}^2$, and the set of corresponding values predicted by a theoretical model as Y^Theo, we write the reduced χ² score as

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red}}^2=\frac{1}{N-k}{\displaystyle \sum _i^N{\left(\frac{{\mathrm{Y}}_i^{\mathrm{Theo}}-{\mathrm{Y}}_i^{\mathrm{Obs}}}{{\sigma }_{{\mathrm{Y}}_i^{\mathrm{Obs}}}}\right)}^2,}\end{array}$$(B.1)

where N is the number of observables and k is the number of free parameters. In general, our MD formulation used in the main part of the paper can be connected to the ${\chi }_{\text{red}}^2$ merit function by ignoring uncertainties on the theoretical predictions of the observables and by ignoring correlations among the observables. Ignoring theoretical uncertainties comes down to ignoring the covariance structures shown in Appendix D. We provide the ${\chi }_{\text{red}}^2$ results here to assess the differences in the solutions to the regression problem by making these simplifications, following the earlier assessment of theoretical uncertainties in Aerts et al. (2018). The modelling procedures applied to the g-mode frequencies of rotating SPB stars in the literature so far consider the mode frequencies as independent observables (Moravveji et al. 2016; Szewczuk & Daszyńska-Daszkiewicz 2018), such that using a ${\chi }_{\text{red}}^2$ is meaningful. Nevertheless, we show the results obtained by using this simplified merit function for two cases, of mode periods and of mode period spacings, where the latter in principle demands using a formalism suitable for correlated data such as the MD. This allows the reader to assess the impact induced by adopting a simplified ${\chi }_{\text{red}}^2$ for applications to correlated data, while also ignoring theoretical uncertainties.

Using ${\chi }_{\text{red}}^2$ as the merit function, the AICc reduces to

$$\begin{array}{c}\hfill \mathrm{AICc}={\chi }_{\mathrm{red}}^2+\frac{2kN}{N-k-1}·\end{array}$$(B.2)

B.2. Observables to fit

When considering χ² as the merit function, there is no variance–covariance matrix whose numerical stability needs to be taken into account through condition numbers. We simply compare the use of different sets of observables and look at the correlation plots and modelling results to assess their use.

The correlation plots show a better constrained region with optimal solutions when selecting the period spacings as observables instead of the periods themselves. This can be seen comparing Figs. B.1 and B.2 for the radiative grid, or Figs. B.3 and B.4 for the Péclet grid.

Fig. B.1. Correlation plot for the radiative grid, using periods in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method. The 50% best models are shown, colour-coded according to the log of their merit function value (at right). The figures on the diagonal show binned parameter distributions, and the larger panel at the top right shows an Hertzsprung–Russell diagram with the 1 and 3σ spectroscopic Teff and astrometric luminosity error boxes.

Fig. B.2. Correlation plot (see Fig. B.1) for the radiative grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.

Fig. B.3. Correlation plot (see Fig. B.1) for the Péclet grid, using periods in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.

Fig. B.4. Correlation plot (see Fig. B.1) for the Péclet grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.

Apart from the correlation structure, the actual χ² values themselves are also lower when fitting period spacings. This is no surprise since the mode periods have much smaller relative errors, of order 0.003% or even smaller in the case of KIC 7760680 (Pápics et al. 2015), than the propagated errors on ΔP, which are about 10 to 20 times larger.

B.3. Theoretical mode period pattern construction

Although the different methods used to construct the theoretical frequency pattern can influence the distribution of the χ² values (Fig. B.5), the best model in the grid and its ${\chi }_{\text{red}}^2$ value remain identical when considering the full grids of models (Table C.5). The optimal solution slightly changes for some of the nested statistical models (Tables C.6–C.8), but this change in most cases only amounts to one grid step-size for one parameter. Additionally, the differences in the correlation plots are minimal, and the same global morphology applies (Figs. B.2, B.6, and B.7).

Fig. B.5. Distribution of the reduced χ2, using ΔP as observables to fit. The inset shows the cumulative distribution function of the 20 best models in the grids.

Fig. B.6. Correlation plot (see Fig. B.1) for the radiative grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the highest amplitude method.

Fig. B.7. Correlation plot (see Fig. B.1) for the radiative grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the highest frequency method.

B.4. Nested statistical regression models

Investigating the best regression models with varying complexity in the CBM of the models reveals a different trend than for the Mahalanobis distance. Comparing the best models of the full grids with six free parameters (Table C.5) to those of the partial grids with five (Tables C.6 and C.7) or four free parameters (Table C.8), reveals that the more complex models are generally preferred over the simpler ones. In contrast to when using the Mahalanobis distance, the increased fit quality plays a dominant role over the penalty for the increased model complexity. The main reason for this is that no uncertainty is propagated for the theoretical predictions when using a χ² merit function.

When comparing the different partial grids we find a clear preference for the models with an exponentially decaying CBM coefficient. This is in agreement with the results of Moravveji et al. (2016). Using a χ² merit function the authors found that an exponentially decaying CBM prescription outperformed a step-like mixing coefficient in the CBM region. We thus recover the previous result under the restricted assumption that the theoretical model predictions are error-free.

B.5. Radiative grid versus Péclet grid

According to the ${\chi }_{\text{red}}^2$ values the radiative grid gives a better fit to the observations than the Péclet grid. Although there is a clear difference between the best models of the two grids (Table C.5), the distribution of their ${\chi }_{\text{red}}^2$ values does not differ too much (Fig. B.5). However, a different conclusion arises when using periods as the observables without considering the period spacings. Specifically, the ${\chi }_{\text{red}}^2$ distributions are again similar between the two grids, but the preferred best model stems from the Péclet grid (Fig. B.8 and Table C.5).

Fig. B.8. Distribution of the reduced χ2, using periods as observables to fit. The inset shows the cumulative distribution function of the 20 best models in the grids.

We note that the best solutions according to the Mahalanobis distance are not among the best solutions delivered by the ${\chi }_{\text{red}}^2$ metric (Table C.1). On the other hand, the best models according to the ${\chi }_{\text{red}}^2$ still score among the better models following the Mahalanobis distance (listed in Table C.5), although they are rated significantly worse than the best ones for that metric.

Even though the conclusions of the analysis of this star with the ${\chi }_{\text{red}}^2$ merit function would not have been influenced by the way we construct the theoretical pulsation pattern, the choice of observables (mode periods versus period spacings) changes the answer to the question of which temperature gradient is preferred when only looking at the best point estimator, without taking into account any precision estimation.

Appendix C: Best models

Tables listing the parameters of the best fit models.

Appendix D: Variance–covariance matrices

Figures showing the variance–covariance matrices for the two grids, the different pattern construction methods, and the different sets of asteroseismic observables.

Fig. D.1. Variance–covariance matrix for the periods of the radiative grid, selected using the longest sequence method.

Fig. D.2. Variance–covariance matrix for the periods of the radiative grid, selected using the highest amplitude method.

Fig. D.3. Variance–covariance matrix for the periods of the radiative grid, selected using the highest frequency method.

Fig. D.4. Variance–covariance matrix for the periods of the Péclet grid, selected using the longest sequence method.

Fig. D.5. Variance–covariance matrix for the periods of the Péclet grid, selected using the highest amplitude method.

Fig. D.6. Variance–covariance matrix for the periods of the Péclet grid, selected using the highest frequency method.

Fig. D.7. Variance–covariance matrix for the period spacings of the radiative grid, selected using the longest sequence method.

Fig. D.8. Variance–covariance matrix for the period spacings of the radiative grid, selected using the highest amplitude method.

Fig. D.9. Variance–covariance matrix for the period spacings of the radiative grid, selected using the highest frequency method.

Fig. D.10. Variance–covariance matrix for the period spacings of the Péclet grid, selected using the longest sequence method.

Fig. D.11. Variance–covariance matrix for the period spacings of the Péclet grid, selected using the highest amplitude method.

Fig. D.12. Variance–covariance matrix for the period spacings of the Péclet grid, selected using the highest frequency method.

All Tables

All Figures

Fig. 1. Radial profiles of a 3 M⊙ star halfway through the main sequence with a central hydrogen content Xc = 0.4. Top panels: temperature gradients. Middle panels: structure (solid) and composition (dashed) components of the Brunt−Väisälä profile, 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). Bottom panels: mode inertia of two g-modes with different radial orders. (a) Near-core region with radiative temperature gradient. (b) Near-core region with temperature gradient according to Péclet number.

In the text

Fig. 2. Correlation plot for the radiative grid, using periods in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the longest sequence method. The 50% best models are shown, colour-coded according to the log of their merit function value (at right). The figures on the diagonal show binned parameter distributions, and the larger panel at the top right shows a Hertzsprung–Russell diagram with the 1 and 3σ spectroscopic Teff and astrometric luminosity error boxes.

In the text

	Fig. 3. Correlation plot (see Fig. 2) for the radiative grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.
In the text

	Fig. 4. Correlation plot (see Fig. 2) for the Péclet grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.
In the text

	Fig. 5. Correlation plot (see Fig. 2) for the Péclet grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the highest amplitude method.
In the text

	Fig. 6. Correlation plot (see Fig. 2) for the Péclet grid, using period spacings in a Mahalanobis distance merit function. The theoretical pulsation pattern was constructed according to the highest frequency method.
In the text

	Fig. 7. Distribution of the Mahalanobis distances using ΔP as observables to fit. The inset shows the cumulative distribution function of the 20 best models in the grids.
In the text

Fig. 8. Period spacing patterns of the observations, of the best models of the two grids, and of the nested grid where no CBM is included in the models, as extracted by the longest sequence method. The formal errors on the observations are smaller than the symbol sizes. The vertical bars in the bottom left corner show the maximum considered uncertainty for the theoretical predictions approximated by the variance–covariance matrix of that particular grid. Middle and bottom panels: 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 Pápics et al. (2015), including the correction factor of four.

In the text

Fig. 9. Parameters of the selected best model with their uncertainty ranges indicated by the black lines. If no uncertainty was found through Bayes’ theorem, its upper limit is indicated by the dashed grey lines. The blue and red symbols indicate the results for the radiative and Péclet grid, respectively, while the labels on the y-axis indicate how the theoretical frequency pattern was constructed. The top row of figures represents the full grids, while the rows underneath represent the nested grids with fewer parameters.

In the text

	Fig. 10. Same as Fig. 9, but the different lines per extraction method are for each of the different perturbed patterns from the bootstrapping; the upper lines are for the unperturbed pattern.
In the text

Fig. B.1. Correlation plot for the radiative grid, using periods in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method. The 50% best models are shown, colour-coded according to the log of their merit function value (at right). The figures on the diagonal show binned parameter distributions, and the larger panel at the top right shows an Hertzsprung–Russell diagram with the 1 and 3σ spectroscopic Teff and astrometric luminosity error boxes.

In the text

	Fig. B.2. Correlation plot (see Fig. B.1) for the radiative grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.
In the text

	Fig. B.3. Correlation plot (see Fig. B.1) for the Péclet grid, using periods in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.
In the text

	Fig. B.4. Correlation plot (see Fig. B.1) for the Péclet grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the longest sequence method.
In the text

	Fig. B.5. Distribution of the reduced χ2, using ΔP as observables to fit. The inset shows the cumulative distribution function of the 20 best models in the grids.
In the text

	Fig. B.6. Correlation plot (see Fig. B.1) for the radiative grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the highest amplitude method.
In the text

	Fig. B.7. Correlation plot (see Fig. B.1) for the radiative grid, using period spacings in a χ2 merit function. The theoretical pulsation pattern was constructed according to the highest frequency method.
In the text

	Fig. B.8. Distribution of the reduced χ2, using periods as observables to fit. The inset shows the cumulative distribution function of the 20 best models in the grids.
In the text

	Fig. D.1. Variance–covariance matrix for the periods of the radiative grid, selected using the longest sequence method.
In the text

	Fig. D.2. Variance–covariance matrix for the periods of the radiative grid, selected using the highest amplitude method.
In the text

	Fig. D.3. Variance–covariance matrix for the periods of the radiative grid, selected using the highest frequency method.
In the text

	Fig. D.4. Variance–covariance matrix for the periods of the Péclet grid, selected using the longest sequence method.
In the text

	Fig. D.5. Variance–covariance matrix for the periods of the Péclet grid, selected using the highest amplitude method.
In the text

	Fig. D.6. Variance–covariance matrix for the periods of the Péclet grid, selected using the highest frequency method.
In the text

	Fig. D.7. Variance–covariance matrix for the period spacings of the radiative grid, selected using the longest sequence method.
In the text

	Fig. D.8. Variance–covariance matrix for the period spacings of the radiative grid, selected using the highest amplitude method.
In the text

	Fig. D.9. Variance–covariance matrix for the period spacings of the radiative grid, selected using the highest frequency method.
In the text

	Fig. D.10. Variance–covariance matrix for the period spacings of the Péclet grid, selected using the longest sequence method.
In the text

	Fig. D.11. Variance–covariance matrix for the period spacings of the Péclet grid, selected using the highest amplitude method.
In the text

	Fig. D.12. Variance–covariance matrix for the period spacings of the Péclet grid, selected using the highest frequency method.
In the text