© T. Steindl et al. 2022

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_⊙ are expected to become unstable to stellar pulsations (e.g. Steindl et al. 2021b) as δ Scuti (e.g. Zwintz et al. 2014) or slowly pulsating B-type stars (e.g. Gruber et al. 2012). Pre-main-sequence δ Scuti stars are expected to have effective temperatures in the range from 6300 to 10 300 K (Steindl et al. 2021b) and pulsate in low-order pressure modes driven by the heat mechanism (Dupret et al. 2005; Steindl et al. 2021b). 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 δ Scuti stars. Such regularities are often shown in echelle diagrams and allow the asteroseismic modelling of δ 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 δ 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 δ 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 Sect. 2. In Sect. 3, we review, discuss, and introduce different approaches for the asteroseismic modelling of pre-main-sequence δ Scuti stars. These approaches were then extensively tested using our model grid, which we describe in Sect. 4. We introduced a simple algorithm based on the results of the first paper of this series (Steindl et al. 2021b) to automatically extract echelle diagrams, which we discuss in Sect. 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 Sect. 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 Sect. 7.

2. Computational set-up

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. 2021b). 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 α_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

$$\begin{array}{c}\hfill D_{\mathit{OV}}=D_{\mathrm{conv},0}exp\left(\frac{-2z}{f{H}_{\mathrm{p}}}\right),\end{array}$$(1)

where D_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₀ describes the position f₀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_above = F, $f_{0,\mathrm{above}}=f_{\mathrm{below}}=\frac{F}{2}$, and $f_{0,\mathrm{below}}=\frac{F}{4}$. The minimum mixing coefficient has values between D_mix = 0 and 100 cm s⁻¹.

The initial composition is set to X_¹H = 1 − X_²H − X_³He − X_⁴He − Z, where X_³He = 85 ppm, X_⁴He = Y + 2 * Z − X_³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_⊙ and 1.5 R_⊙ stellar seeds we used in Steindl et al. (2021b), 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 − τ) relation by Eddington (1926) successfully recreates the instability region for pre-main-sequence δ Scuti stars (Steindl et al. 2021b). 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 Sect. 5.1.3 of the first paper in this series (Steindl et al. 2021b). 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_add) is controlled by a factor β, such that

$$\begin{array}{c}\hfill L_{\mathrm{add}}=\beta ϵ\frac{\mathit{GMM}}{R},\end{array}$$(2)

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

For β = 0, no energy is deposited inside the stellar interior. This case is hence also referred to as cold accretion while 0 <  β ≤ 1 corresponds to hot accretion. With time-dependent accretion rates, β 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 β 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_outer with a linear increase as a function of the mass coordinate m_r (see also Kunitomo et al. 2017) according to

$$\begin{array}{c}\hfill l=\frac{L_{\mathrm{add}}}{M}\max \{0,\frac{2}{M_{\mathrm{outer}}^2}(\frac{m_r}{M}-1+M_{\mathrm{outer}})\}.\end{array}$$(3)

We calculate models with different values of M_outer, ranging between 0.01 and 0.2. Our models are calculated within the constant accretion scenario, meaning that the mass accretion rate ̇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

$$\begin{array}{c}\hfill M=\{\begin{array}{cc}{M}_0\hfill & \mathrm{if}t<t_{\exp },\hfill \\ M_{0}exp\left(\frac{-10(t-t_{\exp })}{\tau }\right)\hfill & \mathrm{else}\hfill \end{array}\end{array}$$(4)

where ̇M₀ is the initial constant accretion rate chosen between 10⁻⁷ and 10⁻⁵ M_⊙ yr⁻¹, t is the age of star since the beginning of the simulation, τ is the timescale for ten e-folding times in the mass-accretion drop chosen between 0.1 and 1 Myr, and t_exp is the time until the exponential drop-off calculated as

$$\begin{array}{c}\hfill t_{\exp }=(M_{\star }-0.01-M_0(1-exp(-10))\tau \times {10}^6)M_0^{-1}.\end{array}$$(5)

Here, M_⋆ 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, …, 8 and the frequency range from 5−100 d⁻¹. 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, δ 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 δ 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 Sect. 5.1.5 of Steindl et al. 2021b). As a consequence, 1309 calculations provide pre-main-sequence δ 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–10 300 K) for pre-main-sequence p-mode pulsations (Steindl et al. 2021b). This results in a total of 182 547 distinct δ 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.

Fig. 1. Parameters of 182 547 distinct pre-main-sequence δ Scuti models. Top left panel: 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.

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-χ²-square statistic:

$$\begin{array}{c}\hfill {\chi }^2=\frac{1}{n}{\displaystyle \sum _{i=1}^n}\frac{{(f_{i,\mathrm{obs}}-f_{i,\mathrm{mod}})}^2}{{\sigma }_f},\end{array}$$(6)

where f_i, obs are the observed frequencies, f_i, mod the model frequencies, and σ_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-χ²-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-χ²-square statistic was evaluated for all sets of model frequencies, we investigated the properties of all models fulfilling χ² ≤ 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σ or 3σ uncertainty box of T_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-χ²-square statistic, such that

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{spec}}^2& =\frac{1}{n+2}({\displaystyle \sum _{i=1}^n}\frac{{(f_{i,\mathrm{obs}}-f_{i,\mathrm{mod}})}^2}{{\sigma }_f}+\frac{{(T_{\mathrm{eff},\mathrm{obs}}-T_{\mathrm{eff},\mathrm{m}od})}^2}{{\sigma }_{T_{\mathrm{eff}}}}\hfill \\ \hfill & +\frac{{(log{g}_{\mathrm{obs}}-log{g}_{\mathrm{m}od})}^2}{{\sigma }_{T_{logg}}}).\hfill \end{array}$$(7)

However, the latter will reject any model that lies outside the spectroscopic error box if we later apply the condition ${\chi }_{\mathrm{spec}}^2\le 1$. A more statistically motivated approach is to reject all models with a pseudo-χ² statistic greater than the expected χ² 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 δ Scuti stars. Hence, the modelling of pre-main-sequence δ 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-χ²-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 Eq. (6), while the other two are based on Eq. (7). The approaches are as follows: (1) all models with χ² ≤ 1; (2) all models with χ² ≤ 1 and that lie within the 3σ spectroscopic error box; (3) all models with χ² ≤ 1 and that lie within the 1σ spectroscopic error box; (4) all models with ${\chi }_{\mathrm{spec}}^2\le 1$; (5) all models with ${\chi }_{\mathrm{spec}}^2<{\chi }_{\mathrm{cut}}^2$, where ${\chi }_{\mathrm{cut}}^2$ 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

We first applied our modelling approaches to one model (#10001). With an effective temperature of 7432 K, mass of 1.85 M_⊙, 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 δ Scuti star (Zwintz et al. 2014; Steindl et al. 2021a,b). The echelle spectrum of this model is shown in Fig. 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.

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

Stellar pulsation modes are represented by spherical harmonics, and higher degree modes suffer from cancellation effects, such that modes with l ≥ 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 Sect. 5). Hence, for this example we only used the radial and dipole modes of model #10001. As is evident from the upper panel of Fig. 3, a range of models with log g ∼ 4.2 has low values for the pseudo-χ² statistics according to Eq. (6). The five models with the lowest pseudo-χ² (≪1) also include, next to model #10001, model #52851 with T_eff = 8137 K, M_⋆ = 2.02 M_⊙, 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-χ² 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-χ² 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.

Fig. 3. Hexbin plot of the resulting pseudo-χ2 statistic according to Eq. (6) (upper panel) and Eq. (7) (lower panel) for the theoretical model #10001 with a Rayleigh limit of $\frac{1}{27}$ 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 Fig. 1. The correct position of model #10001 is marked with the green circle and the assumed spectroscopic error box is drawn in green.

Fig. 4. Pseudo-χ2 statistic distribution of effective temperature, star mass, surface gravity, and age for the frequencies of model #10001. Models that fall within the 1σ spectroscopic box are shown in green, those of the 3σ spectroscopic box are in blue, and all others are in black. The horizontal dashed grey line marks the position of pseudo-χ2 = 1, and the vertical red line shows the actual value for model #10001.

As a comparison, the lower panel of Fig. 3 shows the resulting pseudo-${\chi }_{\mathrm{spec}}^2$ according to Eq. (7). The minimum of the pseudo-${\chi }_{\mathrm{spec}}^2$ 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 }_{\mathrm{spec}}^2$ 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 δ 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 Sect. 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 uncertainty is 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 δ Scuti star). The results are shown in Table B.1 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.

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 δ 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σ of the stellar model’s value. Table B.2 shows the results for 3σ, and Table B.3 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σ 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σ 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 }_{\mathrm{spec}}^2\le 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 B.3) 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.

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.

5. Automatic echelle diagrams from observed pulsation frequencies

Amplitude spectra of δ 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 δ Scuti stars with echelle diagrams similar to the one shown in Fig. 2. From thousands of δ 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 δ 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 δ 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 δ Scuti stars (i.e. Netzel et al. 2022), and we propose a similar technique to automatically find echelle diagrams of pre-main-sequence δ 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 δ Scuti instability region from Steindl et al. (2021b). 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.

Fig. 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).

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-χ² value from Eq. (6) for every model in our model grid. If any model results in χ² ≤ 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 χ² ≤ 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 δ 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 Figs. 5 and C.1, especially). As a consequence, for frequencies with a period of about 0.06 d, almost all values between ∼0.35 and ∼0.8 are included within a boundary for l = 2 modes (see Fig. 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 χ² ≤ 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 d⁻¹. An example output is given in Fig. 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 Sect. 6, but the application to more stars will be the topic of future work.

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

6. Re-discussion of the pre-main-sequence δ 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 Sect. 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_⁴He = 0.29, our values lie in the X_⁴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_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_⊙, 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_⊙, effectively extending our model grid. Figure D.1 shows an updated version of Fig. 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 Myr¹ 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_⊙ 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 Fig. 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_⊙, which is already more than 5 σ from the value reported in Murphy et al. (2021).

Fig. 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).

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

Fig. 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 }_{\mathrm{spec}}^2<22.36$. The 1σ spectroscopic error box for HD 139614 is shown in magenta.

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

$$\begin{array}{cc}& \text{AIC}=2k-2ln{\mathcal{L}}_{\max ,}\hfill \\ \hfill & \text{BIC}=-2ln{\mathcal{L}}_{\max }+klnn,\hfill \end{array}$$

where ℒ_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

$$\begin{array}{cc}& \text{AIC}=nln(\text{RSS}/n)+2k+n,\hfill \\ \hfill & \text{BIC}=nln(\text{RSS}/n)+klnn,\hfill \end{array}$$

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

$$\begin{array}{cc}& {\text{AIC}}_{\mathrm{this}\mathrm{work}}=-106.9{\text{AIC}}_{\mathrm{Murphy}}=-100.5\mathrm{\Delta }\text{AIC}=-6.4,\hfill \\ \hfill & {\text{BIC}}_{\mathrm{this}\mathrm{work}}=-103.3{\text{BIC}}_{\mathrm{Murphy}}=-107.1\mathrm{\Delta }\text{BIC}=2.8.\hfill \end{array}$$

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 Ṁ₀ = 5 × 10⁻⁶,M_⊙ yr⁻¹, e-folding time τ = 0.5 Myr, fractional energy β = 0.1, and fractional mass M_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 D.1) to the ten-dimensional model grid (see Table 4), with only small changes in the expectation values and the uncertainties. The best χ² 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:

$$\begin{array}{cc}& {\text{AIC}}_{6\mathrm{d}}=-114.9{\text{AIC}}_{\mathrm{Murphy}}=-100.5\mathrm{\Delta }\text{AIC}=-14.4,\hfill \\ \hfill & {\text{BIC}}_{6\mathrm{d}}=-117.1{\text{BIC}}_{\mathrm{Murphy}}=-107.1\mathrm{\Delta }\text{BIC}=-10.0.\hfill \end{array}$$

Hence, the AIC provides very strong evidence and the BIC strong evidence (with Δ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 δ 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.

7. Conclusion

In this work, we investigated different approaches in the pulsational modelling of pre-main-sequence δ 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-χ² statistic as a goodness of fit indicator for the pulsation frequencies. While the first approach uses a simple pseudo-χ² 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-χ² 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 δ 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-χ² statistic and sets the pseudo-χ² cut-off depending on the amount of free parameters, is the most reliable. We note that approach 2, which uses the simpler pseudo-χ² statistics and only considers models within 3σ 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_⊙. 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 δ Scuti stars. This algorithm is based on theoretical Petersen diagrams taken from the non adiabatic pre-main-sequence models of Steindl et al. (2021b). 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 δ 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 δ Scuti stars is unavoidable. If, in contrast, simpler model grids are utilised, resulting parameters are often extracted with incorrect expectation values and underestimated uncertainties.

Acknowledgments

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

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

Appendix C: Petersen diagrams

Fig. C.1. Peterson diagram for the first overtone.

Fig. C.2. Peterson diagram for the second overtone.

Fig. C.3. Peterson diagram for the third overtone.

Fig. C.4. Peterson diagram for the fourth overtone.

Appendix D: Re-discussion of HD 139614

Fig. D.1. 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.

All Tables

All Figures

	Fig. 1. Parameters of 182 547 distinct pre-main-sequence δ Scuti models. Top left panel: 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.
In the text

	Fig. 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.
In the text

Fig. 3. Hexbin plot of the resulting pseudo-χ2 statistic according to Eq. (6) (upper panel) and Eq. (7) (lower panel) for the theoretical model #10001 with a Rayleigh limit of $\frac{1}{27}$ 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 Fig. 1. The correct position of model #10001 is marked with the green circle and the assumed spectroscopic error box is drawn in green.

In the text

Fig. 4. Pseudo-χ2 statistic distribution of effective temperature, star mass, surface gravity, and age for the frequencies of model #10001. Models that fall within the 1σ spectroscopic box are shown in green, those of the 3σ spectroscopic box are in blue, and all others are in black. The horizontal dashed grey line marks the position of pseudo-χ2 = 1, and the vertical red line shows the actual value for model #10001.

In the text

	Fig. 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).
In the text

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

In the text

Fig. 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).

In the text

	Fig. 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 }_{\mathrm{spec}}^2<22.36$. The 1σ spectroscopic error box for HD 139614 is shown in magenta.
In the text

	Fig. C.1. Peterson diagram for the first overtone.
In the text

	Fig. C.2. Peterson diagram for the second overtone.
In the text

	Fig. C.3. Peterson diagram for the third overtone.
In the text

	Fig. C.4. Peterson diagram for the fourth overtone.
In the text

	Fig. D.1. 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.
In the text