# Forward Asteroseismic Modeling of Stars with a Convective Core from Gravity-mode Oscillations: Parameter Estimation and Stellar Model Selection

C. Aerts, G. Molenberghs, M. Michielsen, M. Pedersen, R. Björklund, Calum Johnston, J. Mombarg, D. Bowman, B. Buysschaert, P. Pápics, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

C. Aerts, G. Molenberghs, M. Michielsen, M. Pedersen, R. Björklund, et al., Forward Asteroseismic Modeling of Stars with a Convective Core from Gravity-mode Oscillations: Parameter Estimation and Stellar Model Selection. The Astrophysical Journal Supplement, 2018, 237 (1), pp.15. 10.3847/1538-4365/aaccfb. hal-02297145

HAL Id: hal-02297145

https://hal.science/hal-02297145v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Forward Asteroseismic Modeling of Stars with a Convective Core from Gravity-mode Oscillations: Parameter Estimation and Stellar Model Selection

C. Aerts<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, G. Molenberghs<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, M. Michielsen<fnref n="1" />, M. G. Pedersen<fnref n="1" />, R. Björklund<fnref n="1" />, C. Johnston<fnref n="1" />, J. S. G. Mombarg<fnref n="1" />, D. M. Bowman<fnref n="1" />, B. Buysschaert<fnref n="1" />,<fnref n="6" />, P. I. Pápics<fnref n="1" /> ID, S. Sekaran<fnref n="1" />, J. O. Sundqvist<fnref n="1" />,

A. Tkachenko<fnref n="1" />, K. Truyaert<fnref n="1" />, T. Van Reeth<fnref n="1" />, and E. Vermeyen<fnref n="1" />

Received 2018 April 13; revised 2018 May 31; accepted 2018 June 13; published 2018 July 19

## Abstract

We propose a methodological framework to perform forward asteroseismic modeling of stars with a convective core, based on gravity-mode oscillations. These probe the near-core region in the deep stellar interior. The modeling relies on a set of observed high-precision oscillation frequencies of low-degree coherent gravity modes with long lifetimes and their observational uncertainties. Identification of the mode degree and azimuthal order is assumed to be achieved from rotational splitting and/or from period spacing patterns. This paper has two major outcomes. The first is a comprehensive list and discussion of the major uncertainties of theoretically predicted gravity-mode oscillation frequencies based on linear pulsation theory, caused by fixing choices of the input physics for evolutionary models. Guided by a hierarchy among these uncertainties of theoretical frequencies, we subsequently provide a global methodological scheme to achieve forward asteroseismic modeling. We properly take into account correlations among the free parameters included in stellar models. Aside from the stellar mass, metallicity, and age, the major parameters to be estimated are the near-core rotation rate, the amount of convective core overshooting, and the level of chemical mixing in the radiative zones. This modeling scheme allows for maximum likelihood estimation of the stellar parameters for fixed input physics of the equilibrium models, followed by stellar model selection considering various choices of the input physics. Our approach uses the Mahalanobis distance instead of the often-used $\chi^2$ statistic and includes heteroscedasticity. It provides estimation of the unknown variance of the theoretically predicted oscillation frequencies.

**Key words:** asteroseismology – methods: statistical – stars: interiors – stars: massive – stars: oscillations (including pulsations) – stars: rotation

## 1. Introduction

Asteroseismology of low-mass stars with stochastically excited pressure-mode oscillations often relies on simple scaling of solar oscillation frequencies. The so-called scaling relations are easy to apply and provide estimates of the stellar mass, radius, and age, assuming that the stars under consideration behave similarly to the Sun (e.g., Chaplin & Miglio 2013; Hekker & Christensen-Dalsgaard 2017, for reviews). If the aim is to improve the input physics of stellar models, however, one must go beyond application of scaling relations. This requires testing various assumptions and choices for the input physics that enters the computation of the stellar equilibrium models on which one relies for the estimation of theoretical oscillation frequencies. In practice, estimation of various parameters that are left free in the computation of equilibrium models is required, along with proper error assessment. Subsequently, statistical model selection must be applied after the parameter estimation has been achieved, in order to evaluate the capacity of the input physics in explaining the asteroseismic data.

Approaches to achieve parameter estimation followed by evaluation of the input physics from stellar model selection in asteroseismology have been developed and improved, particularly since space photometry became available. An

enlightening astrostatistics tutorial of various statistical methods for asteroseismology is available in Appourchaux (2014). Usually, parameter estimation is done from forward modeling by adopting a grid-based approach and considering millions of stellar structure models of different evolutionary stage. Early applications to low-mass stars with solar-like oscillations were made by, e.g., Miglio & Montalbán (2005), Quirion et al. (2010), and to subdwarfs by Brassard et al. (2001). In those studies, a $\chi^2$ comparison between observed and theoretically predicted frequencies was adopted to achieve the parameter estimation, while selection of the best model input physics among various choices has been limited so far. Our focus is on the latter aspect. Applications in the space asteroseismology era have so far also been focused on low-mass stars with solar-like oscillations (e.g., Gruberbauer et al. 2013; Appourchaux et al. 2014; Bellinger et al. 2016; Deheuvels et al. 2016; Silva Aguirre et al. 2017, to list a few) and on white dwarfs (e.g., Giammichele et al. 2017). Gruberbauer et al. (2012) developed a Bayesian framework that includes systematic uncertainties due to surface effects that occur for solar-like oscillations (e.g., Ball & Gizon 2017; Trampedach et al. 2017). The cause of this effect is unknown, yet it is dominant over other choices of the input physics when fitting oscillation frequencies and must therefore be treated properly. These unknown surface effects

<footnote n="1">¹Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven, Belgium; conny.aerts@kuleuven.be</footnote>

<footnote n="2">²Department of Astrophysics, IMAPP, Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands</footnote>

<footnote n="3">³Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA</footnote>

<footnote n="4">⁴I-BioStat, Universiteit Hasselt, Martelarenlaan 42, B-3500 Hasselt, Belgium</footnote>

<footnote n="5">⁵I-BioStat, KU Leuven, Kapucijnenvoer 35, B-3000 Leuven, Belgium</footnote>

<footnote n="6">⁶LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195, Meudon, France</footnote>
can be avoided by fitting frequency separations and ratios thereof, rather than the frequency values themselves. This has been done in most applications of solar-like oscillations so far.

The potential of the detected oscillations in terms of probing power is totally different for gravity modes than for pressure modes. The latter mainly probe envelope physics, while gravity modes dominantly probe near-core physics, and mixed modes in evolved stars have combined capacity. The exploitation of gravity-mode oscillations in stars born with a convective core and a radiative envelope (i.e., spectral type from early F to O) faces different challenges but offers new opportunities compared to solar-like oscillations. Gravity modes are strongly affected and thus optimally suited to probe three phenomena in the deep stellar interior: convective core overshooting, near-core rotation, and chemical mixing. The aim of the current paper is twofold: (i) to assess the major uncertainties in predicted theoretical frequencies for the case of gravity modes with long lifetimes, leading to the requirement of heteroscedasticity in the statistical treatment; and (ii) to provide a methodological framework for forward modeling, taking into account that gravity-mode oscillations depend nonlinearly on the free parameters of stellar equilibrium models (e.g., Townsend & Teitler 2013) and that these parameters are correlated (e.g., Figure 11 in Pápics et al. 2014).

While thousands of low-mass stars were observed with μmag level precision by *Kepler* (see, e.g., Chaplin & Miglio 2013; Hekker & Christensen-Dalsgaard 2017, for reviews), the satellite observed far fewer intermediate- and high-mass stars with gravity-mode oscillations, because the exoplanet hunting from planetary transits required observing a large sample of cool low-mass stars as potential hosts. For this reason, few stars with mass above $2 M_\odot$ in the field of view of the nominal mission were monitored. Moreover, despite notable new detections of nonradial oscillations in tens of such stars made with the *MOST* (e.g., Walker et al. 2005; Aerts et al. 2006; Saio et al. 2006; Cameron et al. 2008; Zwintz et al. 2013) and *CoRoT* (e.g., Neiner et al. 2009; Poretti et al. 2009; Degroote et al. 2010; Briquet et al. 2011; Pápics et al. 2012; Zwintz et al. 2014) satellites, the time-series photometry assembled by these two missions was limited from weeks to a few months. While this is sufficient to detect gravity modes, it is often insufficient to identify them in terms of their spherical harmonic mode wavenumbers (*l*, *m*). Such mode identification is a necessary prerequisite to model the interior physical properties of stars and to improve the input physics of theoretical models.

The bias toward asteroseismology of low-mass stars is also related to the timescales of the oscillations. Solar-like oscillations of low-mass stars are excited stochastically in their outer convective envelope and have periodicities ranging from minutes for the core hydrogen burning phase to hours for red giants. This stands in sharp contrast to the coherent gravity modes with periods up to several days, leading to complex beating patterns among the modes that may cover several years (e.g., Kurtz et al. 2014; Saio et al. 2015; Van Reeth et al. 2015; Bowman et al. 2016; Pápics et al. 2017). Such modes are self-driven by a heat- or flux-blocking mechanism (e.g., Aerts et al. 2010, Chapters 2 and 3) and have lifetimes of tens to millions of years, i.e., long compared with the duration of the data sets (see Figure 3 discussed below). Such gravity modes occur in stars with masses from roughly 1.4 to $40 M_\odot$. Our long-term aim to improve the evolutionary models of such

intermediate- and high-mass stars using asteroseismology is motivated by the fact that they are the dominant suppliers of the chemicals in galaxies, while their evolution models are far more uncertain than those of low-mass stars.

The evolution of stars born with a well-developed convective core is appreciably affected by the size and properties of their convective core region, the interior rotation, and the transport of chemicals and angular momentum in the radiative envelope. Asteroseismic (hereafter referred to as seismic) modeling of such stars had to await proper and unique mode identification for a sufficient number of detected nonradial oscillations. Early achievements were focused on the identification of a few low-order pressure modes detected in extensive ground-based data of a few bright A- and B-type stars (e.g., Breger et al. 1999; Aerts et al. 2003; De Ridder et al. 2004; Handler et al. 2006; Handler 2009). At best, these studies provided a rough estimate of the core overshooting and core-to-envelope rotation because the pressure modes have limited probing power near the stellar core (e.g., Dupret et al. 2004; Dziembowski & Pamyatnykh 2008; Handler et al. 2009; Briquet et al. 2012). Since the 4 yr nominal *Kepler* mission, seismic modeling based on numerous identified gravity-mode oscillations can be achieved with sufficiently high precision to improve the input physics of stellar evolution models of stars with $M \gtrsim 1.4 M_\odot$. Meanwhile, several stars have been modeled (Kurtz et al. 2014; Moravveji et al. 2015, 2016; Saio et al. 2015; Schmid & Aerts 2016; Van Reeth et al. 2016; Kallinger et al. 2017; Sowicka et al. 2017; Szewczuk & Daszyńska-Daszkiewicz 2018), but at a very different level of depth, relying on a variety of fixed choices for the input physics, and making diverse assumptions on the importance of the free parameters of the stellar models. Now that suitable data for such inference studies are available, we need an appropriate statistical methodology to perform seismic modeling, which goes beyond the methods currently used. This is the topic of the current paper. We focus on core hydrogen burning stars for which the *Kepler* mission provided suitable data, but our methods are readily applicable to stars with a convective core in more advanced stages of stellar evolution.

## 2. The Physical Ingredients for Stellar Oscillation Computations

### 2.1. Preliminaries

The seismic modeling presented here involves comparing the values of observed gravity-mode oscillation frequencies with those predicted from theoretical stellar structure models based on chosen input physics for a given set of stellar parameters. We opt to focus the application on gravity modes because they have the highest mode energy in the near-core regions that we wish to probe (e.g., Figure 2 in Triana et al. 2015; Van Reeth et al. 2016). These modes are far less sensitive to poor descriptions of the physical conditions in the outer envelope of a star, so we are not in need of corrections for surface effects as in the case of high-order pressure modes (e.g., Ball & Gizon 2017; Trampedach et al. 2017).

We take an *observationally driven* approach in the sense that we use all detected oscillation frequencies, irrespective of the excitation mechanism, as long as the mode frequency has unambiguous mode identification in terms of its spherical harmonic wavenumbers *l* and *m*. Together with the mode frequency and the radial and horizontal amplitude, wavenumbers define the displacement vector due to the mode (see
Equation (3.132) in Aerts et al. [2010], and Section 3.3 in that monograph for a general description of linear stellar oscillations). Although we have a measurement of the surface-integrated amplitude of the modes in the line of sight, theories to predict the intrinsic amplitudes for coherent modes remain challenging. While amplitude predictions for radial pulsators based on first principles are available to some extent (e.g., Smolec & Moskalik 2007; Geroux & Deupree 2013), amplitude-saturation mechanisms based on nonlinear mode coupling of nonradial modes have only been developed for pressure modes in main-sequence A-type stars (e.g., Dziembowski & Krolikowska 1985; Buchler et al. 1997; Nowakowski 2005) and not for gravity modes in rotating stars. The predictive power of these available nonlinear mode-coupling computations is insufficient to rely on them in forward modeling. Moreover, although generally successful, details of the excitation mechanisms are insufficiently understood to rely on them for seismic modeling because modes predicted to be excited are often not observed or vice versa. Hence, a secure way to proceed in evaluating the input physics of stellar structure models from detected oscillations is to rely solely on their measured frequency values and visit the mode excitation problem after the forward modeling has delivered the most likely models. The current paper focuses only on the modeling of detected frequencies and parameter estimation; for extensive excitation computations of intermediate- and high-mass pulsators, we refer to the recent paper by Szewczuk & Daszyńska-Daszkiewicz (2017).

As further discussed below, coherent gravity modes have long lifetimes compared to the duration of the data sets. It is hence fully justified to assume that they cause a delta-function convolved with the spectral window function in a Fourier transform and the observational frequency uncertainties are inversely proportional to the total time base of the data. Here, we adopt the notation that $f_i^*$ is the observed cyclic oscillation frequency of mode *i*, which has known error $\epsilon_i^*$, for *i* = 1, ..., *n*.

The typical uncertainty of measured oscillation frequencies derived from data with high duty cycle that does not suffer from aliasing, as is the case for contemporary high-cadence uninterrupted space photometry, is well below the resolving power of the data set. This resolving power is about 2.5 times the Rayleigh limit of the data set (Loumos & Deeming 1978), the latter being equal to the inverse of the total time base of the data. The Rayleigh limit is $0.00068 \text{ day}^{-1}$ (0.0079 µHz) for the 4 yr nominal *Kepler* data (e.g., Bowman 2017) and $0.0066 \text{ day}^{-1}$ (0.077 µHz) for a *CoRoT* 150-day-long run (e.g., Degroote et al. 2009). Recent BRITE photometry has time spans between roughly 100 and 180 days and leads to frequency resolutions similar to *CoRoT* long runs for the brightest stars in the sky (Pablo et al. 2016).

Each oscillation frequency has its own measurement error $\epsilon_i^*$, because it is dependent on the measured mode amplitude and on the noise properties of the data in that particular frequency regime, aside from the total time base (see, e.g., Equation (5.52) in Aerts et al. 2010). For all significant frequencies detected from well-sampled uncorrelated data that do not suffer from systematic uncertainty caused by aliasing, the frequency error is at least a factor of 10 below the Rayleigh limit. An even larger factor applies to oversampled space photometry, but one has to correct for the correlated nature of the data (Degroote et al. 2009). Methods of frequency analysis and prewhitening

applied to heat-driven pulsators are often based on the assumptions of uncorrelated, homoscedastic white Gaussian noise. Degroote et al. (2009) showed that these assumptions are not strictly true for the *CoRoT* asteroseismology data (32 s sampling), but that the deviations thereof are often sufficiently small to still apply the methods, provided that a correction for correlated data is used. However, the prewhitening process inherently implies introduction of uncertainties connected with the limited resolving power of the data set, even for the 4 yr long *Kepler* data, and may result in dependencies among the frequencies and their errors. For these reasons, a good approach is to allow for heteroscedasticity in the measurement errors whenever frequencies deduced from a prewhitening procedure are used in forward modeling.

For *CoRoT/BRITE* and *Kepler* data of core hydrogen burning gravity-mode pulsators, the uncertainties of the measured frequencies that satisfy the commonly used significance criterion of having an amplitude above four times the noise level (Breger et al. 1993) are below 0.001 and $0.0001 \text{ day}^{-1}$, respectively. Any comparison with theoretically predicted frequencies based on stellar models hence requires the frequencies to be computed at such levels. For all the above reasons, we consider $0.001 \text{ day}^{-1}$ to be a typical frequency error derived from a quasi-uninterrupted 0.5 to 1 yr light curve, and we allow for heteroscedasticity in the statistical methodology.

## 2.2. Stellar Structure Models

In order to compute theoretical cyclic frequencies of oscillations, denoted here as $f_i^{\text{th}}$, we rely on stellar structure models in hydrostatic equilibrium and perturb them in a linear approach. In this work, we use spherically symmetric 1D stellar models for single stars, computed with the publicly available MESA code; we refer to the extensive papers by Paxton et al. (2011, 2013, 2015, 2018) for a full description of the suite of routines, the available choices of the input physics, and the list of the numerous free parameters that one can use, as well as the code of conduct adopted by the MESA development team. Of course, our methodology can also be applied to any other stellar evolution code than the one we use here to illustrate it. As a word of caution, we stress that any user must check carefully that the stellar structure models computed numerically comply to a high level of precision with the equations of stellar structure on which they are based in the first place. Performing checks (e.g., on choosing mesh points and tolerances for convergence) to make sure that the models are astrophysically appropriate and that the oscillations are computed with sufficient numerical accuracy is the responsibility of the user (not of the code developers!). This is particularly relevant for applications to gravity modes, whose mode cavities are highly sensitive to discontinuities in the interior profiles of the physical quantities. Typically, gravity-mode predictions require some 5000–10,000 mesh points within the star to meet the observational frequency errors, not necessarily equidistantly distributed. Hence, gravity-mode computations need to be done with a refined mesh compared to computations that only require smooth evolutionary tracks in the Hertzsprung–Russell diagram or predictions of surface abundances (for which a total of ~1000 mesh points is usually sufficient). An appropriate mesh is particularly important in the transition regions between convective and radiative zones, to ensure that the Bruntt–Väisälä (BV) frequency and µ-gradients are well computed, as these determine the propagation cavity of the gravity modes. A
typical MESA inlist based on 5000 mesh points is provided through the links in Appendix A.

Less than 10% of all early F- to O-type stars have a surface magnetic field at the current observational threshold of spectropolarimetry (typically a few gauss; Wade et al. 2016). Moreover, it is currently unknown whether magnetic fields occur at the interface of the convective core and radiative envelope, let alone what their shape and strength are. Theoretical predictions, e.g., as in Braithwaite (2009), are too uncertain to be used as input physics. Rather, they can be evaluated with asteroseismology. In this way, it was recently found that the effect of the Lorentz force on oscillation frequencies is far below that of the Coriolis force in a magnetic gravity-mode pulsator (Buysschaert et al. 2018).

Ignoring rotational effects in equilibrium models is less justified than magnetism, because a large fraction of early F to O stars are fast rotators (Zorec & Royer 2012). Interior rotation or magnetism of stars introduces a multitude of instabilities and accompanying transport of elements and of angular momentum (see Heger et al. 2000, for a comprehensive description). Each of these phenomena has been implemented in MESA, with its own diffusive mixing and angular momentum transport coefficient as a free parameter. However, these ingredients are subject to considerable uncertainties because they remain essentially uncalibrated. The theoretical descriptions of these instabilities give rise to steep and discontinuous local fluctuations in the chemical mixing profile throughout the radiative envelope. While the consequences of these fluctuations remain invisible in evolutionary tracks, they affect the gravity modes appreciably. Because their nature may be numerical rather than physical (e.g., Figures 28 and 29 in Paxton et al. 2013), we performed numerous tests on the effect of each of these instabilities separately, as well as jointly, on gravity-mode predictions (Truyaert 2016). From these tests, we deduce that the current prescriptions of the instabilities induced by rotation or magnetism are not of practical use for forward seismic modeling of gravity modes. Either the discontinuities induced by the instabilities are physical and we conclude from the gravity-mode computations that real stars with such detected oscillations do not behave according to the theoretical predictions, or they are numerical in nature, and thus the MESA models that include them are not applicable. It is therefore better not to rely on the theoretical uncalibrated descriptions for these instabilities, but rather to estimate the global level of macroscopic chemical mixing in the radiative envelope, as of now denoted as $D_{\text{ext}}(r)$, from forward seismic modeling. We therefore assume that the magnetic and centrifugal forces are negligible for the computation of the equilibrium models along the evolutionary track.

Deviations from spherical symmetry due to magnetism or rotation are thus only treated at the level of the nonradial oscillations. This common approach allows us to evaluate the quality of the physical ingredients of the stellar models by exploiting the differences between $f_i^{\text{th}}$ and $f_i^*$, keeping in mind $\epsilon_i^*$. This approach is suitable as long as the departure from spherical symmetry of the equilibrium models is modest and linear pulsation theory is applicable (e.g., Ledoux 1951; Saio 1981; Dziembowski & Goode 1996; Townsend 2003a). Current pulsation computations taking full account of the deformation due to rotation are at best 2D and limited in the model assumptions (see Ballot et al. 2012; Rieutord et al. 2016, for pioneering work). One of our aims is to provide guidance to

future 2D models from application of our asteroseismic methodology. At present, such applications should be restricted to stars rotating at less than about half of their critical velocity.

Equilibrium models of stellar structure for the core hydrogen burning phase require as basic input parameters the birth mass *M* and the mass fractions representing the initial chemical composition: (*X*, *Z*) or (*X*, *Y*), where *X* represents hydrogen, *Y* helium, and *Z* the ensemble of all other chemical elements. Here, we work with (*M*, *X*, *Z*). For these three parameters, evolutionary models are computed as of the zero-age main sequence (ZAMS), for which $X_c = X$, with $X_c$ the hydrogen mass fraction in the convective core. The latter is a proxy for the stellar age during the core hydrogen burning stage, which is the application we focus on here.

### 3. Theoretical Frequency Uncertainties Stemming from Pulsation Theory

We assume that a “sparse” grid of equilibrium models for (*M*, *X*, *Z*) has been computed, with typical grid steps for intermediate-mass stars as in Table 1 of Van Reeth et al. (2016). We stress that the application of forward modeling requires pulsation computations to be done along the full evolutionary track from ZAMS to terminal-age main sequence (TAMS), with a relatively small step in $X_c$ (typically 0.001). Moreover, it is essential to keep in mind that there is a large difference in mode density as a function of $X_c$ when performing modeling. Indeed, as discussed in detail by Buysschaert et al. (2018; see Figures 2 and 3 and the discussion in Section 3.2), unconstrained forward modeling based on frequency fitting without knowing the identification of the modes artificially gives preference to near-TAMS models because there are many more modes for a fixed frequency range at evolved stages ($X_c \le 0.1$) than at earlier stages. For evolutionary stages beyond hydrogen burning, it is even more critical to have sufficiently small time steps given the much shorter timescales. Here, we only treat the case from ZAMS to TAMS.

For chosen time steps along the evolutionary track, or equivalently $X_c$ values, the equilibrium models are perturbed in a linear framework to compute the 3D eigenvectors of the oscillation modes and their accompanying eigenfrequencies. In this way, one obtains the spectrum of linear nonradial oscillations at a particular $X_c$ value. Aside from these four basic parameters (*M*, *X*, *Z*, $X_c$), a multitude of choices for the input physics occurs. In Section 4, we offer a hierarchy among choices in the specific case of stars with a convective core on the main sequence, from the viewpoint of the uncertainties for the theoretical eigenfrequencies of gravity modes.

#### 3.1. *The Traditional Approximation*

In practice, once we have an equilibrium 1D stellar structure model, its nonradial oscillation mode frequencies are computed with the publicly available code GYRE (Townsend & Teitler 2013; Townsend et al. 2018). We adopt an adiabatic framework, the Cowling approximation (Cowling 1941), a static atmosphere model, and the standard inner and outer boundary conditions as outlined in Appendix A of Townsend & Teitler (2013). A typical GYRE input list is given in Appendix A. Furthermore, we use the traditional approximation (hereafter denoted as TA; e.g., Eckart 1960; Chapman & Lindzen 1970; Unno et al. 1989; Bildsten et al. 1996; Townsend 2003a) for the computation of the gravity modes.
Mathematically, the TA is only valid for high-order gravity modes because it requires the mode frequency to remain far below the BV frequency and the tangential eigenvector component to be much larger than the radial component. Only in that case can the radial component of the Coriolis force be neglected with respect to the buoyancy force in the linearized momentum equation. However, the TA also provides a good approximation for low-order gravity modes (e.g., Section 4.2 in Lai 1997).

The TA allows oscillation modes to be computed from the Laplace tidal equations using an analytical approach (see, e.g., Townsend 2003b; Mathis 2013, for comprehensive descriptions) and offers an excellent framework for the computation of gravity-mode oscillations of slowly to moderately rotating stars that are not too deformed from spherical symmetry by rotation. From polytropic models, Ballot et al. (2010, 2012) compared 1D and 2D treatments to predict gravity-mode frequencies and showed the TA to be acceptable in the super-inertial regime: $f_j^* \lesssim 2f_{\text{rot}}$, with $f_{\text{rot}}$ the cyclic rotation frequency of the star ($\Omega_{\text{rot}}/2\pi$, often called the rotation rate). This remains valid even when the rotation rate is a significant fraction of the critical rotation rate (recall that the centrifugal distortion is proportional to $\Omega_{\text{rot}}^2$).

The limits of the TA were recently reevaluated for *Kepler* gravity-mode pulsators of spectral type F by Ouazzani et al. (2017). Formally, these authors found the TA to be appropriate for retrograde gravity modes with frequencies up to at least a quarter of the critical rotation rate; for prograde and zonal modes the TA can be used up to about half the critical rate. However, the TA was even applied to the Be star HD 163868 rotating near critically (Dziembowski et al. 2007) and compared with a second-order perturbative treatment developed by Lee & Baraffe (1995) for several Be stars based on *MOST* photometry (Cameron et al. 2008) and for the F-type star KIC 5608334 from *Kepler* photometry (Saio et al. 2018a). Despite a different behavior in the mode trapping (Figure 6 in Saio et al. 2018a), both treatments were found to be similarly appropriate when compared with the observational uncertainties of the detected frequencies in the 4 yr *Kepler* photometry of the γ Dor star KIC 5608334. Hence, it is very sensible to perform forward seismic modeling in the TA, given its computational convenience.

The TA is not needed for stars that rotate well below their critical rate (Reese et al. 2006; see also Figure 2 below). For those, first-order rotational splitting in the Ledoux approximation (Ledoux 1951) is sufficient to predict the theoretical frequencies. Such low rotation rates occur for a minority of intermediate- and high-mass stars. Nevertheless, pulsators with ultraslow rotation and high amplitudes have been found (e.g., Aerts et al. 2003; Kurtz et al. 2014; Triana et al. 2015). Schmid & Aerts (2016) made a comparison between the theoretical frequency values for a *Kepler* binary with two F-type hybrid pulsators and found the Ledoux splitting to be inadequate for pulsation frequency matching of gravity modes when $f_{\text{rot}} \gtrsim 0.1 \text{ day}^{-1}$. Here, we conduct a more systematic test by providing an overview of the differences between oscillation frequencies computed with and without the Coriolis force, for a range of stellar masses.

We computed theoretical frequencies for zonal dipole modes for a set of 12 baseline models with parameters listed in Table 1. While any forward seismic modeling application to a real star needs to consider the full evolutionary stage from the

<caption>Table 1</caption>

Parameters of 12 Benchmark Models Used to Estimate Theoretical Frequency Uncertainties of Their 660 Zonal Dipole Modes with Radial Order $n_{pg} \in [-50, +5]$, where $n_{pg} < 0$ for Gravity Modes and $n_{pg} > 0$ for Pressure Modes

<table><thead><tr><th>Mass (M<sub>☉</sub>)</th><th>X<sub>c</sub></th><th>Ω<sub>crit</sub> (rad day<sup>-1</sup>)</th><th>Ṁ (10<sup>-9</sup> M<sub>☉</sub> yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1.4 M<sub>☉</sub></td><td>0.7</td><td>22.8</td><td>0.0</td></tr><tr><td>1.4 M<sub>☉</sub></td><td>0.4</td><td>16.6</td><td>0.0</td></tr><tr><td>1.4 M<sub>☉</sub></td><td>0.1</td><td>12.1</td><td>0.0</td></tr><tr><td>2.8 M<sub>☉</sub></td><td>0.7</td><td>19.5</td><td>0.00099</td></tr><tr><td>2.8 M<sub>☉</sub></td><td>0.4</td><td>11.1</td><td>0.00135</td></tr><tr><td>2.8 M<sub>☉</sub></td><td>0.1</td><td>5.4</td><td>0.00066</td></tr><tr><td>14 M<sub>☉</sub></td><td>0.7</td><td>11.0</td><td>7.0</td></tr><tr><td>14 M<sub>☉</sub></td><td>0.4</td><td>6.3</td><td>13.4</td></tr><tr><td>14 M<sub>☉</sub></td><td>0.1</td><td>2.8</td><td>250.4</td></tr><tr><td>28 M<sub>☉</sub></td><td>0.7</td><td>8.7</td><td>303.0</td></tr><tr><td>28 M<sub>☉</sub></td><td>0.4</td><td>4.8</td><td>539.0</td></tr><tr><td>28 M<sub>☉</sub></td><td>0.1</td><td>1.6</td><td>417.0</td></tr></tbody></table>

**Note.** Fixed choices for the input parameters and input physics of these benchmark models are as follows: $\alpha_{\text{MLT}} = 2.0$, exponential core convective overshooting with $f_{\text{ov}} = 0.02$ and the radiative temperature gradient in the overshoot zone (undershooting below the convective outer envelope is not included), $Z = 0.014$, chemical mixture in Przybilla et al. (2013), OPAL opacities, the MESA EOS, the NACRE rates, the basic rate network based on the eight isotopes <sup>1</sup>H, <sup>3</sup>He, <sup>4</sup>He, <sup>12</sup>C, <sup>14</sup>N, <sup>16</sup>O, <sup>20</sup>Ne, and <sup>24</sup>Mg (Paxton et al. 2011), no rotation, no atomic diffusion, and an Eddington gray atmosphere. The oscillation frequencies are computed in the adiabatic approximation while ignoring the Coriolis force; $\Omega_{\text{crit}}$ stands for the angular critical rotation rate. We use the so-called predictive mixing scheme of MESA (available as of version 10000) to determine the position of the convective core boundary (Paxton et al. 2018).

ZAMS to the TAMS, we illustrate the differences between the mode frequencies/periods of these benchmark models and those of models of the same $X_c$ but with different input physics or different pulsational computations for only three snapshots: close to ZAMS ($X_c = 0.7$), close to TAMS ($X_c = 0.1$), and about halfway through the core hydrogen burning ($X_c = 0.4$). We provide a few example period spacing patterns in Figure 1 below, as well as histograms of the frequency differences for all comparisons we made. The masses of the benchmark models were chosen so as to encapsulate roughly the mass range for which high-order gravity modes have been detected so far. The lowest-mass models have $1.4 M_{\odot}$ with a growing convective core and a thin convective outer envelope throughout the main sequence, while the three other masses of 2.8, 14, and $28 M_{\odot}$ have a shrinking core and only develop a convective envelope past the TAMS. The two highest-mass models are subject to a considerable radiation-driven wind throughout their evolution.

Figure 1 shows the period spacing patterns for the dipole zonal modes of the models with $1.4$ and $28 M_{\odot}$, for the three considered evolutionary stages and for the indicated rotation rate. It can be seen that rotation causes a downward-tilted period spacing pattern, a property well known from theory (e.g., Bouabid et al. 2013) and also observed in *Kepler* data of pulsating F and B stars (e.g., Van Reeth et al. 2015; Ouazzani et al. 2017; Pápics et al. 2017). In Figure 2 we compare the dipole zonal frequencies of modes with the same radial order $n_{pg}$, computed without the Coriolis force versus in the TA. It can be seen in Figure 2 that the inclusion of the Coriolis force increases the gravity-mode frequencies for the entire mass
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Gravity-mode period spacing patterns for a 1.4 $M_\odot$ (left) and a 28 $M_\odot$ (right) model at $X_c = 0.7, 0.4, 0.1$, for two indicated rotation rates taken into account in the pulsation computations vs. those for a nonrotating (benchmark) model. Mode trapping occurs for all cases, except the 1.4 $M_\odot$ model at $X_c = 0.7$.</figcaption></figure>

range. The frequency differences increase for increasing radial order. Almost all frequency differences are larger than the typical frequency error expected for a light curve spanning about 0.5 to 1 yr (see Table 2). We conclude that the Coriolis force cannot be neglected in forward seismic modeling of gravity modes and that rotational effects must be included to compute theoretical gravity-mode frequencies, even for stars rotating as slowly as 10% of their critical rotation rate.

Schmid & Aerts (2016) treated two F-type super-inertial gravity-mode pulsators in a nonsynchronized close binary, with near-core rotation rates of 0.09 and $0.14\ \text{day}^{-1}$ and spin parameters in the near-core region ($2\Omega_{\text{rot}}/f_{\text{corotating}}$) ranging from 0.14 to 0.26. They assessed that the impact of ignoring the Coriolis force on the seismically modeled isochrone age is about 5%; the radii are affected by about 3%, while the masses remain unaffected for these rotation rates. These two pulsators can hence be considered as limiting cases for the validity of ignoring the Coriolis force.

Computations of frequencies in the inertial frame of reference based on the TA require estimation of the rotation frequency, $f_{\text{rot}}$, in addition to identification of the mode degree *l* and azimuthal order *m* for each frequency (Van Reeth et al. 2016; Ouazzani et al. 2017). This implies that forward seismic modeling of gravity-mode frequencies is at least a 5D optimization problem, since $f_{\text{rot}}$ must be estimated, along with ($M, X, Z, X_c$). It was shown by Van Reeth et al. (2016) and by Ouazzani et al. (2017) that the slope of period spacing series of gravity modes of consecutive radial order $n_{\text{pg}}$ in F-type pulsators allows for the simultaneous identification of the mode wavenumbers ($l, m$) and the derivation of $f_{\text{rot}}$. This method also works for B-type pulsators (Pápics et al. 2017) and is adopted here as an essential step in our seismic modeling scheme (see Section 6).

## 3.2. The Adiabatic Approximation

Moravveji et al. (2015) performed the first forward seismic modeling of a slowly pulsating B star. It was based on 19 detected dipole zonal modes with radial order $n_{\text{pg}}$ ranging from −14 to −32 for the B8.3V star KIC 10526294 (see also Triana et al. 2015). Moravveji et al. (2015) found that the adiabatic and nonadiabatic approximations to compute the mode frequencies are equivalent for the frequency precision yielded

by a nominal *Kepler* light curve for these gravity modes in that mass regime.

Figure 3 shows the frequency differences for the zonal dipole modes of the 12 benchmark models described in Table 1. Comparison of the range of the abscissa axes of Figures 2 and 3 reveals that the difference between adiabatic and nonadiabatic frequencies is almost always smaller than a typical frequency error for a 0.5 to 1 yr uninterrupted light curve and is much smaller than the effect of rotation. Moreover, various theoretical uncertainties caused by the choice of the input physics, as discussed in the next section, give rise to much larger frequency differences.

The behavior of some modes in Figure 3 attracts attention. It mainly concerns modes of the 14 $M_\odot$ models and some of the 2.8 $M_\odot$ models. We show the period spacing patterns of the converged gravity modes of the model with 14 $M_\odot$ in Figure 4. It can be seen that the modes at $X_c = 0.7$ with periods above some 1.7 days have deviating nonadiabatic behavior. These deviations between the adiabatic and nonadiabatic approximation disappear at later evolutionary stages. Similar results occur, although more sporadically, for the pressure modes of the 2.8 $M_\odot$ model, but the deviations in frequency always remain below $0.012\ \text{day}^{-1}$. Our results are in agreement with those obtained earlier by J. Daszyńska-Daszkiewicz in the framework of the HELAS project.<fnref n="7" /> Her computations revealed that the difference between adiabatic and nonadiabatic frequencies of pressure modes increases as the radial order increases and can be as large as $0.05\ \text{day}^{-1}$ for $n_{\text{pg}} = 6$ or 7. Indeed, high-order pressure modes are much more sensitive to the very outer envelope of the stars, where nonadiabatic effects occur, while the gravity modes are dominantly determined by the physics of the deep stellar interior, where the adiabatic approximation is appropriate.

Our nonadiabatic computations for the excited gravity modes reveal mode lifetimes ranging from about $10^6$ yr for the 1.4 $M_\odot$ models near the ZAMS to a minimum of only 3000 yr near the TAMS. We do point out, however, that MESA does not offer a time-dependent convection treatment, so the mode excitation predictions are not optimal for this mass. We refer to Dupret et al. (2005) for more appropriate computations. The 2.8 $M_\odot$ models reveal only a few excited modes close to the ZAMS

<footnote n="7">⁷ http://helas.astro.uni.wroc.pl/deliverables.php?active=opalmodel&lang=en</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Frequency differences for dipole zonal modes of the same radial order $n_{pg}$, for the 12 benchmark models with parameters listed in Table 1. The frequency differences $\Delta f_i^{th}$ are equal to the frequencies computed without taking the Coriolis force into account minus those where the Coriolis force is treated in the TA, for the case of a rotation rate equal to 10% (left) and 40% (right) of the critical rate. The bottom panels give the global distribution of all $\Delta f_i^{th}$-values adopting 81 frequency bins, while the inset is a zoom-in for 51 bins focusing on an interval corresponding with all $\Delta f_i^{th}$ that fall within a typical frequency error derived from a quasi-uninterrupted 0.5 to 1 yr light curve. The top panel shows the values of $\Delta f_i^{th}$ as a function of the radial order $n_{pg}$, $X_c = 0.7$ as triangles, $X_c = 0.4$ as circles, and $X_c = 0.1$ as crosses. See also Table 2 for the total percentage of the frequencies shown in the inset per mass and per evolutionary stage.</figcaption></figure>

(the $X_c = 0.7$ case), with a relatively short lifetime of some 2500 yr. On the other hand, the 14 and 28 $M_\odot$ models have only a few gravity modes excited near the TAMS ($X_c = 0.1$), with lifetimes of only about 10 and 45,000 yr, respectively. These mode excitation results are in excellent agreement with those by Szewczuk & Daszyńska-Daszkiewicz (2017), who considered models with masses between 2.5 and 19 $M_\odot$ (these authors did not list the mode lifetimes).

The required CPU time for the computation of nonadiabatic frequencies is much larger than for adiabatic computations. Moreover, uncertain physics is involved in the nonadiabatic outer envelope of the star, while the frequencies of detected gravity modes are dominantly determined by the physical properties of the deep interior. For these reasons, we adopt the TA in the adiabatic approach in the forward modeling, because such modeling has the purpose to deduce the global physical properties of the star. However, once the forward modeling results are known, instability computations should be performed for the best models, in order to test the mode excitation and the nonadiabatic frequencies of the detected modes, as a way to improve our knowledge of opacities and assess the properties of possible missing input physics.

# 4. Theoretical Frequency Uncertainties Stemming from the Equilibrium Models

As discussed in Section 3 and illustrated in Moravveji et al. (2015, 2016), the four basic parameters ($M, X, Z, X_c$) to be estimated from forward seismic modeling of stars born with a convective core are insufficient. Equilibrium models contain numerous additional parameters, some of which have a significant impact on the oscillation frequencies, while others do not. This allows us to estimate the parameters that have the largest effect on gravity-mode frequencies, after careful choice of the input physics that leaves the gravity-mode frequencies unchanged.

In order to make appropriate choices for which parameters to estimate and which to keep fixed when describing the macro- and microphysics, we need to know which processes affect the frequency values less than the resolving power of the data set. Here, we provide a hierarchy in the multitude of choices to be made among, e.g., the equation of state (EOS), opacities, nuclear reaction rates, chemical mixtures, the atmosphere model as a boundary condition, the starting model at the ZAMS, etc. Whenever a different yet “fixed choice” of input physics is made without parameters to be estimated, we adopt the terminology that we are dealing with a *different stellar model M(θ, ψ)* with free parameters *θ* and fixed physics *ψ*. For each stellar model, estimation of the optimal value of its free parameters *θ* and the accompanying uncertainties of those parameters is achieved by minimizing the discrepancy between $f_i^{th}$ predicted from $M(\theta, \psi)$ and the observations $f_i^*$. Subsequently, *stellar model selection* is done by comparing the best correspondence between observed and predicted frequencies according to a well-defined criterion, for various choices of the input physics *ψ*, using the optimally estimated parameters *θ* within each stellar model. This stellar model selection takes into account the number of free parameters in each stellar model and applies a penalty accordingly.

In practical applications, the difference between parameter estimation for a fixed type of stellar model $M(\theta, \psi)$ and selection of the most appropriate input physics represented by $\psi$ depends on the type of star. It is therefore essential that modelers have knowledge of the hierarchy among the parameters to be estimated versus the physical properties that can be fixed. For example, forward seismic modeling of pulsating subdwarfs or white dwarfs requires the mass of the thin outer hydrogen layer as a free parameter to be estimated along with the overall stellar mass, because it is vital to evaluate the mode trapping (e.g., Charpinet et al. 2011; Córsico et al. 2012), while rotation can be treated at the level of Ledoux
<table><caption>**Table 2**
Percentage of the 660 Zonal Dipole Modes of the Models in Table 1 That Fall within a Frequency Error Typical of a 0.5 to 1 yr Data Set When Changing One Aspect of the Physics</caption><thead><tr><th>Comparison</th><th>1.4 M<sub>☉</sub></th><th>2.8 M<sub>☉</sub></th><th>14 M<sub>☉</sub></th><th>28 M<sub>☉</sub></th></tr></thead><tbody><tr><td>No Coriolis ↔ TA@10%Ω<sub>crit</sub></td><td>10% (11,9,9)</td><td>12% (16,11,9)</td><td>13% (22,9,9)</td><td>17% (33,9,9)</td></tr><tr><td>No Coriolis ↔ TA@40%Ω<sub>crit</sub></td><td>1% (2,2,0)</td><td>3% (2,7,0)</td><td>5% (4,7,4)</td><td>6% (7,7,4)</td></tr><tr><td>Adiabatic ↔ nonadiabatic</td><td>100% (100,100,100)</td><td>96% (93,96,100)</td><td>83% (100,77,63)</td><td>96% (100,91,89)</td></tr><tr><td>Ledoux ↔ Schwarzschild</td><td>86% (67,91,100)</td><td>80% (80,71,89)</td><td>100% (100,100,100)</td><td>99% (96,100,100)</td></tr><tr><td>α<sub>MLT</sub> = 2.0 ↔ α<sub>MLT</sub> = 1.0</td><td>10% (8,6,15)</td><td>7% (2,6,13)</td><td>73% (55,73,91)</td><td>81% (59,82,100)</td></tr><tr><td>Gray ↔ Paczynski</td><td>59% (44,42,91)</td><td>90% (89,82,100)</td><td>93% (89,91,100)</td><td>97% (91,100,100)</td></tr><tr><td>Gray ↔ FASTWIND</td><td>...</td><td>...</td><td>94% (93,94,94)</td><td>95% (95,94,96)</td></tr><tr><td>Precomputed ZAMS ↔ from Hayashi</td><td>7% (0,0,22)</td><td>27% (0,2,80)</td><td>33% (6,8,82)</td><td>45% (20,27,82)</td></tr><tr><td>Standard network ↔ Fe, Ni added</td><td>98% (98,96,100)</td><td>90% (82,87,100)</td><td>93% (89,91,100)</td><td>100% (100,100,100)</td></tr><tr><td>OPAL ↔ OP</td><td>41% (6,51,67)</td><td>0% (0,0,0)</td><td>37% (20,41,51)</td><td>53% (26,56,74)</td></tr><tr><td>Przybilla ↔ Asplund</td><td>2% (7,0,0)</td><td>0% (0,0,0)</td><td>33% (42,51,7)</td><td>49% (56,69,22)</td></tr><tr><td>No atomic diffusion ↔ settling & levitation</td><td>...</td><td>18%</td><td>19%</td><td>43%</td></tr></tbody></table>

**Note.** The first three rows concern changes in the computation of the pulsation modes, while all other rows are results for changes in the input physics of the equilibrium models. The frequency differences are computed for models that differ less than 0.001 in $X_c$. These percentages are visually represented in histograms in the top right panels of Figures 2–10. The higher the percentage, the more justified it is to fix this phenomenon in the input physics of the stellar models for forward seismic modeling. The three values in brackets denote the subpercentages for $X_c$ = (0.1, 0.4, 0.7), to reveal the dependence on evolutionary stage. In the case of the models with atomic diffusion, the reaction network was extended with <sup>56</sup>Fe and <sup>58</sup>Ni. For the bottom row, the combined effect of gravitational settling, concentration diffusion, thermal diffusion, and radiative levitation was considered; these comparisons were limited to $X_c$ = 0.7 for CPU reasons.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Same as Figure 2, but for adiabatic vs. nonadiabatic frequencies for modes of the same radial order $n_{pg}$, showing that almost all frequency differences fall within the frequency precision typical for an uninterrupted 1 yr light curve.</figcaption></figure>

splitting for such stars and does not require use of the TA (e.g., Giammichele et al. 2017; Hermes et al. 2017). Similarly, mass loss is considered not to be of importance for intermediate-mass stars during core hydrogen burning. For high-mass stars burning hydrogen in the core, modelers should make sure that ignoring mass loss or fixing it to a specified value is appropriate by considering a second stellar model, without adding the complexity of having $\dot{M}$ as an extra parameter to estimate (see Figure 6 below). Hence, two models $M(\theta, \psi)$, one with and one without mass loss, could be considered to check its effect on the seismic parameter estimation, relative to the measurement error of the frequencies of the detected modes. For stars beyond core hydrogen burning, however, $\dot{M}$ could be included as a parameter to estimate in $\theta$.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Adiabatic and nonadiabatic gravity-mode period spacing patterns for a 14 $M_{\odot}$ model at $X_c$ = 0.7, 0.4, 0.1. When invisible, there is no difference between the black and the blue circles.</figcaption></figure>

Many additional cases to those above can be discussed in the context of considering different stellar models, the choices depending on the type of star to be modeled. Below we focus on the most important ingredients to consider as different models for the specific case of core hydrogen burning stars with a convective core, whose most important parameters for the macrophysics to be estimated will be discussed first.

### 4.1. Frequency Errors due to Choices of Macrophysics

Aside from the four basic parameters needed to compute the equilibrium models during the main-sequence life of a star, ($M$, $X$, $Z$, $X_c$), and an estimate of the near-core rotation frequency $f_{rot}$ as input for the TA oscillation frequency computations, a multitude of parameters due to macrophysical phenomena are included in state-of-the-art evolution codes. We evaluate the most important ones, from the viewpoint of the uncertainty they cause in the computation of theoretical gravity-mode oscillation frequencies.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Same as Figure 2, but for models with the Schwarzschild instead of Ledoux convection criterion (left) and for a mixing-length parameter $\alpha_{MLT} = 1.0$ instead of 2.0 (right).</figcaption></figure>

### 4.1.1. Convection and Mixing Parameters

In this work, we consider stars with $M \gtrsim 1.4 M_\odot$ at birth. The central nuclear burning for such stars is dominated by the CNO cycle, which is highly temperature sensitive. Such stars have a convective core, and its mass is of major importance for the evolution of the star. For this reason, the core mass is also a dominant parameter to estimate along with $(M, X, Z, X_c, f_{rot})$. The timescale of the convection in the core is much shorter than the nuclear timescale, so we can assume full and instantaneous mixing there. Hence, the core has uniform chemical composition as long as hydrogen burning is active. We adopt the Ledoux criterion and use MESA's predictive mixing scheme to locate the boundary of the convective core for the equilibrium models in Table 1. We show in the left panel of Figure 5 the difference in the oscillation frequencies if we instead adopt the Schwarzschild criterion. It can be seen that ignoring the $\mu$-gradient in the stability criterion somewhat impacts the gravity modes for the 1.4 and 2.8 $M_\odot$ models, and more so at evolved rather than at early evolutionary stages. The frequency differences always remain below $10^{-1}$ day$^{-1}$. Paxton et al. (2018) provide a thorough discussion on the subtle differences that may occur in the inhomogeneous region beyond the core boundary owing to the choice of the convection criterion, depending on the stellar mass and evolutionary stage. However, as becomes clear from Table 2, the implications for the modes due to these small differences are inferior to those of several other effects connected with choices for the input physics.

The observations of gravity-mode period spacings in a sample of F stars (e.g., Van Reeth et al. 2015) and B stars (e.g., Pápics et al. 2017) reveal detailed structure. Indeed, the period spacing morphology appreciably deviates from a constant value, even for stars that hardly rotate (Degroote et al. 2010; Pápics et al. 2014). Estimation of the core mass from forward seismic modeling of core hydrogen burning stars based on gravity modes has so far only been done for a few stars, and all of them have demanded the inclusion of convective core overshooting and chemical mixing in the radiative envelope

(Moravveji et al. 2015, 2016; Schmid & Aerts 2016). Each of these two necessary ingredients is described by at least one parameter. The overall mass in the core is dependent on the treatment of the overshooting, while the chemical profiles throughout the star are determined by the element transport in the radiative envelope. We describe the physical meaning of these two required ingredients here.

As in many stellar evolution codes, convective energy transport in MESA relies on the time-independent mixing-length theory, while time-dependent convection (turbulent pressure) is ignored. Mixing-length theory is characterized by one free parameter, $\alpha_{MLT}$, representing the characteristic length scale over which convective fluid elements travel before they dissipate; it is expressed in local pressure scale heights and takes typical values between 1 and 2. For low-mass stars, $\alpha_{MLT}$ determines the convective motions in the outer convective envelope, and one often fixes it to the solar value. However, the intermediate-mass stars we consider here have very thin and shallow convective envelopes. Moreover, along with the high-mass stars, they have a fully mixed convective core, where physical circumstances are very different than in Sun-like stars. For this reason, $\alpha_{MLT}$ is a free parameter. We show in the right panel of Figure 5 the effect of changing $\alpha_{MLT}$ from 2.0 to 1.0. It can be seen from Figure 5 and deduced from Table 2 that this change has a major impact for the two lowest masses considered here, particularly at evolved stages. This is in line with observational results for $\delta$ Sct stars (Bowman & Kurtz 2018). The choice for $\alpha_{MLT}$ is less important for the two high-mass stars, particularly close to the ZAMS. The results in Figure 5 are as expected since the two lower-mass models have shallow outer convection zones whose extent is determined by $\alpha_{MLT}$, while the high-mass stars have radiative envelopes early in their lives but develop convective envelopes near the TAMS. Given the results in Table 2, one should either estimate $\alpha_{MLT}$, particularly for intermediate-mass stars, or consider at least two different fixed values to assess the impact on the forward modeling.

Complications to the above picture arise in the determination of the core boundary because convective fluid elements
“overshoot” the boundary between the convective core and radiative envelope owing to their inertia. The extent and properties of this overshoot region, including the level of mixing of the chemical elements inside it and its thermal structure, are of major importance for the future evolution of the star (Kippenhahn et al. 2012, Section 30.4, for a thorough discussion). Gravity modes allow us to seismically probe the detailed properties of this near-core overshoot region. The parameters describing the aspects of convective core overshooting are among the most important ones for modeling of the stars in the mass range we consider (see Section 3 in Paxton et al. 2018, for a thorough discussion).

In MESA, convection is implemented as a diffusive process with an associated mixing coefficient. In such a framework, the convective core has a large diffusive mixing coefficient $D_{\text{mix}}$, typically between $10^{11}$ and $10^{17}$ cm$^2$ s$^{-1}$ during core hydrogen burning for the mass regime we consider here. As such, it is natural to describe the mixing due to core overshooting as a diffusive process as well, with coefficient $D_{\text{ov}}$ containing at least one dimensionless parameter (Equation (9) in Paxton et al. 2013). Several options are available to define the level and shape of overshoot mixing, each characterized by a set of free parameters—Paxton et al. (2013, Appendix B.7) and Pedersen et al. (2018). Among them is the option to consider an exponentially decaying prescription described by the parameter $f_{\text{ov}}$, and the simpler step function overshoot for which the core is extended into a fully mixed overshoot region from its boundary by an amount $\alpha_{\text{ov}}$.

The temperature gradient in the overshoot region is unknown (see Salaris & Cassisi 2017, for a thorough discussion on chemical mixing and temperature gradients). In the standard version of MESA, it is the radiative gradient, but the code is sufficiently flexible to change this into the adiabatic one (so-called convective penetration) or any transition between these two. These choices each constitute a different stellar model $M(\theta, \psi)$. In regions where the temperature gradient takes a value between the adiabatic one and the adiabatic one increased with a positive $\mu$ gradient, semiconductive mixing occurs. This phenomenon describes the growing oscillatory motion of a convective fluid element in that region, which implies partial chemical mixing in that zone. The effect of semiconvection on the temperature gradient is unknown. In MESA, semiconductive mixing in the vibrationally unstable region is again implemented as a diffusive process, characterized by the dimensionless coefficient $\alpha_{\text{sc}} \in [0, 1]$. The larger the coefficient, the faster the matter mixes in the semiconductive zone.<fnref n="8" />

As already discussed above, “extra” macroscopic mixing $D_{\text{ext}}$ occurs in the zone above the core overshoot zone in the radiative envelope of intermediate- and high-mass stars. Indeed, such extra mixing was required to bring theoretical frequencies in agreement with observational data (Moravveji et al. 2015, 2016). Coherent oscillation modes (Aerts et al. 2014) and damped modes (Rogers & McElwaine 2017) may also induce mixing, but these options have so far hardly been taken into account in stellar evolution theory. In order not to introduce bias on the importance of the various phenomena

giving rise to element transport, we take a pragmatic approach and make use of one diffusive mixing coefficient $D_{\text{ext}}$ as a free parameter in the forward seismic modeling. This parameter represents the average overall effect of (rotationally and/or magnetically and/or pulsationally induced) mixing in the radiative zones of a star as probed by the gravity modes. Only if this free parameter can be estimated with good precision from seismic data for a considerable sample of stars can we start to unravel the role of the various individual causes of mixing at play in the stellar envelope. As said, the first steps in this direction were taken for a few stars (Moravveji et al. 2015, 2016; Schmid & Aerts 2016), which showed the need to include $D_{\text{ext}} \neq 0$ just adjacent to the core overshoot region in the radiative envelope to explain the structure in the observed period spacings. Meanwhile, Pedersen et al. (2018) assessed whether gravity modes can distinguish between a constant $D_{\text{ext}}$ and an actual profile $D_{\text{ext}}(r)$ stemming from 2D simulations of internal gravity waves. These authors found that a distinction can only be made on the basis of combined observational restrictions from gravity modes and surface abundances. In the first instance, it is therefore justified to consider $D_{\text{ext}}$ to be a single parameter, describing the mixing adjacent to the core overshoot zone to be estimated.

As a side remark, we point out that extra turbulent mixing in the surface layers of intermediate-mass stars is sometimes also considered in the context of asteroseismology or galactic archeology. Such extra turbulent mixing in the near-surface layers is usually assumed to have a density dependence. Aside from this, envelope undershooting may also occur. While both these phenomena are important for chemical tagging at advanced evolutionary stages of low-mass stars (e.g., Dotter et al. 2017), they hardly affect the gravity-mode period spacings (Pedersen et al. 2018). Indeed, mixing in the zone adjacent to the core overshoot zone is much more important for the fitting of detected gravity-mode frequencies, because these modes mainly probe the near-core regions. At first instance, such surface mixing can therefore be ignored in forward seismic modeling of gravity modes but can be evaluated from high-precision surface abundances if available.

The transport of angular momentum during a star’s evolution is also not calibrated by observations, just as the element transport. Moreover, the level of its inhibition by $\mu$-gradients is not known. Even for the simple case of low-mass stars born with a radiative core, where the macroscopic mixing is thought to be absent, current theory of angular momentum transport fails (e.g., Eggenberger et al. 2017). The importance of a convective core and near-core properties during core hydrogen burning for appropriate understanding and interpretation of angular momentum of red giants was recently highlighted (Tayar & Pinsonneault 2013; Eggenberger et al. 2017). Seismic modeling of intermediate-mass stars from gravity modes during the core hydrogen burning is highly relevant to understand the angular momentum of secondary clump red giants, as well as more massive blue supergiants.

Given the major uncertainties in the transport of chemicals and angular momentum in current stellar evolution theory, an appropriate approach in terms of forward seismic modeling is to estimate the chemical mixing by seismic inference from gravity modes, assuming rigid rotation near the core as the first step (see Van Reeth et al. 2016; Aerts et al. 2017; Ouazzani et al. 2017). Subsequently, the rotation profile and the accompanying angular momentum transport can be evaluated

<footnote n="8">⁸ In this work, we ignore thermohaline mixing (“fingering convection”) due to a negative $\mu$-gradient value in the models. This phenomenon of slow mixing typically occurs in evolved low-mass stars with shell burning. It can also occur during the hydrogen core burning, but the mixing profile it would cause is not well known (Deal et al. 2017). Should a profile become available, this phenomenon could be included separately with yet another diffusive coefficient $\alpha_{\text{th}}$, and it could easily be included in our statistical formalism.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Same as Figure 2, but for a gray vs. Paczynski static atmosphere model as a boundary condition for MESA (left) and for a gray atmosphere vs. FASTWIND dynamical atmosphere model for the two highest masses (right).</figcaption></figure>

by considering $f_{\text{rot}}(r)$ instead of constant $f_{\text{rot}}$ for the best seismic models, once their mixing properties have been derived seismically. This approach has been applied to the above-mentioned $3.25 M_{\odot}$ star KIC 10526294 by Moravveji et al. (2015) and Triana et al. (2015), pointing to the need of angular momentum transport. The required transport could meanwhile be explained by internal gravity waves (Rogers 2015). Such types of analyses need to be done for stars with various levels of near-core rotation to deduce a broad range of $f_{\text{rot}}(r)$ covering the entire hydrogen-burning phase, in order to be able to predict angular momentum transport for the more evolved phases (Aerts et al. 2017; Ouazzani et al. 2018).

In summary, near-core rotation, convective core overshooting, and mixing in the radiative envelope adjacent to the core overshoot zone must all be included in the forward seismic modeling of gravity-mode pulsators with $M \gtrsim 1.4 M_{\odot}$, implying that the optimization problem is at least seven-dimensional (+7D), with $\theta = (M, X, Z, X_c, f_{\text{rot}}, D_{\text{ov}}, D_{\text{ext}})$. For intermediate-mass stars, one should additionally test the impact of a different mixing-length value.

### 4.1.2. Choice of a Static Atmosphere as a Boundary Condition

All stars experience mass loss during their life. The occurrence and duration of the episodes of mass loss are vastly different for stars with different masses. Low- and intermediate-mass stars born with $M \lesssim 9 M_{\odot}$ experience substantial mass loss due to a dust-driven wind on the asymptotic giant branch, while more massive stars experience a radiation-driven wind during and/or beyond core hydrogen burning until they explode as supernovae. Two important parameters characterizing the wind are the mass-loss rate $\dot{M}$ and the wind velocity profile (Kudritzki & Puls 2000).

Stellar evolution computations treat the stellar wind in a simplistic way in the sense that a mass-loss rate from a recipe is adopted and the appropriate mass is peeled off the outer layers of a star at each time step in the evolution, taking into account the accompanying angular momentum loss and assuming that hydrostatic equilibrium holds at every time step. This is a crude

approximation, because the driving force of the wind is ignored in the equation of motion and a static rather than dynamical atmosphere is considered as a boundary condition to solve the equations of stellar structure. MESA and GYRE offer various static atmosphere models from precomputed tables as choices for these boundary conditions. The effect of different choices of the atmosphere model on the frequency values of gravity mode is illustrated in the left panel of Figure 6, where we show the frequency differences obtained for a gray versus Paczynski atmosphere model. It can be seen that frequencies are only modestly affected and that the majority of frequency differences are within the resolving power of a typical data set, except for the lowest considered mass $M = 1.4 M_{\odot}$ at evolved stages. The latter models have a thin convective outer envelope, and the interplay among convection, pulsation, and atmospheric properties becomes important for the modes. For stars without a (thin) convective outer envelope, the range of the frequency differences is always below $0.01 \text{ day}^{-1}$ (Figure 6), which justifies that we neglect this effect compared to other more important ones.

### 4.1.3. Effect of Mass Loss due to a Radiation-driven Wind

So far, the mass-loss rate has not been considered as a parameter to estimate in forward seismic modeling for core hydrogen burning stars, given that one is already dealing with a + 7D optimization problem. Cluster ensemble asteroseismology from *Kepler* low-mass red giants did allow a rough evaluation of the upper limit of mass-loss rate $\dot{M}$ due to a dust-driven wind on the red giant branch (Miglio et al. 2012). However, the potential to estimate $\dot{M}$ seismically and take into account angular momentum loss during stellar evolution is appealing in view of the core rotation rates found for intermediate-mass stars (Aerts et al. 2017).

We test here how the parameter $\dot{M}$ affects the gravity-mode frequency predictions. As shown in the left panel of Figure 6, forward seismic modeling of gravity modes is not affected by the choice of the atmosphere in the case of static models to compute the boundary conditions, because such modes are
determined by the properties of the deep stellar interior. Nevertheless, the mode cavities are affected by the stellar wind, and for a dynamical atmosphere, energy wave leakage and angular momentum loss are predicted to occur (Townsend 2000a, 2000b). We considered a simplified test by comparing the difference in adiabatic gravity-mode frequencies for a static gray atmosphere versus a dynamical atmosphere for particular values of the mass-loss rate as indicated in Table 1. For each of the parameters of those benchmark models in Table 1, we took the outcome of the MESA model for the luminosity, radius, and effective temperature in the case of a gray atmosphere. Along with the mass and mass-loss rate, these values were then used for the calculation of a radiation-driven stellar wind structure with the non-LTE (NLTE) radiative transfer code FASTWIND (Puls et al. 2005). In this code the wind structure is derived from the mass continuity equation and a predescribed velocity field derived from standard radiation-line-driven wind theory, $v(r) = v_\infty(1 - R_*/r)^{0.8}$, with the stellar radius $R_*$ and the terminal wind speed $v_\infty$ as given by MESA. NLTE Hopf functions are used for the temperature stratification using an exact $T(\tau_R)$ relation from a converged plane-parallel NLTE model (with $\tau_R$ the Rosseland optical depth) and taking into account sphericity effects (see Santolaya-Rey et al. 1997, for an extensive description). This atmosphere and wind structure was then stitched to the interior structure from MESA. The point of connection was chosen as the one where the opacities of both models match best at $\tau_R = 2/3$. A smooth and continuous connection was enforced by linear interpolation of the various quantities in log space around the point of connection.

The comparison between the frequencies of the modes for a static gray atmosphere versus a dynamical FASTWIND atmosphere is shown in the right panel of Figure 6 for the highest two considered masses (for the lower masses, computations with FASTWIND are not meaningful). We find it justified to ignore the wind for gravity modes, but not for pressure modes. However, we note that our simple test here only accounts for the wind outflow indirectly through a modified atmospheric structure, i.e., we do not explicitly account for potential wave leakage when solving the pulsation equations within GYRE. As such, the effects seen in Figure 6 quite likely constitute lower limits.

Any radiation-driven wind effect has so far been ignored in seismic modeling of β Cep stars based on ground-based data. Given that stars born with a mass above some $24 M_\odot$ experience strong mass loss throughout their entire life and stars with $M \gtrsim 3 M_\odot$ as of the blue supergiant phase, $\dot{M}$ should be considered as an additional parameter to estimate from forward seismic modeling of such stars. Unfortunately, seismic data sets to achieve this are not yet available for stars with mass above some $24 M_\odot$. Major progress on this front is within reach of the TESS continuous viewing zones (Ricker et al. 2016) and the long pointings of the PLATO mission (Rauer et al. 2014).

### 4.1.4. Treatment of the Pre-main-sequence Evolution

Solving the stellar structure equations requires not only boundary conditions at the center and the surface of the star but also initial conditions. In essence, an initial “birth model” is required to “start” the evolution from core hydrogen burning in full equilibrium, and the age of the star is determined as of that so-called ZAMS model. However, partial nuclear burning already occurs during the contraction phase of the protostar, so there is no one-to-one equivalence between the stellar age and

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as Figure 2, but comparing frequencies for models computed from the Hayashi track vs. models based on an interpolation within a precomputed grid of ZAMS models.</figcaption></figure>

$X_c$ at the ZAMS. Moreover, the physical processes during star formation, including rotation, magnetism, and accretion during contraction from the Hayashi track toward the start of core hydrogen burning in equilibrium at the ZAMS, remain largely unknown. In view of this, one tends to circumvent pre-main-sequence evolutionary computations. MESA comes with a library of precomputed ZAMS models for $Z = 0.02$ and with masses between 0.08 and $100 M_\odot$ (Paxton et al. 2011). An option to start the evolutionary computations offered to users is to interpolate between these precomputed ZAMS models, in order to achieve a ZAMS model for the chosen mass. We have tested this procedure against the more time-consuming computation of contraction models for the chosen $(M, X, Z)$ from the Hayashi track to the ZAMS to end up with a self-consistent ZAMS model according to the chosen $(M, X, Z)$ and frozen input physics as initial conditions. Comparing the frequencies between those two cases is nontrivial, as one in fact evolves two slightly different stars as of the ZAMS. For this comparison, Figure 7 shows the difference in the oscillation frequencies of models at equal values of $X_c$ to within 0.001. This figure demonstrates that taking an interpolated ZAMS model from MESA’s precomputed libraries versus a ZAMS model resulting from contraction starting on the Hayashi track with consistent input physics for the chosen $(M, X, Z)$ changes the frequencies more than the typical frequency error derived from a 0.5 to 1 yr light curve for the majority of modes. It is therefore essential to compute models from the Hayashi track to end up with proper initial conditions in forward seismic modeling.

Unlike the interior rotation, core overshooting, and envelope mixing, we cannot treat the “global characteristics” of the ZAMS model as a simple free parameter to estimate during the seismic modeling. This aspect is rather a typical case where various choices for the input physics during the pre-main-sequence phase each constitute a different stellar model $M(\theta, \psi)$. The dependence of the frequency values on the pre-main-sequence structure is good news for future forward seismic modeling of pre-main-sequence pulsators. While
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Same as Figure 2, but comparing frequencies for models computed with the basic network vs. a network with <sup>56</sup>Fe and <sup>58</sup>Ni added.</figcaption></figure>

applications of seismic modeling to this early phase of stellar evolution must await suitable photometry from space, the same methodology as described here can be used for these early evolutionary stages once sufficiently identified modes can be measured (e.g., Zwintz et al. 2014, 2017).

## 4.2. Errors due to Choices of Microphysics

Aside from the aspects of the macrophysics discussed above, a myriad of choices for the microphysics can also be made to compute the stellar models. It concerns, e.g., the EOS, thermonuclear and weak reactions, nuclear reaction networks, opacities, etc. Here, we focus on the core hydrogen burning phase for stars that do not undergo electron degeneracy. Such cases are the best-understood ones in terms of the EOS (an ideal gas with radiation is appropriate). Moreover, the choices of the nuclear reactions and of the nuclear reaction networks among a variety of sets affect the seismic properties least during this simplest evolutionary phase, but they may appreciably alter the more advanced phases of stellar evolution. For this work, we used the standard MESA EOS, the NACRE rates (Angulo et al. 1999), and the MESA basic rate network based on eight isotopes: <sup>1</sup>H, <sup>3</sup>He, <sup>4</sup>He, <sup>12</sup>C, <sup>14</sup>N, <sup>16</sup>O, <sup>20</sup>Ne, and <sup>24</sup>Mg (Paxton et al. 2011). As an illustration of the effect of such a choice, we also computed the oscillation frequencies in the case where the <sup>56</sup>Fe and <sup>58</sup>Ni isotopes were added to the network and provide the comparison in Table 2 and Figure 8. It can be seen that, as expected, this only has a minor effect compared to other aspects of changing the input physics.

We now turn our attention to some of the more critical aspects of the microphysics that affect the oscillation frequencies appreciably.

### 4.2.1. Opacities and Chemical Mixtures

It is well known that the choice of opacities affects the oscillation properties of B-type pulsators appreciably (e.g., Daszyńska-Daszkiewicz & Walczak 2010). In their extensive study of KIC 10526294, Moravěji et al. (2015) quantified the effect of choosing different opacity tables and mixtures for the equilibrium models on the resulting frequency predictions for

gravity modes (see their Figure 3). Two different sets of opacities—OPAL and OP—were combined with three different chemical mixtures to end up with six different models $M(\theta, \psi)$, adopting the MESA EOS. The six corresponding dipole zonal frequencies for the mode with radial order $n_{pg} = -32$ differ up to $\simeq 0.004 \text{ day}^{-1}$, while the frequency differences for the mode with radial order $n_{pg} = -14$ amount to $\simeq 0.008 \text{ day}^{-1}$. This already demonstrates that the theoretical frequency uncertainty due to different choices of opacity tables cannot be ignored.

Figure 9 shows the deviations between the theoretical frequencies for the benchmark models in Table 1 due to a different choice of opacities (left) and chemical mixture (right). Both these choices have a major effect on the frequency values in the sense that the differences are larger than the frequency error of a typical 0.5 to 1 yr light curve. Both the pressure-mode and gravity-mode frequencies get shifted with a systematic effect that cannot be ignored. We conclude that, with the current poor knowledge of opacities, asteroseismology in the mass regime considered here should be based on various stellar models $M(\theta, \psi)$ relying on different opacity tables and chemical mixtures.

### 4.2.2. Microscopic Atomic Diffusion

In addition to the aforementioned parameter $D_{\text{ext}}$ to be estimated in the radiative envelope, composition changes may also occur as a result of transport processes following microscopic atomic diffusion in the radiative layers of the star. The computation of such microscopic processes requires detailed evaluation of frequency-dependent absorption coefficients and reevaluation of opacities from all the relevant ions involved in the various layers of the star, for each time step along the evolution. This obviously constitutes an immense challenge for the evolutionary computations (CPU-wise) and is usually omitted for massive stars, as these evolve on a timescale of millions of years. For cool stars with an extensive convective envelope, it is necessary to consider models with gravitational settling for evolutionary computations (e.g., Michaud et al. 2004) and for proper chemical tagging for galactic archeology (Dotter et al. 2017). Here, we investigate the effect of atomic diffusion for intermediate- and high-mass stars, which also requires the inclusion of radiative levitation (e.g., Deal et al. 2016, 2017).

Within MESA, four components of atomic diffusion can be included, each for various chemical elements: gravitational settling, concentration diffusion, thermal diffusion, and radiative levitation (Hu et al. 2011; Paxton et al. 2015). Including atomic diffusion in MESA requires any of these phenomena to be turned “on” or “off,” i.e., the activation of any of these particular microscopic effects can be included or not. Moreover, each of the phenomena can be included with or without a parameter to be estimated. In our methodological setup, including atomic diffusion for each different chemical element would imply for each a different stellar model $M(\theta, \psi)$. For low-mass stars (e.g., Verma et al. 2017) and for white dwarfs (e.g., Romero et al. 2017), gravitational settling is an important ingredient for asteroseismic modeling. On the other hand, radiative levitation acting on Fe and Ni is critical to take into account when estimating oscillation frequencies for hot subdwarfs in core helium burning (e.g., Charpinet et al. 2011; Bloemen et al. 2014).
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Same as Figure 2, but comparing frequencies for models based on OPAL vs. OP opacities (left) and for models based on the solar mixture by Asplund et al. (2009) vs. the standard mixture of OB stars in the solar neighborhood by Przybilla et al. (2013), adopting OPAL opacities for both mixtures (right).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Same as Figure 2, but comparing models without and with atomic diffusion, taking into account gravitational settling, concentration diffusion, thermal diffusion, and radiative levitation. This comparison is limited to the case of $X_c = 0.7$ for computational reasons. For all these cases, $^{56}$Fe and $^{58}$Ni were added to the network.</figcaption></figure>

For the stars under study here, we show in Figure 10 the comparison between the frequencies due to models without and with atomic diffusion. For this comparison, we used the treatment by Hu et al. (2011), as implemented in MESA, taking into account gravitational settling, concentration diffusion, thermal diffusion, and radiative levitation for the elements H, He, C, N, O, Ne, Mg, Fe, and Ni. For this comparison we limited to the models from the ZAMS until $X_c = 0.7$ for CPU reasons, where we adopted the most recent MESA implementation (Paxton et al. 2018). In this comparison, we omit the models with $M = 1.4 M_\odot$, because settling is the dominant effect for those, and this requires some level of turbulent mixing (combined or not with envelope undershooting) to be included in the model computations, in order to end up with

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Change in the gravity-mode period spacing patterns for the 2.8 and 14 $M_\odot$ models at $X_c = 0.7$ when settling alone and levitation in addition are considered.</figcaption></figure>

appropriate surface abundances of metals (see, e.g., Turcotte et al. 1998, for pioneering work). One usually takes such ad hoc turbulent mixing to be a function of density (e.g., Richer et al. 2000; Deal et al. 2016; Dotter et al. 2017). Since our baseline models do not have such density-dependent turbulence, nor undershooting, the comparison with the baseline models indicated in Table 2 is not meaningful for this mass.

Figure 10 and Table 2 reveal that atomic diffusion has a large effect on the frequencies for all the modes. We find that this is mainly due to the radiative levitation for the three masses we considered here. This is illustrated for the period spacing patterns of the 2.8 and 14 $M_\odot$ models in Figure 11, where it is seen that the mode trapping changes when radiative levitation is taken into account, because it modifies the $\mu$-gradients and hence the convection zones. The importance of radiative levitation was already stressed by Deal et al. (2017) in their study of the solar-like oscillations of the F-type exoplanet host star 94 Ceti. Radiative levitation has so far hardly been included in forward seismic modeling of gravity modes, mainly because it is too demanding in computation time. Our
results show it to be worthwhile to undertake future studies on
how to include the combined effect of settling and levitation in
the forward modeling of intermediate- and high-mass stars.

As a side remark in this section, we stress that atomic
diffusion may be of importance locally in and near the driving
layer in the envelope for mode excitation of opacity-driven
modes during the main sequence, since more modes are
detected than predicted to be excited in massive stars (e.g.,
Moravveji et al. 2016; Daszyńska-Daszkiewicz et al. 2017).
For this reason, mode excitation is not part of our forward
seismic modeling scheme but is treated as a posteriori
evaluation of how good the best seismic models are in terms
of explaining overall stellar properties, along with other
observables such as a spectroscopic T<sub>eff</sub> and log g, surface
abundances, and luminosities deduced from Gaia astrometry.

In general, our statistical framework outlined in the
following sections is only based on gravity-mode frequencies
and does *not* rely on surface properties from spectroscopy, nor
on astrometry or interferometry. The reason is that the gravity
modes provide direct observational information on the deep
stellar interior, at the level of εᵢ* / fᵢ* ≪ 0.1%, while the relative
precisions for T<sub>eff</sub>, log g, and abundances from spectroscopy of
stars in the mass range considered here are far worse, often
above 10% (e.g., Morel et al. 2008; Martins et al. 2012a,
2012b; Van Reeth et al. 2015). Moreover, our aim is to tune the
near-core interior properties of the stars. After deducing the
best seismic models, we will confront their predictions with
measured surface properties, to assess the physical conditions
in the outer envelopes of the stars. In this respect, our setup and
application are quite different than the case of low-mass stars
with pressure modes. Of course, whenever additional spectro-
scopic, interferometric, or astrometric data become available at
a precise enough level for our scientific aims, they can also be
used as *input* in the forward modeling, together with the
oscillation modes. This is not yet the case for the mass regime
we treat here.

**5. Statistical Model Formulation and Parameter Estimation**

We now place the principles outlined in the previous
Sections 2–4 into a formal and appropriate statistical frame-
work. We assume that we have various stellar model grids
based on different stellar models M(θ, ψ), along with their
pulsational frequencies, at our disposal. We want to reach two
major goals:

1. To fit as closely as possible the observed and identified gravity-mode oscillation frequencies $f_i^*$ ($i=1, ..., n$) by model (theoretical) values $f_i^{\text{th}}$. As outlined above, we are in the situation that some of the uncertainties on the theoretical frequencies are smaller than the errors of the observed frequencies, while others are larger or of the same order of magnitude. The overall theoretical errors (for $f_i^{\text{th}}$) are unknown, because additional factors to those shown in Figures 2–10 are at play, such as the limitation to 1D stellar structure models, the absence of certain forces, missing input physics, etc. Hence, heteroscedasticity must be included in the statistical model formulation.

2. To select the most likely physical model $\mathcal{M}(\theta, \psi)$ for a single star or for an ensemble of stars, without introducing a priori bias about the unknown input physics, and in particular keeping in mind that we are

dealing with a different physical regime compared to
Sun-like stars with a radiative core: here, the Coriolis
force, core overshooting, and transport of elements in the
radiative envelope are dominant phenomena that cannot
be ignored.

This setting is different from many statistical problems of
parameter estimation, where a single parameter vector is
estimated from a set of repetitions. Here, we are *not dealing with replicates* of one and the same mode frequency, but rather
with a single measurement for a set of observed frequencies $f_i^*$ ($i = 1, ..., n$) belonging to modes that each have their own
mode cavity and probing power for the interior physics. Each
of the measurements, $f_i^*$, is accompanied by its own single
error measurement $\epsilon_i^*$ for which we assume $\epsilon_i^* \sim N(0, \lambda_i)$.

We now introduce convenient vectorial notations to outline
the above aims in suitable statistical frameworks. An observa-
tion takes the form of a vector **Y**\*, consisting of the observed
frequency values, **Y**\*<sub>i</sub> = **f**\*<sub>i</sub> with *i* = 1, ..., *n*, where the length
of the vector, *n*, is star specific (the number of identified gravity
modes). Let us denote the measured frequency precisions as Λ,
a diagonal matrix with elements λᵢ. Given a physical model
*M*(θ, ψ), a vector θ of stellar parameters, and a vector of fixed
properties ψ that is assumed known, the model predicts the
observation: **Y**(θ, ψ), where **Y** is a vector of *n* predicted
frequencies **f**\*<sub>i</sub> from model *M*(θ, ψ) for the identified
frequencies corresponding with **Y**\*<sub>i</sub>. Denote by *p* the length of
vector θ. For convenience, ψ will sometimes be dropped from
notation: **Y**(θ) ≡ **Y**(θ, ψ). The inclusion of a parameter ψ in
the statistical model allows us to separate key parameters to be
optimized (grouped in θ) from the ones that are considered
fixed for the fit. This is necessary when the overall parameter
vector (grouping θ and ψ) is high-dimensional, as in our
problem of forward seismic modeling.

In our application, the number of detected gravity modes
with well-determined frequency, *n*, is typically between 10 and
40 when dealing with nominal *Kepler* data but smaller for other
data sets, while *p* is at least 7: aside from the four standard
input parameters for the computation of stellar evolution,
(*M*, *X*, *Z*, *X*<sub>c</sub>), we need to estimate as well the near-core
rotation with period *P*<sub>rot</sub> = 1/*f*<sub>rot</sub> = 2π/*Ω*<sub>rot</sub>, at least one
parameter connected with the core overshooting, *D*<sub>ov</sub>, and
one parameter describing the mixing in the radiative envelope
responsible for the mode trapping. Ideally, however, *p* is larger,
as one also wants to estimate the *profiles* of the core overshoot,
envelope mixing, and interior rotation. As a concrete example
of this approach in the literature, θ was taken to be
(*M*, *X*, *Z*, *X*<sub>c</sub>, *D*<sub>ov</sub>, *D*<sub>ext</sub>) (hence *p* = 6), while ψ stood for fixed
values of α<sub>MLT</sub> = 1.8, α<sub>sc</sub> = 0.01, M̄ = 0.0, along with the
fixed choices of the MESA EOS, the chemical mixture by
Przybilla et al. (2013), and OPAL opacities, in “Model 4” of
the study of the ultraslow rotator KIC 10526294 by Moravveji
et al. (2015). The latter authors tested the performance of this
Model 4 with respect to a simpler stellar model having no extra
diffusive mixing in the envelope, i.e., a stellar model with
θ = (*M*, *X*, *Z*, *X*<sub>c</sub>, *D*<sub>ov</sub>) (*p* = 5) and ψ defined by α<sub>MLT</sub> = 1.8,
*D*<sub>ext</sub> = 0.0, α<sub>sc</sub> = 0.01, M̄ = 0.0 with the same fixed choices
of the MESA EOS, chemical mixture, and opacities (“Model 1”
in Moravveji et al. 2015). The result was that Model 4 gave a
better representation of the detected zonal dipole gravity-mode
frequencies than Model 1, hence the conclusion that this star
undergoes envelope mixing in addition to core overshooting. In
such an experiment, Models 1 and 4 are called *nested models*, because Model 4 has one more parameter to estimate than Model 1 but is otherwise the same.

Commonly, the estimation of $\theta$ is evaluated in all points of a grid. This grid can result from a Cartesian product in the components of $\theta$, but this is not a necessary requirement. Also, in practical applications, the grid will be finer in one component than in another. Let us index the grid points by a single index $j = 1, ..., q$. This leads to a collection of $q$ theoretically generated vectors in each grid point $j$, denoted as $\mathbf{Y}_j \equiv \mathbf{Y}_j(\theta_j, \psi)$. These vectors $\mathbf{Y}_j$ each have length $n$. They result from the theoretical model $\mathcal{M}(\theta, \psi)$ that lies at the basis of the grid. We consider and distinguish between the following four problems to solve:

Problem 1: Given one set of observed identified frequencies $\mathbf{Y}^*$ and a theoretical model $\mathcal{M}(\theta, \psi)$ used to compute a grid with $j = 1, ..., q$ grid points, find the most plausible value for $\theta$ from the grid, i.e., select $\theta_j$ for which $\mathbf{Y}_j(\theta_j, \psi)$ comes closest to $\mathbf{Y}^*$. By extension, a small set of candidate $\theta_j$ can be selected. Assess the errors of the selected $\theta_j$.

Problem 2: We wish to compare the capacity of various theoretical models to describe the data. In order to tackle this, it is convenient to compare the estimates among two theoretical models, $\mathcal{M}(\theta^{(m)}, \psi^{(m)})$ and $\mathcal{N}(\theta^{(n)}, \psi^{(n)})$. The problem is thus for a given value of the vector $\theta^{(m)}$ from a benchmark model of grid ($m$), to find candidates $\theta^{(n)}$ from grid ($n$) that provide the best match.

Problem 3: Search for a statistical model to replace theoretical model $\mathcal{M}(\theta, \psi)$, delivering an outcome vector sufficiently close to $\mathbf{Y}$ but much faster to evaluate and lending itself for error estimation of $\theta$ and for statistical model selection.

Problem 4: Find the most adequate theory $\mathcal{M}(\theta, \psi)$, among a collection of competing theories, using a set of $N$ stars with $\mathbf{Y}_i^*$ (where $t = 1, ..., N$).

We provide methodology to treat each of these four problems in the following subsections, adopting vectorial notation as is common in statistics.

### 5.1. Problem 1: A Mahalanobis-distance-based Solution to Match the Identified Gravity Modes of One Star

The problem is to find the value for $\theta$ that best predicts the stellar observation $\mathbf{Y}^*$, say, $\theta_0$, with corresponding value $\mathbf{Y}_0 \equiv \mathbf{Y}_0(\theta_0, \psi)$. We need to select $\theta_0$ such that the distance between $\mathbf{Y}^*$ and $\mathbf{Y}_0$ is smallest, keeping in mind the measured frequency precisions $\Lambda$. One should accommodate two features. The first is that the variability in one component of $\theta$ may be larger than in another; in addition, the units may vary from component to component. This implies that standardization is necessary. The second feature is that correlations occur between the components of $\theta$ (see Figure 5 in Moravveji et al. 2015; Pedersen et al. 2018). The stronger a component correlates with the other components, the less independent information it contributes. An appropriate solution to this problem is offered by the Mahalanobis distance (Johnson & Wichern 2002):

$$ \theta_0 = \arg \min \{(\mathbf{Y}(\theta) - \mathbf{Y}^*)^\top \mathbf{V}^{-1} (\mathbf{Y}(\theta) - \mathbf{Y}^*)\}, \quad (1) $$

where $\mathbf{V} = \text{var}(\mathbf{Y})$, the variance-covariance matrix of the vector $\mathbf{Y}(\theta, \psi)$. The Mahalanobis distance naturally takes into account correlations among the components of $\theta$ and reduces to

the Euclidean distance in absence of any such correlations, i.e., when $\mathbf{V}$ is a diagonal matrix.

To evaluate Equation (1), we need an expression for $\mathbf{V}$. This matrix can be seen as the intrinsic variability in $\mathbf{Y}$ across a relevant range of $\mathbf{Y}$. The grid based on $\mathcal{M}(\theta, \psi)$, especially when sufficiently fine, will provide an approximation of this quantity when computed over a reasonable range of $\theta$. This implies that it can be estimated from the grid. Considering all grid points gives

$$ \hat{\mathbf{V}} = \frac{1}{q-1} \sum_{j=1}^{q} (\mathbf{Y}_j - \bar{\mathbf{Y}})(\mathbf{Y}_j - \bar{\mathbf{Y}})^\top, \quad (2) $$

where the average is

$$ \bar{\mathbf{Y}} = \frac{1}{q} \sum_{j=1}^{q} \mathbf{Y}_j. $$

In practice, we will typically evaluate Equation (1) using all grid points:

$$ \theta_0 = \arg \min_{j=1}^{q} \{(\mathbf{Y}_j - \mathbf{Y}^*)^\top \mathbf{V}^{-1} (\mathbf{Y}_j - \mathbf{Y}^*)\}. \quad (3) $$

Taking the measurement errors of the frequencies $\mathbf{Y}^*$ into account in the computation of the Mahalanobis distance is done by replacing $\hat{\mathbf{V}}$ by $\hat{\mathbf{V}} + \Lambda$ in Equation (3). Note that the inverse $\mathbf{V}^{-1}$ may need to be replaced with the Moore-Penrose inverse (Penrose 1955) for models that are not of full rank (see the statistical models treated in Problem 3). Expression (3) provides a more appropriate measure than the often-used $\chi^2$ based on a Euclidian distance, which does not take into account the correlated nature of the components of the vector $\theta$.

As for a first estimate of the error of $\theta_0$, we consider the Mahalanobis distance to comply with a normal distribution. In that case, it follows a $\chi^2$ distribution. Let the critical point be $c_p$, i.e., $P(\chi_p^2 \ge c_p) = 0.05$. We then consider the empirical prediction ellipsoid across a relevant part of the grid:

$$ C = \{\theta_j | j = 1, ..., q; (\mathbf{Y}_j(\theta_j, \psi) - \mathbf{Y}^*)^\top \hat{\mathbf{V}}^{-1} (\mathbf{Y}_j(\theta_j, \psi) - \mathbf{Y}^*) \le c_p\}. \quad (4) $$

The theoretical ellipsoid would range over all $\theta$ in continuous space, which is impractical, of course. As a consequence, the prediction ellipsoid depends on the coarseness of the grid, so it is essential to have it sufficiently fine yet covering an appropriately broad range for $\theta$. We can circumvent this below when treating Problem 3, where $\theta$ is determined from a statistical model rather than a theoretical model. In the case in which this can be done successfully, we improve the error estimation from Equation (4) by relying on the errors of the parameters of the statistical model to determine the errors of $\theta$.

As an illustration of the benefit of $\theta_0$ with respect to a $\chi^2$ statistic, we rely on two very dense seismic model grids computed to explain the 19 identified dipole zonal modes of the slowly pulsating B star KIC 10526294 by Moravveji et al. (see Table 1 in Moravveji et al. 2015, for the frequencies, and Table 3 for the grids and selection of the best models 4 and 10 from a $\chi^2$ statistic). Both these grids have $p = 5$, and their input physics only differs in the treatment of the core overshooting, with one grid using exponential overshoot described by the parameter $f_{ov}$ and the other grid relying on a
<table><caption>**Table 3**
Parameter Estimation Based on a $\chi^2$ Statistic (Moravveji et al. 2015) and on the Mahalanobis Distance $\theta_0$ from the 19 Dipole Zonal Modes of KIC 10526294</caption><thead><tr><th>Statistic</th><th><i>M</i> (<i>M</i><sub>☉</sub>)</th><th><i>Z</i></th><th><i>X</i><sub>c</sub></th><th><i>D</i><sub>ov</sub></th><th><i>D</i><sub>mix</sub> (cm<sup>2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>$\chi^2$</td><td>3.25</td><td>0.014</td><td>0.627</td><td><i>f</i><sub>ov</sub> = 0.017</td><td>56.2</td></tr><tr><td>$\theta_0$</td><td>3.24</td><td>0.017</td><td>0.644</td><td><i>f</i><sub>ov</sub> = 0.018</td><td>177.8</td></tr><tr><td>$\chi^2$</td><td>3.19</td><td>0.019</td><td>0.628</td><td>$\alpha_{ov} = 0.21$</td><td>56.2</td></tr><tr><td>$\theta_0$</td><td>3.23</td><td>0.017</td><td>0.640</td><td>$\alpha_{ov} = 0.20$</td><td>177.8</td></tr></tbody></table>

step overshoot with parameter $\alpha_{ov}$ (both have the radiative temperature gradient in the overshoot zone). This star rotates at less than $1\%v_{crit}$, and it revealed zonal modes with unambiguous identification from rotational splitting, so rotation can be ignored for the matching of the mode frequencies.

With our description based on the Mahalanobis distance in Equation (3), we improve the modeling compared to Moravveji et al. (2015) by taking into account the correlations among the parameters $\theta = (M, Z, X_c, D_{ov}, D_{ext})$, as well as allowing nonlinear dependence of the observed frequencies on those five parameters. The parameter estimation from the two grids based on $\chi^2$ and $\theta_0$ is listed in Table 3. The period spacing patterns of the best solutions are compared with the observed ones in Figure 12, while the two variance-covariance matrices for the two grids, $\hat{V} + \Lambda$ (hence taking into account the measured frequency errors), are shown in Figure 13. The observational errors decrease from $0.000057 \text{ day}^{-1}$ for $f_1 = 0.472220 \text{ day}^{-1}$ to $0.000036 \text{ day}^{-1}$ for $f_{19} = 1.013415 \text{ day}^{-1}$. Hence, almost all the variance represented in Figure 13 comes from the variability in the model grid predictions of the frequencies.

It is found that the estimated parameter values for $\theta$ are not exactly the same yet very similar in value. The $\chi^2$ selection by Moravveji et al. (2015) focused on the best agreement for the lowest-order modes, while the Mahalanobis distance selects models that give the overall best agreement. From the middle two panels of the figure, we deduce that the largest discrepancies between observed and theoretically predicted modes for the $\chi^2$ prediction occur for the highest-order modes and reach to about 900 s, while the Mahalanobis prediction has the poorest agreement for the lowest-order modes at a level of about 600 s. All differences between the frequencies of the four considered models and KIC 10526294 translate into $|f_i^* - f_i^{th}|/f_i^* < 1\%$. Moreover, the discrepancies between the model and observed frequencies in the bottom panel of Figure 12 are smaller than most of the theoretical uncertainties found in Sections 3 and 4 (see upper left insets of Figures 2–10). Hence, with the current knowledge of stellar interiors, this is an excellent overall agreement (we have stretched the ordinate axes in Figure 12 for optimal visibility of the small discrepancies).

Given that these two dense grids resulted from a zoom-in based on an initial screening of a much sparser one, optimally tuned to comply with the data, we expect the error estimation according to Equation (4) to encapsulate these entire dense grids. Their ranges are given in Table 2 of Moravveji et al. (2015, Mixing Grid and Step Overshoot Grid), and the distribution of the Mahalanobis distances is shown in Figure 14. The dashed vertical line is the $\chi_5^2 = 11.07$ cutoff value. About half of the models in the two grids have a lower distance, and as expected, the two prediction ellipsoids for $\theta$ span the entire grids in both cases. As a final note, we point out that the best models selected by the Mahalanobis distance do not change the inversion results for the rotation profile of the

star determined in Triana et al. (2015). Indeed, it was already shown in that paper that this profile was hardly dependent on the best forward model selection as long as it provided a reasonable match to the zonal mode frequencies.

## 5.2. Problem 2: A Mahalanobis-distance-based Solution to Compare Two Stellar Models

A pertinent question in seismic modeling of stars with a convective core is whether the core overshooting is better described by a step function or by an exponentially decaying one. Pedersen et al. (2018) assessed the capacity of the nominal Kepler data to achieve such a comparison based on gravity modes. Another case of comparing the capacity of two stellar models was already mentioned above, namely, with and without chemical transport in the radiative envelope (Moravveji et al. 2015; Pedersen et al. 2018). Such assessments require comparison between two theoretical models, $M(\theta^{(m)}, \psi^{(m)})$ and $N(\theta^{(n)}, \psi^{(n)})$. We assume, of course, that both models produce observation vectors $\mathbf{Y}$ and $\mathbf{Z}$, respectively, of the same length $n$ and having the same meaning for all of their components. That said, the length of the vector and the meaning of the components of $\theta^{(m)}$ and $\theta^{(n)}$ may be different.

A Mahalanobis-based solution, similar to that of Problem 1, is considered. It should be clear that there is some asymmetry to the problem, as each of the two grids has its own metric. Denote the grids by the following:

1. For $M(\theta^{(m)}, \psi^{(m)})$, we have resulting grid points $j = 1, \dots, q$, with $\mathbf{Y}_j(\theta_j^{(m)}, \psi^{(m)})$.

2. For $N(\theta^{(n)}, \psi^{(n)})$, we have resulting grid points $s = 1, \dots, r$, with $\mathbf{Z}_s(\theta_s^{(n)}, \psi^{(n)})$.

Let the variance-covariance matrix $V_m$, like in Equation (2), correspond to the first grid, with the equivalent quantity $V_n$ corresponding to the second grid. Suppose that a benchmark model computed from the first theory $M(\theta_0^{(m)}, \psi^{(m)})$ with parameters $\theta_0^{(m)}$ results in frequencies $\mathbf{Y}_0$; then we need to find

$$ \theta_0^{(n)} = \arg\min_{s=1}^{r} \{(\mathbf{Z}_s - \mathbf{Y}_0)^{\top} \mathbf{V}_n^{-1} (\mathbf{Z}_s - \mathbf{Y}_0)\}. \quad (5) $$

In the other direction, we pick a benchmark model computed from the second theory $N(\theta_0^{(n)}, \psi^{(n)})$ with parameters $\theta_0^{(n)}$ resulting in frequencies $\mathbf{Z}_0$, and we need to solve

$$ \theta_0^{(m)} = \arg\min_{j=1}^{q} \{(\mathbf{Y}_j - \mathbf{Z}_0)^{\top} \mathbf{V}_m^{-1} (\mathbf{Y}_j - \mathbf{Z}_0)\}. \quad (6) $$

As an illustration of this Problem 2, we rely on the 4D grids of stellar models computed by Pedersen et al. (2018, Table 2) and their high-order dipole modes with periods ranging from 0.8 to 3 days. For computational reasons, these grids were centered around two evolutionary stages ($X_c$ near 0.5 and 0.1). We consider the $n = 34$ frequencies of the dipole prograde gravity modes as $\mathbf{Y}$ to be those of the benchmark model from the step overshoot grid with $\theta_0^{(m)} = (M, X, X_c, \alpha_{ov}) = (3.25, 0.71, 0.5, 0.20)$, and we search for the best model in the grid computed for exponential overshooting with $f_{ov}$ instead of $\alpha_{ov}$ to ascertain which of the exponential overshoot models explains the one with the step overshooting of the benchmark model best. We refer to Pedersen et al. (2018) for the details on the fixed input physics $\psi$, which is the same for both grids except for the mixing efficiency in the overshoot region.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison between observed and theoretical period spacing patterns (top panel) and their difference (second panel), the corresponding mode period (third panel), and mode frequency (bottom panel) differences between the models and observations. The gray dashed lines indicate the observed frequencies listed in Table 1 of Moravveji et al. (2015). The errors of the measured mode periods are represented in Figure 13.</figcaption></figure>

We show the variance-covariance matrix $\hat{V}$ in the left panel of Figure 15, revealing once again that the dominant variability among the models in the grid occurs for modes in the high-frequency domain. The Mahalanobis distance leads to the parameter estimation given by $\theta_0^{(n)} = (M, X, X_c, f_{ov}) = (3.25, 0.71, 0.505, 0.015)$ from Equation (5). The distribution of the best solutions for the parameters $\theta$ is shown, along with the Mahalanobis distance distribution, in Figure 16 (red). The parameter distributions were constructed by giving every grid point $s = 1, \dots, r$ a weight according to the minimum Mahalanobis distance of the entire grid divided by the Mahalanobis distance at this grid point $s$.

Conversely, we select the dipole prograde gravity-mode frequencies $Z$ from a benchmark model with exponential

overshoot $\theta_0^{(n)} = (M, X, X_c, f_{ov}) = (3.25, 0.71, 0.5, 0.015)$ and search for the best model in the grid computed for step overshooting with parameter $\alpha_{ov}$ to find $\theta_0^{(m)} = (M, X, X_c, \alpha_{ov}) = (3.25, 0.71, 0.495, 0.20)$ from Equation (6) as represented in Figure 16 (blue), with the variance-covariance matrix $\hat{V}$ given in the right panel of Figure 15. It is visually clear that this problem set is asymmetric, as demonstrated by the asymmetry of the distributions of the parameters $\theta^{(n)}$ (red) and $\theta^{(m)}$ (blue). This reflects the fact that the parameters $\alpha_{ov}$ and $f_{ov}$ correlate differently with $(M, X, X_c)$. Such a different correlation structure was also hinted at from modeling of close binary companions with a convective core (Claret & Torres 2017) but was not yet taken up coherently in treating determinations of core masses in such objects from isochrone
<figure><img src="image_16.png" /><figcaption>Figure 13. Variance-covariance matrix $\hat{V} + \Lambda$ of dimension 19 × 19, representing the 19 detected zonal mode frequencies of KIC 10526294, for two considered dense grids of models, including the observational errors (left: for the exponential overshoot grid, whose best model is represented by the red circles in Figure 12; right: for the step overshoot grid, whose best model is indicated in blue in Figure 12).</figcaption></figure>

fitting. The bottom right panel of Figure 16 shows the power of the Mahalanobis distance to perform stellar aging from gravity modes, as there is no confusion between young and old models. Finally, the mass and initial hydrogen distributions in Figure 16 show the opportunity to determine stellar masses with high precision from gravity modes, typically better than 5%, while the chemical composition is harder to derive, as it suffers more from the degeneracies among the parameters.

### 5.3. Problem 3: Fitting a Statistical Model to Frequencies Generated from a Physical Model

Given that evaluating theoretical models and their frequencies is computationally expensive, one might aim to replace the $Y$ derived from a theoretical model $M(\theta, \psi)$ in the $q$ grid points with a statistical model, sufficiently close in outcome but easier to evaluate. Such a situation typically occurs when grids of stellar evolution and pulsation computations turn out to be insufficiently fine for forward seismic modeling and one wants to avoid having to recompute models with a too-extensive number of extra grid points (see Moravveji et al. 2016). Therefore, we write such a statistical model as

$$Y_j \sim N(\mu_j = \mu(\beta; \theta_j, \psi), \Sigma). \quad (7)$$

In other words, we treat the vectors with $n$ frequencies resulting from a theoretical model $M(\theta, \psi)$ as outcome variables $Y_j$, and we want to describe these outcome variables by a function $\mu_j = \mu(\beta, \theta_j, \psi)$. The latter can be linear or nonlinear. The statistical model will typically contain unknown parameters to be estimated, grouped in a vector $\beta$. Because a perfect match between the frequencies resulting from the theoretical stellar model $M(\theta, \psi)$ and from the statistical approximate model never occurs, we introduce an error term for each grid point, $\epsilon_j$, with variance-covariance matrix $\Sigma$. Hence, we rewrite the statistical model (7) as

$$Y_j = \mu(\beta; \theta_j, \psi) + \epsilon_j, \quad (8)$$

where $\epsilon_j$ is the error of the estimated frequency vector in the grid point $j$ resulting from the statistical model. At the level of

the individual frequencies labeled with $i=1, \dots, n$, Equation (8) can be written as

$$Y_{ji} = \mu_i(\beta_i; \theta_j, \psi) + \epsilon_{ji}, \quad (9)$$

with $\text{var}(\epsilon_{ji}) = \sigma_{ii}$, the $i$th diagonal element of $\Sigma$. Also, $\text{cov}(\epsilon_{ji}, \epsilon_{jk}) = \sigma_{ik}$ is the covariance, and $(i, k)$ is the off-diagonal element of $\Sigma$.<fnref n="9" /> So, in general, each regression component $i$ is allowed its own variance $\sigma_{ii}$ for $\epsilon_{ji}$, implying heteroscedasticity, as required following Sections 3 and 4.

With this setup, we assume that the normal distribution offers a sensible statistical working model to represent the theoretical models $M(\theta, \psi)$. These theoretical models lead to uncertainties for the theoretically predicted frequencies as illustrated in Figures 2–10. Several of these have symmetric distributions, while others have skewed distributions due to systematic uncertainties. Nevertheless, there are several arguments supporting the choice of a normal statistical theory for $Y_j$ in Equation (8). First, apart from computational convenience, this rests on the fact that, for inferences to be valid, first- and second-order moment assumptions suffice, even if normality is not satisfied (Johnson & Wichern 2002). This also explains why normal regression and linear regression lead to the same results. Second, White (1982) provides a theory for misspecified statistical models, indicating that fitting misspecified models with a normal distribution still leads to the best possible description given a class of statistical models considered. Third, transformations can be applied, if deemed appropriate, to improve the quality of the normal approximation. Fourth, the fact that $\mu(\beta; \theta_j, \psi)$ is allowed to be highly nonlinear brings appropriate flexibility to the model formulation. For all these reasons, the normal distribution is a very defensible model hypothesis to make (see also Gruberbauer et al. 2013; Appourchaux 2014, for similar arguments in the case of solar-like oscillations).

<footnote n="9">⁹ Note that we use the notation $\sigma_{ik}$ for the generic element of a variance-covariance matrix $\Sigma$, as is commonly done in the statistical literature. The elements correspond to variances if $i=k$ and to covariances if $i \neq k$ (Johnson & Wichern 2000). The corresponding standard deviations are then $\sqrt{\sigma_{ii}}$.</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 14. Histograms of the Mahalanobis distances where the overall range was divided over 1000 bins, for each of the two considered dense grids of models for KIC 10526294 (left: exponential overshoot; right: step overshoot). The vertical dashed line indicates a value of 11.07 ($P(\chi_5^2 \ge 11.07) = 0.05$), and the inset shows the distance over that interval.</figcaption></figure>

<figure><img src="image_18.png" /></figure>

### 5.3.1. Maximum Likelihood Estimation

One convenient way to estimate parameters is via maximum likelihood estimation (MLE; Welsh 1996). When each component of $Y_{ji}$ has its own regression parameters, the overall parameter vector partitions as $\beta = (\beta_1^\top, ..., \beta_n^\top)$, and we can separately solve each of the n problems. In this case, we can write $Y_{ji} = \mu_{ji}(\beta_i; \theta_j, \psi) + \varepsilon_{ji}$, where the length of the entire vector $\beta$ then partitions as $u = u_1 + ... + u_n$, with $u_i$ the length of the subvector pertaining to the regression for the ith frequency component. The normal likelihood for component i is then written as

$$L(\beta_i, \sigma_{ii}) = \prod_{j=1}^{q} \frac{1}{\sqrt{\sigma_{ii}} \sqrt{2\pi}} \exp \left\{ -\frac{1}{2} \frac{(Y_{ji} - \mu_{ji})^2}{\sigma_{ii}} \right\}, \quad (10)$$

where $\mu_{ji} \equiv \mu_{ji}(\beta_i; \theta_j, \psi)$ is used for shorthand. One proceeds by writing the corresponding log-likelihood:

$$\ell(\beta_i, \sigma_{ii}) = -\frac{q}{2} \ln(\sigma_{ii}) - \frac{q}{2} \ln(2\pi) - \frac{1}{2} \sum_{j=1}^{q} \frac{(y_{ji} - \mu_{ji})^2}{\sigma_{ii}},$$

where $y_{ji}$ stands for the realized value of the random variable $Y_{ji}$. In practice, it is sufficient to work with the kernel of the log-likelihood:

$$\tilde{\ell}(\beta_i, \sigma_{ii}) = -\frac{1}{2} \ln(\sigma_{ii}) - \frac{1}{2} \sum_{j=1}^{q} \frac{(y_{ji} - \mu_{ji})^2}{\sigma_{ii}}.$$

Its first derivative w.r.t. $\beta_i$, the so-called score function for $\beta_i$, is

$$S(\beta_i) = \frac{\partial \ell}{\partial \beta_i} = \frac{1}{\sigma_{ii}} \sum_{j=1}^{q} (y_{ji} - \mu_{ji}) \cdot \frac{\partial \mu_{ji}}{\partial \beta_i}. \quad (11)$$

The score equation to solve hence is

$$S(\beta_i) = \frac{\partial \ell}{\partial \beta_i} = \frac{1}{\sigma_{ii}} \sum_{j=1}^{q} (y_{ji} - \mu_{ji}) \frac{\partial \mu_{ji}}{\partial \beta_i} = 0.$$

For nonlinear models for $\mu_{ji}$, this first derivative will typically not allow for an explicit expression and may well be cumbersome (Seber & Wild 2003). Note that $S(\beta_i) = 0$ if and only if

$$\sum_{j=1}^{q} (y_{ji} - \mu_{ji}) \frac{\partial \mu_{ji}}{\partial \beta_i} = 0. \quad (12)$$

Evidently, the log-likelihood function also depends on the argument $\sigma_{ii}$. The score equation for this parameter is

$$S(\sigma_{ii}) = -\frac{1}{2} \frac{1}{\sigma_{ii}} + \frac{1}{2} \frac{1}{\sigma_{ii}^2} \sum_{j=1}^{q} (y_{ji} - \mu_{ji})^2 = 0. \quad (13)$$

Once Equation (12) has been solved, the solution to Equation (13) follows as

$$\hat{\sigma}_{ii} = \frac{1}{q} \sum_{j=1}^{q} [y_{ji} - \mu_{ji}(\hat{\beta}_i)]^2, \quad (14)$$

where $\hat{\beta}_i$ is the solution to Equation (12), the MLE for $\beta_i$.

Because of small-sample bias in the MLE for $\sigma_{ii}$, one often uses the least-squares version:

$$\hat{\sigma}_{ii} = \frac{1}{q - u_i} \sum_{j=1}^{q} [y_{ji} - \mu_{ji}(\hat{\beta}_i)]^2. \quad (15)$$

Clearly, when the number of grid points q is large relative to the number of parameters considered for optimization, $u_i$, as is the case in any practical application of forward statistical modeling, the bias becomes negligible.

To estimate the standard errors of $\theta$ and $\sigma^2$, we proceed by calculating the second derivatives of the log-likelihood
<figure><img src="image_19.png" /><figcaption>Figure 15. Variance-covariance matrix $\hat{V}$ of dimension 34 × 34, representing 34 prograde dipole mode frequencies of benchmark model $\theta_0^{(m)} = (M, X, X_c, \alpha_{ov}) = (3.25, 0.71, 0.5, 0.20)$ for a dense grid of models with exponential overshooting (left) and of benchmark model $\theta_0^{(n)} = (M, X, X_c, f_{ov}) = (3.25, 0.71, 0.5, 0.015)$ for a dense grid of models with step overshooting (right).</figcaption></figure>

function, the so-called Hessian matrix, denoted by $H(\cdot)$. The argument is left unspecified because we can calculate it for both $\theta$ and $\sigma^2$. The negative of the Hessian is the information matrix, denoted $I(\cdot)$. The inverse of this produces the variance-covariance matrix of the parameter estimators, denoted by $V(\cdot)$. We often work with the expected values of both $H(\cdot)$ and $I(\cdot)$, for ease of computation and because they are asymptotically equivalent. We will denote these expected values by $\mathcal{H}(\cdot)$ and $\mathcal{I}(\cdot)$. At this point, it is useful to note that, in fact, the entire variance-covariance matrix $\Sigma$ can be estimated, by vectorizing Equation (14) or Equation (15):

$$ \hat{\Sigma} = \frac{1}{q} \sum_{j=1}^{q} [\mathbf{Y}_j - \boldsymbol{\mu}_j(\hat{\boldsymbol{\beta}})][\mathbf{Y}_j - \boldsymbol{\mu}_j(\hat{\boldsymbol{\beta}})]^\top \quad (16) $$

and

$$ \tilde{\Sigma} = \frac{1}{q-p} \sum_{j=1}^{q} [\mathbf{Y}_j - \boldsymbol{\mu}_j(\hat{\boldsymbol{\beta}})][\mathbf{Y}_j - \boldsymbol{\mu}_j(\hat{\boldsymbol{\beta}})]^\top. \quad (17) $$

The Hessian for $\beta_i$ is

$$ H(\beta_i) = \frac{\partial^2 \ell}{\partial \beta_i \partial \beta_i^2} = \frac{\partial S(\beta_i)}{\partial \beta_i} = -\frac{1}{\sigma_{ii}} \sum_{j=1}^{q} \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right) \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right)^\top + \frac{1}{\sigma_{ii}} \sum_{j=1}^{q} (y_{ji} - \mu_{ji}) \frac{\partial^2 \mu_{ji}}{\partial \beta_i \partial \beta_i^\top}. $$

Hence,

$$ \mathcal{H}(\beta_i) = -\frac{1}{\sigma_{ii}} \sum_{j=1}^{q} \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right) \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right)^{\top} $$

because $E[y_{ji} - \mu_{ji}] = E[y_{ji}] - \mu_{ji} = \mu_{ji} - \mu_{ji} = 0$. Further, the variance-covariance matrix of $\beta_i$ takes the form

$$ V(\beta_i) = -\mathcal{H}(\beta_i)^{-1} = \left[ \frac{1}{\sigma^2} \sum_{j=1}^{q} \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right) \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right)^{\top} \right]^{-1}. $$

Turning to $\sigma_{ii}$,

$$ \begin{aligned} H(\sigma_{ii}) &= \frac{q}{2} \frac{1}{\sigma_{ii}^2} - \frac{1}{\sigma_{ii}^3} \sum_{j=1}^{q} (y_{ji} - \mu_{ji})^2, \\ \mathcal{H}(\sigma_{ii}) &= \frac{q}{2} \frac{1}{\sigma_{ii}^2} - \frac{1}{\sigma_{ii}^3} q \sigma_{ii} = -\frac{q}{2} \frac{1}{\sigma_{ii}^2}, \\ V(\hat{\sigma}_{ii}) &= \frac{2 \sigma_{ii}^2}{q}. \end{aligned} $$

From the variance-covariance matrices (note that for $\sigma_{ii}$ it is a scalar), we derive the standard errors. For $\beta_i$, we take the square root of the diagonal of $V(\beta_i)$. For $\sigma_{ii}$, we find s.e.($\sigma_{ii}$) = $\sqrt{2/q} \sigma_{ii}$. Strictly speaking, we should also consider the covariance between the estimated $\beta_i$ and $\sigma_{ii}$. However, this is asymptotically zero. Indeed,

$$ \frac{\partial S(\beta_i)}{\partial \sigma_{ii}} = \frac{\partial S(\sigma_{ii})}{\partial \beta_i} = H(\beta_i, \sigma_{ii}) = -\frac{1}{\sigma_{ii}^2} \sum_{j=1}^{q} (y_{ji} - \mu_{ji}) \frac{\partial \mu_{ji}}{\partial \beta_i}. $$

Hence,

$$ \mathcal{H}(\beta_i, \sigma_{ii}) = -\frac{1}{\sigma_{ii}^2} \sum_{j=1}^{q} E(y_{ji} - \mu_{ji}) \frac{\partial \mu_{ji}}{\partial \beta_i} = 0. $$

### Special Case: Linear Statistical Models

In linear models, all of the above estimators and precision estimators have closed forms. Write the linear model as $\mu_{ji} = \mathbf{x}_{ji}^\top \beta_i$, where $\mathbf{x}_{ji}$ is a vector of predictors, consisting of the components of $\theta_j$ and functions thereof, and create a convenient design matrix and outcome vector:

$$ X_i = \begin{pmatrix} \mathbf{x}_{1i}^\top \\ \mathbf{x}_{2i}^\top \\ \vdots \\ \mathbf{x}_{qi}^\top \end{pmatrix}, \quad Y_i = \begin{pmatrix} y_{1i}^\top \\ y_{2i}^\top \\ \vdots \\ y_{qi}^\top \end{pmatrix}, $$
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Distributions of the Mahalanobis distance and of the four parameter predictions for $\theta$ resulting from the comparison of 34 zonal dipole mode frequencies of a benchmark model with $(M, X, X_c, \alpha_{ov}) = (3.25, 0.71, 0.5, 0.20)$ as found from a grid of models with exponential overshoot (red) and of a benchmark model with $(M, X, X_c, f_{ov}) = (3.25, 0.71, 0.5, 0.015)$ compared to a grid of step overshoot models (blue).</figcaption></figure>

$$ 
\begin{aligned}
S(\beta_i) = 0 &\iff \sum_{j=1}^{q} (y_{ji} - \mathbf{x}_{ji}^{\top} \beta_i) \mathbf{x}_{ji} = 0 \\
&\iff \sum_{j=1}^{q} \mathbf{x}_{ji} y_{ji} = \left( \sum_{j=1}^{q} \mathbf{x}_{ji} \mathbf{x}_{ji}^{\top} \right) \beta_i \\
&\iff \mathbf{x}_{i}^{\top} Y_i = (\mathbf{x}_{i}^{\top} \mathbf{X}_i) \beta_i.
\end{aligned}
$$

This leads to the classical ordinary least-squares (OLS) solution for linear regression (Neter et al. 1990):

$$ \beta_i = (X_i^\top X_i)^{-1} X_i^\top Y_i. \qquad (18) $$

The expected Hessian becomes

$$ \mathcal{H}(\beta_i) = -\frac{1}{\sigma_{ii}} \sum_{j=1}^{q} \mathbf{x}_{ji} \mathbf{x}_{ji}^{\top} = -\frac{1}{\sigma_{ii}} X_i^{\top} X_i, $$

and hence

$$ V(\beta_i) = \sigma_{ii} (X_i^\top X_i)^{-1}. \qquad (19) $$

For $\sigma_{ii}$, all remains the same. Note that the space spanned by $\mu(\hat{\beta}; \theta_j, \psi)$ is at most of dimension length $(\theta_j) = t$, even when $t < n$.

### 5.3.2. Including a Measurement Error Component

In the above, Model (8) is heteroscedastic in the sense that it allows for a different variance $\sigma_{ii}$ for each of the components $Y_{ji}$ of the vector $Y_j$. In concrete applications, one would like to assess how the difference between the $Y$ vector derived from the statistical model and that derived from the “true” physical model $\mathcal{M}(\theta, \psi)$ compares to the measured frequency precision that applies to $Y^*$. Indeed, in practice, the various grids will have been computed for a particular star, whose gravity-mode oscillations and their errors have been measured for its $i = 1, \dots, n$ frequencies.

We recall that the measured frequency precisions are denoted as $\Lambda_j$ and that its diagonal matrix elements are $\lambda_{j,ii}$. Taking measurement errors into account is done by replacing the variance–covariance matrix $\Sigma$ by $\Sigma + \Lambda_j$. The regression equation then becomes

$$ Y_{ji} = \mu_{ji}(\beta_i) + \varepsilon_{ji} + \varepsilon_{ji}^*, \qquad (20) $$

with $\varepsilon_{ji} \sim N(0, \sigma_{ii})$, $\varepsilon_j^* \sim N(0, \lambda_{ji})$, and $\lambda_{ji}$ known. Note that, for generality, we allow the measurement error to depend on the component within the vector $Y$, as well as the actual grid point $j$. Alternatively, $j$ could refer to a star when the data set is made up of a collection of stars (see Problem 4). The likelihood
then becomes

$$L(\beta_i, \sigma_{ii}) = \prod_{j=1}^{q} \frac{1}{\sqrt{\sigma_{ii} + \lambda_{ji} \sqrt{2\pi}}} \exp \left\{ -\frac{1}{2} \frac{(y_{ji} - \mu_{ji})^2}{\sigma_{ii} + \lambda_{ji}} \right\}. \quad (21)$$

The kernel of the log-likelihood becomes

$$\tilde{\ell}(\beta_i, \sigma_{ii}) = -\frac{1}{2} \sum_{j=1}^{q} \ln(\sigma_{ii} + \lambda_{ji}) - \frac{1}{2} \sum_{j=1}^{q} \frac{(y_{ji} - \mu_{ji})^2}{\sigma_{ii} + \lambda_{ji}}.$$

The corresponding score equations are

$$S(\beta_i) = \sum_{j=1}^{q} \frac{(y_{ji} - \mu_{ji})}{\sigma_{ii} + \lambda_{ji}} \cdot \frac{\partial \mu_{ji}}{\partial \beta_i}, \quad (22)$$

$$S(\sigma_{ii}) = -\frac{1}{2} \sum_{j=1}^{q} \frac{1}{\sigma_{ii} + \lambda_{ji}} + \frac{1}{2} \sum_{j=1}^{q} \frac{(y_{ji} - \mu_{ji})^2}{(\sigma_{ii} + \lambda_{ji})^2}. \quad (23)$$

Expressions for Hessians, expected Hessians, and variance-covariance matrices are

$$\begin{align*}
H(\beta_i) &= -\sum_{j=1}^{q} \frac{1}{\sigma_{ii} + \lambda_{ji}} \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right) \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right)^{\top} \\
&\quad + \sum_{j=1}^{q} \frac{(y_{ji} - \mu_{ji})}{\sigma_{ii} + \lambda_{ji}} \frac{\partial^2 \mu_{ji}}{\partial \beta_i \partial \beta_i^{\top}}, \\
\mathcal{H}(\beta_i) &= -\sum_{j=1}^{q} \frac{1}{\sigma_{ii} + \lambda_{ji}} \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right) \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right)^{\top}, \\
V(\beta_i) &= \left[ \sum_{j=1}^{q} \frac{1}{\sigma_{ii} + \lambda_{ji}} \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right) \left( \frac{\partial \mu_{ji}}{\partial \beta_i} \right)^{\top} \right]^{-1}, \\
H(\sigma_{ii}) &= \frac{1}{2} \sum_{j=1}^{q} \frac{1}{(\sigma_{ii} + \lambda_{ji})^2} - \sum_{j=1}^{q} \frac{(y_{ji} - \mu_{ji})^2}{(\sigma_{ii} + \lambda_{ji})^3}, \\
\mathcal{H}(\sigma_{ii}) &= \frac{1}{2} \sum_{j=1}^{q} \frac{1}{(\sigma_{ii} + \lambda_{ji})^2} - \sum_{j=1}^{q} \frac{1}{(\sigma_{ii} + \lambda_{ji})^2} \\
&= -\frac{1}{2} \sum_{j=1}^{q} \frac{1}{(\sigma_{ii} + \lambda_{ji})^2}, \\
V(\sigma_{ii}) &= \left[ \frac{1}{2} \sum_{j=1}^{q} \frac{1}{(\sigma_{ii} + \lambda_{ji})^2} \right]^{-1}.
\end{align*}$$

Also here, the asymptotic covariance between $\beta_i$ and $\sigma_{ii}$ is equal to zero.

In this case, the score equations corresponding to Equations (22)–(23) cannot be solved separately. One can therefore iterate between the two, i.e., first solve $S(\beta_i) = 0$ with $\sigma_{ii}$ fixed and then $S(\sigma_{ii}) = 0$ with $\beta_i$ fixed. These two steps are then iterated until convergence. To achieve this, one can use the iterated profile likelihood method, in which case $\sigma_{ii}$ can be considered just another component of the parameter vector: $(\beta_{i1}, \dots, \beta_{iu}, \sigma_{ii})^\top$.

### Special Case: Linear Model

Also here, simplification is possible for linear models. We assume the same linear mean structure as in Section 5.3.1, with the same design matrices $X_i$ and $Y_i$, and we additionally create a diagonal weight matrix $W_i$, with elements $w_{ji} = (\sigma_{ii} + \lambda_{ji})^{-1}$

along the diagonal. The score equation for $\beta_i$ takes the form

$$\begin{align*}
S(\beta_i) = 0 &\iff \sum_{j=1}^{q} w_{ji} (y_{ji} - \mathbf{x}_{ji}^{\top} \beta_i) \mathbf{x}_{ji} = 0 \\
&\iff \sum_{j=1}^{q} w_{ji} \mathbf{x}_{ji} y_{ji} = \left( \sum_{j=1}^{q} \mathbf{x}_{ji} w_{ji} \mathbf{x}_{ji}^{\top} \right) \beta_i \\
&\iff \mathbf{x}_{i}^{\top} W_i Y_i = (\mathbf{x}_{i}^{\top} W_i X_i) \beta_i \\
&\iff \beta_i = (\mathbf{x}_{i}^{\top} W_i X_i)^{-1} \mathbf{x}_{i}^{\top} W_i Y_i.
\end{align*}$$

Here, $\beta_i$ depends on $\sigma_{ii}$ through the weights $w_{ji}$. This means that the above expression produces the MLE estimator only after the MLE for $\sigma_{ii}$ is found. Both optimizations can be coupled in the following iterative scheme:

1. Update the variance: $\sigma_{ii}^{(\text{new})} = \sigma_{ii}^{(\text{old})} - \mathcal{H}^{-1}(\sigma_{ii}^{(\text{old})}) S(\sigma_{ii}^{(\text{old})})$.

2. With this new $\sigma_{ii}^{(\text{new})}$, the weights can be recomputed.

3. These weights are then used to update $\beta_i$ using $\beta_i^{(\text{new})} = (\mathbf{x}_{i}^{\top} W_i X_i)^{-1} \mathbf{x}_{i}^{\top} W_i Y_i$.

The above triplet of steps is repeated until convergence, at which time the variance can be estimated from

$$\begin{align*}
V(\beta_i) &= (\mathbf{x}_{i}^{\top} W_i X_i)^{-1}, \\
V(\sigma_{ii}) &= 2 \left( \sum_{j=1}^{q} w_{ji}^2 \right)^{-1}.
\end{align*}$$

In the case where linear statistical models do not offer a proper fit to mimic the theoretical models $\mathcal{M}(\theta, \psi)$, Appendix B offers a practical guide to try nonlinear models by means of fractional polynomials.

### 5.3.3. Statistical Model Selection

For the overall set of $n$ frequencies per star, we may want to choose between various candidate statistical models, where each $\beta_i$ is of dimension $u_i$ for $i = 1, \dots, n$. Let us assume we have two nested statistical models and we want to compare them: a “null” model, $H_0$, where each of the estimated parameters $\beta_{0,i}$ has dimension $u_{0,i}$ and a more complex alternative model described by $\beta_i$, each with dimension $u_i > u_{0,i}$. Hence, the parameters $\beta_{0,i}$ constitute a subspace of that of the more complex model described by $\beta_i$. In the case of such nested models, the comparison of their appropriateness can be done using the classical likelihood-based test statistics: likelihood ratio, score, and Wald test statistics with expressions given by, respectively:

$$T_{LR} = 2 \left[ \sum_{i=1}^{n} \ell(\hat{\beta}_i) - \sum_{i=1}^{n} \ell(\beta_{0,i}) \right], \quad (24)$$

$$T_S = \sum_{i=1}^{n} \left[ \frac{\partial \ell(\beta_i)}{\partial \beta_i} \bigg|_{\beta_i = \beta_{0,i}} \right]^{\top} \left[ - \frac{\partial^2 \ell(\beta_i)}{\partial \beta_i \partial \beta_i} \bigg|_{\beta_i = \beta_{0,i}} \right]^{-1} \\
\times \left[ \frac{\partial \ell(\beta_i)}{\partial \beta_i} \bigg|_{\beta_i = \beta_{0,i}} \right], \quad (25)$$

$$T_W = \sum_{i=1}^{n} [\hat{\beta}_i - \beta_{0,i}]^{\top} \left[ - \frac{\partial^2 \ell(\beta_i)}{\partial \beta_i \partial \beta_i} \bigg|_{\beta_i = \hat{\beta}_i} \right] [\hat{\beta}_i - \beta_{0,i}]. \quad (26)$$
Classical likelihood theory implies that, under $H_0$, $T_{LR}$, $T_S$, and $T_W$ are asymptotically equivalent and $\chi^2$ distributed, with $z = \sum_{i=1}^n(u_i - u_{0,i})$ the overall difference in the total number of degrees of freedom of the model $H_0$ and the alternative more complex model (Cox & Hinkley 1990). To reduce complexity, the second-derivative matrix occurring in Equations (25) and (26) can be replaced by $-\mathcal{H}(\beta_i)$.

A model selection procedure such as that based on the likelihood ratio test statistic suffers from two limitations. First, it can be used only with nested models (i.e., when one model is a submodel of the other). Second, there is the known tendency to favor the more complex models more than should be the case. Various information criteria, such as Akaike's and Bayesian information criteria (AIC and BIC, respectively, but there are several others), have been formulated to this effect. Several such criteria are in themselves functions of the likelihood ratio test statistic, often being equal to it, but augmented with a term that depends on the difference in the number of model parameters between both models being compared. The AIC is defined as

$$ \text{AIC} \equiv 2 \sum_{i=1}^{n} u_i - 2 \sum_{i=1}^{n} \ell(\hat{\beta}_i) \quad (27) $$

(Akaike 1987), and the BIC as

$$ \text{BIC} \equiv \left( \sum_{i=1}^{n} u_i \right) \ln(n) - 2 \sum_{i=1}^{n} \ell(\hat{\beta}_i) \quad (28) $$

(Schwarz 1978). The latter is a particular type of approximation to Bayesian posterior probabilities for statistical model selection, which penalizes stronger for complexity than the AIC (Claeskens & Hjort 2008). Both the AIC and BIC are only relevant for relative comparisons, where a lower value implies a better statistical model. The AIC and BIC offer an alternative to the likelihood ratio, score, and Wald test statistics discussed above, whether the statistical models are nested or not.

A drawback for the use of AIC and BIC is that there is no formal (distribution) theory that can be used to select a model. Therefore, one often resorts to rules of thumb, such as that based on the difference between the AIC or BIC values of two models A and B. A guide for values of $\Delta$AIC $\equiv$ AIC(A)-AIC (B) or $\Delta$BIC $\equiv$ BIC(A)-BIC(B) is as follows: [0, 2] for weak, [2, 6] for positive, [6, 10] for strong, and [10, $+\infty$] for very strong evidence in favor of model B (Molenberghs & Verbeke 2005, Section 4). A related method, Bayes factors, was described by Kass & Raftery (1995).

### 5.3.4. Error Estimation Based on Statistical Models

Once the optimal statistical model has been found to replace the “true” model $M(\theta, \psi)$, we need to assess the error of the corresponding $\theta$ prediction based on the Mahalanobis distance as in Problem 1. In this case, the uncertainty of $\beta$ needs to be taken into account. Assume that $\beta$ and its variance-covariance matrix have been estimated from MLE, say, $V \equiv V(\beta)$. It is easiest to proceed in a Monte Carlo way. First, sample $R$ values $\beta_r \sim N(\hat{\beta}, V)$, $r = 1, ..., R$. Then, apply Equation (4) for the statistical model with parameter $\beta_r$ separately, leading to $C_r$.

Finally, consider

$$ C = \bigcup_{r=1}^{R} C_r. $$

Generating a copy $\beta_r \sim N(\hat{\beta}, V)$ can easily be done in the following way:

1. Generate $\gamma_r \sim N(0_p, I_p)$, where $0_p$ is a column vector of zeros of length $p$ and $I_p$ is the $p \times p$ identity matrix.

2. Decompose $V = LL^\top$, where $L$ is the (lower triangular) Cholesky decomposition of $V$ (Nash 1990).

3. Define $\beta_r = \hat{\beta} + L\gamma_r$.

### 5.3.5. Application

As an illustration of Problem 3, we once more rely on the 4D grids of stellar models by Pedersen et al. (2018, Table 2), already used in the previous section (see Figure 16). We consider the $n = 34$ frequencies of the dipole prograde gravity-mode frequencies of the benchmark model with parameters $(M, X, X_c, f_{ov}) = (3.25, 0.71, 0.5, 0.015)$ as an “observed” star with $Y^*$, and we search for the best model in the grid computed for step overshoot described by $\alpha_{ov}$ that minimizes the difference between $Y$ and $Y^*$, where we now search for statistical models to replace the values of the 34 frequencies (omitting this time the grid models with $X_c \sim 0.1$). For this illustration, we constructed two linear nested models according to Equation (9) for $i = 1, ..., 34$ with the following:

1. Statistical model 1:

$$ \theta = (\theta_1, \theta_2, \theta_3, \theta_4) = (M, X, X_c, \alpha_{ov}); $$

2. Statistical model 2:

$$ \theta = (\theta_1, \theta_2, \theta_3, \theta_4, \theta_5, \theta_6, \theta_7, \theta_8, \theta_9, \theta_{10}) = (M, X, X_c, \alpha_{ov}, M \cdot X, M \cdot X_c, M \cdot \alpha_{ov}, X \cdot X_c, X \cdot \alpha_{ov}, X_c \cdot \alpha_{ov}), $$

where we perform OLS to estimate $\beta$ via Equations (18) and (19). For each of the 34 frequencies, this leads to a $\beta$-vector of dimension 5 in the case of Statistical model 1 (an intercept estimate and four regression coefficients for each of the four components of $\theta$) and of dimension 11 for Statistical model 2. The choice for Statistical model 2 was taken on the basis of correlations between the parameters of the vector $\theta$ in previous seismic modeling applications (e.g., Moravveji et al. 2015, Figure 5). The results from the step overshoot grid itself and from the two linear statistical models are shown in Figures 17 and 18. The accompanying parameter estimates obtained from the two statistical models are as follows:

1. Statistical model 1: $\hat{M} = 3.20_{3.15}^{0.25}M_\odot$, $\hat{X} = 0.69_{0.68}^{0.70}$, $\hat{X}_c = 0.490_{0.485}^{0.515}$, $\hat{\alpha}_{ov} = 0.14_{0.12}^{0.21}$;

2. Statistical model 2: $\hat{M} = 3.25_{3.15}^{0.30}M_\odot$, $\hat{X} = 0.72_{0.68}^{0.73}$, $\hat{X}_c = 0.512_{0.490}^{0.525}$, $\hat{\alpha}_{ov} = 0.14_{0.10}^{0.22}$

where the super- and subscripts denote the ranges of the parameters by performing the error determination using the Monte Carlo method discussed in Section 5.3.4. It can be seen that both statistical models represent the 34 “observed” frequencies from the benchmark model very well and in both cases better than the corresponding best set of 34 frequencies retrieved from the grid points of the step overshoot grid. The standard deviations of the 34 frequencies estimated from Equation (14) are in the range $\sqrt{\sigma_{ii}} \in [0.0022, 0.0281] \text{ day}^{-1}$ for model 1 and $[0.0022, 0.0264] \text{ day}^{-1}$ for model 2. This application shows the benefit of replacing a model grid by
<figure><img src="image_21.png" /><figcaption>**Figure 17.** Comparison between “observed” period spacing patterns from a benchmark model with ($M, X, X_c, f_{ov}$) = (3.25, 0.71, 0.5, 0.015) and their estimated values (top panel) from a grid based on step overshoot models of similar age and from two linear statistical models deduced from this step overshoot grid. The differences of the period spacing (second panel), of the corresponding mode periods (third panel), and of the mode frequencies (bottom panel) are also shown.</figcaption></figure>

statistical models based on such a grid and application of MLE, because it allows a refined parameter estimation without having to compute an infinitely fine grid.

Even though our formalism is suitable for treating nonlinear statistical models, applications to gravity-mode oscillations of main-sequence stars will most of the time work well for linear statistical models whenever the dimensions of the physical quantities in the vector $\theta$ are appropriately chosen. In practice, one works best with $M$ expressed in solar masses, $X$, $X_c$, and $Z$ in percentages, $f_{ov}$ or $\alpha_{ov}$ in local pressure scale height, $\log D_{\text{ext}}$ with $D_{\text{ext}}$ expressed in cm$^2$ s$^{-1}$, and the frequencies of rotation and of the oscillations in day$^{-1}$. In that case, it is well known that deviations from linearity of the dependence of the frequencies on the stellar parameters $\theta$ remain modest, except for the more rapid phases near the TAMS. Early illustrations of

these dependencies can be found in Figure 2 of Ausseloos et al. (2004) for pressure modes and in Figure 5 of Moravveji et al. (2015) for gravity modes of B-type pulsators.

As for the statistical model selection, we compare the two nested linear models by means of the various tests described above. For the likelihood ratio test, we get $T_{LR} = 0.24$. Comparing this with the 5% cutoff of a $\chi^2_{204}$ (given that $z = 34 \cdot (11 - 5)$) points out that Statistical model 2 is better than Statistical model 1, in line with Figure 17. This is not surprising, as the second model has many more degrees of freedom. Given that fact, it is better to make the model comparison by means of the AIC and BIC, according to Equations (27) and (28), because these penalize for the complexity of the model. We obtain the following values: AIC(model 1) = -6.29 and BIC(model 1) = 253.20, while
<figure><img src="image_22.png" /><figcaption>Figure 18. Distributions of the Mahalanobis distance and of the four parameter estimations for $\theta$ resulting from the comparison of 34 prograde dipole mode frequencies of a benchmark model with ($M$, $X$, $X_c$, $f_{ov}$) = (3.25, 0.71, 0.5, 0.015) from a grid based on step overshoot models of similar age and for two linear statistical models deduced from this step overshoot grid.</figcaption></figure>

AIC(model 2) = 401.47 and BIC(model 2) = 972.33. Hence, $\Delta$AIC = 407.73 and $\Delta$BIC = 719.13. We thus conclude to have found very strong evidence that model 1 outperforms model 2 according to both the AIC and BIC, i.e., the improvement for the estimation of the frequencies in the grid points achieved by the second model is insufficient to justify the more complex model with the extra degrees of freedom.

The astrophysical conclusion of this application (see Figure 18) is that the mass and age of a star can be well estimated from MLE based on model grids at the level of variance considered in this example, even if the shape of the core overshooting is unknown. The resulting uncertainty of the overshoot parameter is relatively large owing to its correlated nature with the global stellar mass and age and with other aspects of the input physics. Note that the latter correlation structure is not properly taken into account when estimating overshooting from binary or cluster modeling in the literature so far.

## 5.4. Problem 4: Stellar Model Selection from an Ensemble of Pulsators

Suppose that we have not one star as assumed so far, but rather an ensemble of $N$ stars, as well as a collection of candidate theories:

1. Observations $Y_t^*$, $t = 1, ..., N$;

2. Theories $M(\theta^{(m)}, \psi^{(m)})$, $m = 1, ..., M$.

The goal is, using all stars simultaneously, to select the most appropriate theory among the $M$ candidate theories. Using the solution to Problem 1, i.e., applying Equation (3) to every star $t$ for every grid $m$, we find the pair {$Y_{mt}, \theta_{mt}$}, with $\theta_{mt}$ the fitted stellar properties for star $t$ from theory $m$, and $Y_{mt}$ the predicted observation for star $t$ from theory $m$.

The discrepancy between observed values and values predicted from theory $m$ is

$$D_m = \sum_{t=1}^{N} D_{mt} = \sum_{t=1}^{N} \{ (Y_t^* - Y_{mt})^\top V_m^{-1} (Y_t^* - Y_{mt}) \},$$

with $V_m$ the dispersion matrix, calculated as in Equation (2), for the $m$th theory. Here, the stellar models can take the form of an astrophysical model, as well as that of a statistical model built according to the principles laid out in Problem 3.

As argued above, it is again sensible to take the working view that $Y_t^* \sim N(Y_{mt}, V_m)$. Let $t_m$ be the length of vector $\theta^{(m)}$. Then $D_{mt}$ follows a $\chi^2_{n-t_m}$ distribution, and hence $D_m$ as a whole follows a $\chi^2_{N(n-t_m)}$ distribution. Using these distributions, the tail probability $p(M) = P(\chi^2_{N(n-t_m)} \ge D_m)$ can be calculated. The theory with the largest $p(M)$ fits the observed data of the ensemble of stars best.

We are not yet able to include an example of such an application, since we first must apply the above set of Problems 1, 2, and 3 to a sufficient number of observed gravity-mode pulsators. This will be the subject of subsequent
papers with applications of the current methodological
framework.

5.5. Bayesian Approach for Parameter Estimation

Statistical inference in asteroseismology is nowadays often
done with a Bayesian approach, where the likelihood function
based on data is combined with prior knowledge to update the
information on the parameters θ and to revise the probability
associated with the parameters (posterior probability). In that
sense, an MLE is related to Bayesian analysis with a so-called
noninformative (or flat) prior. We refer to Gruberbauer et al.
(2012) and Appourchaux (2014) for thorough discussions in
the case of solar-like oscillations in low-mass stars, where the
input physics of the model M(θ, ψ) is usually taken to be very
similar to the one of solar models calibrated from helioseis-
mology. Indeed, for the “easy” case of low-mass stars with
sun-like oscillations, one can make the reasonable assumption
that such stars have input physics very similar to the Sun,
because they are slow rotators and have a convective envelope.

In the context of gravity modes and the results obtained in
Sections 2–4, it is essential to start any Bayesian analysis with a
noninformative prior about the input physics and parameters to
be estimated. Given the results on theoretical uncertainties of
gravity-mode frequencies of rotating stars with a convective
core, as summarized in Table 2, knowledge of the choice of θ
to be estimated and of the input physics in the model M(θ, ψ)
to be fixed for the stellar models has yet to be assembled for
large samples of stars with a variety in mass, rotation,
metallicity, and binary properties. As long as we do not know
the level and shape of the interior profiles for core over-
shooting, rotation, and chemical mixing, one should not impose
insecure assumptions about these quantities, because they have
a large effect on the gravity modes. For this reason, a Bayesian
approach is not obvious yet for gravity-mode pulsators with a
convective core, because we essentially do not have any
objective prior knowledge on the choice and properties of
M(θ, ψ). In this case, one should not introduce biases about
the ingredients of the stellar models (see Table 2), but rather
complete Problem 4 for samples of tens to hundreds of
stars first.

However, using prior knowledge in the forward seismic
modeling of gravity-mode pulsators may be applied at the level
of identification of the degree *l* and/or the azimuthal order *m* of
the gravity modes. In all of the above, we have assumed that
each observed frequency *Y*<sub>*i*</sub><sup>*</sup> = *f*<sub>*i*</sub><sup>*</sup> has known (*l*, *m*). The
occurrence of rotational splitting indeed implies unambiguous
labeling of the mode wavenumbers (*l*, *m*) (e.g., Kurtz
et al. 2014; Moravveji et al. 2015; Saio et al. 2015), but this
may not necessarily be the case when only period spacing
patterns are available for mode identification (e.g., Van Reeth
et al. 2016; Ouazzani et al. 2017; Pápics et al. 2017; Saio
et al. 2018b). In case of ambiguity (e.g., Degroote et al. 2010;
Pápics et al. 2012; Buysschaert et al. 2018), several options for
(*l*, *m*) should be kept open, just as for the case of solar-like
oscillations with doubtful identifications (see Appourchaux
2014, for an enlightening discussion). In this sense, model
selection is ubiquitous in asteroseismology, even for the “well-
known” case of solar-like oscillations, not only due to lack of
secure mode identification (see the F-type pulsator HD 49333
as in Benomar et al. 2009; Gruberbauer et al. 2009) but also
given problems in mode detection (e.g., Deheuvels et al. 2010;
Corsaro & De Ridder 2014; Davies et al. 2016). For heat-

driven coherent modes, identification of *l* and/or *m* can
sometimes be achieved from time-resolved long-term multicolor
photometry and/or line-profile variations (see Aerts et al. 2010,
Chapter 6, for an extensive discussion and methodology to
achieve this). In case of successful empirical identification of
(*l*, *m*), this information can be used as an informative prior. Lack
of identification may also occur at the level of a unique value for
the radial orders *n*<sub>pg,i</sub> of each *f*<sub>*i*</sub><sup>*</sup>. Indeed, this labeling may be
ambiguous for the regime of high-order gravity modes, given
the density of the mode spectrum at low frequencies. One may
then want to consider a Bayesian approach instead of MLE, so
as to incorporate prior probabilities for the *θ* considering
different priors of (*l*, *m*, *n*) for each of the modes. The ensuing
estimation method is known as *maximum a posteriori
estimation*.

Given that the posterior involves the likelihood function as
well, its evaluation is commonly approached in a simulation-based
fashion, using Markov chain Monte Carlo (MCMC) methods
(Louis & Carlin 2009). Care has to be taken in order not to miss
the best solution of the score Equations (22) and (23). An optimal
approach is to compare MLE and Bayesian treatments. Such
comparisons in asteroseismology have been done quite exten-
sively before (e.g., Appourchaux et al. 2012a, 2012b), but only in
the case of solar-like oscillations, where M(θ, ψ) is far better
known and the effect of rotation can be ignored or treated in the
Ledoux approximation (e.g., Appourchaux et al. 2008, 2014;
Benomar et al. 2009; Deheuvels et al. 2012; Gruberbauer
et al. 2013; Lagarde et al. 2016; Bossini et al. 2017; Handberg
et al. 2017; Rodrigues et al. 2017).

Investigating MLE versus Bayesian MCMC has not been
done in the context of forward seismic modeling of gravity
modes. For applications to intermediate-mass or high-mass
stars as treated here, interior rotation and convective core
overshooting are crucial and cannot be ignored in the parameter
estimation because the gravity-mode frequency values depend
heavily on it. The same is true, although to a lesser extent, for
chemical mixing in the µ-gradient zone near the receding core,
as this determines the mode trapping in that region. This
fortunate circumstance in principle allows the properties of the
core overshooting and of the near-core physics from gravity
modes to be estimated, in a similar way that mode trapping
properties allow the chemical profiles of white dwarf stars to be
deduced (Giammichele et al. 2018). However, it does imply
that we must first model an ensemble of stars without biasing
prior knowledge on these uncalibrated physical phenomena.

An extensive grid-based MLE and stellar model selection
procedure based on a χ²-type analysis for gravity modes and
with estimation of the core overshooting parameter was applied
to the B-type gravity-mode oscillator KIC 10526294 by
Moravveji et al. (2015). This MLE and model selection
approach was then followed by rotation profile estimation from
inversion by Triana et al. (2015), who used the AIC to deduce
the most likely rotation profile. This combination of MLE for
forward seismic modeling and stellar model selection from a
grid-based approach, followed by statistical model selection to
deduce the optimal rotation profile, is the only application of its
kind so far, given that Kurtz et al. (2014), Saio et al. (2015),
Kallinger et al. (2017), and Szewczuk & Daszyńska-Daszkiewicz
(2018) performed the modeling with simplistic parameter
estimation, if at all, without considering convective overshooting
as a free parameter.
**6. Overall Forward Modeling Scheme for Gravity Modes in Rotating Stars**

Forward seismic modeling requires observed and identified pulsation mode frequencies and fits them with theoretically predicted frequencies deduced from stellar models. As shown in Figure 2, it is essential to take into account the Coriolis force in the computations of the theoretical gravity-mode frequencies when performing forward seismic modeling of intermediate- and high-mass stars, even for moderate rotators. We summarize here our overall modeling scheme, including the Coriolis force, which has so far been ignored in the modeling, with the exception of the studies in Moravveji et al. (2016) and Schmid & Aerts (2016).

While mode frequencies are derived from space photometry, the patterns that lead to the mode identification of gravity modes as the start of the modeling in absence of rotational splitting rather concern the period spacings of the modes, $P_i^* = 2\pi/f_i^*$ (see Figure 12). Once mode identification has been achieved, we perform parameter estimation from the observed frequencies, because these are the measured quantities from the data. Moreover, this has the major advantage that we can easily deal with hybrid pulsators, which reveal both pressure and gravity modes, and can involve the task of searching for regularities in frequency when identifying pressure modes in addition to gravity-mode period spacing patterns (e.g., Bowman & Kurtz 2018).

As shown by Van Reeth et al. (2016), the shape of the observed period spacings of gravito-inertial modes allows us to break down the +7D estimation problem in a multistep approach, where the rotation frequency in the near-core region, $f_{\text{rot}}^{\text{core}}$, is estimated first by marginalizing over the other six parameters. Indeed, the near-core rotation of the star sets the “slope” of the period spacing pattern, and the value of this slope turns out to be sufficiently independent of ($M, X, Z, X_c, D_{\text{ov}}$) to apply a two-stage approach. Any “dips” in the tilted period spacing pattern are caused by mode trapping in the near-core region due to the value and shape of $D_{\text{ov}}(r)$ and $D_{\text{ext}}(r)$, and these can be estimated after $f_{\text{rot}}^{\text{core}}$ has been derived.

Van Reeth et al. (2016) and Ouazzani et al. (2017) developed methods to estimate $f_{\text{rot}}^{\text{core}}$ from the measured slope of the period spacings of gravity modes for intermediate-mass stars; the method by Van Reeth et al. (2016) was meanwhile also successfully applied to B-type pulsators (Pápics et al. 2017). With $f_{\text{rot}}^{\text{core}}$ estimated, the modeling process can be continued with exploitation of the *value* and *morphology* of the observed period spacing patterns, including the structure of the dips due to mode trapping in the near-core region. In summary, we propose the following forward modeling scheme:

1. Deduce the frequencies $f_i^*$ of pressure and gravity modes and their measurement error, $\epsilon_i^*$, in both the super-inertial and sub-inertial regime, from space photometry. Select the modes revealing rotational splitting and the $P_i^* = 1/f_i^*$ of $i = 1, ..., n$ that constitute period spacing patterns.

2. Identify the degree $l$ and azimuthal order $m$ of rotationally split modes, whenever detected. Use the slope of the observed period spacing patterns ($\Delta P_i, P_i^*$) to identify the degree $l$ and azimuthal order $m$ of the gravity modes $P_i^*$, by estimating the near-core rotation frequency $f_{\text{rot}}^{\text{core}}$, e.g., with the methods in Van Reeth et al. (2016) or Ouazzani et al. (2017). If pressure modes are identified from

measured rotational splitting, deduce the rotation frequency in the stellar envelope, $f_{\text{rot}}^{\text{env}}$, from the computation of the Ledoux constant (Equations (3.354) and (3.361) in Aerts et al. 2010). The estimation of $f_{\text{rot}}^{\text{core}}$ and/or $f_{\text{rot}}^{\text{env}}$ requires a sparse grid of stellar models $\mathcal{M}(\theta, \psi)$ for an appropriate range in mass, core overshoot, and metallicity and the computation of their pressure- and gravity-mode frequencies in an inertial frame of reference. This step will deliver a rough range for the mass and metallicity, as well as a narrow range for $f_{\text{rot}}^{\text{core}}$ (and of $f_{\text{rot}}^{\text{env}}$ if rotational splitting of pressure modes is detected).

3. Compute several dedicated +7D fine grids of stellar models with free parameters mass, age, and initial chemical composition, for the range of the estimated near-core rotation $f_{\text{rot}}^{\text{core}}$ and its error, and for a wide range of core overshooting $D_{\text{ov}}$ and envelope mixing $D_{\text{ext}}$. For the latter three quantities, assume particular shapes for their profiles, $f_{\text{rot}}(r), D_{\text{ov}}(r), D_{\text{ext}}(r)$, and also for the temperature gradient in the overshoot zone—radiative or adiabatic or a gradual transition between these two— along with fixed choices of the input physics according to Table 2, and treat those as different stellar models $\mathcal{M}(\theta, \psi)$. For all the model grids, compute their gravity-mode frequencies and identify the radial orders $n_{\text{pg}}$; of the detected modes $f_i^*$. For each model grid $\mathcal{M}(\theta, \psi)$, predict and estimate the parameters $\theta$ in +7D space and compute the uncertainties of $\theta$ from MLE, following the methodology developed for Problem 1 and possibly Problem 3 in Section 5. The user can get a feel for the effect of different input physics on the mode frequencies and on the forward modeling in the context of correlated parameters $\theta$ from attacking Problem 2 for various model grids $\mathcal{M}(\theta, \psi)$.

4. Select the most likely equilibrium models $\mathcal{M}(\theta, \psi)$ following Problem 4 for each star separately ($N=1$), and iterate over the scheme for an ensemble of $N$ gravity-mode pulsators. Select the overall best input physics and interpret its shortcomings vis-à-vis the ensemble of pulsators from the estimated values of the variance-covariance matrices $V(\sigma_{ii})$. The latter offer a guide to improve the input physics of stellar models.

**7. Conclusions**

In this work, we have provided a proper methodological framework to perform forward seismic modeling based on detected gravity-mode oscillation frequencies that constitute a period spacing pattern. We allow for correlations among the parameters of the stellar models, in order to properly take into account parameter degeneracies. The method relies on the computation of the Mahalanobis distance in a grid-based approach. We focused on stars rotating up to about half their critical rotation, such that the traditional approximation provides a valid approach to compute the theoretical oscillation frequencies. The method requires gravity-mode oscillation frequencies that have been measured with a precision better than 0.001 day⁻¹ and that can be identified in terms of their spherical degree $l$ and azimuthal order $m$. The parameter estimation is done in a +7D parameter space for each of the fixed choices of the input physics. As a minimum, it delivers estimation of the seven parameters ($M, X, Z, X_c, f_{\text{rot}}^{\text{core}}, D_{\text{ov}}, D_{\text{ext}}$) and their uncertainties, as well as an estimate of the
variance of each of the theoretical models *M*(*θ*, *ψ*)) and its uncertainty, as a guide to improve the input physics of future stellar structure and evolution models. A major aim of our future work is to achieve an accurate mass of the helium core near the TAMS for a representative sample of gravity-mode pulsators covering a large mass range.

The methodology presented here also offers a way to replace
physical model predictions of theoretical frequencies by
statistical model estimations, such that cumbersome computa-
tions of ultrafine stellar model grids can be avoided and that
proper error estimation can be achieved from the paradigm of
MLE estimation. A hierarchy in the importance of the physical
ingredients of the stellar structure models in terms of gravity-
mode frequency uncertainty was offered as a guide of the
modeling (Table 2).

We stress that the choice of the model parameters θ to
estimate depends on the type of star, its rotational properties,
and the type of mode frequencies detected in the data. Before
the CPU-intensive procedure of forward seismic modeling is
started, the user must assemble a good understanding of the key
physical ingredients of stellar structure and evolution models
and their computations for the star under study. Only after
proper understanding of stellar structure theory can one make a
sensible choice of the free parameters to estimate from forward
modeling. Our Figures 2–10, along with Table 2, offer a
practical guideline for this choice.

In this work, we did not cover the case of close binary stars whose equilibrium models and oscillation modes are subject to tidal forces. Tides introduce deformations and/or instabilities in the interior structure of the binary components, which may be accompanied by transport of elements and of angular momentum. Moreover, tidal forces affect the solutions to the pulsation equations. This is mainly the case for gravity modes, where (near) resonances between the frequencies of free oscillations and (multiples of) the orbit may occur. Just as for rotation and magnetism, the equilibrium models of close binaries due to tides are subject to major uncertainties. Also in this case, it is best to start forward modeling under the assumption that the tides can be treated as a secondary effect. An iterative scheme, in which a light curve and time-resolved high-precision spectroscopy covering the orbit are simultaneously solved for binarity and oscillations, will allow the oscillation frequencies to be derived. Once this has been achieved, one can perform the forward modeling as outlined here, where strong(er) (model-independent) constraints on the star can be included compared to the case of single stars. In this way, deviations between the observed oscillation frequencies and those of the optimal model by means of σᵢᵢ can be used as a guideline on how to improve the evolutionary models in terms of the missing tidal input physics.

Our formalism presented here will be applied in the near
future to numerous gravity-mode pulsators found in the *Kepler*
data, such as those in Van Reeth et al. (2016), Bowman &
Kurtz (2018), and Pápics et al. (2017). Such future applications
to tens of individual stars will undoubtedly improve our
knowledge of the input physics of stellar evolution models, as
well as pave the way to more automated applications to gravito-
inertial asteroseismology. The *TESS* mission will offer suitable
data for hundreds of gravity-mode pulsators in the mass range
1.4–40 *M*⊙ from its continuous viewing zone (Ricker
et al. 2016), coupled to spectroscopy (e.g., to be assembled
with SDSS-V; Kollmeier et al. 2017). On the longer term, the

PLATO mission will provide thousands of suitable targets for
gravity-mode asteroseismology (Rauer et al. 2014).

We express our sincere appreciation to the developer teams
of the MESA and GYRE codes, for making their software
available to the astronomical community. In particular, we
thank Aaron Dotter for his advice on how to include atomic
diffusion in MESA in an optimal way for the mass range
considered here. Further, we appreciate the numerous sugges-
tions from the referee to improve the presentation of our results
and to include more details of our analyses. C.A., G.M., and T.
V.R. are grateful for the kind hospitality at the Kavli Institute of
Theoretical Physics, University of California at Santa Barbara,
USA, where this work was initiated in April 2017. The authors
are also grateful to Christel Faes of the University of Hasselt
for her inspiring tutorial on Bayesian Inference in the
framework of the Scientific Research Network “Turning
Images into Value through Statistical Parameter Estimation”
funded by the Research Foundation Flanders under FWO grant
WO.010.16N. The research leading to these results has
received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and
innovation program (grant agreement no. 670519: MAMSIE),
from Interuniversitary Attraction Pole research Network P7/06
of the Belgian Government (Belgian Science Policy, Belspo),
and from the National Science Foundation of the United States
under grant NSF PHY-1748958.

*Software:* MESA (Paxton et al. 2011, 2013, 2015, 2018), GYRE (Townsend & Teitler 2013; Townsend et al. 2018), FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005).

Appendix A
MESA and GYRE Inlists

Example MESA and GYRE inlists used for this work are
available from the MESA Inlists section of the MESA Marketplace:
[cococubed.asu.edu/mesa_market/inlists.html](http://cococubed.asu.edu/mesa_market/inlists.html).

In addition, a MESA/GYRE tutorial on gravity-mode
asteroseismology can be retrieved from the Education section,
MESA summer school year 2016: cococubed.asu.edu/
mesa_market/education.html.

Appendix B
Solving the Score Equations in the Case of Nonlinear
Statistical Models with Fractional Polynomials

The aim in Problem 3 is to solve the score equations and
derive optimal values for β. The general case of nonlinear
statistical models to do so may be quite challenging, although
necessary in the case that linear models do not offer a proper fit.
Although classical polynomial predictors are still very
customary and relatively easy to use, they are often inadequate
because quickly very high degree polynomials are needed,
resulting in poor predictive properties. This can, in principle, be
overcome by using predictor functions that are nonlinear in the
parameters. However, such models, in turn, are often ridden
with computational challenges. Therefore, we would like to
point to a model family that shares with polynomial methods
the ease of fitting because it still is part of the generalized linear
models family, on the one hand, and that shares with nonlinear
models the flexibility of describing a wide array of functions.
This family is called “fractional polynomials” and was
advocated by Royston & Altman (1994).
A broad overview of the use and applications of this method in a variety of statistical areas can be found in Royston & Sauerbrei (2008). Several uses in hierarchical data structures are described in Verbeke & Molenberghs (2000) and Molenberghs & Verbeke (2005).

For a given degree *m* and a univariate argument *x* > 0, fractional polynomials are defined as

$$\beta_0 + \sum_{j=1}^{m} \beta_j x^{[p_j]},$$

where the $\beta_j$ are regression parameters, $x^{[p]} = x^p$ if $p \neq 0$ and $x^{[0]} = \ln(x)$. The powers $p_1 < ... < p_m$ are either positive or negative integers, or fractions. Royston & Altman (1994) argue that polynomials with degree higher than 2 are rarely required in practice and further restrict the powers of *x* to a small predefined set of noninteger values:

$$\Pi = \{-2, -1, -1/2, 0, 1/2, 1, 2, \dots, \max(3, m)\}.$$

The full definition includes possible “repeated powers” that involve multiplication with $\ln(x)$. For example, a fractional polynomial of degree *m* = 3 with powers $(-1, -1, 2)$ is of the form $\beta_0 + \beta_1x^{-1} + \beta_2x^{-1}\ln(x) + \beta_3x^2$ (Royston & Altman 1994; Sauerbrei & Royston 1999). Setting *m* = 2, for example, will generate the following:

(1) four “quadratics” in powers of *x*, represented by

(i) $\beta_0 + \beta_1 1/x + \beta_2 1/x^2$,

(ii) $\beta_0 + \beta_1 1/\sqrt{x} + \beta_2 1/x$,

(iii) $\beta_0 + \beta_1 \sqrt{x} + \beta_2 x$,

(iv) $\beta_0 + \beta_1 x + \beta_2 x^2$;

(2) a quadratic in $\ln(x)$: $\beta_0 + \beta_1 \ln(x) + \beta_2 \ln^2(x)$; and

(3) several other curves with shapes different from those of low-degree polynomials.

For given *m*, we consider as the best set of transformations the one producing the highest log-likelihood. For example, the best first-degree fractional polynomial is the one with the highest log-likelihood among the eight models with one regressor ($x^{-2}, x^{-1}, \dots, x^3$). As with conventional polynomials, the degree *m* is selected either informally on a priori grounds or by increasing *m* until no worthwhile improvement in the fit of the best-fitting fractional polynomial occurs. In the above discussion, it is assumed that *x* is strictly positive. If *x* can take zero values, a preliminary transformation of *x* is needed to ensure positivity (e.g., $x + 1$).

Evidently, the *x* in our case would be a component of $\theta$. Given that there are typically several components in $\theta$, the ideas above can be applied to all components simultaneously. This may give rise to a large set of possible predictors.

## ORCID iDs

C. Aerts https://orcid.org/0000-0003-1822-7126

P. I. Pápics https://orcid.org/0000-0002-3697-1664

## References

Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Berlin: Springer)

Aerts, C., De Cat, P., Kuschnig, R., et al. 2006, ApJL, 642, L165

Aerts, C., Molenberghs, G., Kenward, M. G., & Neiner, C. 2014, ApJ, 781, 88

Aerts, C., Thoul, A., Daszyńska, J., et al. 2003, Sci, 300, 1926

Aerts, C., Van Reeth, T., & Tkachenko, A. 2017, ApJL, 847, L7

Akaike, H. 1987, Psychometrika, 52, 317

Angulo, C., Arnould, M., Rayet, M., et al. 1999, NuPhA, 656, 3

Appourchaux, T. 2014, in Asteroseismology, 22nd Canary Islands Winter School of Astrophysics, ed. P. L. Palle & C. Esteban (Cambridge: Cambridge Univ. Press), 123

Appourchaux, T., Antia, H. M., Benomar, O., et al. 2014, A&A, 566, A20

Appourchaux, T., Benomar, O., Gruberbauer, M., et al. 2012a, A&A, 537, A134

Appourchaux, T., Chaplin, W. J., García, R. A., et al. 2012b, A&A, 543, A54

Appourchaux, T., Michel, E., Auvergne, M., et al. 2008, A&A, 488, 705

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Ausseloos, M., Scuflaire, R., Thoul, A., & Aerts, C. 2004, MNRAS, 355, 352

Ball, W. H., & Gizon, L. 2017, A&A, 600, A128

Ballot, J., Lignières, F., Prat, V., Reese, D. R., & Rieutord, M. 2012, in ASP Conf. Ser. 462, Progress in Solar/Stellar Physics with Helio- and Asteroseismology, ed. H. Shibahashi, M. Takata, & A. E. Lynas-Gray (San Francisco, CA: ASP), 389

Ballot, J., Lignières, F., Reese, D. R., & Rieutord, M. 2010, A&A, 518, A30

Bellinger, E. P., Angelou, G. C., Hekker, S., et al. 2016, ApJ, 830, 31

Benomar, O., Appourchaux, T., & Baudin, F. 2009, A&A, 506, 15

Bildsten, L., Ushomirsky, G., & Cutler, C. 1996, ApJ, 460, 827

Bloemen, S., Hu, H., Aerts, C., et al. 2014, A&A, 569, A123

Bossini, D., Miglio, A., Salaris, M., et al. 2017, MNRAS, 469, 4718

Bouabid, M.-P., Dupret, M.-A., Salmon, S., et al. 2013, MNRAS, 429, 2500

Bowman, D. M. 2017, PhD thesis, Univ. Central Lancashire

Bowman, D. M., & Kurtz, D. W. 2018, MNRAS, 476, 3169

Bowman, D. M., Kurtz, D. W., Breger, M., Murphy, S. J., & Holdsworth, D. L. 2016, MNRAS, 460, 1970

Braithwaite, J. 2009, MNRAS, 397, 763

Brassard, P., Fontaine, G., Billères, M., et al. 2001, ApJ, 563, 1013

Breger, M., Pamyatnykh, A. A., Pikall, H., & Garrido, R. 1999, A&A, 341, 151

Breger, M., Stich, J., Garrido, R., et al. 1993, A&A, 271, 482

Briquet, M., Aerts, C., Baglin, A., et al. 2011, A&A, 527, A112

Briquet, M., Neiner, C., Aerts, C., et al. 2012, MNRAS, 427, 483

Buchler, J. R., Goupil, M.-J., & Hansen, C. J. 1997, A&A, 321, 159

Buysschaert, B., Aerts, C., Bowman, D. M., et al. 2018, A&A, in press (arXiv:1805.00802)

Cameron, C., Saio, H., Kuschnig, R., et al. 2008, ApJ, 685, 489

Chaplin, W. J., & Miglio, A. 2013, ARA&A, 51, 353

Chapman, S., & Lindzen, R. 1970, Atmospheric. Thermal and Gravitational (Dordrecht: Reidel)

Charpinet, S., Van Grootel, V., Fontaine, G., et al. 2011, A&A, 530, A3

Claeskens, G., & Hjort, N. L. 2008, Model Selection and Model Averaging (Cambridge: Cambridge Univ. Press)

Caret, A., & Torres, G. 2017, ApJ, 849, 18

Corsaro, E., & De Ridder, J. 2014, A&A, 571, A71

Córiscal, A. H., Romero, A. D., Althaus, L. G., & Hermes, J. J. 2012, A&A, 547, A96

Cowling, T. G. 1941, MNRAS, 101, 367

Cox, D., & Hinkley, D. 1990, Theoretical Statistics (Boca Raton, FL, London: CRC Press, Chapman and Hall)

Daszyńska-Daszkiewicz, J., Pamyatnykh, A. A., Walczak, P., et al. 2017, MNRAS, 466, 2284

Daszyńska-Daszkiewicz, J., & Walczak, P. 2010, MNRAS, 403, 496

Davies, G. R., Silva Aguirre, V., Bedding, T. R., et al. 2016, MNRAS, 456, 2183

Deal, M., Escobar, M. E., Vauclair, S., et al. 2017, A&A, 601, A127

Deal, M., Richard, O., & Vauclair, S. 2016, A&A, 589, A140

Degroote, P., Aerts, C., Baglin, A., et al. 2010, Natur, 464, 259

Degroote, P., Briquet, M., Catala, C., et al. 2009, A&A, 506, 111

Deheuvels, S., Brandão, I., Silva Aguirre, V., et al. 2016, A&A, 589, A93

Deheuvels, S., Bruntt, H., Michel, E., et al. 2010, A&A, 515, A87

Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19

De Ridder, J., Telting, J. H., Balona, L. A., et al. 2004, MNRAS, 351, 324

Dotter, A., Conroy, C., Cargile, P., & Asplund, M. 2017, ApJ, 840, 99

Dupret, M.-A., Grigahcène, A., Garrido, R., Gabriel, M., & Scuflaire, R. 2005, A&A, 435, 927

Dupret, M.-A., Thoul, A., Scuflaire, R., et al. 2004, A&A, 415, 251

Dziembowski, W., & Krolikowska, M. 1985, AcA, 35, 5

Dziembowski, W. A., Daszyńska-Daszkiewicz, J., & Pamyatnykh, A. A. 2007, MNRAS, 374, 248

Dziembowski, W. A., & Goode, P. R. 1996, ApJ, 458, 338

Dziembowski, W. A., & Pamyatnykh, A. A. 2008, MNRAS, 385, 2061

Eckart, C. 1960, PhFl, 3, 421

Eggenberger, P., Lagarde, N., Miglio, A., et al. 2017, A&A, 599, A18
Geroux, C. M., & Deupree, R. G. 2013, ApJ, 771, 113

Giammichele, N., Charpinet, S., Brassard, P., & Fontaine, G. 2017, A&A, 598, A109

Giammichele, N., Charpinet, S., Fontaine, G., et al. 2018, Natur, 554, 73

Gruberbauer, M., Guenther, D. B., & Kallinger, T. 2012, ApJ, 749, 109

Gruberbauer, M., Guenther, D. B., MacLeod, K., & Kallinger, T. 2013, MNRAS, 435, 242

Gruberbauer, M., Kallinger, T., Weiss, W. W., & Guenther, D. B. 2009, A&A, 506, 1043

Handberg, R., Brogaard, K., Miglio, A., et al. 2017, MNRAS, 472, 979

Handler, G. 2009, MNRAS, 398, 1339

Handler, G., Jerzykiewicz, M., Rodríguez, E., et al. 2006, MNRAS, 365, 327

Handler, G., Matthews, J. M., Eaton, J. A., et al. 2009, ApJL, 698, L56

Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368

Hekker, S., & Christensen-Dalsgaard, J. 2017, A&ARv, 25, 1

Hermes, J. J., Kawaler, S. D., Romero, A. D., et al. 2017, ApJL, 841, L2

Hu, H., Tout, C. A., Glebbbeek, E., & Dupret, M.-A. 2011, MNRAS, 418, 195

Johnson, R., & Wichern, D. 2002, Applied Multivariate Statistical Analysis (Upper Saddle River, NJ: Prentice-Hall)

Johnson, R. A., & Wichern, D. W. 2000, Applied Multivariate Statistical Analysis (4th ed.; Englewood Cliffs: Prentice-Hall)

Kallinger, T., Weiss, W. W., Beck, P. G., et al. 2017, A&A, 603, A13

Kass, R., & Raftery, A. 1995, J. Am. Stat. Assoc., 90, 773

Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (Berlin: Springer)

Kollmeier, J. A., Zasowski, G., Rix, H.-W., et al. 2017, arXiv:1711.03234

Kudritzki, R.-P., & Puls, J. 2000, ARA&A, 38, 613

Kurtz, D. W., Saio, H., Takata, M., et al. 2014, MNRAS, 444, 102

Lagarde, N., Bossini, D., Miglio, A., Vrard, M., & Mosser, B. 2016, MNRAS, 457, L59

Lai, D. 1997, ApJ, 490, 847

Ledoux, P. 1951, ApJ, 114, 373

Lee, U., & Baraffe, I. 1995, A&A, 301, 419

Louis, T. L., & Carlin, B. P. 2009, Bayesian Methods for Data Analysis (London: Chapman and Hall)

Loumos, G. L., & Deeming, T. J. 1978, Ap&SS, 56, 285

Martins, F., Escolano, C., Wade, G. A., et al. 2012a, A&A, 538, A29

Martins, F., Mahy, L., Hillier, D. J., & Rauw, G. 2012b, A&A, 538, A39

Mathis, S. 2013, in Lecture Notes in Physics, Vol. 865 ed. M. Goupil et al. (Berlin: Springer), 23

Michaud, G., Richard, O., Richer, J., & VandenBerg, D. A. 2004, ApJ, 606, 452

Miglio, A., Brogaard, K., Stello, D., et al. 2012, MNRAS, 419, 2077

Miglio, A., & Montalbán, J. 2005, A&A, 441, 615

Molenberghs, G., & Verbeke, G. 2005, Models for Discrete Longitudinal Data (New York: Springer)

Moravveji, E., Aerts, C., Pápics, P. I., Triana, S. A., & Vandoren, B. 2015, A&A, 580, A27

Moravveji, E., Townsend, R. H. D., Aerts, C., & Mathis, S. 2016, ApJ, 823, L13

Morel, T., Hubrig, S., & Briquet, M. 2008, A&A, 481, 453

Nash, J. C. 1990, The Choleski Decomposition, in: Compact Numerical Methods for Computers: Linear Algebra and Function Minimisation (2nd ed.; Bristol: Hilger)

Neiner, C., Gutiérrez-Soto, J., Baudin, F., et al. 2009, A&A, 506, 143

Neter, J., Wasserman, W., & Kutner, M. 1990, Applied Linear Statistical Models, Regression, Analysis of Variance and Experimental Designs (3rd ed.; Homewood, IL: Richard D. Irwin, Inc.)

Nowakowski, R. M. 2005, AcA, 55, 1

Ouazzani, R.-M., Marques, J. P., Goupil, M.-J., et al. 2018, A&A, submitted

Ouazzani, R.-M., Salmon, S. J. A. J., Antoci, V., et al. 2017, MNRAS, 465, 2294

Pablo, H., Whittaker, G. N., Popowicz, A., et al. 2016, PASP, 128, 125001

Pápics, P. I., Briquet, M., Baglin, A., et al. 2012, A&A, 542, A55

Pápics, P. I., Moravveji, E., Aerts, C., et al. 2014, A&A, 570, A8

Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, A&A, 598, A74

Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3

Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4

Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15

Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34

Pedersen, M. G., Aerts, C., Pápics, P. I., & Rogers, T. M. 2018, A&A, 614, 128

Penrose, R. 1955, PCPS, 51, 406

Poretti, E., Michel, E., Garrido, R., et al. 2009, A&A, 506, 85

Przybilla, N., Nieva, M. F., Irrgang, A., & Butler, K. 2013, EAS, 63, 13

Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669

Quirion, P.-O., Christensen-Dalsgaard, J., & Arentoft, T. 2010, ApJ, 725, 2176

Rauer, H., Catala, C., Aerts, C., et al. 2014, ExA, 38, 249

Reese, D., Lignières, F., & Rieutord, M. 2006, A&A, 455, 621

Richer, J., Michaud, G., & Turcotte, S. 2000, ApJ, 529, 338

Ricker, G. R., Vanderspek, R., Winn, J., et al. 2016, Proc. SPIE, 9904, 99042B

Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, JCoPh, 318, 277

Rodrigues, T. S., Bossini, D., Miglio, A., et al. 2017, MNRAS, 467, 1433

Rogers, T. M. 2015, ApJL, 815, L30

Rogers, T. M., & McElwaine, J. N. 2017, ApJL, 848, L1

Romero, A. D., Córsico, A. H., Castanheira, B. G., et al. 2017, ApJ, 851, 60

Royston, P., & Altman, D. 1994, Applied Statistics, 43, 429

Royston, P., & Sauerbrei, W. 2008, Multivariable Model Building: A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables (Chichester: Wiley)

Saio, H. 1981, ApJ, 244, 299

Saio, H., Bedding, T. R., Kurtz, D. W., et al. 2018a, MNRAS, 477, 2183

Saio, H., Kurtz, D. W., Murphy, S. J., Antoci, V. L., & Lee, U. 2018b, MNRAS, 474, 2774

Saio, H., Kurtz, D. W., Takata, M., et al. 2015, MNRAS, 447, 3264

Saio, H., Kuschnig, R., Gautschy, A., et al. 2006, ApJ, 650, 1111

Salaris, M., & Cassisi, S. 2017, RSOS, 4, 170192

Santolaya-Rey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488

Sauerbrei, W., & Royston, P. 1999, Journal of the Royal Statistical Society, Series A, 162, 71

Schmid, V. S., & Aerts, C. 2016, A&A, 592, A116

Schwarz, G. 1978, AnSta, 6, 461

Seber, G., & Wild, C. 2003, Nonlinear Regression (New York: Wiley)

Silva Aguirre, V., Lund, M. N., Antia, H. M., et al. 2017, ApJ, 835, 173

Smolec, R., & Moskalik, P. 2007, MNRAS, 377, 645

Sowicka, P., Handler, G., Dębski, B., et al. 2017, MNRAS, 467, 4663

Szewczuk, W., & Dązyńska-Dązkiewicz, J. 2017, MNRAS, 469, 13

Szewczuk, W., & Dązyńska-Dązkiewicz, J. 2018, MNRAS, arXiv:1805.07100

Tayar, J., & Pinsonneault, M. H. 2013, ApJL, 775, L1

Townsend, R. H. D. 2000a, MNRAS, 318, 1

Townsend, R. H. D. 2000b, MNRAS, 319, 289

Townsend, R. H. D. 2003a, MNRAS, 340, 1020

Townsend, R. H. D. 2003b, MNRAS, 343, 125

Townsend, R. H. D., Goldstein, J., & Zweibel, E. G. 2018, MNRAS, 475, 879

Townsend, R. H. D., & Teitler, S. A. 2013, MNRAS, 435, 3406

Trampedach, R., Aarslev, M. J., Houdek, G., et al. 2017, MNRAS, 466, L43

Triana, S. A., Moravveji, E., Pápics, P. I., et al. 2015, ApJ, 810, 16

Truyaert, K. 2016, Master thesis, KU Leuven

Turcotte, S., Richer, J., & Michaud, G. 1998, ApJ, 504, 559

Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: Univ. of Tokyo)

Van Reeth, T., Tkachenko, A., Aerts, C., et al. 2015, ApJS, 218, 27

Van Reeth, T., Tkachenko, A., & Aerts, C. 2016, A&A, 593, A120

Verbeke, G., & Molenberghs, G. 2000, Linear Mixed Models for Longitudinal Data (New York: Springer)

Verma, K., Raodeo, K., Antia, H. M., et al. 2017, ApJ, 837, 47

Wade, G. A., Neiner, C., Alecian, E., et al. 2016, MNRAS, 456, 2

Walker, G. A. H., Kuschnig, R., Matthews, J. M., et al. 2005, ApJL, 635, L77

Welsh, A. 1996, Aspects of Statistical Inference (New York: Wiley)

White, H. 1982, Econometrica, 50, 1

Zorec, J., & Royer, F. 2012, A&A, 537, A120

Zwintz, K., Fossati, L., Guenther, D. B., et al. 2013, A&A, 552, A68

Zwintz, K., Fossati, L., Ryabchikova, T., et al. 2014, Sci, 345, 550

Zwintz, K., Moravveji, E., Pápics, P. I., et al. 2017, A&A, 601, A101