Fast and automated oscillation frequency extraction using Bayesian multi-modality

TL;DR

This paper introduces a rapid Bayesian method for extracting individual oscillation mode frequencies in stars, achieving high precision and accuracy, which is crucial for analyzing large datasets from space missions like TESS and PLATO.

Contribution

A novel, fast Bayesian approach building on previous methods to accurately determine stellar oscillation frequencies from space-based observations.

Findings

Achieves frequency precision of about 0.1%

Attains accuracy of about 0.01% compared to literature

Effective for red giant branch stars

Abstract

Since the advent of CoRoT, and NASA Kepler and K2, the number of low- and intermediate-mass stars classified as pulsators has increased very rapidly with time, now accounting for several $10^4$ targets. With the recent launch of NASA TESS space mission, we have confirmed our entrance to the era of all-sky observations of oscillating stars. TESS is currently releasing good quality datasets that already allow for the characterization and identification of individual oscillation modes even from single 27-days shots on some stars. When ESA PLATO will become operative by the next decade, we will face the observation of several more hundred thousands stars where identifying individual oscillation modes will be possible. However, estimating the individual frequency, amplitude, and lifetime of the oscillation modes is not an easy task. This is because solar-like oscillations and especially their evolved version, the red giant branch (RGB) oscillations, can vary significantly from one star to another depending on its specific stage of the evolution, mass, effective temperature, metallicity, as well as on its level of rotation and magnetism. In this perspective I will present a novel, fast, and powerful way to derive individual oscillation mode frequencies by building on previous results obtained with \diamonds. I will show that the oscillation frequencies obtained with this new approach can reach precisions of about 0.1 % and accuracies of about 0.01 % when compared to published literature values for the RGB star KIC~12008916.

Equations

$$P_{\mathrm{isla}}\left(\nu ;\Gamma \right)=\frac{A^2/(\pi \Gamma )}{1+4{\left(\frac{\nu -{\nu }_0}{\Gamma }\right)}^2}$$

Full text

\helveticabold

1 Keywords:

asteroseismology, oscillation mode analysis, stars: oscillations, stars: late-type, methods: statistical, methods: numerical, methods: data analysis, Bayesian

2 Introduction

Despite both CoRoT (Baglin et al., 2006) and NASA Kepler and K2 missions (Borucki et al., 2010; Koch et al., 2010; Howell et al., 2014) have ended, we keep gathering new datasets of oscillating stars from recently launched NASA TESS (Ricker et al., 2014), which is offering the opportunity to extract individual oscillation mode frequencies even with 27-days long observations (see e.g. Huber et al. submitted). The rapidly growing amount of stars with visible stellar oscillation modes, expected to encounter a significant step increase with the advent of the future ESA PLATO space mission (Rauer et al., 2014), poses the problem of facing an adequate type of analysis to fully characterize stellar oscillations on a large sample basis. The full characterization of the oscillation modes is a high priority task because it delivers the largest amount of most valuable information that we can extrapolate from each star, essential to understand stellar structure and evolution in detail (Aerts et al., 2010).

First works focusing on extracting and exploiting oscillation mode frequencies, amplitudes, and linewidths were already carried out (see e.g. Corsaro et al., 2015a, b; Benomar et al., 2015; Davies et al., 2016; Lund et al., 2017; Corsaro et al., 2017a; Handberg et al., 2017; García Saravia Ortiz de Montellano et al., 2018; Vrard et al., 2018), showing that an effort to tackle the problem is currently ongoing. Nevertheless the analysis proves to be challenging because of the large number of oscillation modes involved, although it is still possible to automatize it in some cases.

My view on the problem is that we require a simple approach that can be at the same time powerful, flexible to be adapted to different conditions, and fast in performing a peak bagging analysis for each star. This approach should also be accessible by the global asteroseismic community. For this purpose, I present a novel methodology, based on the public code Diamonds  (Corsaro and De Ridder, 2014, hereafter C14), to extract individual oscillation mode frequencies that is at the same time fast, easy to use, and accurate. For demonstrating its working principle and performance, I apply it to the RGB star KIC 12008916, for which a detailed peak bagging analysis is available (Corsaro et al., 2015a, hereafter C15).

3 Fast oscillation frequency extraction with Diamonds

An efficient and robust approach to the fitting of individual oscillation modes exploits the Bayesian inference based on the nested sampling Monte Carlo algorithm (Skilling, 2004; Sivia and Skilling, 2006, C14). This is because the nested sampling has the advantages of: 1) requiring a factor 100 less sampling points than standard Markov chain Monte Carlo methods to accurately perform a statistical inference; 2) computing the Bayesian evidence of the model, a key element to test models in the light of both their complexity and fit quality; 3) efficiently sampling likelihood distributions that exhibit phase changes, i.e. multiple local maxima.

With the development of Diamonds111https://github.com/EnricoCorsaro/DIAMONDS, a multi-platform publicly available Bayesian inference tool in C++11 (see also C14; C15; Corsaro, 2018, for general guidelines), the characterization of the multitude of oscillation modes typical of the low- and intermediate-mass stars has been made accessible through the nested sampling and the flexible code implementation. In C14 two completely different approaches to peak bagging for a main sequence star were presented, which I summarize below.

The first approach, uni-modal and high-dimensional (hereafter Approach 1), requires priors set up for each of the individual mode properties of frequency centroid, amplitude and linewidth. The resulting sampling of the parameter space will exhibit a single global maximum of the likelihood distribution. Approach 1 was also used extensively to perform peak bagging analysis in about 90 red giant stars (C15; Corsaro et al., 2017a). This procedure yields the most precise and accurate set of asteroseismic measurements, once appropriate prior distributions are obtained. However, the fit is necessarily high-dimensional (with three free parameters for each oscillation mode being included in the fit), and although Diamonds  can still afford a peak bagging fit up to about $k\sim 40$ dimensions without particular problems (see the new calibration for peak bagging by Corsaro et al., 2018), the computational efficiency as well as the computational time required to carry out an individual fit increases quite rapidly with increasing number of oscillation modes being fitted, with a time dependency that follows a $k^{2.4}$ power law222See the Diamonds  2018 User Guide Manual available in the same GitHub repository of Diamonds  for more details.. While this aspect can be overcome by dividing the power spectrum of the star into multiple chunks, which are considered and analyzed independently from one another, one still has to face the problem of setting up priors for each oscillation peak that should be fitted, which in practice represents the main limitation. 2. 2.

The second approach, originally proposed by C14, is multi-modal and low-dimensional (hereafter Approach 2). The initial recipe accounted for a model with three Lorentzian profiles trying to reproduce the $\ell=2,0,1$ triplet typical of main sequence solar-like pulsators. This was done through the fitting of the small frequency spacings between adjacent $\ell=2,0$ modes and adjacent $\ell=0,1$ modes and of a frequency centroid $\nu_{1}$ of the dipolar mode used to locate the triplet of peaks in the spectrum. With this peak bagging model, the resulting sampling exhibits a series of clustered regions in the parameter space, each one corresponding to a local maximum (island) of the likelihood (see Figure 13 of Corsaro and De Ridder, 2014, for a visualization). The approach proved to be able to recover the oscillation mode properties of 27 consecutive peaks of a main sequence star using only nine free parameters, with an average accuracy on the oscillation frequencies of about 0.2 % as compared to estimates from Approach 1 for the same star. However, setting up Diamonds  for this configuration appeared not to be straightforward, and even the computational time required to perform one run (about 1 hour in a 2.6 GHz CPU core for a Kepler  short cadence power spectrum accounting for 27 oscillation peaks) did not yield any significant gain with respect to performing Approach 1 to the same frequency region.

In the following sections, I describe how to revise the original Approach 2 presented by C14 to improve significantly its speed, simplicity in use, and also its reliability in producing accurate and precise estimates of individual oscillation frequencies.

3.1 Islands peak bagging model

The choice of the model to fit the so-called power spectral density (PSD), namely the power spectrum normalized by the frequency resolution of the dataset, requires that the output sampling from Diamonds  has to resemble that of a multi-modal likelihood distribution. In order to create a model that is at the same time as simple as possible so that it can be easily configured, and also effective in recognizing and fitting the actual oscillation peaks that are present in the stellar power spectra, I consider a peak bagging model consisting of a single Lorentzian profile

[TABLE]

which I refer to as the islands model, where $A$ is the amplitude of the oscillation peak, $\Gamma$ its FWHM, and $\nu_{0}$ its frequency centroid. The model is superimposed on a fixed background level that can be fitted in a previous step (see e.g. Corsaro et al., 2017b, for a formulation). It is important to note that in the islands model only $A$ and $\nu_{0}$ are considered as free parameters since $\Gamma$ is fixed to a value that is related to the actual $\nu_{\mathrm{max}}$ of the star. By fixing the FWHM of the peak we can stabilize the fit and obtain a better resolving power of the actual oscillation peak structures that are present in the dataset. The choice of $\Gamma$ proves not to be particularly critical, but the general idea behind is that it should be set to a value equal to or below the minimum linewidth of an oscillation peak among those of a given star. I refer to Section 4 for more discussion about the choice of $\Gamma$ for the application presented in this work. With this islands model, we therefore account for only two free parameters, clearly expecting a considerable gain in speed when performing the fit. In addition, because of the adoption of just one Lorentzian profile and the only need for a rough assumption about the linewidth and amplitude of the observed peak structures, the sampling from Diamonds  proves to be significantly more stable with respect to the original Approach 2 presented by C14, thus producing more reliable results.

3.2 Multi-modal sampling at high threshold

Once the islands model is constructed, one can perform the actual fit to the PSD by providing as uniform prior distribution for the frequency centroid $\nu_{0}$ an input range that spans the actual region of the PSD that we intend to analyze, and a uniform prior range for the amplitude that resembles the average level of the amplitude of the peaks in the spectrum. This is the way we practically make the islands model a peak bagging model producing a multi-modal likelihood distribution. In the test performed, an input range for $\nu_{0}$ comparable to the large frequency separation $\Delta\nu$ offers an optimal choice in terms of resulting frequency resolution of the sampling for detecting all the oscillation peaks reported by C15.

Before performing the fit, it is essential to provide some considerations on the parameters that configure the nested sampling algorithm of Diamonds  (see C14; C15; Corsaro, 2018; Corsaro et al., 2018, for more details). In Corsaro et al. (2018) in particular, it was proved that using a number of $N_{\mathrm{live}}=500$ live points is adequate for Approach 1 to converge to an accurate solution even if a large number of dimensions is involved (up to $k\sim 60$). This relatively low number of live points can sustain the fit because of the uni-modality nature of the problem. In a multi-modal application instead, as indicated by C14, we require a high resolution of live points in the parameter space to be able to locate the individual local maxima of the likelihood distribution. This is because we know a priori that the likelihood will contain several, if not many, different local maxima. In this regard, C14 adopted a value of $N_{\mathrm{live}}=2000$ and performed the fit until the termination condition of $\delta_{\mathrm{final}}=0.01$ was reached (see also Keeton, 2011, fore more discussion). While the choice of the number of live points seems adequate also for the new application, that of adopting a low threshold $\delta_{\mathrm{final}}$ constitutes a severe problem. The result is that not only the fit takes a long time to be performed in its entirety (many nested iterations are required), but that it is also very likely that the computation will fail before reaching the end because the sampling has to take place from all or most of the local maxima at each attempt to draw a new sampling point. Clearly, the more we approach to the threshold, the more difficult it will be to find a new point that satisfies the new higher likelihood constraint.

For the reasons just presented, I have implemented another termination condition in Diamonds  that can take into account a fixed number of nested iterations. This allows to decide when to stop the computation depending on how many sampling points have been drawn. In this way, one is independent of the actual evidence estimated from the sampling, which can change significantly from one fit to another. Given that each nested iteration progressively corresponds to an increasing likelihood value, the peaks in the PSD that have a low height compared to the level of the background, quickly disappear after popping up in the sampling, while the peaks that have a larger height in the PSD are sampled for a longer time because they yield a larger likelihood value to the fit. This can be visualized in the left panel of Figure 1, for the application discussed in Section 3.4, where the sampling from Diamonds  terminates for five different peaks after about 2000 iterations since being sampled333The nested iteration number reported in Figure 1 does not correspond to the real nested iteration number used by Diamonds  during its sampling. The reader can refer to the caption of Figure 1 for more details on how this final sampling is obtained., while continuing until the end for the other six peaks. This effect can be summarized as follows: 1) if we stop too late in the nested sampling process we tend to oversample those peaks that are bigger, thus killing out all the structures that have a small height in the PSD; 2) if we stop too early, we risk to visualize too many structures that belong to the noise realization of the dataset, thus crowding with noise peaks the actual set of oscillation peaks that we can extrapolate from the sampling. Adopting $N_{\mathrm{live}}=2000$ and a high threshold stopping point, set to $N_{\mathrm{nest}}=6000$, produces a very stable and reliable sampling of the structures observed in the PSD for the application presented in this work. All the other configuring parameters of Diamonds  remain unchanged and can be set according to Corsaro et al. (2018).

3.3 Counts histogram for the extraction of individual oscillation frequencies

By using the sampling cloud obtained from Diamonds  for the free parameter $\nu_{0}$, we can proceed with the construction of a corresponding counts histogram. I adopt a resolution per bin in $\nu_{0}$ given by 1/100 of the length of the parameter range, where each bin in frequency will contain the number of sampling points falling in the given bin. This proves to be adequate to detect at least 11 different peaks within a frequency range of length $\Delta\nu$ and to get rid of spurious structures arising from the noise. As a result, when the sampling matches an actual oscillation peak in the PSD, a corresponding peak pops up in the counts histogram. The subsequent step is to consider a detection threshold in counts in order to pick up only those peaks of the histogram that exceed the threshold. A value for the threshold of 3% of the maximum number of counts found in the histogram appear to provide an optimal condition for the test presented in this work in order to detect all the oscillation peaks reported by C15. Once the threshold is defined, one can rely on a relatively simple hill climbing algorithm to progressively gather the local maxima present in the histogram by starting from the left edge of the frequency range. The result for this application is shown in the right panel of Figure 1. The adoption of a counts histogram coupled with a hill climbing algorithm is an important step because it provides a stable visualization and identification of the peak structures sampled from the PSD. This in turn makes the process of extracting the oscillation frequencies fast, automated, and reliable, in contrast to the original Approach 2 presented by C14 where the actual peak structures had to be manually identified from the sampling cloud obtained by Diamonds.

3.4 Comparison with Approach 1 for RGB oscillations

The example shown in Figure 1 is illustrating the extraction process for a frequency range of length $\Delta\nu$ in the low-luminosity low-mass RGB star KIC 12008916, observed nearly continuously by Kepler  for more than 4 years and using the same PSD as in C15. This star has been adopted by C15 as a reference star to present the results of a peak bagging analysis based on Approach 1. Here I used the same range adopted in Figure 5 of C15 to demonstrate that the new Approach 2 is able to select all the actual oscillation frequency peaks identified by C15 (see Figure 1 top right panel), with the exception of the only peak that was not deemed significant by C15 according to the peak significance test. For deriving the final frequency estimates and corresponding uncertainties from the counts histogram I have adopted four different definitions, described in Figure 2. The different definitions have the purpose of testing the impact of: 1) including a weight in the sample average to improve the frequency accuracy, with a weight given by the number of the nested iteration, meaning that those sampling points that correspond to a higher likelihood value have more weight; 2) adopting a symmetric (or asymmetric) frequency range of the oscillation peak with respect to its highest point (in sampling counts) in order to test how the sampling at low likelihood values can bias the estimate of the oscillation frequency. In Figure 2 I compare 22 extracted frequencies in the interval $146.5$-$172.0$ $\mu$Hz (length 2$\Delta\nu$) by using the new Approach 2, to those published by C15. The agreement is remarkable for all definitions presented but it is the best when adopting case (b) of Figure 2 (hence weighting by the nested iteration and adopting the symmetric frequency range), where the new frequencies deviate on average by only 0.01 % and have precisions that are on average of about 0.1 %.

It is important to consider that the newly presented Approach 2 has the advantage that can be applied to the complex pattern of the RGB oscillations, as demonstrated in this work, contrary to the original Approach 2 that was limited to main sequence stars only. A direct comparison between the new and old Approach 2 is therefore not possible for this dataset, but one can in general expect to achieve precisions on the frequency estimates that are of the same order of magnitude given that both approaches are based on the sampling of a multi-modal posterior distribution and involve the use of a low number of free parameters.

Lastly, I note that the extraction of the oscillation frequencies using the new Approach 2 takes about 2 s on a single 2.6 GHz CPU core for the PSD chunk presented in Figure 1. This implies that the entire star, containing roughly 70 different oscillation modes, can be processed in about 15 s. By comparison with Approach 1, where the frequency extraction done by C15 took about 27 min for the entire star KIC 12008916, we gain an useful factor $\sim 100$ in speed.

4 Discussion

The methodology presented in Section 3, building on a previous result by C14, was able not only to find all the actual oscillation frequency peaks in a range of 25 $\mu$Hz containing 22 different oscillation modes of KIC 12008916 – both pressure and mixed dipolar modes and some of which extremely narrow in width and small in height as compared to the background level – but also to deliver frequencies that are as accurate as 0.01% for the entire set analyzed, and precise to a level of 0.1 %. This implies that these frequencies could be used for modeling purposes, although precisions do not reach the level of those extracted from a peak bagging analysis using Approach 1, which ranges from $10^{-2}$ to $10^{-3}$ % for the range of frequencies used in this work, as obtained by C15. This new Approach 2 has two important advantages: 1) it is extremely convenient in terms of computational speed, yielding a factor of about 100 improvement with respect to Approach 1; 2) it does not require any assumption on the location of the frequency peaks, but only a simple average estimate of the oscillation amplitude in a chunk and an estimate of a minimum linewidth for each star, making it an adequate approach also for the complex oscillation patterns found in evolved stars. For the test presented here, I have adopted a FWHM of three times the frequency resolution of the dataset, because this was shown by C15 to be the typical FWHM of a fully unresolved oscillation peak, which is the narrowest observed. Therefore Approach 2 is automated because it does not require any human supervision, with the only requirement of using the global oscillation properties of $\nu_{\mathrm{max}}$ and $\Delta\nu$ to set up the model and priors, proving that it is perfectly suitable for the analysis of a large sample of targets.

Approach 2 can be easily coupled to Approach 1 in an automated sequence where Approach 2 is used to rapidly define the number of peaks, locate their frequency position and build relatively narrow prior ranges that are then feed in Approach 1 to perform a full peak bagging analysis. As a result, also the peak significance test performed with Approach 1 can be fully automated. This is possible because, thanks to the high accuracy of the $\nu_{0}$ values extracted with Approach 2, the frequency priors can be set up as the boundaries given by the extracted uncertainties on the frequency centroids from Approach 2, thus eliminating the problem of locating the peaks before performing the actual fit. Future work will aim at investigating the performances of Approach 2, as well as assessing its calibration, in the light of different $\nu_{\mathrm{max}}$ values, level of background photon noise, actual linewidth of the peaks (e.g. in the presence of blending as for F-type stars), and of different frequency resolutions typical of the new observations from K2 and TESS.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Author Contributions

E.C. has developed the new features implemented in Diamonds, refined Approach 2 from C14, performed its calibration, testing, and the application to the RGB star.

Funding

E.C. is funded by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 664931.

Acknowledgments

I thank Dr. Joris De Ridder for helpful discussions.

Data Availability Statement

The dataset analyzed for this study can be found in MAST, https://archive.stsci.edu/kepler/.