Systematics in Asteroseismic Modelling: Application of a Correlated Noise Model for Oscillation Frequencies

Understanding the systematics that affect the oscillation frequencies in a stellar grid is the key to finding proper functional forms with which to describe the model systematics. As seen in Figure 1, the model systematics can be described as a combination of two components. The first is the frequency-dependent offset between the best-fitting model and the observed frequencies. Instead of calling it the surface term, we refer to it as the model systematic error ($\mathcal{E}$) to represent the systematics caused by all improper and missing physics in a theoretical stellar model. We use a smooth functional form, similar to a surface correction formulae, to model $\mathcal{E}$.

The second component, as mentioned in Section 1, is the systematic uncertainty ($\mathcal{U}$) caused by the frequency resolution of the model grid. The detailed modelling uses global parameters (effective temperature, luminosity, metallicity, etc.) and individual mode frequencies to characterise stars. Comparing the models in Figure 1 we observe that, at a given mean density, changing the mass by a typical mass interval in a grid shifts the mode frequencies horizontally by an amount which increases smoothly with frequency. Changing one of the other model inputs like metallicity, helium fraction, and mixing-length parameter shifts mode frequencies in a similar way. There is also a secondary term to be considered in the systematic uncertainty related to the signature of rapid structural variations in the oscillation frequencies (known as the helium ‘glitch’ signature). The structural variation can be seen in the first adiabatic index ($\Gamma_{1}$). Gough (1990) and Houdek & Gough (2007) assumed that the helium glitch signature arises from the second helium ionisation zone. Later studies indicate that this signature is from the $\Gamma_{1}$ peak between the first and second helium ionisation zones (see Houdayer et al., 2021, for details on modelling of the ionisation region). The helium glitch strongly correlates with the helium fraction in the convective envelope. In Figure 2, we illustrate the signature of the helium glitches extracted from theoretical models with approximately the solar mean density. The glitch signatures follow the sinusoidal function with decaying amplitudes and the signature parameters (amplitude, period, and phase) change with the fundamental inputs (Houdek & Gough, 2007; Verma et al., 2019). For the models presented here, the average frequency shift is at a level of $\sim$0.1$\mu$Hz, which are comparable or larger than the observed frequency uncertainties on some Kepler stars (e.g., Li et al., 2020a) thus should not be ignored. Hence, the systematic uncertainty can be described with a two-term functional form. The primary term ($\mathcal{U}_{1}$) which is a very smooth function of frequency and the secondary term ($\mathcal{U}_{2}$) which is a fast-varying function of frequency.

As a result, we describe the model systematics as the combination of the systematic errors ($\mathcal{E}$) and two systematic uncertainty terms ($\mathcal{U}_{1}$ and $\mathcal{U}_{2}$).

As demonstrated above, proper descriptions of systematic errors and uncertainties are smooth functions of frequency. This is expected because mode frequencies are highly correlated following the so-called asymptotic relation. Thus, a white noise model for the systematic noise term is inherently not a good model. Here we suggest a noise model with a smooth function form which is able to consider the correlation between mode frequencies, and we refer to it as correlated noise model (CNM).

The Gaussian Process (GP) kernel, which is used to generate a covariance matrix for a multivariate Normal distribution, is particularly suitable for building up the CNM. Rasmussen & Williams (2006) consider a GP in a functional form domain whereby the conditional GP acts to marginalize over all possible functional forms weighted by the prior (i.e. the kernel) and the data. Here we chose a Squared Exponential (SE) kernel as it follows a smooth functional form and is able to consider the correlation between oscillation frequencies. The SE kernel is described by the covariance matrix

The kernel $k(\nu,\nu^{\prime})$ models the joint variability of the Gaussian Process random variables. It returns the modelled covariance between each pair of frequencies $\nu$ and $\nu^{\prime}$. When using the kernel to represent the systematic uncertainties, $\nu$ are the computed mode frequencies in the model grid and $\nu^{\prime}$ are the predicted frequencies based on $\nu$. There are two free parameters in the function: the lengthscale $l$ and the noise variance $\sigma$. The lengthscale describes how smooth a generated function is (length of the ‘wiggles’), and the noise variance specifies the average distance of the function away from its mean. Tuning these two parameters allows us to generate functional noise forms with different smoothness and varying ranges. In sampling process, the kernel randomly generates frequency noise ($\nu$ - $\nu^{\prime}$) for a given $\nu$ following the multivariate normal distribution specified by the lengthscale and the noise variance.

In Figure 3, we demonstrate some randomly generated noise realizations using Eq 1 based on a set of computed model frequencies ($\nu$). We use a lengthscale of 5$\Delta\nu$ and a variance of 1.0$\mu$Hz. As seen, the predicted frequency sets ($\nu^{\prime}$) change in a highly correlated way. We also show a set of random samples from the white noise model (WNM) for comparison. The WNM follows a Gaussian distribution and has the same variance value ($\sigma$ = 1.0$\mu$Hz). The predicted frequencies based on the WNM include spiky features, which obviously do not well reflect the true curvature of radial mode frequencies. We suggest here that a CNM based on the SE kernel is the better representation for model systematics because it is more consistent with our prior belief in the nature of the oscillation frequencies.

Here we discuss the determination of the mean function ($\mu_{\mathcal{E}}$), error kernel lengthscale ($l_{\mathcal{E}}$), and error kernel variance ($\sigma_{\mathcal{E}}$). The systematic error kernel is normally determined from frequency differences between observations and the best-fitting model. However, the best-fitting model is unknown to us a priori. Instead, we use specific knowledge of model errors of similar stars to characterise the systematic error and determine three adjusted parameters. In what follows, we justify our choices but note that the method itself could incorporate alternative choices.

Previous studies of the surface correction on stars similar to the Sun provide some useful references for model errors of the Sun. Compton et al. (2018) estimated the surface terms of 66 Kepler main-sequence stars. For stars within a parameter range of $T_{\text{eff}}$ = 5777 $\pm$ 250K and $\log g$ = 4.44 $\pm$ 0.5dex, the relative surface corrections at $\nu_{\rm max}$ ($\delta\nu(\nu_{\rm max})$) vary from -0.0015$\nu_{\rm max}$ to -0.0045$\nu_{\rm max}$. We use this prior information in the form of a weight $w_{\rm\nu_{max}}$ given by a super Gaussian function with the exponent raised to a power of 10 (flat-top Gaussian function)

Here, $\mu_{\nu_{\rm max}}$ and $\sigma_{\nu_{\rm max}}$ are -0.003$\nu_{\rm max}$ and 0.0015$\nu_{\rm max}$, respectively. This likelihood function returns a weight that is flat when $\delta\nu(\nu_{\rm max})$ is in the $\mu_{\nu_{\rm max}}\pm\sigma_{\nu_{\rm max}}$ range, and quickly drops toward zero at $\mu_{\nu_{\rm max}}\pm 1.5\sigma_{\nu_{\rm max}}$. Moreover, previous studies (Ball & Gizon, 2014; Compton et al., 2018; Li et al., 2020b) have also inferred that the model frequency errors at the low-frequency range, below $\sim$0.7$\nu_{\rm max}$, are not as significant as those in high-frequency range. In sun-like stars, the low-frequency-range errors are on average approximately zero with a small amount of 0.001$\nu_{\rm max}$ spread. This can be used as another constraint and we describe it as a Gaussian function as

where $\widetilde{\delta\nu(\nu_{\rm obs})}$ represents the median of frequency offsets for given observed frequencies, $\sigma_{\rm low-\nu}$ is adjusted and defines for how much the $\widetilde{\delta\nu(\nu)}$ deviating from zero is sensible. According to earlier studies (Ball & Gizon, 2014; Compton et al., 2018), we adopt $\sigma_{\rm low-\nu}$ = 0.001$\nu_{\rm max}$.

Now we demonstrate how these two functions are used to characterise the model errors of the Sun. Firstly, we compare observed and model frequencies to calculate frequency differences. We then use Eq. 5 and 6 to obtain a joint weight ($w_{\rm low-\nu}$$\cdot$$w_{\rm\nu_{max}}$) for each model. As a simple example, we use stellar models on an 1M_⊙ evolutionary track to demonstrate this estimation in the left panel of Figure 4. When the joint weight of each model is obtained, we find a couple of potential good-fitting models (with weight value larger than 0.01) on this track and plot their frequency differences as a function of the frequency. The weight distribution indicates the model error is about -5$\mu$Hz at 3,000 $\mu$Hz and goes down to about -10$\mu$Hz at 4,000$\mu$Hz. Now we could fit a cubic polynomial function to all frequency differences in Figure 4 against the frequency weighted by the joint weight to estimate the mean of frequency offsets (the solid line). Moreover, we calculate the weighted standard deviations of frequency differences (as illustrated by the grey shade). The fitted polynomial function represents the mean estimate of systematic error in these theoretical models, and the standard deviation infer the range for it to vary. We thus use them as the mean function and noise variance in the systematic error kernel. In our CNM, the lenghtscale determines the degree of correlation between mode frequency noise. By inspecting the surface correction results obtained on the Kepler LEGACY sample (Compton et al., 2018), we notice that the frequency offsets between model and observations mostly follow slowly-varying smooth functions (higthly-correlated). This is to say the error kernel lengthscale ought to be much larger than the large separation. We test a selection of lengthscale values and compare the kernel predictions with the surface correction results. We find a suitable lengthscale of 5$\Delta\nu$ because the kernel with this lengthscale value well reproduces both the very smooth surface terms expected for stars with approximately one solar mass and the slightly curved ones, expected for the more massive F-type stars(see Fig. 8 in Compton et al., 2018). Using the mean function, variance, and lengthscale obtained above, we illustrate some random draws from the systematic error kernel in the right panel of Figure 4.

Now we determine the lengthscales and variances of the two uncertainty kernels. The two lengthscales and variances are different because the primary kernel describes the general frequency change (a very smooth function of frequency) and the secondary kernel corresponds to the helium glitch signatures (a fast-varying function of frequency). The uncertainty kernels represent the systematic uncertainty at a grid point and their free parameters depend on the local frequency changes between the point and its neighbouring points. Thus, the free parameters in the uncertainty kernels vary for different grid points.

Here we use a simple method to determine free parameters for the grid point at $M$ = 1M_⊙, [Fe/H] = 0.0, $Y_{\rm init}$ = 0.26, $\alpha_{\rm MLT}$ = 2.1, and $\rho$ = $\rho_{\odot}$. A model grid always contains multiple input dimensions. For the grid in this work, the independent inputs are mass, initial metallicity, initial helium abundance, and initial mixing-length parameters. We therefore inspect the frequency changes for each fundamental input (see details in Figure 8). We find that the primary uncertainty kernel needs a large lengthscale to make it act as a very smooth function. We test different lengthscale from 10$\Delta\nu$ to 30$\Delta\nu$ and adopt $l_{\mathcal{U},1}$ =20$\Delta\nu$ because it best recovers the general shape of frequency changes in Figure 8. On the other hand, we use a small lengthscale for the secondary uncertainty kernel to match the quick variation of the glitch signature. By inspecting the second-order variations in Figure 8, we notice that the lengths of the ‘wiggles’ are mostly between 2–3 radial orders. The choice of lengthscale needs to be in a similar range to reproduce the fast variations. We hence test different lengthscale values from 1$\Delta\nu$ to 4$\Delta\nu$ and find the kernel with $l_{\mathcal{U},2}$ = 2$\Delta\nu$ best recovers the test case. Considering the frequency resolution for all input dimensions, we find an average varying range is $\sim$0.75$\mu$Hz at 0.5$\nu_{\rm max}$, $\sim$1.5$\mu$Hz at $\nu_{\rm max}$, and $\sim$2.0$\mu$Hz at 1.5$\nu_{\rm max}$. The variance is hence frequency dependent and we use $\sigma_{\mathcal{U},1}$ = $1.5\nu_{\rm obs}/\nu_{\rm max}$. The secondary kernel variance is also frequency-dependent given that the signature of the helium glitch is a damped sine wave. As shown in Figure 2, the secondary uncertainty is up to $\sim$0.4$\mu$Hz at 0.5$\nu_{\rm max}$ and gradually reduces to $\sim$ 0 at 1.5$\nu_{\rm max}$. We therefore set up the secondary variance as $\sigma_{\mathcal{U},2}$ = 0.1($\nu_{\rm max}/\nu_{\rm obs}$)². Note that the method of measuring frequency differences and estimating the kernel variation is fairly rough. Firstly, the kernel variation is not uniform through a model grid. Secondly, we measure frequency changes in each input dimension independently, but the four fundamental inputs highly degenerate. The grid resolution of oscillation frequency and the free parameters of the uncertainty kernels should be solved as multiple-dimesions problems. We leave these in future studies.

In Figure 5, we use the two uncertainty kernels to generate some frequency noise as a function of frequency. As seen, the primary uncertainty kernel generates very smooth and highly correlated noise, and what the secondary uncertainty kernel provides is similar to the damped sine wave. The combination of the two kernels gives reasonable predictions of frequency noise between stellar model grids.