Fast and automated oscillation frequency extraction using Bayesian multi-modality

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.


INTRODUCTION
Despite both CoRoT (Baglin et al., 2006) and NASA Kepler and K2 missions 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).

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 DIAMONDS 1 , 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.
1. 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 ∼ 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 law 2 . 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. 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 = 2, 0, 1 triplet typical of main sequence solar-like pulsators. This was done through the fitting of the small frequency spacings between adjacent = 2, 0 modes and adjacent = 0, 1 modes and of a frequency centroid ν 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.

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 which I refer to as the islands model, where A is the amplitude of the oscillation peak, Γ its FWHM, and ν 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 ν 0 are considered as free parameters since Γ is fixed to a value that is related to the actual ν 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 Γ 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 3 for more discussion about the choice of Γ 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.

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 ν 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 ν 0 comparable to the large frequency separation ∆ν 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 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 ∼ 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 live = 2000 and performed the fit until the termination condition of δ 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 δ 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 2.4, where the sampling from DIAMONDS terminates for five different peaks after about 2000 iterations since being sampled 3 , 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 live = 2000 and a high threshold stopping point, set to N 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).

Figure 1. Left panel:
Sampling evolution from DIAMONDS of the frequency centroid ν 0 (orange dots) as a function of the nested iteration, for the red giant KIC 12008916. Each nested iteration corresponds to one sampling point. The sampling presented is built as follows: 1) DIAMONDS performed a total of N nest = 6000 iterations (hence 6000 sampling points) using 2000 live points; 2) the run is completed by adding up the set of 2000 live points that is left at the end of the nested sampling, thus resulting in a total of 8000 effective nested iterations (or equivalently sampling points); 3) the sampling from the first 1500 nested iterations was removed to improve the detection of the peak structures in the PSD. This results in a final set of 6500 nested iterations, hence of sampling points used for the analysis. Right panels: (top) The PSD of the corresponding chunk analyzed (in gray) with a 10 frequency-bin smoothing overlaid (blue curve) to highlight the oscillation peak structures. (bottom). The counts histogram from the sampling, with local maxima corresponding to the extracted oscillation frequencies indicated by vertical dotted red lines, as identified by a hill climbing algorithm. The different regions of the histogram, denoted with varying colors, represent the ranges of frequency automatically selected to compute the final estimates of the oscillation frequencies and their corresponding uncertainties, according to the method (b) presented in Figure 2. The estimated frequencies and corresponding 1-σ uncertainties from method (b) are indicated by the red circles and error bars, respectively. The horizontal dashed line represents the threshold level for the selection of the peaks, set to 3 % of the maximum counts in the histogram.

Counts histogram for the extraction of individual oscillation frequencies
By using the sampling cloud obtained from DIAMONDS for the free parameter ν 0 , we can proceed with the construction of a corresponding counts histogram. I adopt a resolution per bin in ν 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 ∆ν 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.

Comparison with Approach 1 for RGB oscillations
The example shown in Figure 1 is illustrating the extraction process for a frequency range of length ∆ν 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  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 ∼ 100 in speed.

DISCUSSION
The methodology presented in Section 2, building on a previous result by C14, was able not only to find all the actual oscillation frequency peaks in a range of 25 µ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 Figure 2. Comparison between the oscillation frequencies extracted from DIAMONDS using the multimodal approach (ν 0 ) and those published by C15 (ν Corsaro15 ). The values are reported in percentage of deviation with respect to the literature values. The uncertainty error bars obtained with the new approach are also overlaid for each oscillation frequency, after they were rescaled by ν Corsaro15 . The horizontal lines at (ν 0 − ν Corsaro15 )/ν Corsaro15 = 0 % denote the perfect matching. The panels from left to right depict the result obtained by considering four different methods to extract the frequencies from the counts histogram: (a) using a symmetric frequency range for each local maximum, with an extent on each side of the maximum that is the minimum found between the two distances going from the maximum to the adjacent midpoints from the next and previous maxima (or the edge of the total frequency range) on the right and left side respectively. The final frequencies are computed as the sampling mean from the sampling falling in the selected range of each local maximum, while the corresponding uncertainty is the standard deviation from the same sampling used to compute the mean; (b) using the same symmetric ranges as (a) but this time performing a weighted mean and weighted standard deviation using the nested iteration of each sampling point as a weight; (c) using an asymmetric frequency range, built with the distances between the maximum and the adjacent midpoints from the next and previous maxima (or edge of the total frequency range) on the right and left side respectively. The mean and standard deviation of the selected sampling of the range are used as in (a) to compute the final frequency and uncertainty estimates from each local maximum; (d) same as the case (c) but this time using a weighted mean and standard deviation as done for the case (b). 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 ν max and ∆ν 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 ν 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 ν 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.