11institutetext: Institut für Astro- und Teilchenphysik, Universität
Innsbruck, Technikerstraße 25, A-6020 Innsbruck, Austria
11email<EMAIL_ADDRESS>

# Pulsational instability of pre-main-sequence models from accreting
protostars

II. Modelling echelle diagrams of $\delta$ Scuti stars without rotational
splitting

T. Steindl 11 K. Zwintz 11 M. Müllner 1111

(submitted 01.02.2022; accepted 25.05.2022)

###### Abstract

Context. The physics of early stellar evolution (e.g. accretion processes) is
often not properly included in the calculations of pre-main-sequence models,
leading to insufficient model grids and hence systematic errors in the
results.

Aims. We aim to investigate current and improved approaches for the
asteroseismic modelling of pre-main-sequence $\delta$ Scuti stars.

Methods. We calculated an extensive grid of pre-main-sequence models including
the early accretion phase and used the resulting equilibrium models as input
to calculate theoretical frequency spectra. These spectra were used to
investigate different approaches in modelling echelle diagrams to find the
most reliable methods. By applying Petersen diagrams, we present a simple
algorithm to extract echelle diagrams from observed pulsation frequencies.

Results. We show that model grids with insufficient input physics and
imperfect modelling approaches lead to underestimated uncertainties and
systematic errors in the extracted stellar parameters. Our re-discussion of HD
139614 leads to different stellar parameters than the ones derived by Murphy
et al. (2021). We performed a model comparison between this previous
investigation and our results by applying the Akaike and Bayesian information
criteria. While the results with regard to our ten-dimensional model are
inconclusive, they show (very) strong evidence of a six-dimensional model with
fixed accretion parameters (leading to almost identical stellar parameters to
those of the ten-dimensional model) to be preferred over the model applied by
Murphy et al. (2021). In general, our modelling approach can provide narrow
constraints on the stellar parameters (i.e. $\Delta R\sim 0.05\,R_{\odot}$,
$\Delta\log\,g\lesssim 0.01$, and $\Delta M_{\star}\sim 0.1\,M_{\odot}$).

Conclusions. The extensively tested modelling approaches and automatic
extraction of echelle diagrams should allow us to study many more pre-main-
sequence $\delta$ Scuti stars in the future and lead to reliable stellar
parameters.

###### Key Words.:

asteroseismology – stars: oscillations (including pulsations) – accretion,
accretion disks – stars: variables: delta Scuti – stars: pre-main-sequence –
stars: protostars

##  1 Introduction

The evolution of intermediate-mass stars is often truncated to their main-
sequence, red-giant, and asymptotic-giant-branch phases, and their final fate
as white dwarfs. However, before stars are able to burn significant amounts of
hydrogen in their cores, they undergo the phase of pre-main-sequence
evolution. During this phase, which starts from the formation of the second
Larson’s core (or second hydrastic core, see Larson 1969 and Masunaga &
Inutsuka 2000), a multitude of physical processes shape the stars and, hence,
their future fate. While these processes play a major role in stellar
evolution (e.g. for the formation of magnetic fields, angular momentum
transport, and planet formation), they are often ignored in model
calculations. Many studies of early stellar evolution still rely on the
pioneering works of Henyey et al. (1955), Hayashi (1961), and Iben (1965),
albeit ignoring first and foremost any accretion processes. While the outdated
view of the pre-main-sequence evolution just assumes the initial models with
huge radii, the addition of accretion processes in stellar modelling allows
for a much more realistic view. This series of papers aims to combine such
more realistic evolution models with asteroseismology, which is the study of
pulsating stars.

Pre-main-sequence stars with masses between $1.5$ and $3.5\,M_{\odot}$ are
expected to become unstable to stellar pulsations (e.g. Steindl et al., 2021a)
as $\delta$ Scuti (e.g. Zwintz et al., 2014) or slowly pulsating B-type stars
(e.g. Gruber et al., 2012). Pre-main-sequence $\delta$ Scuti stars are
expected to have effective temperatures in the range from $6\,300$ to
$10\,300$ K (Steindl et al., 2021a) and pulsate in low-order pressure modes
driven by the heat mechanism (Dupret et al., 2005; Steindl et al., 2021a).
Zwintz et al. (2014) show that their pulsation characteristics are dependent
on their relative evolutionary status, with stars further evolved towards the
ZAMS pulsating at higher frequencies. Bedding et al. (2020) present
regularities in the pulsation spectra of mainly young main-sequence $\delta$
Scuti stars. Such regularities are often shown in echelle diagrams and allow
the asteroseismic modelling of $\delta$ Scuti stars (see e.g. Bernabei et al.,
2009; Chen & Li, 2019; Murphy et al., 2021; Kim et al., 2021).

In slowly rotating stars, echelle diagrams are expected to be much simpler
compared to their fast-rotating counterparts, since pulsation modes are split
according to the angular velocity of the stars (Aerts et al., 2010). It is
possible to automatically infer an echelle diagram and, hence, identify the
pulsation modes from the detected pulsation frequencies of $\delta$ Scuti
stars. From that point on, it is a rather simple exercise to produce a grid of
theoretical stellar evolution models and pulsation spectra and apply
asteroseismic model approaches. Here, we present the limitations of this
method if the underlying stellar grid insufficiently includes physics of early
stellar evolution and if additional free parameters are ignored. In addition,
the uncertainties might be underestimated and the stellar parameters might be
extracted incorrectly.

In this work, we calculated an extensive grid of pre-main-sequence models in
the mass range expected for $\delta$ Scuti stars, starting from the second
hydrostatic core and including the accretion phase. In addition, we allowed
for free parameters in the initial composition and mixing parameters needed to
accurately describe possible evolutionary pathways of early stars. The
computational set-up and calculation of the grid are described in Section 2.
In Section 3, we review, discuss, and introduce different approaches for the
asteroseismic modelling of pre-main-sequence $\delta$ Scuti stars. These
approaches were then extensively tested using our model grid, which we
describe in Section 4. We introduced a simple algorithm based on the results
of the first paper of this series (Steindl et al., 2021a) to automatically
extract echelle diagrams, which we discuss in Section 5. We then performed
asteroseismic modelling using the theoretical pulsation frequencies of our
state-of-the-art pre-main-sequence grid on the star HD 139614, which was
previously discussed by Murphy et al. (2021), and we offer a model comparison
in Section 6. We chose to re-discuss HD 139614 since it is very slowly
rotating (Murphy et al., 2021) and, hence, a perfect specimen for the
methodology described in this work. We conclude our results in Section 7.

##  2 Computational set-up

Table 1: Summary of the covered input physics.

param. | description | value range | step size
---|---|---|---
Z | metallicity | 0.006 - 0.024 | 0.002
Y | helium abundance | 0.204 - 0.234 | 0.01
$\alpha_{\rm MLT}$ | mixing length | 1.8 - 2.4 | 0.2
$M_{\star}$ | final mass | 1.5 - 2.5 | continuous
$F$ | overshooting | 0-0.02 | continuous
$D_{\rm mix}$ | envelope mixing | 0-100 ${\rm cm}{\rm s}^{-1}$ | continuous
$\beta$ | fractional energy | 0.05-0.2 | continuous
$M_{\rm outer}$ | fractional mass | 0.01-0.2 | continuous
$\dot{M}_{0}$ | accretion rate | $10^{-7}-10^{-5}\,M_{\odot}/{\rm yr}$ | continuous
$\tau$ | e-folding time | 0.1 -1 Myr | continuous

111As each combination of the parameters Z, Y, and $\alpha_{\rm MLT}$ call for
an initial model, we limited these parameters to periodic grids.

Figure 1: Parameters of 182 547 distinct pre-main-sequence $\delta$ Scuti
models. The top left panel shows a Kiel diagram with all models as black
crosses and a few random evolutionary tracks as magenta lines. All other
panels show the histograms of the occurence distribution of 15 different
parameters, where the y-axis shows the number of models in a given bin.

###  2.1 Stellar evolution models

We used version v-12778 of Modules for Experiments in Stellar Astrophysics
(MESA, Paxton et al., 2011, 2013, 2015, 2018, 2019) to calculate the stellar
evolution models. The software instrument MESA solves the fully coupled
structure and composition equations simultaneously for a one-dimensional,
spherically symmetric stellar model (Paxton et al., 2011). We refer the reader
to the MESA instrument papers for a full description of the numerical methods.
More details about the inputs for the microphysics are given in Appendix A.

The physics for the stellar evolution models is similar to the set-up used in
the first paper of this series (Steindl et al., 2021a). We used the OPAL
opacity tables (Seaton, 2005) and the standard MESA equation of state (Paxton
et al., 2011) while ignoring stellar rotation and magnetic fields. Convection
is treated via the mixing length theory developed by Cox & Giuli (1968).
Convective boundaries are found with the Ledoux criterion and the convective
premixing approach (Paxton et al., 2018) while not allowing for semi-
convection. We computed stellar evolution models with mixing lengths
$\alpha_{\rm MLT}=1.8,2.0,2.2,2.4$. Beyond the convective boundary, we
employed exponential overshooting (Herwig, 2000; Paxton et al., 2011)
according to

$\displaystyle
D_{OV}=D_{\mathrm{conv,}0}\,\exp\left(\frac{-2z}{fH_{p}}\right)\,,$ (1)

where $D_{\mathrm{conv,}0}$ is the mixing coefficient near the Schwartzschild
boundary, $H_{p}$ is the pressure scale height, $z$ is the distance to the
radiative layer, and $f$ is a free parameter (Herwig, 2000) set differently
for the overshooting towards the surface and the overshooting towards the core
of the star. An additional free parameter, $f_{0}$ describes the position
$f_{0}H_{p}$ into the convective zone where the switch between convection and
convective overshooting occurs. The resulting four parameters for the
overshooting above and below a convective zone in our models are controlled by
a parameter $F$ between $0$ and $0.2,$ which results in $f_{\rm above}=F$,
$f_{0,{\rm above}}=f_{\rm below}=\frac{F}{2}$, and $f_{0,{\rm
below}}=\frac{F}{4}$. The minimum mixing coefficient has values between
$D_{\rm mix}=0$ and $100\,{\rm cm}s^{-1}$.

The initial composition is set to
$X_{{}^{1}\mathrm{H}}=1-X_{{}^{2}\mathrm{H}}-X_{{}^{3}\mathrm{He}}-X_{{}^{4}\mathrm{He}}-Z$,
where $X_{{}^{3}\mathrm{He}}=85$ ppm,
$X_{{}^{4}\mathrm{He}}=Y+2*Z-X_{{}^{3}\mathrm{He}}$, where
$Y=0.204,0.214,0.224,$ or $0.234$ based on the present day surface abundance
of the sun and an initial metallicity of $Z,$ which is a value between $0.006$
and $0.024$ spaced out by $0.002$. The mass fractions of metals were taken
according to Asplund et al. (2009) updated based on Nieva & Przybilla (2012)
and Przybilla et al. (2013). The calculation of models is stopped once the
central hydrogen fraction has dropped by $0.01$ relatively to the starting
value. That is, the evolution includes the earliest part of the main-sequence
lifetime.

The initial models are the same $0.01\,M_{\odot}$ and $1.5\,R_{\odot}$ stellar
seeds we used in Steindl et al. (2021a), where we followed the approach of
Kunitomo et al. (2017). Hence, this is similar to the approach of Stahler et
al. (1980) and Hosokawa et al. (2011). These stellar seeds correspond to
slightly evolved protostars (Kunitomo et al., 2017) with a high value for the
initial entropy of the star (Baraffe et al., 2012).

For atmospheric boundary conditions, we found that only the temperature-
opacity ($T{-}\tau$) relation by Eddington (1926) successfully recreates the
instability region for pre-main-sequence $\delta$ Scuti stars (Steindl et al.,
2021a). Hence, we used the Eddington atmosphere for our calculations in this
work.

####  2.1.1 Treatment of accretion

We describe the energy budget of accretion in Section 5.1.3 of the first paper
in this series (Steindl et al., 2021a). Here, we only briefly discuss the
details of the treatment of accretion, which are, in part, different to our
previous work.

The energy of the accreted material that is added to the star ($L_{\rm add}$)
is controlled by a factor $\beta,$ such that

$\displaystyle L_{\rm add}=\beta\epsilon\frac{GM\dot{M}}{R}\,,$ (2)

where $\epsilon\leq 0.5$ describes the accretion from a thin disc at the
object’s equator (Hartmann et al., 1997; Baraffe et al., 2009), and
$\epsilon=0.5$ is a common choice (Baraffe et al., 2009; Jensen & Haugbølle,
2018; Steindl et al., 2021a). Furthermore, $G$ is the gravitational constant,
$M$ is the current mass of the star, $\dot{M}$ the mass accretion rate, and
$R$ the current radius of the star.

For $\beta=0$, no energy is deposited inside the stellar interior. This case
is hence also referred to as cold accretion while $0<\beta\leq 1$ corresponds
to hot accretion. With time-dependent accretion rates, $\beta$ is also
expected to vary (e.g. Baraffe et al., 2012; Vorobyov et al., 2017; Jensen &
Haugbølle, 2018; Elbakyan et al., 2019, and references therein). Since the
accretion rates used in this work are constant, besides an exponential drop
towards the end of the accretion phase, we used a fixed value of $\beta$
between $0.05$ and $0.2$.

In MESA, the internal energy of the accreted material is added as extra heat
to the star. The distribution of this extra heat in the stellar model is
implemented differently in the literature. As a first approach, Baraffe &
Chabrier (2010) used a uniform distribution that is most likely too simple to
mirror the processes in real stars. As in the first paper of this series, we
distributed the heat in an outer region of fractional mass $M_{\rm outer}$
with a linear increase as a function of the mass coordinate $m_{r}$ (see also
Kunitomo et al., 2017) according to

$\displaystyle l=\frac{L_{\rm add}}{M}{\rm max}\left\\{0,\frac{2}{M_{\rm
outer}^{2}}\left(\frac{m_{r}}{M}-1+M_{\rm outer}\right)\right\\}.$ (3)

We calculate models with different values of $M_{\rm outer}$, ranging between
$0.01$ and $0.2$. Our models are calculated within the constant accretion
scenario, meaning that the mass accretion rate $\dot{M}$ is constant in the
beginning of the evolution. In contrast with the first paper of the series,
where we used an abrupt change in mass accretion rate, we used an exponential
drop with different timescales for the models calculated in this work. Hence,
the mass accretion rate is given by

$\displaystyle\dot{M}=\cases{\dot{\hfil}}{M}_{0}$ ift<texp, (4) (5)

where $\dot{M}_{0}$ is the initial constant accretion rate chosen between
$10^{-7}$ and $10^{-5}\,M_{\odot}/{\rm yr}$, $t$ is the age of star since the
beginning of the simulation, $\tau$ is the timescale for ten e-folding times
in the mass-accretion drop chosen between $0.1$ and $1$ Myr, and $t_{\rm exp}$
is the time until the exponential drop-off calculated as

$\displaystyle t_{\rm
exp}=\left(M_{\star}-0.01-\dot{M}_{0}\,\left(1-\exp(-10)\right)\,\tau\times
10^{6}\,\right)\dot{M}_{0}^{-1}.$ (7)

Here, $M_{\star}$ is the final mass of the star.

###  2.2 Stellar pulsations

The stellar oscillation code GYRE (Townsend & Teitler, 2013; Townsend et al.,
2018; Goldstein & Townsend, 2020) uses the equilibrium model resulting from
MESA computations and calculates the linear eigenfrequencies for normal
oscillation modes. In this work, we used version 6.0.1 of GYRE and refer the
reader to the instrument papers (Townsend & Teitler, 2013; Townsend et al.,
2018; Goldstein & Townsend, 2020) for numerical details.

We calculated modes with harmonic degrees $l=0,1,$ and $2$ (radial, dipole,
and quadrupole modes) with azimuthal order $m=0$ for p modes for radial orders
$n=1,\cdots,8$ and the frequency range from $5-100$ $\mathrm{d}^{-1}$. We
chose regularity-enforcing inner boundary conditions and outer boundary
conditions according to the formulation of Dziembowski (1971). In contrast to
the first paper, we performed the GYRE calculations within the adiabatic
framework, since we are only interested in the pulsation frequencies and the
calculation of non-adiabatic frequencies is computationally much more
expensive.

In general, $\delta$ Scuti stars often show very rich amplitude spectra since
the pulsation frequencies of rotating stars are split according to their
azimuthal order and rotational properties. In this work, we ignored the
splitting of modes and focused on the modes with an azimuthal order of $m=0$.
The inclusion of rotational effects will be the subject of future work.

###  2.3 Computational sample

To produce a computational sample of pre-main-sequence $\delta$ Scuti stars,
we performed a total of $1500$ evolutionary calculations with a random
combination of the input physics described above. A summary of the covered
parameter space is given in Table 1. The convergence of these calculations is
dependent on the choice of input physics (see the discussion in Section 5.1.5
of Steindl et al. (2021a)). As a consequence, $1309$ calculations provide pre-
main-sequence $\delta$ Scuti models, while the remaining $191$ failed to
converge, mostly very early in the evolution. During the calculation, we saved
a model for the subsequent GYRE calculation if the effective temperature was
within the expected instability region ($6300$ K - $10300$ K) for pre-main-
sequence p-mode pulsations (Steindl et al., 2021a). This results in a total of
182 547 distinct $\delta$ Scuti models. The parameters listed in Table 1 are
direct input parameters for the models. Every model additionally includes a
star age, effective temperature, surface gravity, and luminosity. The latter
are a direct consequence of the models saved for the GYRE calculations. Figure
1 provides an overview over all models calculated for the sample.

Figure 2: Echelle diagram of theoretical model #10001. The red, green, and
blue symbols correspond to the frequencies of the fundamental and first seven
overtones for the radial, dipole, and quadrupole modes, respectively. The
echelle diagram is repeated to the right of the dashed black line for better
visibility.

##  3 Modelling approach

In this work, we explored different modelling approaches. For the first ones,
we followed previous studies (e.g. Bernabei et al., 2009; Chen & Li, 2019;
Murphy et al., 2021; Kim et al., 2021) with the main idea of using a
pseudo-$\chi^{2}$-square statistic:

$\displaystyle\chi^{2}=\frac{1}{n}\sum_{i=1}^{n}\frac{\left(f_{i,{\rm
obs}}-f_{i,{\rm mod}}\right)^{2}}{\sigma_{f}},$ (8)

where $f_{i,{\rm obs}}$ are the observed frequencies, $f_{i,{\rm mod}}$ the
model frequencies, and $\sigma_{f}$ the uncertainty on the observation. We
used the Rayleigh limit ($\frac{1}{T}$, where $T$ is the timespan of
observations) as the uncertainty of all observed frequencies. In the
calculation of the pseudo-$\chi^{2}$-square statistic, only frequencies with
the same mode identification (radial order $n$ and mode degree $l$) are
considered, essentially rendering this a modelling of the echelle diagram
including the fundamental frequency and the first seven overtones. After the
pseudo-$\chi^{2}$-square statistic was evaluated for all sets of model
frequencies, we investigated the properties of all models fulfilling
$\chi^{2}\leq 1$. Expectation values and uncertainties of stellar parameters
are then obtained by additionally applying a spectroscopic observation box
(i.e. ignoring all models that do not fall within the $1\sigma$ or $3\sigma$
uncertainty box of $T_{\rm eff}$ and $\log\,g$) and calculating the mean or
median as well as the standard deviation of the remaining models.

A different option is to include the spectroscopic observations within the
pseudo-$\chi^{2}$-square statistic, such that

$\displaystyle\begin{split}\chi^{2}_{\rm
spec}=\frac{1}{n+2}&\left(\sum_{i=1}^{n}\frac{\left(f_{i,{\rm obs}}-f_{i,{\rm
mod}}\right)^{2}}{\sigma_{f}}+\frac{\left(T_{{\rm eff},{\rm obs}}-T_{{\rm
eff},\rm mod}\right)^{2}}{\sigma_{T_{\rm eff}}}+\right.\\\
&\quad\quad\left.+\frac{\left(\log\,g_{{\rm obs}}-\log\,g_{\rm
mod}\right)^{2}}{\sigma_{T_{\log\,g}}}\right).\end{split}$ (9)

However, the latter will reject any model that lies outside the spectroscopic
error box if we later apply the condition $\chi^{2}_{\rm spec}\leq 1$. A more
statistically motivated approach is to reject all models with a
pseudo-$\chi^{2}$ statistic greater than the expected $\chi^{2}$ value for a
pre-defined p value, for instance $28.87$ for a p value of $0.05$ and 18
degrees of freedom.

This modelling approach has significant downsides, as discussed by Aerts et
al. (2018) who proposed the use of a Mahalanobis distance for the forward
modelling of gravity-mode pulsators. The calculation of the Mahalanobis
distance, however, is in need of a variance-covariance matrix obtained from
the model grid (Aerts et al., 2018). The model grid calculated in this work is
not dense enough to provide a good estimate for the variance-covariance
matrix. Hence, we are not able to use the Mahalanobis distance as a model
statistic. The work of Aerts et al. (2018) focused on the in-depth modelling
of g-mode pulsators for which usual observation time spans of years lead to
very low Rayleigh limits. This consequently allows for a much more in-depth
study, than the observation time spans of about $30$ to $100$ days usually
found for pre-main-sequence $\delta$ Scuti stars. Hence, the modelling of pre-
main-sequence $\delta$ Scuti stars pursues lower goals than the forward
modelling of gravity mode pulsators. The aim of this work is to investigate
whether this simple pseudo-$\chi^{2}$-square statistic is able to correctly
infer the stellar parameters from an identified echelle diagram and to what
extent the resulting expectation values and uncertainties can be trusted. As
such, we applied five distinct modelling approaches, three of which are based
on Equation 8, while the other two are based on Equation 9. The approaches are
as follows: (1) all models with $\chi^{2}\leq 1$; (2) all models with
$\chi^{2}\leq 1$ and that lie within the $3\sigma$ spectroscopic error box;
(3) all models with $\chi^{2}\leq 1$ and that lie within the $1\sigma$
spectroscopic error box; (4) all models with $\chi^{2}_{\rm spec}\leq 1$; (5)
all models with $\chi^{2}_{\rm spec}\leq\chi^{2}_{\rm cut}$, where
$\chi^{2}_{\rm cut}$ corresponds to a p value of $0.05$ for the number of free
parameters given by the number of frequencies plus two. Approaches 2 and 3 are
usually applied in the literature. This study aims to test the validity of
these approaches.

##  4 Estimating the stellar parameters of theoretical echelle diagrams

###  4.1 Result for a single model

Figure 3: Hexbin plot of the resulting pseudo-$\chi^{2}$ statistic according
to Equation (8) (upper panel) and Equation (9) (lower panel) for the
theoretical model #10001 with a Rayleigh limit of $\frac{1}{27}$
$\mathrm{d}^{-1}$. The colours of the bins show the minimal value of all
models that fall within the bin. Black lines show the same evolutionary tracks
as in Figure 1. The correct position of model #10001 is marked with the green
circle and the assumed spectroscopic error box is drawn in green.  Figure 4:
Pseudo-$\chi^{2}$ statistic distribution of effective temperature, star mass,
surface gravity, and age for the frequencies of model #10001. Models that fall
within the $1\sigma$ spectroscopic box are shown in green, those of the
$3\sigma$ spectroscopic box are in blue, and all others are in black. The
horizontal dashed grey line marks the position of pseudo-$\chi^{2}=1,$ and the
vertical red line shows the actual value for model #10001.

We first applied our modelling approaches to one model (#10001). With an
effective temperature of $7432$ K, mass of $1.85\,M_{\odot}$, and logarithmic
surface gravity of $4.208$ at an age of $8.12$ Myr, this model represents
typical values of a bona fide pre-main-sequence $\delta$ Scuti star (Zwintz et
al., 2014; Steindl et al., 2021b, a). The echelle spectrum of this model is
shown in Figure 2. To model these theoretical frequencies, we chose a time
span of $27$ days, corresponding to one-sector observations of the Transiting
Exoplanet Survey Satellite (TESS, Ricker et al., 2015), and perturbed all
frequencies with a uniformly chosen random number between $-\frac{1}{T}$ and
$+\frac{1}{T}$ where $\frac{1}{T}$ is the resulting Rayleigh limit. In
addition, we assumed a typical uncertainty of $200$ K on the effective
temperature and $0.1$ on the surface gravity for the spectroscopic error box.

Stellar pulsation modes are represented by spherical harmonics, and higher
degree modes suffer from cancellation effects, such that modes with $l\geq 2$
are expected to have lower amplitudes in the light curves than radial and
dipole modes (Aerts et al., 2010). In addition, the frequencies of quadrupole
modes make a detection of a convincing echelle ridge in observed frequencies
much more difficult (see discussion in Section 5). Hence, for this example we
only used the radial and dipole modes of model #10001. As is evident from the
upper panel of Figure 3, a range of models with $\log\,g\sim 4.2$ has low
values for the pseudo-$\chi^{2}$ statistics according to Equation (8). The
five models with the lowest pseudo-$\chi^{2}$ ($\ll 1$) also include, next to
model #10001, model #52851 with $T_{\rm eff}=8137$ K,
$M_{\star}=2.02\,M_{\odot}$, $\log\,g=4.22$, and an age of $7.35$ Myr. The
stellar parameters of model #52851 differ significantly from the the ones of
model #10001. The stellar parameters of many more models, with
pseudo-$\chi^{2}$ statistics well below one, differ to even higher degrees.
Thus, it is already evident that the choice of how to enforce the
spectroscopic observations plays a major role in the extracted stellar
parameters. Figure 4 shows the pseudo-$\chi^{2}$ distribution for these four
parameters, with models falling into the spectroscopic error box coloured
differently. This further proves that the inclusion of the spectroscopic
parameters is a delicate matter. The choice to do so will influence the
resulting expectation value and uncertainties.

As a comparison, the lower panel of Figure 3 shows the resulting
pseudo-$\chi^{2}_{\rm spec}$ according to Equation 9. The minimum of the
pseudo-$\chi^{2}_{\rm spec}$ is clearly centred around the position of model
#10001 in the Kiel diagram. The latter, however, is expected since it directly
influences the pseudo-$\chi^{2}_{\rm spec}$ value. To this point, we assumed
that the spectroscopic input parameters for the modelling approach were
correct. When applied to real data, the spectroscopic input parameters will
most likely not exactly be the correct physical values but (hopefully) within
the given uncertainty.

To find a preferable method of modelling pre-main sequence $\delta$ Scuti
stars, we need to compare their performance in reproducing the stellar
parameters and, in an optimal case, as much input physics as possible. For a
method to be applicable, it needs to produce the correct parameters even if
the spectroscopic parameters are not equivalent to the physical properties of
the star. Hence, we also draw the expectation value from a Gaussian
distribution around the theoretical model’s value with a standard deviation
$200$ K.

We applied the five approaches discussed in Section 3 to model #10001. The
resulting stellar parameters are given in Table 2. If we assume the
expectation value of the spectroscopic uncertainty box to be correct, only
approach 5 exclusively gives correct stellar parameters; that is, the value of
the model is within the uncertainty associated with the value extracted by the
method. All other approaches disagree for at least one parameter, with
approach 2 performing second best and approach 1 performing worst. However,
the expectation value of approach 5 is in general not closer to the model
parameter, but the associated is uncertainty larger, rendering the results
more stable for this model.

The expectation value of the spectroscopic parameters are most likely not
correct in practice. Hence, we investigated the influence of the spectroscopic
uncertainty box by drawing a temperature and surface gravity offset from a
normal distribution centred around zero with standard deviations of $200$ K
and $0.1$ and applying the methods again. The results, given in Table 2, draw
a similar picture. Methods 2 and 5 perform the best, obtaining the correct
value for 86 and 91% of the stellar parameters. The other approaches fall
behind with 29, 63, and 34% for methods 1, 3, and 5, respectively. We repeated
the same analysis for a Rayleigh limit corresponding to a light curve duration
of $108$ days (four sectors of TESS observations, which are currently the
maximum time base available for a pre-main-sequence $\delta$ Scuti star). The
results are shown in Table 5 and are very similar. Again, methods 2 and 5 far
outperform methods 1 and 4. Method 3 is the only one to perform better for the
assumed $108-$ day precision than the $27-$ day precision. Nevertheless, it
trails the performance of methods 2 and 5.

Table 2: Extracted stellar parameters of model #10001 for the different
approaches (App.) with a Rayleigh limit of a 27-day light curve.

| model | App. 1 | App. 2 | App. 3 | App. 4 | App. 5
---|---|---|---|---|---|---
| #10001 | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}_{\rm spec}$ | $\chi^{2}_{\rm spec}$
|  |  | $3\sigma$ | $1\sigma$ | $\leq 1$ | $\leq 28.87$
spectroscopic parameters assumed correct
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.09(20) | 1.83(8) | 1.78(6) | 1.77(6) | 1.75(10)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7650(210) | 7470(90) | 7430(25) | 7460(140)
$\log\,g$ | 4.208 | 4.224(14) | 4.205(6) | 4.202(5) | 4.201(6) | 4.198(11)
$R$ ($R_{\odot}$) | 1.772 | 1.85(6) | 1.77(2) | 1.75(2) | 1.75(2) | 1.75(4)
$L$ ($L_{\odot}$) | 8.6 | 20(8) | 9.7(1.2) | 8.5(5) | 8.4(2) | 8.5(8)
age (Myr) | 8.12 | 6.3(2.5) | 8.4(1.4) | 8.9(1.3) | 8.9(9) | 8.8(1.5)
Z | 0.018 | 0.018(4) | 0.018(3) | 0.018(3) | 0.018(3) | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7488$ K and $\log\,g=4.275$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.09(20) | 1.83(8) | 1.79(6) | 1.78(5) | 1.77(10)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7700(230) | 7490(110) | 7480(30) | 7500(150)
$\log\,g$ | 4.208 | 4.224(14) | 4.206(6) | 4.202(5) | 4.202(4) | 4.199(11)
$R$ ($R_{\odot}$) | 1.772 | 1.85(6) | 1.77(3) | 1.75(2) | 1.75(2) | 1.75(4)
$L$ ($L_{\odot}$) | 8.6 | 20(8) | 10(1.4) | 8.7(6) | 8.6(2) | 8.7(9)
age (Myr) | 8.12 | 6.3(2.5) | 8.3(1.5) | 8.7(1.3) | 8.7(9) | 8.7(1.4)
Z | 0.018 | 0.018(4) | 0.018(3) | 0.018(3) | 0.018(3) | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7553$ K and $\log\,g=4.127$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.09(20) | 1.84(9) | 1.80(7) | 1.80(5) | 1.78(10)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7740(250) | 7550(110) | 7560(30) | 7560(150)
$\log\,g$ | 4.208 | 4.224(14) | 4.206(6) | 4.203(6) | 4.203(4) | 4.20(1)
$R$ ($R_{\odot}$) | 1.772 | 1.85(6) | 1.77(3) | 1.76(2) | 1.76(2) | 1.75(4)
$L$ ($L_{\odot}$) | 8.6 | 20(8) | 10.2(1.5) | 9.0(7) | 9.1(2) | 9.0(9)
age (Myr) | 8.12 | 6.3(2.5) | 8.1(1.5) | 8.5(1.2) | 8.2(1.0) | 8.5(1.4)
Z | 0.018 | 0.018(4) | 0.018(3) | 0.018(3) | 0.018(3) | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7261$ K and $\log\,g=4.145$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.09(20) | 1.80(7) | 1.77(6) | 1.76(5) | 1.74(9)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7580(160) | 7390(60) | 7280(30) | 7330(130)
$\log\,g$ | 4.208 | 4.224(14) | 4.203(6) | 4.201(6) | 4.201(4) | 4.196(11)
$R$ ($R_{\odot}$) | 1.772 | 1.85(6) | 1.76(2) | 1.75(2) | 1.745(15) | 1.74(4)
$L$ ($L_{\odot}$) | 8.6 | 20(8) | 9.2(9) | 8.1(3) | 7.6(2) | 7.9(7)
age (Myr) | 8.12 | 6.3(2.5) | 8.5(1.5) | 9.2(1.3) | 10.0(1.2) | 9.2(1.5)
Z | 0.018 | 0.018(4) | 0.018(3) | 0.018(3) | 0.019(2) | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7368$ K and $\log\,g=4.17$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.09(20) | 1.82(7) | 1.78(6) | 1.77(5) | 1.75(9)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7640(200) | 7430(80) | 7380(30) | 7410(140)
$\log\,g$ | 4.208 | 4.224(14) | 4.204(6) | 4.202(5) | 4.201(5) | 4.197(11)
$R$ ($R_{\odot}$) | 1.772 | 1.85(6) | 1.76(2) | 1.75(2) | 1.745(18) | 1.74(4)
$L$ ($L_{\odot}$) | 8.6 | 20(8) | 9.6(1.2) | 8.4(4) | 8.1(2) | 8.2(7)
age (Myr) | 8.12 | 6.3(2.5) | 8.4(1.5) | 9.0(1.2) | 9.4(1.4) | 9.0(1.5)
Z | 0.018 | 0.018(4) | 0.018(3) | 0.018(3) | 0.02(3) | 0.018(4)
correct percentage
% | - | 29 | 86 | 63 | 34 | 91

222The expectation value of the extracted parameters adopted as the mean of
all chosen models. Only for the metallicity, Z, we adopted the median as the
expectation value since the grid is evenly spaced. Values in parentheses give
the uncertainties on the extracted parameters adopted from the standard
deviation of all models. Parameters that agree with the model are given in
bold.

###  4.2 Result for the model grid

Approaches 2 and 5 performed best for the modelling of model #10001. We
continued to apply all approaches to more different models. To do so, we drew
75 random echelle diagrams of $l=0$ and $l=1$ modes from our model grid,
perturbed the frequencies with a Rayleigh limit corresponding to light curves
of $27$ and $108$ days, drew random spectroscopic offsets as we did before,
and applied all five modelling approaches. Since we are interested in pre-
main-sequence $\delta$ Scuti stars, we limited the age of the models to be
below $20$ Myr. Table 3 presents the results of this exercise as the
percentage of stellar parameters that were extracted within $1\sigma$ of the
stellar model’s value. Table 6 shows the results for $3\sigma,$ and Table 7
presents the usual uncertainty for each approach as the median of all 75
models. Approaches 2 and 5 again perform best, providing the correct stellar
parameter within $1\sigma$ in 75% (76%) and 74% (72%) of the cases for a
Rayleigh limit corresponding to a light curve of 27 (108) days. These
approaches also yield the correct stellar parameters within $3\sigma$ in
almost all cases. Interestingly, approach 1, which performed worst in the case
of model #10001, does not fall far behind. However, we want to stress that the
large amount of models often results in a very large uncertainty on the
extracted stellar parameters. Hence, in general we discourage the use of this
approach if spectroscopic parameters are available. If no spectroscopic
parameters are available, as is often the case for stars in not-well-studied
clusters or field stars, this result indicates that meaningful stellar
parameters can still be extracted. The reasonable success of approach 3 for
the model #10001 can in general not be reproduced. Only about 50% of the
extracted stellar parameters coincide with the model’s parameters. This is,
however, not a consequence of the small uncertainties, but rather the
stringent spectroscopic constraints being applied. The same argument holds for
approach 4, since it requires a $\chi^{2}_{\rm spec}\leq 1$ and hence an even
stronger constraint on the spectroscopic parameters.

Our results for the model grid verify the results of model #10001 in the sense
that approaches 2 and 5 are to be preferred. The resulting uncertainties (see
Table 7) are mostly comparable, with the exception that the effective
temperature and luminosity are slightly better constrained for approach 5
while the surface gravity and age are slightly better constrained for approach
2. Since the cut-off value for approach 5 is dependent on the number of
extracted pulsation frequencies and hence keeps the number of free parameters
in mind, this seems to be the most reliable method. However, the
implementation of any of these methods is similar enough that we encourage
both methods, that is, approaches 2 and 5, be applied.

Table 3: Percentage rate of extracted stellar parameters that agree with the model’s values within $1\sigma$. |  $\chi^{2}$ |  $\chi^{2}$ |  $\chi^{2}$ |  $\chi^{2}_{\rm spec}$ | $\chi^{2}_{\rm spec}$
---|---|---|---|---|---
|  |  $3\sigma$ |  $1\sigma$ |  $\leq 1$ | $\leq 28.87$
|  % |  % |  % |  % | %
Rayleigh limit corresponding to 27 day light curve
$M_{\star}$ |  63 |  76 |  49 |  69 | 76
$T_{\rm eff}$ |  73 |  83 |  37 |  16 | 64
$\log\,g$ |  63 |  70 |  46 |  63 | 81
$R$ |  69 |  71 |  50 |  64 | 81
$L$ |  74 |  81 |  43 |  31 | 70
age |  71 |  77 |  57 |  61 | 84
Z |  64 |  71 |  59 |  71 | 64
total |  68 |  75 |  49 |  54 | 74
Rayleigh limit corresponding to 108 day light curve
$M_{\star}$ ($M_{\odot}$) |  71 |  77 |  49 |  47 | 74
$T_{\rm eff}$ (K) |  71 |  81 |  39 |  13 | 67
$\log\,g$ |  69 |  67 |  43 |  46 | 79
$R$ ($R_{\odot}$) |  69 |  73 |  44 |  40 | 79
$L$ ($L_{\odot}$) |  74 |  81 |  43 |  16 | 67
age (Myr) |  74 |  79 |  46 |  40 | 71
Z |  70 |  77 |  56 |  64 | 64
total |  71 |  76 |  48 |  38 | 72

###  4.3 Constraining the input physics of accretion and other free parameters

The models calculated in this work include free parameters that govern the
input physics of accretion. These include the e-folding timescale, the amount
of injected heat, the region of injected heat, and the mass-accretion rate.
According to our results, echelle diagrams seem to be unresponsive to these
parameters; that is, we cannot use the pulsation frequencies to constrain, for
example, the mass-accretion rate, since in most cases multiple models with
widely different mass accretion rates produce similar pulsation spectra.
However, the different input physics of accretion leads to a better estimate
of the early evolution of the star in question, and hence contributes to the
distribution of stellar models that provide similar theoretical pulsation
spectra.

A similar discussion is possible for the other free parameters included in our
model grid: the mixing length of convection, envelope mixing, and convective
overshoot parameter. Our results again show that, for the Rayleigh limits
studied in this work, echelle diagrams are not able to differentiate between
the models. However, similarly to the input physics of accretion, these
parameters will add to the distribution of models agreeing with the observed
frequencies and hence provide a more complete picture of possible evolutionary
paths.

Figure 5: Peterson diagram for the fundamental mode. Left panel: Radial modes
($l=0$). Middle panel: Dipolar modes ($l=1$). Right panel: Quadrupole modes
($l=2$).  Figure 6: Example output for automatic echelle diagram. The small
grey crosses indicate all frequencies given to the algorithm as input. The
open circles (triangles) are the radial (dipole) frequencies found by the
algorithm, and the coloured circles (triangles) indicate the pulsation modes
of one model that agrees with the echelle diagram. Again, the echelle diagram
is repeated to the right of the dashed black line for better visibility.

##  5 Automatic echelle diagrams from observed pulsation frequencies

Amplitude spectra of $\delta$ Scuti stars often show plenty of significant
frequencies, such that the extraction of echelle diagrams is mostly
impossible. Bedding et al. (2020) showed regularities in young main-sequence
$\delta$ Scuti stars with echelle diagrams similar to the one shown in Figure
2. From thousands of $\delta$ Scuti stars observed with TESS and Kepler, they
were able to find regularities in the amplitude spectra for only $60$ of them.
The dense amplitude spectra of many $\delta$ Scuti stars, most probably due to
rotational splitting of non-radial pulsation modes, make it difficult to find
regularities in the form of echelle diagrams. Here, we present an automatic
approach, based on the period ratios of the pulsation modes.

It is widely accepted that the ratio between the first overtone and the
fundamental mode in radial pulsation modes of $\delta$ Scuti stars show a
ratio of $0.77$ (Breger, 1979). The ratios of many more combinations of
different pulsation modes show similar distinct values, to the extent that
diagrams showing such ratios have become known as Peterson diagrams (Petersen,
1973; Andreasen et al., 1983; Andreasen, 1983; Aerts et al., 2010), which show
the period ratio of two modes on the y-axis and the period of the pulsation
mode with lower radial order on the x-axis. The latter have been used to
identify the modes of $\delta$ Scuti stars (i.e. Netzel et al., 2022), and we
propose a similar technique to automatically find echelle diagrams of pre-
main-sequence $\delta$ Scuti stars. The principle idea is rather simple. We
take two significant frequencies and compare their position in the Peterson
diagram with the location of period ratios from theoretical models. For the
latter, we used all unstable modes in the $\delta$ Scuti instability region
from Steindl et al. (2021a). Figure 5 shows the resulting boundaries for the
fundamental radial, dipole, and quadrupole modes. We used this to search
echelle ridges including the fundamental mode (F) and first seven overtones
(O1, O2, O3, O4, O5, O6, O7). Our approach to finding ridges is based on the
following steps.

We start the search for possible ridges by choosing all frequencies that fall
within the period range expected for the radial fundamental mode
($0.0425-0.18$ d). Then, we calculate the period ratios with all other
frequencies and check whether the position in the Petersen diagram lies within
a boundary of the first seven overtones. If so, the frequency is assigned to
the mode of the boundary and added to the ridge. If at least four different
frequencies corresponding to different modes are found, we keep the ridge for
further analysis. In the case that multiple frequencies have been assigned to
the same overtone mode, we split the ridges in such a way that there are two
distinct ridges without multiple frequencies assigned to an overtone mode.
This procedure is repeated for dipole and quadrupole modes, and for searches
starting at the first, second, third, and fourth overtone. This first step
results in a list of possible ridges with mode degree and mode identification.

In a second step, we look to verify these ridges by checking every frequency
combination within the ridge, whether the resulting period ratio is within the
expected boundary in the Petersen diagram. For example, if the ridge contains
F, O1, O3, and O4, we check the combinations O1/O3, O1/O4, O3/O4. If more than
70% (or 85% for l=2 modes) of the total period ratios fall within the
boundary, the ridge is kept for the next step, otherwise removed from the
analysis.

Thirdly, we perform a model search by calculating the pseudo-$\chi^{2}$ value
from Equation 8 for every model in our model grid. If any model results in
$\chi^{2}\leq 1$, we keep the ridge. Otherwise, we remove the ridge from the
list.

Finally, we combine each ridge with any other ridges with different mode
degree and look for possible models explaining both ridges, that is
$\chi^{2}\leq 1$ when combining all frequencies of both ridges.

This method finds echelle diagrams from a list of significant frequencies that
agree with predictions from pre-main-sequence $\delta$ Scuti models. However,
care must be taken when interpreting the results, since it is possible that
ridges are built by chance, especially for very dense amplitude spectra. A few
other things have to be kept in mind. It is possible that what has been found
is not the full ridge; that is, a frequency might be missing because it fell
just outside the boundaries in the Petersen diagram. However, by comparing the
incomplete ridge with the corresponding model frequencies, it is often simple
to identify the missing frequencies. The constraints on the period ratios for
degree l=2 modes is less confined than for degree l=0 and l=1 modes (see
Figures 5 and 9, especially). As a consequence, for frequencies with a period
of about $0.06$ d, almost all values between $\sim 0.35$ and $\sim 0.8$ are
included within a boundary for $l=2$ modes (see Figure 5). Therefore, one
finds many more $l=2$ ridges; however, most of them do not correspond to real
echelle diagrams. Additionally, the amplitudes of $l=2$ modes are expected to
be smaller than radial and dipole modes anyway (Aerts et al., 2010), such that
finding an echelle ridge should be less probable. Splitting the frequencies in
step 1 can lead to an overwhelming number of different ridges, especially in
dense amplitude spectra. In such cases, it is best to only use frequencies
with higher significance, that is, a signal-to-noise ratio above six instead
of four, as is usual in asteroseismology (Breger et al., 1993; Kuschnig et
al., 1997).

Our approach to automatically finding echelle diagrams from the observed
frequencies needs to be verified in each case. The simplest approach is to
compare any spectroscopic parameters available for the star, with the
positions of the models with $\chi^{2}\leq 1$ in the Kiel diagram. We verified
that this approach works by applying it to frequency lists from our model
grid. Here, we only show the application to frequencies of model #10001, but
the results were similar for many more models. We applied the algorithm
multiple times, each time taking F, O1, O3, O4, and O5 from the radial modes
and O1, O2, O3, O4, and O6 from the dipole modes and perturbing them with a
Rayleigh limit corresponding to a 27-day light curve. We added twenty randomly
chosen frequencies in the range from $20$ to $50$ $\mathrm{d}^{-1}$. An
example output is given in Figure 6, clearly showing that the algorithm
extracted the correct echelle diagram. It remains to be seen how this approach
operates for the pulsation frequencies of real stars. We present an example in
Section 6, but the application to more stars will be the topic of future work.

Figure 7: Dependence of stellar parameters of HD 139614 on the choice of
initial helium abundance. The black points show the resulting parameters for
different values of initial helium abundance when we used the surface gravity
as a second classical constraint. The grey counterpart corresponds to the
luminosity as second classical constraint. Error bars give the resulting
standard deviation. The green area shows our model result, where we adopted
the values from approach 5 and the red symbols correspond to the result of
Murphy et al. (2021).  Figure 8: Possible evolutionary history for HD 139614.
The figure shows the evolutionary tracks (black lines) of all evolutionary
models of the extended grid that provided at least one pulsation model with
$\chi^{2}_{\rm spec}<22.36$. The $1\sigma$ spectroscopic error box for HD
139614 is shown in magenta.

##  6 Re-discussion of the pre-main-sequence $\delta$ Scuti star HD 139614

###  6.1 Parameter estimation of HD 139614

Murphy et al. (2021) discussed the pre-main-sequence star HD 139614 in detail.
They extracted frequencies from the TESS light curve and essentially applied
modelling approach 3 to infer the stellar parameters of HD 139614. A major
difference between the analysis by Murphy et al. (2021) and our analysis
presented in the following are the adopted stellar models. Murphy et al.
(2021) used the classical description of pre-main-sequence models. Their model
grid included different values for the initial mass, mixing length, and
metallicity. In this work, we employed more advanced models of the pre-main-
sequence evolution, starting from accreting protostellar seeds that include
many more free parameters (see Section 2), which leads to a better estimate of
the possible evolutionary past of HD 139614. There is also a difference in the
adopted value for the initial helium abundance. While Murphy et al. (2021)
adopted an inital helium abundance of $X_{{}^{4}\mathrm{He}}=0.29$, our values
lie in the $X_{{}^{4}\mathrm{He}}=0.216-0.282$ range. An additional difference
is that the authors used the stellar luminosity to constrain the spectroscopic
parameters, while we used the surface gravity reported by Fairlamb et al.
(2015). Hence, the spectroscopic parameters for HD 139614 used in our study
are $\log\,g=4.31(12)$(Fairlamb et al., 2015) and $T_{\rm eff}=7650(200)$
(Murphy et al., 2021).

We used the frequencies given by Murphy et al. (2021) and our method to
automatically find echelle diagrams. We indeed find the same mode
identification as Murphy et al. (2021) including all other possible
identifications mentioned. Our automatic approach additionally yielded their
F14 as a fundamental radial mode, a possibility discussed by the authors but
ultimately declined as they do not seem to have had corresponding stellar
models. We used the same identification as Murphy et al. (2021) for the
subsequent modelling. This ensures that differences in resulting stellar
parameters are independent of any mode identification ambiguity but are the
consequence of an improved stellar modelling. In addition, we later include
this fundamental radial mode and discuss the influence on the resulting
stellar parameters.

Since Murphy et al. (2021) reported a mass of $1.520(18)\,M_{\odot}$, we
extended our model grid to lower stellar masses. We calculated an additional
300 models, with final masses between $1.2$ and $1.5\,M_{\odot}$, effectively
extending our model grid. Figure 13 shows an updated version of Figure 1. We
performed the modelling with approaches 2, 3, and 5, since approaches 2 and 5
have proven to be the most reliable, and approach 3 is comparable to that
applied by Murphy et al. (2021). The evolutionary calculations of our model
grid include the earliest parts of the main sequence. Some of these show
theoretical pulsation frequencies in agreement with the echelle diagram shown
by HD 139614, but they are too old to represent the evolutionary stage of the
star. The few too far evolved models will influence the calculation of the
mean and standard deviation needed for the extraction of the stellar
parameters, with the effect being strongest for the age itself. As such, we
decided to apply an age cut-off of 20 Myr333We experimented with different age
cut-offs, looking both at the resulting stellar parameters and the
age-$\chi^{2}$ distribution. While lower cut-offs change the resulting
parameters only to a small degree, and within the uncertainties, higher values
lead to an asymmetric age-$\chi^{2}$ distribution and hence to changes in the
extracted stellar parameters and their uncertainties. Hence, we believe 20 Myr
to be the best choice. In addition, 20 Myr seems to be a reasonable upper
limit since HD 139614 is located in the Upper Centaurus Lupus association. and
only use models with ages below that.

Table 4 shows the resulting stellar parameters extracted from the different
modelling approaches. As is evident, the values differ significantly from the
values reported by Murphy et al. (2021), but they are essentially the same
independently of whether or not we include the radial fundamental mode. In
addition, it is evident that the extracted uncertainties are larger than the
ones reported by Murphy et al. (2021). We obtain the largest difference for
the stellar mass, where our modelling yields a value about $0.15\,M_{\odot}$
larger, depending on the modelling approach.

###  6.2 Model comparison

We then conducted a comparison between our models and the models presented by
Murphy et al. (2021) with the aim of illustrating the differences. First, we
investigated the effect of using a fixed helium abundance in the models. For
this, we calculated a grid of classical pre-main-sequence models, followed the
modelling approach by Murphy et al. (2021), and applied different values for
the initial helium abundance.

As is evident in Figure 7, the stellar mass and stellar age are most sensitive
to this particular choice of input physics. Changing the helium abundance from
$Y=0.29$ to $Y=0.25$ leads to a stellar mass for HD 139614 of around
$1.62\,M_{\odot},$ which is already more than $5\,\sigma$ from the value
reported in Murphy et al. (2021).

It is important to note that the initial helium abundance is not constrained
in the echelle diagram; that is, the models result in a similar goodness of
fit. In the same manner, the other input physics (mixing length, e-folding
timescale, amount of injected heat, region of injected heat, mass accretion
rate, envelope mixing, and overshooting) are not constrained but lead to
different stellar evolutions, all producing pulsation frequencies that fit the
observed echelle diagram. This result proves that next to the mass and
metallicity, the initial helium abundance plays a major role in the
determination of stellar parameters. Our model grid additionally allows for
variations in envelope mixing and overshooting coefficients. The results of
our model grid furthermore indicate that overshooting is an additional
important parameter, while the influence of envelope mixing is significantly
smaller. In comparison, given the current observational constraints on the
pulsation frequencies of HD 139614, the accretion parameters seem to have
little importance on the estimated stellar parameters, but they are necessary
if additional information about possible previous evolutionary pathways is
desired. The inclusion of these parameters allowed us to produce Figure 8,
which shows the possible evolutionary past of HD 139614 in the constant
accretion scenario. Evidently, this provides a much more realistic picture
than the classical evolutionary concept would let us assume. With HD 139614
still holding onto a protoplanetary disc, this allows for a possible
investigation using additional constraints from disc parameters, which shall
be the focus of future work.

To further compare the models applied to HD 139614 by Murphy et al. (2021) and
in this work, we applied Akaike’s information criterion (AIC) and the Bayesian
information criterion (BIC). Both are information criteria commonly used for
model comparison if the calculation of evidence in the form of likelihood
ratios is not possible (Trotta, 2008; Aerts et al., 2018), as in the case of
asteroseismic modelling (Aerts et al., 2018). The two information criteria are
defined as

$\displaystyle\text{AIC}=2k-2\ln\mathcal{L}_{max,}$
$\displaystyle\text{BIC}=-2\ln\mathcal{L}_{max}+k\ln n\,,$

where $\mathcal{L}_{max}$ is a maximum likelihood estimator, $n$ is the
effective number of observations, and $k$ is the number of free parameters in
the model. Both are only relevant in comparison, with the basic picture being
’the lower the better’ (Aerts et al., 2018).

Following Degroote et al. (2009), the AIC and BIC, under the assumption of
Gaussian white noise, reduce to

$\displaystyle\text{AIC}=n\ln(\text{RSS}/n)+2k+n,$
$\displaystyle\text{BIC}=n\ln(\text{RSS}/n)+k\ln n\,,$

where RSS is the residual sum of squares, which in our particular case can
easily be calculated from the $\chi^{2}$ value of the best model. We derived
the AIC and BIC for our modelling (using the best $\chi^{2}$ value of $0.107$)
and the modelling by Murphy et al. (2021) (using their best $\chi^{2}$ value
of $0.3864$) to be

$\displaystyle\text{AIC}_{\rm this\,work}=-106.9\quad\quad\text{AIC}_{\rm
Murphy}=-100.5\quad\quad\Delta\text{AIC}=-6.4,$ $\displaystyle\text{BIC}_{\rm
this\,work}=-103.3\quad\quad\text{BIC}_{\rm
Murphy}=-107.1\quad\quad\Delta\text{BIC}=2.8\,.$

Hence, the AIC provides strong evidence that our model should be preferred,
while the BIC provides positive evidence that the analysis of Murphy et al.
(2021) should be preferred (Aerts et al., 2018). It is important to state that
the BIC penalises models with more free parameters stronger than the AIC
(Trotta, 2008). Additionally, other than, for example, Bayesian evidence, the
information criteria apply the penalty to extra parameters regardless of
whether they are constrained (Trotta, 2008), which is not the case for our
ten-dimensional model when applied to HD 139614. As a consequence, the BIC
might be overly harsh with regard to our model.

We investigated this further and calculated an additional model grid, for
which we fixed the accretion parameters to these values: accretion rate
$\dot{M}_{0}=5\cdot 10^{-6},M_{\odot}/{\rm yr}$, e-folding time $\tau=0.5$
Myr, fractional energy $\beta$= $0.1$, and fractional mass $M_{\rm outer}$ =
$0.1$. As discussed above, these parameters are not constrained by the echelle
diagram in combination with the observational uncertainty, and hence they
impose a strong penalty on the calculation of the AIC and especially the BIC.
For this model grid, we calculated a total of 500 stellar evolution models,
leading to roughly 50 000 theoretical pulsation spectra. Applying this six-
dimensional model grid to HD 139614 with modelling approaches 2, 3, and 5
leads to very similar stellar parameters (see Table 8) to the ten-dimensional
model grid (see Table 4), with only small changes in the expectation values
and the uncertainties. The best $\chi^{2}$ of $0.107$ is also very close to
the value obtained from the ten-dimensional model grid. Comparison of the
information criteria now results in the following:

$\displaystyle\text{AIC}_{\rm 6\,d}=-114.9\quad\quad\text{AIC}_{\rm
Murphy}=-100.5\quad\quad\Delta\text{AIC}=-14.4,$ $\displaystyle\text{BIC}_{\rm
6\,d}=-117.1\quad\quad\text{BIC}_{\rm
Murphy}=-107.1\quad\quad\Delta\text{BIC}=-10.0\,.$

Hence, the AIC provides very strong evidence and the BIC strong evidence (with
$\Delta\text{BIC}=-10$ as the limit between strong and very strong evidence
(Aerts et al., 2018)) for this model to be preferred over the model by Murphy
et al. (2021).

We illustrate the influence of missing essential ingredients in the model
grids used to perform asteroseismic modelling of pre-main-sequence $\delta$
Scuti stars. When a grid consists of too many fixed parameters, the resulting
fundamental stellar properties will have underestimated uncertainties and
possibly inaccurate expectation values. Asteroseismology is a powerful tool,
for example, in the description of exoplanets and their host stars. In this
regard, it is crucial that the parameters estimated for the host stars be
reliable, as any false information will propagate through to the estimated
parameters of the planets.

Table 4: Extracted stellar parameters of HD 139614 for the different
approaches and the two discussed versions of the echelle diagram.

| Murphy et al. (2021) | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}_{\rm spec}$
---|---|---|---|---
| modelling | $3\sigma$ | $1\sigma$ | $\leq 22.36$
same mode identification as in Murphy et al. (2021)
$M_{\star}$ ($M_{\odot}$) | 1.520(18) | 1.7(1) | 1.66(6) | 1.66(7)
$T_{\rm eff}$ (K) | 7615(64) | 7830(250) | 7670(110) | 7650(130)
$\log\,g$ | - | 4.288(8) | 4.285(5) | 4.284(8)
$R$ ($R_{\odot}$) | - | 1.55(3) | 1.53(2) | 1.54(3)
$L$ ($L_{\odot}$) | 6.7(3) | 8.2(1.3) | 7.4(5) | 7.3(6)
age (Myr) | 10.75(77) | 12(3) | 12(3) | 12(3)
Z | 0.010(1) | 0.014(5) | 0.012(4) | 0.014(5)
| Murphy et al. (2021) | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}_{\rm spec}$
| modelling | $3\sigma$ | $1\sigma$ | $\leq 23.68$
additionally include the fundamental radial mode
$M_{\star}$ ($M_{\odot}$) | 1.520(18) | 1.71(9) | 1.66(6) | 1.66(7)
$T_{\rm eff}$ (K) | 7615(64) | 7830(250) | 7680(110) | 7660(130)
$\log\,g$ | - | 4.289(7) | 4.286(5) | 4.284(8)
$R$ ($R_{\odot}$) | - | 1.55(3) | 1.54(2) | 1.54(3)
$L$ ($L_{\odot}$) | 6.7(3) | 8.2(1.3) | 7.4(5) | 7.3(6)
age (Myr) | 10.75(77) | 12(3) | 12(3) | 12(3))
Z | 0.010(1) | 0.014(4) | 0.012(3) | 0.014(5)

444The expectation value of the extracted parameter adopted as the mean of all
chosen models. Only for the metallicity, Z, did we adopt the median as the
expectation value since the grid is evenly spaced. Values in parentheses give
the uncertainty on the extracted parameter adopted from the standard deviation
of all models.

##  7 Conclusion

In this work, we investigated different approaches in the pulsational
modelling of pre-main-sequence $\delta$ Scuti stars. For this, we calculated a
grid of stellar evolution models starting from the protostellar accretion
phase, including various free parameters such as the final mass, initial
composition, mixing parameter, and details of the accretion process. The
resulting equilibrium models of pre-main-sequence stars are used as input for
adiabatic stellar pulsation calculations to deliver a grid of pulsation
spectra. We introduced a total of five different modelling approaches all
based on a pseudo-$\chi^{2}$ statistic as a goodness of fit indicator for the
pulsation frequencies. While the first approach uses a simple
pseudo-$\chi^{2}$ statistic, the other approaches additionally make use of
spectroscopic constraints. Approaches 2 and 3 only consider models that fall
within a given spectroscopic error box, while approaches 4 and 5 include the
spectroscopic constraints within the pseudo-$\chi^{2}$ statistic itself. All
different modelling techniques are first applied to a single theoretical model
and then to a set of 75 randomly chosen models from our model grid. During
this endeavour, we include frequency shifts according to a Rayleigh limit
corresponding to light curves with time bases of 27 and 108 days. In addition,
we draw a random offset for the the expectation value of the spectroscopic
parameters from a Gaussian distribution to insure the inclusion of measurement
uncertainties in these aspects as well.

Our results clearly indicate that two approaches (i.e. 2 and 5) are to be
preferred in the modelling of pre-main-sequence $\delta$ Scuti stars. Other
approaches often result in unreliable stellar parameters (approaches 3 and 4)
or in very large uncertainties (approach 1). According to our results, we
conclude that modelling approach 5, which included the spectroscopic
parameters in the calculation of the pseudo-$\chi^{2}$ statistic and sets the
pseudo-$\chi^{2}$ cut-off depending on the amount of free parameters, is the
most reliable. We note that approach 2, which uses the simpler
pseudo-$\chi^{2}$ statistics and only considers models within $3\sigma$ of the
spectroscopic constraints, does not perform worse. The work presented here
demonstrates that typical uncertainties for the stellar radius are in the
range of $0.05\,R_{\odot}$. For lots of stars, we expect resulting
uncertainties to be smaller than this since this is the median value within
all studied models.

In addition, we introduced a simple algorithm for the extraction of echelle
diagrams in the observed pulsation spectra of pre-main-sequence $\delta$ Scuti
stars. This algorithm is based on theoretical Petersen diagrams taken from the
non adiabatic pre-main-sequence models of Steindl et al. (2021a). We show that
the algorithm correctly identifies the echelle ridge for perturbed frequencies
of our modelling grid.

Finally, we combine the presented work and apply it to rediscuss the
asteroseismic modelling of the pre-main-sequence $\delta$ Scuti star HD
139614. This star was recently modelled by Murphy et al. (2021), but the
applied model grid used an outdated description of pre-main-sequence models
and a less dimensional model grid. Our results indicate that this simplified
model grid is insufficient to properly describe the complete picture of the
evolutionary pathways of HD 139614. Hence, our study results in different
stellar parameters, which most notably have higher uncertainties. At this
point, it is important to stress that these uncertainties stem from an
improved grid of stellar models and an improved modelling approach and should
hence depict the current limitations of asteroseismic modelling of pre-main-
sequence stars. Additionally, we repeated the analysis of Murphy et al. (2021)
with classical model grids to show the dependence of the results on the
initial helium composition. Most notably, this explains the difference in
extracted stellar mass for HD 139614 and results in far greater uncertainty
with regard to stellar age than mentioned by Murphy et al. (2021).We performed
a model comparison by applying commonly used information criteria to
illustrate the differences between our approach and the one by Murphy et al.
(2021).

This work clearly shows that producing an extensive grid of stellar models to
perform this kind of asteroseismic modelling of pre-main-sequence $\delta$
Scuti stars is unavoidable. If, in contrast, simpler model grids are utilised,
resulting parameters are often extracted with incorrect expectation values and
underestimated uncertainties.

###### Acknowledgements.

We are grateful to Bill Paxton and his collaborators for their valuable work
on the stellar evolution code MESA. This research has made use of matplotlib,
a Python library for publication quality graphics (Hunter, 2007); SciPy
(Virtanen et al., 2020); NumPy (Van Der Walt et al., 2011); MESA SDK for Linux
(Version 20.3.1) (Townsend, 2020).

## References

  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology
  * Aerts et al. (2018) Aerts, C., Molenberghs, G., Michielsen, M., et al. 2018, ApJS, 237, 15
  * Andreasen (1983) Andreasen, G. K. 1983, A&A, 121, 250
  * Andreasen et al. (1983) Andreasen, G. K., Hejlesen, P. M., & Petersen, J. O. 1983, A&A, 121, 241
  * Angulo et al. (1999) Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A, 656, 3
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Baraffe & Chabrier (2010) Baraffe, I. & Chabrier, G. 2010, A&A, 521, A44
  * Baraffe et al. (2009) Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJ, 702, L27
  * Baraffe et al. (2012) Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118
  * Bedding et al. (2020) Bedding, T. R., Murphy, S. J., Hey, D. R., et al. 2020, Nature, 581, 147
  * Bernabei et al. (2009) Bernabei, S., Ripepi, V., Ruoppo, A., et al. 2009, A&A, 501, 279
  * Breger (1979) Breger, M. 1979, PASP, 91, 5
  * Breger et al. (1993) Breger, M., Stich, J., Garrido, R., et al. 1993, A&A, 271, 482
  * Buchler & Yueh (1976) Buchler, J. R. & Yueh, W. R. 1976, ApJ, 210, 440
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094
  * Chen & Li (2019) Chen, X. & Li, Y. 2019, ApJ, 872, 156
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028
  * Cox & Giuli (1968) Cox, J. P. & Giuli, R. T. 1968, Principles of stellar structure
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240
  * Degroote et al. (2009) Degroote, P., Briquet, M., Catala, C., et al. 2009, A&A, 506, 111
  * Dupret et al. (2005) Dupret, M. A., Grigahcène, A., Garrido, R., Gabriel, M., & Scuflaire, R. 2005, A&A, 435, 927
  * Dziembowski (1971) Dziembowski, W. A. 1971, Acta Astron., 21, 289
  * Eddington (1926) Eddington, A. S. 1926, The Internal Constitution of the Stars
  * Elbakyan et al. (2019) Elbakyan, V. G., Vorobyov, E. I., Rab, C., et al. 2019, MNRAS, 484, 146
  * Fairlamb et al. (2015) Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1
  * Goldstein & Townsend (2020) Goldstein, J. & Townsend, R. H. D. 2020, ApJ, 899, 116
  * Gruber et al. (2012) Gruber, D., Saio, H., Kuschnig, R., et al. 2012, MNRAS, 420, 291
  * Hartmann et al. (1997) Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, ApJ, 475, 770
  * Hayashi (1961) Hayashi, C. 1961, PASJ, 13, 450
  * Henyey et al. (1955) Henyey, L. G., Lelevier, R., & Levée, R. D. 1955, PASP, 67, 154
  * Herwig (2000) Herwig, F. 2000, A&A, 360, 952
  * Hosokawa et al. (2011) Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Iben (1965) Iben, Icko, J. 1965, ApJ, 141, 993
  * Iglesias & Rogers (1993) Iglesias, C. A. & Rogers, F. J. 1993, ApJ, 412, 752
  * Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411
  * Jensen & Haugbølle (2018) Jensen, S. S. & Haugbølle, T. 2018, MNRAS, 474, 1176
  * Kim et al. (2021) Kim, S.-L., Lee, J. W., Lee, C.-U., et al. 2021, AJ, 162, 212
  * Kunitomo et al. (2017) Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A, 599, A49
  * Kuschnig et al. (1997) Kuschnig, R., Weiss, W. W., Gruber, R., Bely, P. Y., & Jenkner, H. 1997, A&A, 328, 544
  * Langanke & Martínez-Pinedo (2000) Langanke, K. & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481
  * Larson (1969) Larson, R. B. 1969, MNRAS, 145, 271
  * Masunaga & Inutsuka (2000) Masunaga, H. & Inutsuka, S.-i. 2000, ApJ, 531, 350
  * Murphy et al. (2021) Murphy, S. J., Joyce, M., Bedding, T. R., White, T. R., & Kama, M. 2021, MNRAS, 502, 1633
  * Netzel et al. (2022) Netzel, H., Pietrukowicz, P., Soszyński, I., & Wrona, M. 2022, MNRAS, 510, 1748
  * Nieva & Przybilla (2012) Nieva, M. F. & Przybilla, N. 2012, A&A, 539, A143
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Petersen (1973) Petersen, J. O. 1973, A&A, 27, 89
  * Pols et al. (1995) Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964
  * Potekhin & Chabrier (2010) Potekhin, A. Y. & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82
  * Przybilla et al. (2013) Przybilla, N., Nieva, M. F., Irrgang, A., & Butler, K. 2013, in EAS Publications Series, Vol. 63, EAS Publications Series, ed. G. Alecian, Y. Lebreton, O. Richard, & G. Vauclair, 13–23
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Rogers & Nayfonov (2002) Rogers, F. J. & Nayfonov, A. 2002, ApJ, 576, 1064
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713
  * Seaton (2005) Seaton, M. J. 2005, MNRAS, 362, L1
  * Stahler et al. (1980) Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637
  * Steindl et al. (2021a) Steindl, T., Zwintz, K., Barnes, T. G., Müllner, M., & Vorobyov, E. I. 2021a, A&A, 654, A36
  * Steindl et al. (2021b) Steindl, T., Zwintz, K., & Bowman, D. M. 2021b, A&A, 645, A119
  * Timmes & Swesty (2000) Timmes, F. X. & Swesty, F. D. 2000, ApJS, 126, 501
  * Townsend (2020) Townsend, R. 2020, MESA SDK for Linux
  * Townsend et al. (2018) Townsend, R. H. D., Goldstein, J., & Zweibel, E. G. 2018, MNRAS, 475, 879
  * Townsend & Teitler (2013) Townsend, R. H. D. & Teitler, S. A. 2013, MNRAS, 435, 3406
  * Trotta (2008) Trotta, R. 2008, Contemporary Physics, 49, 71
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Vorobyov et al. (2017) Vorobyov, E. I., Elbakyan, V., Hosokawa, T., et al. 2017, A&A, 605, A77
  * Zwintz et al. (2014) Zwintz, K., Fossati, L., Ryabchikova, T., et al. 2014, Science, 345, 550

##  Appendix A MESA microphysics

The MESA EOS is a blend of the OPAL Rogers & Nayfonov (2002), SCVH Saumon et
al. (1995), PTEH Pols et al. (1995), HELM Timmes & Swesty (2000), and PC
Potekhin & Chabrier (2010) EOSes. Radiative opacities are primarily from OPAL
(Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al.
(2005) and the high-temperature, Compton-scattering-dominated regime by
Buchler & Yueh (1976). Electron conduction opacities are from Cassisi et al.
(2007).

Nuclear reaction rates are a combination of rates from NACRE (Angulo et al.
1999), JINA REACLIB (Cyburt et al. 2010), and additional tabulated weak
reaction rates. (Fuller et al. (1985); Oda et al. (1994); Langanke & Martínez-
Pinedo (2000). Screening is included via the prescription of Chugunov et al.
(2007). Thermal neutrino loss rates are from Itoh et al. (1996).

##  Appendix B Estimating the stellar parameter of theoretical echelle
diagrams

Table 5: Extracted stellar parameters of model #10001 for the different
approaches with a Rayleigh limit of a 108-day light curve.

| #10001 | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}_{\rm spec}$ | $\chi^{2}_{\rm spec}$
---|---|---|---|---|---|---
|  |  | $3\sigma$ | $1\sigma$ | $\leq 1$ | $\leq 28.87$
spectroscopic parameters assumed correct
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.11(18) | 1.85(5) | 1.81(3) | 1.81(2) | 1.76(9)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7680(200) | 7460(80) | 7430(30) | 7480(130)
$\log\,g$ | 4.208 | 4.225(12) | 4.208(4) | 4.205(3) | 4.205(2) | 4.200(8)
$R$ ($R_{\odot}$) | 1.772 | 1.85(5) | 1.77(2) | 1.76(1) | 1.759(6) | 1.75(3)
$L$ ($L_{\odot}$) | 8.6 | 19(7) | 9.9(1.2) | 8.6(4) | 8.5(2) | 8.6(7)
age (Myr) | 8.12 | 5.9(1.6) | 7.8(1.4) | 8.1(8) | 8.4(7) | 8.8(1.4)
Z | 0.018 | 0.018(3) | 0.018(2) | 0.018(2) | 0.020(1) | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7488$ K and $\log\,g=4.275$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.10(18) | 1.86(6) | 1.82(3) | 1.79(4) | 1.77(10)
$T_{\rm eff}$ (K) | 7431 | 8700(700) | 7720(230) | 7500(110) | 7480(30) | 7510(140)
$\log\,g$ | 4.208 | 4.225(12) | 4.208(4) | 4.205(3) | 4.203(3) | 4.200(9)
$R$ ($R_{\odot}$) | 1.772 | 1.85(5) | 1.78(2) | 1.762(11) | 1.75(2) | 1.75(3)
$L$ ($L_{\odot}$) | 8.6 | 20(8) | 10.1(1.3) | 8.9(6) | 8.66(13) | 8.8(8)
age (Myr) | 8.12 | 6.1(1.6) | 7.8(1.4) | 8.2(1.1) | 8.4(4) | 8.7(1.4)
Z | 0.018 | 0.018(3) | 0.018(2) | 0.018(2) | 0.018(1) | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7553$ K and $\log\,g=4.127$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.13(18) | 1.88(8) | 1.83(5) | - | 1.78(10)
$T_{\rm eff}$ (K) | 7431 | 8900(800) | 7770(250) | 7560(120) | - | 7560(140)
$\log\,g$ | 4.208 | 4.227(12) | 4.209(6) | 4.206(4) | - | 4.20(9)
$R$ ($R_{\odot}$) | 1.772 | 1.86(5) | 1.78(2) | 1.767(15) | - | 1.75(3)
$L$ ($L_{\odot}$) | 8.6 | 21(8) | 10.5(1.6) | 9.2(7) | - | 9.0(9)
age (Myr) | 8.12 | 5.6(1.6) | 7.5(1.3) | 8.1(1.1) | - | 8.5(1.4)
Z | 0.018 | 0.018(3) | 0.018(2) | 0.018(1) | - | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7261$ K and $\log\,g=4.145$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.11(18) | 1.83(5) | 1.80(4) | - | 1.74(9)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7580(160) | 7400(40) | - | 7370(110)
$\log\,g$ | 4.208 | 4.225(12) | 4.206(4) | 4.204(3) | - | 4.198(9)
$R$ ($R_{\odot}$) | 1.772 | 1.85(6) | 1.77(2) | 1.76(1) | - | 1.74(3)
$L$ ($L_{\odot}$) | 8.6 | 19(7) | 9.3(9) | 8.3(2) | - | 8.0(6)
age (Myr) | 8.12 | 6.0(1.6) | 8.1(1.4) | 8.6(8) | - | 9.3(1.5)
Z | 0.018 | 0.018(3) | 0.018(2) | 0.020(1) | - | 0.018(4)
spectroscopic parameters assumed $T_{\rm eff}=7368$ K and $\log\,g=4.17$
$M_{\star}$ ($M_{\odot}$) | 1.849 | 2.11(18) | 1.84(5) | 1.80(4) | 1.78(4) | 1.76(9)
$T_{\rm eff}$ (K) | 7431 | 8800(800) | 7640(190) | 7430(60) | 7380(30) | 7440(120)
$\log\,g$ | 4.208 | 4.225(12) | 4.207(4) | 4.204(4) | 4.202(4) | 4.199(8)
$R$ ($R_{\odot}$) | 1.772 | 1.85(6) | 1.77(2) | 1.754(13) | 1.750(13) | 1.74(4)
$L$ ($L_{\odot}$) | 8.6 | 20(8) | 9.7(1.1) | 8.4(3) | 8.2(2) | 8.4(7)
age (Myr) | 8.12 | 6.0(1.6) | 8.1(1.4) | 8.5(7) | 8.6(8) | 9.0(1.4)
Z | 0.018 | 0.018(3) | 0.018(2) | 0.02(2) | 0.02(2) | 0.018(4)
correct percentage
% | - | 20 | 86 | 74 | 20 | 89

555The expectation value of the extracted parameters adopted as the mean of
all chosen models. Only for the metallicity, Z, did we adopt the median as the
expectation value since the grid is evenly spaced. Values in parentheses give
the uncertainties on the extracted parameters adopted from the standard
deviation of all models. Parameters that agree with the model are given in
bold. Values are not given if there are fewer than five models.

Table 6: Percentage rate of extracted stellar parameters that agree with the model’s values within $3\sigma$. |  $\chi^{2}$ |  $\chi^{2}$ |  $\chi^{2}$ |  $\chi^{2}_{\rm spec}$ | $\chi^{2}_{\rm spec}$
---|---|---|---|---|---
|  |  $3\sigma$ |  $1\sigma$ |  $\leq 1$ | $\leq 28.87$
|  % |  % |  % |  % | %
Rayleigh limit corresponding to 27 day light curve
$M_{\star}$ |  100 |  100 |  78 |  94 | 99
$T_{\rm eff}$ |  100 |  100 |  70 |  40 | 97
$\log\,g$ |  100 |  100 |  73 |  93 | 99
$R$ |  100 |  100 |  79 |  93 | 100
$L$ |  100 |  100 |  73 |  63 | 97
age |  97 |  99 |  76 |  94 | 99
Z |  100 |  100 |  77 |  97 | 100
total |  100 |  100 |  75 |  82 | 99
Rayleigh limit corresponding to 108 day light curve
$M_{\star}$ |  100 |  100 |  71 |  83 | 99
$T_{\rm eff}$ |  100 |  100 |  71 |  40 | 97
$\log\,g$ |  100 |  100 |  69 |  79 | 99
$R$ |  100 |  100 |  73 |  81 | 100
$L$ |  100 |  100 |  69 |  51 | 94
age |  100 |  100 |  76 |  80 | 99
Z |  100 |  100 |  77 |  84 | 100
total |  100 |  100 |  72 |  71 | 98
Table 7: Median $1\sigma$ uncertainty on the extracted parameters.

|  $\chi^{2}$ |  $\chi^{2}$ |  $\chi^{2}$ |  $\chi^{2}_{\rm spec}$ | $\chi^{2}_{\rm spec}$
---|---|---|---|---|---
|  |  $3\sigma$ |  $1\sigma$ |  $\leq 1$ | $\leq 28.87$
|  % |  % |  % |  % | %
Rayleigh limit corresponding to 27 day light curve
$M_{\star}$ ($M_{\odot}$) |  0.19 |  0.11 |  0.08 |  0.07 | 0.13
$T_{\rm eff}$ (K) |  700 |  300 |  110 |  28 | 150
$\log\,g$ |  0.014 |  0.008 |  0.006 |  0.005 | 0.015
$R$ ($R_{\odot}$) |  0.06 |  0.04 |  0.03 |  0.0325 | 0.05
$L$ ($L_{\odot}$) |  8 |  3.4 |  1.4 |  0.6 | 2.0
age (Myr) |  1.7 |  1.3 |  0.85 |  0.75 | 1.7
Z |  0.005 |  0.004 |  0.004 |  0.004 | 0.004
Rayleigh limit corresponding to 108 day light curve
$M_{\star}$ ($M_{\odot}$) |  0.13 |  0.08 |  0.06 |  0.05 | 0.11
$T_{\rm eff}$ (K) |  600 |  300 |  110 |  25 | 150
$\log\,g$ |  0.010 |  0.006 |  0.005 |  0.004 | 0.009
$R$ ($R_{\odot}$) |  0.04 |  0.03 |  0.02 |  0.017 | 0.04
$L$ ($L_{\odot}$) |  6.5 |  3.1 |  1.2 |  0.4 | 1.8
age (Myr) |  1.2 |  0.9 |  0.7 |  0.5 | 1.6
Z |  0.004 |  0.003 |  0.003 |  0.003 | 0.005
666We note that for approaches 3 and 4 it is possible that not a single model
fulfils both the stringent spectroscopic and the asteroseismic demands. The
values given here are the medians of the models for which at least two
theoretical pulsation spectra from different evolutionary calculations fit
both constraints.

##  Appendix C Petersen diagrams

Figure 9: Peterson diagram for the first overtone.  Figure 10: Peterson
diagram for the second overtone.  Figure 11: Peterson diagram for the third
overtone.  Figure 12: Peterson diagram for the fourth overtone.

##  Appendix D Re-discussion of HD 139614

Figure 13: Similar to Figure 1, but including the extension for the modelling
of HD 139614. The additional models calculated are added to the Kiel diagram
and stacked upon the histograms in pale green.  Table 8: Extracted stellar
parameters of HD 139614 when applying the six-dimensional model grid.

| Murphy et al. (2021) | $\chi^{2}$ | $\chi^{2}$ | $\chi^{2}_{\rm spec}$
---|---|---|---|---
| modelling | $3\sigma$ | $1\sigma$ | $\leq 22.36$
same mode identification as in Murphy et al. (2021)
$M_{\star}$ ($M_{\odot}$) | 1.520(18) | 1.71(11) | 1.67(9) | 1.66(10)
$T_{\rm eff}$ (K) | 7615(64) | 7840(260) | 7670(120) | 7660(130)
$\log\,g$ | - | 4.289(9) | 4.286(8) | 4.283(9)
$R$ ($R_{\odot}$) | - | 1.55(3) | 1.54(3) | 1.54(4)
$L$ ($L_{\odot}$) | 6.7(3) | 8.2(1.3) | 7.4(6) | 7.4(7)
age (Myr) | 10.75(77) | 12(3) | 13(3) | 12(3)
Z | 0.010(1) | 0.018(5) | 0.016(4) | 0.014(5)

777The expectation value of the extracted parameter adopted as the mean of all
chosen models. Only for the metallicity, Z, did we adopt the median as the
expectation value since the grid is evenly spaced. Values in parentheses give
the uncertainty on the extracted parameter adopted from the standard deviation
of all models. We note that this model grid consists of about a factor of
three fewer model calculations than than the ten-dimensional model grid.