Asteroseismology of $δ$ Scuti stars: emulating model grids using a
  neural network
Young $\delta$ Scuti stars have proven to be valuable asteroseismic targets
but obtaining robust uncertainties on their inferred properties is challenging.
We aim to quantify the random uncertainties in grid-based modelling of $\delta$
Sct stars. We apply Bayesian inference using nested sampling and a neural
network emulator of stellar models, testing our method on both simulated and
real stars. Based on results from simulated stars we demonstrate that our
method can recover plausible posterior probability density estimates while
accounting for both the random uncertainty from the observations and neural
network emulation. We find that the posterior distributions of the fundamental
parameters can be significantly non-Gaussian, multi-modal, and have strong
covariance. We conclude that our method reliably estimates the random
uncertainty in the modelling of $\delta$ Sct stars and paves the way for the
investigation and quantification of the systematic uncertainty.
Owen J. Scutt, Simon J. Murphy, Martin B. Nielsen, Guy R. Davies, Timothy R. Bedding, Alexander J. Lyttle
# Asteroseismology of \(\delta\) Scuti stars: emulating model grids using a neural network

###### Abstract

Young \(\delta\) Scuti stars have proven to be valuable asteroseismic targets but obtaining robust uncertainties on their inferred properties is challenging. We aim to quantify the random uncertainties in grid-based modelling of \(\delta\) Sct stars. We apply Bayesian inference using nested sampling and a neural network emulator of stellar models, testing our method on both simulated and real stars. Based on results from simulated stars we demonstrate that our method can recover plausible posterior probability density estimates while accounting for both the random uncertainty from the observations and neural network emulation. We find that the posterior distributions of the fundamental parameters can be significantly non-Gaussian, multi-modal, and have strong covariance. We conclude that our method reliably estimates the random uncertainty in the modelling of \(\delta\) Sct stars and paves the way for the investigation and quantification of the systematic uncertainty.

keywords: asteroseismology - stars: variables: Scuti - stars: fundamental parameters - methods: data analysis - methods: statistical

## 1 Introduction

Stellar ages for individual stars are notoriously difficult to measure (Soderblom, 2010). One method is to model a cluster with isochrones, which is particularly sensitive to high-mass stars at the main-sequence (MS) turn-off (Lipatov et al., 2022). Other techniques, such as the lithium depletion boundary (e.g. Galindo-Guil et al., 2022) or kinematics (Miret-Roig et al., 2022; Zerjal et al., 2023), are able to use low-mass stars, which are much more abundant. However, methods that utilize intermediate-mass stars for measuring stellar ages have been lacking.

Asteroseismology - the study of stellar oscillations - is highly sensitive to age and has long held promise as an independent method for age determination (e.g., Aerts, 2015). Like other techniques, asteroseismology is model-dependent, but the physics of those models is generally different from the high- and low-mass stars (Soderblom, 2010), hence the techniques are highly complementary (Kerr et al., 2022, 2022). Until recently, however, asteroseismology of intermediate-mass stars (the so-called \(\delta\) Scuti variables) has been hampered by the difficulties in identifying which modes are excited. The discovery of regular patterns in the pulsation mode frequencies of some \(\delta\) Sct stars (Bedding et al., 2020) has opened up a pathway to determine their masses, ages, and metallicities, without the requirement that the star resides in a cluster or association.

In recent years, oscillations in large numbers of \(\delta\) Sct stars have been measured using white-light photometry from space telescopes such as _CoRoT_(e.g., Paparo et al., 2013; Michel et al., 2017; Barcelo Forteza et al., 2018), _Kepler_(e.g., Uytterhoeven et al., 2011; Balona et al., 2015; Garcia Hernandez et al., 2017; Bowman and Kurtz, 2018; Guzik, 2021) and _TESS_(e.g., Antoci et al., 2019; Hasanzadeh et al., 2021; Barac et al., 2022; Chen et al., 2022). Observed oscillation frequencies can be compared against grids of model frequencies to find a best-fitting set of parameters (Murphy et al., 2021, 2022). It is somewhat more challenging to understand the resulting uncertainties, which are not uniquely determined by the spacing of the model grid (Pedersen, 2020), and instead depend more strongly on the underlying physics (Steindl et al., 2021). Part of the challenge is that models can be computationally expensive and calculating new evolutionary tracks on-the-fly for Monte Carlo sampling is prohibitive.

In order to treat the uncertainties more robustly, we aim to convert a discrete grid of stellar models into a continuous function. We use a neural network to emulate a grid of stellar models that has been pre-computed over the range of expected stellar parameters. We combine the trained neural network with a Bayesian sampler to formally treat random uncertainties in the observables. This yields estimates for the posterior probability density of the fundamental properties which quantifies their uncertainties. It also allows us to infer viable frequencies for modes that were not detected, but which might exist in the data at low signal-to-noise.

In the following section we describe the grid of stellar models on which the neural network is trained, and in Sec. 3 we discuss the details of the network architecture and training method. In Sec. 4 we present the method used to perform the Bayesian inference, and show results for a selection of simulated and real sets of observations (Sec. 5).
## 2 The Stellar Model Grid

We used the model grid of Murphy et al. (in prep), consisting of evolutionary tracks computed with MESA (r15140; Paxton et al., 2011, 2013, 2015, 2018, 2019) and pulsation models calculated with GYRE (v6.0.1; Townsend & Teitler, 2013). Provisional versions of this grid have already been used to model the pulsations of \(\delta\) Sct stars (Murphy et al., 2022; Kerr et al., 2022a,b; Currie et al., 2022), and the physics of the models are described in Murphy et al. (2022).

A well-sampled grid was needed to train the neural network emulator. Here, evolutionary tracks were spaced by \(0.02\,\mathrm{M}_{\odot}\) in mass \(M\) and \(0.001\) in initial metallicity \(Z_{\mathrm{in}}\). For \(Z_{\mathrm{in}}>0.010\), the spacing was increased to \(0.002\). The grid is shown in mass-metallicity space in Fig. 1. A common problem in MESA is that pre-MS models sometimes fail to converge and the evolution is terminated (see, e.g., Steindl et al., 2021). In such cases, we attempted to re-calculate the track with a slightly increased mass (\(M\)+= 0.001) up to five times before abandoning that track. Abandoned tracks appear as gaps in the grid in Fig. 1.

It is also important to ensure the tracks are sampled well in age. Computational errors are minimised by keeping the time interval small throughout the evolution, even if not all time steps are saved as outputs. The internal sampling is described in Murphy et al. (in prep). For outputs, we saved evolutionary and pulsation models every \(0.05\,\mathrm{Myr}\) from \(2\,\mathrm{Myr}\) until \(10.5\,\mathrm{Myr}\), in order to adequately sample the rapid evolutionary changes that occur on the pre-MS. After this the evolution is somewhat slower, and sampling of \(3\,\mathrm{Myr}\) was deemed adequate up to ages of \(40\,\mathrm{Myr}\). Beyond this, the tracks were instead sampled according to changes in position on the HR diagram (limits of \(\Delta\log T_{\mathrm{eff}}=0.0006\) and \(\Delta\log L=0.002\)), with an upper limit of \(100\,\mathrm{Myr}\) between samples. Where large gaps occurred in the grid, or when the specific \(M\)-\(Z_{\mathrm{in}}\) combination demanded it, we manually recalculated tracks with finer sampling. This explains the variations in the number of samples per track in Fig. 1.

For each pulsation model, we computed the frequencies of radial modes (spherical degree \(\ell=0\)) having radial orders \(n\) from \(1\) to \(11\), and dipole (\(\ell=1\)) modes having \(n\sim 1\)-\(10\). This encompasses the range of radial orders observed for real stars (e.g. Bedding et al., 2022). We calculated the mean frequency separation between radial orders (\(\Delta\nu\)) using the radial modes having \(n=5\)-\(9\), by fitting a straight line to the mode frequencies as a function of \(n\)(see White et al., 2011), using

\[\nu=\Delta\nu(n+\ell/2+\epsilon). \tag{1}\]

The variable \(\epsilon\) is the intercept of that line with the y-axis, and describes the distance of the radial mode ridge from the y-axis in an echelle diagram. In addition to the individual mode frequencies, we stored the values of \(\Delta\nu\) and \(\epsilon\) for each model in the grid, since these asteroseismic quantities relate to astrophysical quantities (Murphy et al. in prep.).

To reduce the effect of the strong covariance between stellar age \(\tau\) and mass \(M\), and ease the training of the neural network, we used the assumption that the MS lifetime is approximately proportional to \(M^{-3.2}\) and defined the scaled age (e.g., Davies & Miglio, 2016)

\[\mathcal{K}=10^{-4}\,\tau\,(M/\mathrm{M}_{\odot})^{3.2}. \tag{2}\]

This scaled age serves as an estimate of the fractional MS age of our models.

## 3 Constructing the Neural Network

To overcome the discretely sampled nature of the model grid, we used a neural network consisting of a series of fully connected dense layers in place of standard interpolation for continuous stellar model emulation. The network was trained on the model grid, learning to predict observable parameters given stellar model input parameters. This way, the network learned the map from observables to model parameters and could be used for likelihood estimation during inference. To this end, we used the fundamental parameters \(M\), \(Z_{\mathrm{in}}\) and scaled age (\(\mathcal{K}\)) as inputs for parameter augmentation and network training. Outputs consist of the classical observables (\(L\) and \(T_{\mathrm{eff}}\)); asteroseismic quantities (\(\Delta\nu\) and \(\epsilon\)); and \(11\) radial and \(10\) dipole mode frequencies. We refer to these \(25\) outputs collectively as the 'observable parameters'.

Once the input and output parameters were defined, we carried out dataset-wide parameter augmentation to improve the training of the network. We converted all parameters (excluding \(\epsilon\)) to the decimal logarithm and applied a Z-score standardisation to all parameters (including \(\epsilon\)). Both of these operations restricted all parameters to similar ranges, to avoid the neural network assigning erroneously high importance to parameters spanning several orders of magnitude during training. We found the combination of the two operations to be optimal for this investigation.

To further simplify the training process, we performed principal component analysis on the observable parameters in the model grid, as follows. For all models, we calculated the covariance matrix of all \(25\) observable parameters. The eigenvectors of the resulting covariance matrix, or 'principal components', were ranked in order of descending eigenvalue, returning a list of principal components explaining the most to the least variance in the observable parameters. We determined how many principal components to include using the explained variance ratio, which describes the percentage of the variance of the observable space present in just the chosen principal components. We found that \(9\) principal components explained all but \(10^{-4}\) per cent of the total variance. This sufficiently explained the covariance of the \(25\) parameters in the full observable space.

The use of principal components presents the neural network with a simpler map to learn--replacing the fundamental parameters by the reduced dimensions of the 'latent parameters'--and also removes covariance information from the observables, which is redundant for the neural network. We then added a custom non-trainable layer to the

Figure 1: The model grid used in this work, where each symbol represents an evolutionary track with a particular metallicity (\(Z_{\mathrm{in}}\)) and mass. Colour-coding indicates the number of models along each track for which pulsation frequencies were calculated (see Sec. 2).
neural network, which projects the latent parameters back into the full observable parameter space before the network outputs predictions.

Finally, we split the model grid into a training and a testing set for the neural network. The training set was randomly selected to comprise 80 per cent of the model grid, to be seen by the network during training. The test set, composed of the remaining 20 per cent of the grid, was unseen by the network during the training process and was used solely for evaluation of network prediction performance. This served as a check that the network is capable of model grid interpolation -- the training set became a sparser representation of the original model grid, with the test set providing models guaranteed to hold combinations of parameters previously unknown to the network.

In addition to the data augmentation prior to training, the hyperparameters of a neural network can be tuned to promote faster and more stable learning. To quantify network performance for comparison between different hyperparameter permutations, we compared their validation loss profiles over multiple network training sessions. We adopted a 'grid search' method for testing potential combinations of network hyperparameters. This involved creating a grid of potential values for the number of fully connected dense layers (ranging from 3 to 8 in steps of 1), activation functions, optimizers, learning rates, loss functions, and batch sizes. We populated a grid with these hyperparameters, and then tested the resulting network at each position in the hyperparameter grid for successful validation loss minimisation.

We found the optimal network consisted of 6 fully-connected dense layers of 64 neurons, each using an exponential linear unit activation function (Clevert et al., 2015), followed by the custom layer for projection from latent to observable parameters, and a final dense output layer with linear activation function. We used the Adam optimizer (Kingma and Ba, 2014) with a learning rate of \(10^{-4}\), and the mean-squared-error loss function. A batch size of \(6\times 10^{4}\) models provided a good compromise between speed and training stability. We used a validation split of 25 per cent of the training set -- where the test set is used to evaluate neural network success after training, a 'validation set' is used for evaluation of neural network success during training. Once primary training was complete with the learning rate above, we saved and recompiled the network with a slower but less volatile learning rate of \(7\times 10^{-5}\), and restarted training until no validation loss reduction was observed for \(10^{4}\) training epochs. The network and custom latent-to-observable projection layer were constructed using the TensorFlow sequential API (Abadi et al., 2015).

Once the optimal network from the grid search was trained, we evaluated the network performance across the full observable parameter space. Using the test set previously removed from the model grid, we plotted distributions of the decimal logarithm prediction residuals for each parameter. This allowed us to visualize any bias and uncertainty inherent in the neural network predictions. We used the median absolute deviation of these prediction residual distributions, shown in Fig 2, to quantify network prediction uncertainty for observable parameters. We found an uncertainty in network predictions of \(8\times 10^{-4}\) dex and \(2\times 10^{-4}\) dex for log \(L\) and \(\log T_{\text{eff}}\), respectively, and a mean uncertainty of \(\sim 3\times 10^{-4}\) dex for the mode frequencies.

## 4 Inferring the fundamental stellar parameters

To perform the Bayesian inference on the input model parameters, \(\theta=(M,\log Z_{\text{in}},\log\mathcal{K})\), for a given observed set of mode frequencies, we sampled the posterior distribution

\[P(\theta|D)=\frac{\mathcal{P}(\theta)\mathcal{L}(D|\theta)}{\mathcal{E}(D)}. \tag{3}\]

Here, \(\mathcal{P}(\theta)\) is the prior on the input model parameters, and \(\mathcal{L}(D|\theta)\) is the likelihood of observing a set of parameters (\(D\)) given the model parameters. The evidence, \(\mathcal{E}(D)\), is calculated at each step during the sampling.

In addition to the input model parameters, we also included a variable offset term, \(\Delta n\), as input to account for the possible ambiguity in assigning the radial orders of the observed modes. This ambiguity arises because the radial orders of a set of modes cannot always be determined from the observed mode frequencies alone, and are typically decided by comparison to stellar models. Including \(\Delta n\) in the sampling allows us to marginalize over this uncertainty when estimating the posterior distribution of the input model parameters. We expect this uncertainty to only lead to an error of \(\pm 1\) radial order, and so we chose the prior on \(\Delta n\) to be a set of \(\delta\)-functions at \(\Delta n=-1,0\) and \(1\).

Table 1 lists a summary of the prior functions. For the priors on \(M\), log \(Z_{\text{in}}\) and log \(\mathcal{K}\), we chose to use \(\beta\) distributions, since they can be bounded to match the limits of the stellar model grid. In addition, the shape parameters of the \(\beta\) distributions may be chosen such that the priors reflect our expectation of the distribution of real observations of \(M\), \(Z_{\text{in}}\) and \(\mathcal{K}\).

Our choice of range and shape of the prior on \(\mathcal{K}\) is motivated by the age range and distribution we expect to target, and also be able to observe. Mode identification for \(\delta\) Sct stars is currently possible up to approximately one third of the MS age, after which the coupling between the buoyancy dominated and acoustic modes spoils the regular mode frequency patterns. Furthermore, at older ages the physics of mixing and overshooting become more important, and those were not treated as variables in the model grid. Hence, the prior on age extends to approximately one third of the expected MS lifetime. The lower limit on the prior on the scaled age was chosen because stars in our mass range of interest (see below) do not evolve to cross the \(\delta\) Sct instability strip until ages \(\geq 2\) Myr. Due to the motion of stars through the \(\delta\) Scitnstability strip, the age prior is biased toward lower ages, with a fall-off in the age prior distribution toward older stars.

The prior on \(Z_{\text{in}}\) ranges from approximately 0.07 to 1.5 times the solar metal mass fraction of 1.42 per cent used in the models (Asplund et al., 2009), which covers the metallicity distribution of stars forming within approximately 1 kpc of the Sun at the current age of the Galaxy (Hayden et al., 2020). The existence of metal-poor \(\delta\) Sct stars in modern star-forming regions (e.g. HD 139614 in Upper Centaurus Lupus; Murphy et al., 2021) suggests that slightly sub-solar metallicities are more common than slightly super-solar metallicities in young \(\delta\) Sct stars. We therefore skewed the prior probability density toward sub-solar values.

Finally, the mass range was chosen to ensure that models exist both within and on either side of the instability strip (Dupret et al., 2004; Murphy et al., 2019). Our slight skew towards lower masses accounts for the similar skew present in the stellar initial mass function (Krumholz, 2014).

Figure 3 shows samples drawn from these prior density distributions, both in terms of the sampled variables and those transformed to \(M\), \(Z_{\text{in}}\) and age. These priors are applied to the inference performed for all targets (see below). Additional priors may be applied on a target-by-target basis if, for example, the mass can be constrained by other sources such as orbiting companions, or limits can be placed on the metallicity by spectroscopy.

For each of the samples drawn from the prior distributions, the
neural network produces the following outputs: a set of mode frequencies, the effective temperature, and the luminosity. Given a set of outputs we then evaluated the likelihood of the observations by

\[\log\mathcal{L}\left(D|\theta\right)=\log\mathcal{L}(D_{\mathrm{S}}|\theta)+\log \mathcal{L}(D_{\mathrm{C}}|\theta). \tag{4}\]

We separate the log-likelihood into the seismic variables, \(D_{\mathrm{S}}\), and the classical (non-seismic) variables, \(D_{\mathrm{C}}\). The contribution to the likelihood of the mode frequencies is given by

\[\log\mathcal{L}(D_{\mathrm{S}}|\theta)=\sum_{i}\log\mathcal{N}\left\langle v_{ i}^{\mathrm{obs}},\sqrt{\sigma_{v_{i}^{\mathrm{obs}}}^{2}+\sigma_{v_{i}^{ \mathrm{SN}}}^{2}}\right\rangle, \tag{5}\]

and that of the classical observables is given by

\[\log\mathcal{L}(D_{\mathrm{C}}|\theta)= \log\mathcal{N}\left(\log L^{\mathrm{obs}},\sqrt{\sigma_{L^{ \mathrm{obs}}}^{2}+\sigma_{L^{\mathrm{SN}}}^{2}}\right)+\] \[\log\mathcal{N}\left(T_{\mathrm{eff}}^{\mathrm{obs}},\sqrt{ \sigma_{T^{\mathrm{obs}}}^{2}+\sigma_{T^{\mathrm{SN}}}^{2}}\right). \tag{6}\]

The width of the probability densities used in the inference is given by two terms that specify the observational uncertainty (superscript 'obs'), and the noise due to the precision of the neural network's ability to emulate the model grid (superscript 'NN'). Based on the spread of the residuals presented in Fig. 2, this emulation uncertainty is approximately \(4\times 10^{-4}\) dex, which equates to a relative uncertainty of \(\approx 0.1\) per cent on the output parameters. This additional uncertainty was added in quadrature to the uncertainty of the observed mode frequencies.

In the following we will use simulated frequencies corresponding to those obtained from 2 sectors of data from the _TESS_ mission (Ricker et al., 2015). We therefore adopted an uncertainty on the mode frequencies of \(0.02\mathrm{d}^{-1}\), which is the frequency resolution of the resulting power spectra. The uncertainties on \(\log L\) and \(\log T_{\mathrm{eff}}\) depend on the target in question, but for the simulations shown below these were fixed to \(\sigma_{L}=0.05\) dex and \(\sigma_{T}=200\)K.

The neural network residuals also showed a bias of \(\sim 10^{-5}\) dex, which equates to an offset of \(0.01\) per cent on each of the output parameters. This offset is small compared to the combination of the observed and neural network uncertainties, and so we did not consider it in the analysis. However, if either of these sources of uncertainty were decreased by, for example, improving the estimates of the observed mode frequencies, the importance of this bias would need to be re-evaluated.

We performed the sampling using the nested sampling method from the Dynesty Python package (Skilling, 2004; Speagle, 2020). Nested sampling determines iso-likelihood contours in the input parameter space, which were iteratively redefined until samples were consistently drawn around the global likelihood maximum. In the Dynesty package, this process is terminated when the change in model log-evidence \(\Delta\log\mathcal{E}\) is less than a predefined value chosen according to the Dynesty documentation1. The method presented above is not restricted to using Dynesty, and other sampling methods

\begin{table}
\begin{tabular}{c c} \hline Parameter & Prior function \\ \hline \(M\left[M_{\odot}\right]\) & \(\beta_{\mathrm{S}}^{2}(1.3,2.3)\) \\ \(\log Z_{\mathrm{in}}\) & \(\beta_{\mathrm{S}}^{2}(-3.1,-1.6)\) \\ \(\log\mathcal{K}\) & \(\beta_{\mathrm{L,2}}^{2}(-3,-0.3)\) \\ \(\Delta n\) & \(\epsilon_{R}^{2}\left\{-1,0,1\right\}\) \\ \hline \end{tabular}
\end{table}
Table 1: Prior density functions used in Eq. 3. The priors on \(M\), \(\log Z_{\mathrm{in}}\) and \(\mathcal{K}\) are given by \(\beta_{\mathrm{b}}^{\alpha}\), where \(a\) and \(b\) are the shape parameters of the \(\beta\)-distributions, and the prior on \(\Delta n\) is a series of \(\delta\)-functions at integer values. In all cases the arguments to the distribution functions denote lower and upper limits.

Figure 2: Residual distributions of predictions from the neural network. The residuals are that of the decimal logarithms of \(T_{\mathrm{eff}}\), \(L\) (blue), and mode frequencies (\(n\), \(\ell\)) (orange). The boxes show a central line at the median value of the distribution, with edges at the lower and upper quartiles. Whiskers extend to the 5th and 95th percentile range. The dashed line indicates complete agreement between the network predictions and model grid values.
may be used, such as MultiNest(Buchner et al., 2014) or EMCEE (Foreman-Mackey et al., 2013).

## 5 Results

### Simulated stars

In order to test our methodology and validate the accuracy of the neural network emulator, we have performed tests based on 25 simulated stars in a 'hare-and-hounds' exercise. To produce these simulated stars, we proceeded as follows. Values of stellar mass and initial metallicity were selected to lie in between values in the grid, but still within the defined parameter range of the grid. We calculated stellar models and pulsation frequencies using MESA and GYRE, using the same settings as for the grid. We selected ages from the newly calculated tracks, which then defined the truth values for mass, metallicity and age for our simulated stars and their associated 'true' observables. We only selected a subset of the calculated modes, to better reflect typical observations of \(\delta\) Sct stars. We selected modes at four consecutive radial orders for each degree, within the bounds of \(n=1\)-\(8\).

To simulate noisy observations, we added noise to the observable parameters of the simulated stars. These random offsets were drawn from a normal distribution, with mean of zero and a standard deviation of: \(0.02\,\mathrm{d}^{-1}\) for the mode frequencies, \(200\,\mathrm{K}\) for the effective temperatures, and \(0.05\,\mathrm{dex}\) for the log-luminosity.

#### 5.1.1 Exemplar simulated star

Figure 4 shows the posterior probability estimates for a simulated star with one of the best results. The figure demonstrates our ability to quantify random uncertainty on our inferred properties from the posterior and shows that the true properties of this simulated star lie comfortably with the posterior distribution. In this case our method is performing as required but it is worth noting that, even for this exemplar, the posterior still contains significant correlation between the parameters.

We see that the posterior distribution is not well described by a series of separable 1-D normal distributions. Instead, there are strong covariances between inferred parameters, which are to be expected from stellar evolution theory. However, the 1-D marginal posterior distributions show evidence of not being normally distributed and, in the case of the stellar age, even somewhat multi-modal. To examine the degree of accuracy and precision of the results we will study the summary statistics of the posterior distribution. This does not capture all the detail that is of value, but is nonetheless useful as a test of our methods.

#### 5.1.2 Results for 25 simulated stars

While examining a single simulated star is useful, it is hard to draw conclusions on the validity of our approach because we are looking at a single realisation of noise on the observables. We now consider all 25 simulated stars, including the exemplar above (simulated star index 5), to look at the statistics of our posterior probability distributions when compared to the truth values of the input properties. As part of our method, we fitted a parameter to account for our uncertainty in our assumption of the radial order label \(n\). In our tests on simulated stars, we recovered the correct radial order in all cases with no meaningful uncertainty on the posterior of the radial order labels.

For each parameter of each star, we computed the difference between the truth value and the inferred value (the mean of the posterior samples for that parameter) divided by the uncertainty (the standard deviation of the posterior samples for that parameter). If our inference is perfect, and our posterior distributions are well behaved, then this metric should be drawn from a normal distribution with zero mean and unit variance. However, as observed above in our exemplar, multi modality, non-Gaussianity, and other pathological behaviour in the posterior can bias our metrics away from our assumed normal distribution.

Figure 5 shows the metric for each simulated star and each input property of the star. It is clear that the majority of our simulated star results are consistent with the truth value given the uncertainty. And broadly, the numbers of metrics at the 1 and 2 sigma levels are consistent with expectations. There are however some outliers or results which we will discuss. The index for the most significant outlier in terms of metallicity is 14 and an interesting behaviour in the age distribution is observed in simulated star index 21.

#### 5.1.3 Further tests of simulated star 14

Simulated star 14 appears as an outlier in metallicity by \(\sim 2.5\sigma\). To examine this behaviour we have produced 10 more realisations of this simulated star. That is, we have taken the same truth values as inputs, but redrawn the simulated noise on each observable parameter using the same noise distributions.

Figure 6 shows the posterior samples for the original and subsequent runs for simulated star 14. Firstly, it is clear that the posterior for \(Z_{\mathrm{in}}\)is multi-modal and contains significant covariance. Secondly, there is a bias of the posterior distributions away from the truth value in both metallicity and age that cannot be explained simply as a result of the realisation noise on the observables.

We have checked for the possible origins of this bias. We examined the prior probability distribution, but found it to be smooth and nearly flat over the region of the posterior. We have performed multiple realisations of the noise and still observed this bias and therefore also exclude the noise or likelihood as the source of the bias.

A possible source of error is in the neural network emulation producing differences in the predicted mode frequencies. While these errors are typically small, of order \(3\,\times\,10^{-4}\) dex, the noise from the neural network is not random noise that would be expected to reduce with more realisations of the observables. Instead, the error is systematic and will produce a bias. The systematic error will always be present in emulation and this will lead to a bias, but it is the magnitude of this bias that is interesting. For this simulated star, the bias is similar to the reported uncertainty, which is about \(1.5\,\mathrm{Myr}\). However, this error can be reduced by extending the training time of the neural network, or increasing the grid search density around the optimal neural network architecture.

#### 5.1.4 Further tests of simulated star 21

Simulated star 21 shows an interesting behaviour in the age posterior. Figure 7 shows the posteriors for the original simulated star 21 and for 10 more realisations, as we did for simulated star 14 above. No meaningful bias is observed in the posteriors for mass or metallicity, given the priors we apply.

The true age of the simulated star was \(12.96\,\mathrm{Myr}\), which corresponds to the pre-MS evolution stage. The age posterior is clearly bimodal, with solutions around \(10\,\mathrm{Myr}\) and \(160\,\mathrm{Myr}\). This behaviour is observed in all the realisations, lending confidence that this is not a
result of the noise being added. In fact, this bi-modality is consistent with our understanding of the evolution of these stars and illustrates the difficulty of distinguishing the phase of the MS evolution where the track crosses its pre-MS evolution in the HR diagram.

### Application to HD 99506

We applied our methods to HD 99506, which is one of the high-frequency \(\delta\) Sct stars discussed by Bedding et al. (2020). We used the following inputs: \(T_{\rm eff}=7970\pm 250\) K and \(L\)/L\({}_{\odot}=7.58\pm 0.37\) (taken from Table 1 of Bedding et al., 2020), and the mode frequencies that we have measured and listed in Table 2. We chose only the mode frequencies that were obvious, leaving out tentatively identified modes such as the \(n=4\) and \(n=9\) radial modes, and the \(n=1\) and \(n=8\) dipole modes. The identified modes span two radial orders more than any of the simulation and so, despite the gaps at some orders, they provide tighter constraints. The resulting posteriors on \(M\), \(Z_{\rm in}\) and age are unimodal and indicate a percent-level random uncertainty (Fig. 8). The inferred age (\(9.71\pm 0.31\) Myr) corresponds to the pre-MS phase, before the onset of pp-chain H-burning but after the temporary pre-MS CNO burning phase.

The calculation of well-sampled posteriors for uncertainty estimates is a marked improvement on what is possible using discrete grid points and \(\chi^{2}\) minimisation (e.g. Kerr et al., 2022), where an arbitrary threshold in \(\chi^{2}\) needs to be adopted. It is especially useful that the posteriors are marginalized, given the aforementioned correlation in astrophysical parameters demonstrated with the simulated stars.

The neural network is also able to generate posterior predictions for

Figure 4: Samples drawn from the posterior distribution of simulated star 5. For clarity we transform the initial metallicity to % and \(\mathcal{K}\) to stellar age \(\tau\). The model input values used to generate the simulated star 5 data are shown in blue.

Figure 3: Left: Samples drawn from the one-dimensional prior distributions of the sampled parameters. The diagonal frames show the marginalized distributions (black) and the functions used to draw the samples (orange). The prior used for \(\Delta n\) are \(\delta\)-functions at \(\Delta n=-1,0\), and \(1\). The off-diagonal frames show the two-dimensional distributions of the input parameters. Right: The samples from the left frames (black) transformed to the same units as the output from the stellar model grid (blue).
each mode frequency, using the posterior samples as inputs. This can be useful for estimating the validity of uncertain mode identifications. In Fig. 9, we see that the leftmost (lower frequency) of two close peaks at the \(n=4\) radial mode is a good match and could perhaps be identified. The weak peak at \(n=9\) would also have fitted well. Inclusion of these would have resulted in tighter posteriors. On the other hand, none of the missing dipole mode frequencies, nor the \(n=1\) or 2 radial modes, would have been good additions. If we had supplied those modes as input, the posteriors would have broadened markedly.

## 6 Conclusions

We have presented a method for performing Bayesian inference on fundamental stellar properties of \(\delta\) Sct stars using a neural network. This method emulates the stellar model and oscillation codes, MESA and GYRE, by learning from a grid of models that encompasses the physical properties of stars in or near the \(\delta\) Sct instability strip. We used a nested sampling method to estimate the posterior distribution of the fundamental stellar properties, given a set of mode frequencies and classical observables. The resulting posterior distribution reflects the random observational uncertainty as well as the uncertainty of the neural network. By providing samples from the posterior probability density, which might be multimodal, non-Gaussian, and strongly covariate, we formally quantified the statistical uncertainty in the fundamental stellar properties. This improves our ability to investigate the systematic uncertainty in the stellar models.

We used a test set that was initially unseen by the training algorithm to evaluate the performance of the trained neural network. We found that the neural network is capable of reaching an average frequency precision of \(\approx 3\times 10^{-4}\) dex, with an offset of \(\approx 5\times 10^{-5}\) dex. These performance metrics may improve if the network were retrained with, for example, additional grid points or with the aim to reach a lower target loss. However, the flexibility of neural networks

Figure 5: Difference between the inferred and model values of \(\mathbf{M}\), \(Z_{\rm in}\) and age relative to the uncertainty of the inference, for a set of 25 simulated stars. The uncertainty is taken as the standard deviation of the marginalized posterior distributions of each of the parameters (see Sec. 5.1.2 for details).

Figure 6: Samples drawn from the posterior distribution of simulated star 14 are shown in red. Subsequent runs using the same truth values (M=1.696 M\({}_{\odot}\), \(Z_{\rm in}\)=0.0136, \(r\)=8.76 Myr) but different noise realisations are shown in black. The truth values are shown in blue (see Sec. 5.1.3 for details).

Figure 7: Samples drawn from the posterior distribution of simulated star 21 are shown in red. Subsequent runs using the same truth values (M=1.777 M\({}_{\odot}\), \(Z_{\rm in}\)=0.0182, \(r\)=12.97 Myr) but different noise realisations are shown in black. The truth values are shown in blue (see Sec. 5.1.4 for details).
allows for the extension of the model grid to include additional variables, such as initial helium abundance or convective overshoot, by increasing the number of neurons in the network architecture (see, e.g., Hendriks & Aerts, 2019; Lyttle et al., 2021).

We applied our method to 25 simulated stars to quantify our accuracy and precision in the recovery of our input stellar properties. We showed the method to be capable of faithfully reproducing the true input parameters with the exception of simulated star 14. On investigation, we found a bias in the reported metallicity and age of this simulated star to be a result of the error in prediction by the neural network. The bias is of order 1.5 Myr in age for this star, but was confirmed to be systematic in nature. Further improvements in the accuracy of the neural network emulation would reduce the size of this effect.

Finally, we have applied our method to observations of a real \(\delta\) Sct star, HD 99506. We found this star to be in the pre-MS stage of evolution and report a random uncertainty of only 3 per cent in age. Systematic uncertainty, such as that arising from missing or imperfect physics, has not been accounted for in this number but our methods pave the way for the quantification of the systematic uncertainty in future work. Similarly, the inclusion of additional physics such as rotation, and additional modes such as those of higher degree or different azimuthal order, would be trivial extensions to this framework in future.

## Acknowledgements

SJM was supported by the Australian Research Council (ARC) through Future Fellowship FT210100485. MBN and GRD acknowledge support from the UK Space Agency. TRB acknowledges support from Australian Research Council through Laureate Fellowship FL220100117. OJS and AJL acknowledge the support of the Science and Technology Facilities Council. This paper has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (CartographY GA. 804752). This paper includes data collected by the _TESS_ mission. Funding for the _TESS_ mission is provided by the NASA's Science Mission Directorate.

## Data Availability

The data will be supplied upon reasonable request.

\begin{table}
\begin{tabular}{c c c} \hline frequency & \(n\) & \(\ell\) \\ \([\mathrm{d}^{-1}]\) & & \\ \hline
33.48997 & 3 & 0 \\
46.46549 & 5 & 0 \\
53.51968 & 6 & 0 \\
60.60165 & 7 & 0 \\
67.65639 & 8 & 0 \\
35.99870 & 3 & 1 \\
50.04504 & 5 & 1 \\
57.18150 & 6 & 1 \\
64.19957 & 7 & 1 \\ \hline \end{tabular}
\end{table}
Table 2: Identified modes for HD 99506 used as inputs in the modelling.

Figure 8: The posterior distributions in \(M\), \(Z_{\mathrm{in}}\) and age for HD 99506, based on _TESS_ observations.

Figure 9: The posterior predicted frequencies overlaid on observed mode frequencies for HD 99506. The greyscale is the observed amplitude spectrum smoothed by a Gaussian of width 4 times the frequency resolution.
## Software

Below we include additional software used in this work which has not explicitly been mentioned above.

* Python Van Rossum and Drake Jr (1995)
* matplotlib Hunter (2007)
* Numpy Harris et al. (2020)
* Scipy Virtanen et al. (2020)
* Pandas Reback et al. (2020)
* corner Foreman-Mackey (2016)
* lightkurve Lightkurve Collaboration et al. (2018)
* echelle Hey and Ball (2020)
