# PBjam: A Python Package for Automating
Asteroseismology of Solar-like Oscillators

M. B. Nielsen, G. R. Davies, W. H. Ball, A. J. Lyttle, T. Li, O. J. Hall, W. J.

Chaplin, P. Gaulme, L. Carboneau, J. M. J. Ong, et al.

► To cite this version:

M. B. Nielsen, G. R. Davies, W. H. Ball, A. J. Lyttle, T. Li, et al., PBjam: A Python Package for Automating Asteroseismology of Solar-like Oscillators. The Astronomical Journal, 2021, 161, 10.3847/1538-3881/abcd39. insu-03745718

HAL Id: insu-03745718

https://insu.hal.science/insu-03745718v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# PBjam: A Python Package for Automating Asteroseismology of Solar-like Oscillators*

M. B. Nielsen¹,²,³ Ⓓ, G. R. Davies¹ Ⓓ, W. H. Ball¹,² Ⓓ, A. J. Lyttle¹,² Ⓓ, T. Li (李坦达)¹,² Ⓓ, O. J. Hall¹,² Ⓓ, W. J. Chaplin¹,² Ⓓ, P. Gaulme⁴ Ⓓ, L. Carboneau¹,² Ⓓ, J. M. J. Ong (王加冕)⁵ Ⓓ, R. A. García⁶,⁷ Ⓓ, B. Mosser⁸ Ⓓ, I. W. Roxburgh¹,⁹ Ⓓ, E. Corsaro¹⁰ Ⓓ, O. Benomar³,¹¹ Ⓓ, A. Moya¹,¹² Ⓓ, and M. N. Lund² Ⓓ

¹ School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, UK; m.b.nielsen.1@bham.ac.uk

² Stellar Astrophysics Centre (SAC), Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark

³ Center for Space Science, NYUAD Institute, New York University Abu Dhabi, P.O. Box 129188, Abu Dhabi, UAE

⁴ Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, D-37077, Göttingen, Germany

⁵ Department of Astronomy, Yale University, 52 Hillhouse Avenue, New Haven, CT 06511, USA

⁶ IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

⁷ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

⁸ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 92195 Meudon, France

⁹ Astronomy Unit, School of Physics and Astronomy, Queen Mary University of London, London E1 4NS, UK

¹⁰ INAF-Osservatorio Astrofisico di Catania, via S. Sofia 78, I-95123 Catania, Italy

¹¹ Solar Science Observatory, NAOJ and Department of Astronomical Science, Sokendai (GUAS), Mitaka, Tokyo, Japan

¹² Electrical Engineering, Electronics, Automation and Applied Physics Department, E.T.S.I.D.I, Polytechnic University of Madrid (UPM), Madrid E-28012, Spain
*Received 2020 July 7; revised 2020 October 16; accepted 2020 November 22; published 2021 January 12*

## Abstract

Asteroseismology is an exceptional tool for studying stars using the properties of observed modes of oscillation. So far the process of performing an asteroseismic analysis of a star has remained somewhat esoteric and inaccessible to nonexperts. In this software paper we describe PBjam, an open-source Python package for analyzing the frequency spectra of solar-like oscillators in a simple but principled and automated way. The aim of PBjam is to provide a set of easy-to-use tools to extract information about the radial and quadropole oscillations in stars that oscillate like the Sun, which may then be used to infer bulk properties such as stellar mass, radius, age, or even structure. Asteroseismology and its data analysis methods are becoming increasingly important as space-based photometric observatories are producing a wealth of new data, allowing asteroseismology to be applied in a wide range of contexts such as exoplanet, stellar structure and evolution, and Galactic population studies.

*Unified Astronomy Thesaurus concepts: Astronomy software (1855); Publicly available software (1864); Asteroseismology (73); Stellar oscillations (1617); Astronomy data analysis (1858)*

## 1. Introduction

In the past few decades asteroseismology<fnref n="13" /> has become an important tool for characterizing stars. The frequencies of modes in which many different stars oscillate are sensitive probes of their physical properties (Aerts et al. 2010; Chaplin & Miglio 2013). These mode frequencies are used as observational constraints in stellar modeling (see, e.g., Appourchaux et al. 2012; Davies et al. 2016; Lund et al. 2017, for recent examples from Kepler), where they routinely allow for estimates of the physical properties of stars down to the percent level (see, e.g., Metcalfe et al. 2012; Lebreton & Goupil 2014; Silva Aguirre et al. 2017). The types of oscillations that a star can exhibit vary widely depending on its physical properties (Handler 2013), but in the following we will focus on stars that oscillate like the Sun, namely those with convective layers at the surface.

Asteroseismology uses time series from two main sources: either radial velocity or photometric intensity. These observations stem from what is now a multitude of different ground- and space-based observatories and projects, including the ground-based Global Oscillation Network Group (GONG; Harvey et al. 1996), Birmingham Solar-Oscillations Network (BiSON; Hale et al. 2016), and Stellar Observations Network Group (SONG; Grundahl et al. 2017) projects, and the Solar and Heliospheric Observatory (SOHO; Domingo et al. 1995),

CoRoT (Baglin et al. 2009), Kepler (Borucki et al. 2010), and Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) spacecraft. Analysis of the observations from this multitude of telescopes is in principle possible in the time domain (Foreman-Mackey et al. 2017), but this remains computationally expensive. Asteroseismology of stars that oscillate like the Sun is therefore typically performed in frequency space, where the oscillations modes are often readily visible.

In cool main-sequence (MS) stars like the Sun, as well as in subgiant (SG) and red giant (RG) stars, the oscillations that are most clearly visible are acoustic *p*-modes (García & Ballot 2019). These are stochastically excited, damped harmonic oscillations, which take the form of a set of semi-regularly spaced near-Lorentzian peaks in the power spectrum (Tassoul 1980; Anderson et al. 1990). The central frequencies of these peaks correspond to the resonance frequencies of the star, which can be compared to the frequencies of stellar structure models and so allow constraints to be placed on the properties of the star (e.g., Brown et al. 1994; Silva Aguirre et al. 2015; Angelou et al. 2017).

The measurement of the oscillation mode frequencies, also known as “peakbagging,” is the overall aim of PBjam.<fnref n="14" /> This is typically done by fitting a parametric model to the power spectrum of the photometric time series, for example, for each individual star (Appourchaux 2003; Davies et al. 2014; Corsaro

<footnote>* Release 1.0.0 Zenodo, doi:10.5281/zenodo.4300079.</footnote>

<footnote n="13">¹³ Helioseismology in the case of the Sun.</footnote>

<footnote n="14">¹⁴ Documentation and usage instructions are available at https://pbjam.readthedocs.io.</footnote>
& De Ridder 2014). This process involves two main parts: the mode identification and the model fitting, where the mode identification step informs the choice of the parametric model of the spectrum that is then fit to the observed oscillation spectrum.

Identification of the oscillation frequencies is often done manually and requires some knowledge of how the modes may appear at a given evolutionary stage of the star. The modes in most oscillating stars are described using spherical harmonic functions of angular degree *l* and azimuthal order *m*, each of which have a number of overtones of radial order *n*. The mode identification requires assigning unique labels (*n*, *l*, and *m*) to the modes that are visible in the spectrum. This is particularly difficult in SG and RG stars where the modes start coupling to the internal gravity modes (*g*-modes; see, e.g., Mosser et al. 2012) and thereby start rapidly varying in frequency as the star evolves. For F-type stars the mode widths increase substantially compared to cooler stars (Appourchaux et al. 2014), blending the *l* = 2, 0 modes in the spectrum and making them appear almost identical to the *l* = 1 modes, thus any distinction is difficult (Appourchaux et al. 2008; Benomar et al. 2009; White et al. 2012). Both these cases as well as simpler ones, like Sun-like stars, also become further complicated when the signal-to-noise ratio (S/N) of the observations is low.

Following the mode identification, the next step is fitting the chosen model to the spectrum of the oscillation frequencies. Finding the model that has the highest probability of explaining the observed spectrum can often be computationally expensive. This is especially true for observations made by the Kepler spacecraft, which observed more than 20,000 oscillating stars for several years with a cadence of ≈30 minutes (see, e.g., Yu et al. 2018) and several hundred with a shorter cadence of ≈1 minute (Lund et al. 2017; Serenelli et al. 2017). The spectra of many of these stars still need to be analyzed in detail. Furthermore, at the time of writing TESS is continuously producing new time series of stars across almost the entire sky, and with the future launch of the PLAnetary Transits and Oscillations of stars (PLATO; Rauer et al. 2014) mission, a fast and automated method of peakbagging is becoming increasingly important.

With PBjam we focus on solving these two issues for solar-like oscillators in an automated fashion. These stars are the most numerous type of oscillator, and include MS stars with masses $\lesssim 1.6 M_\odot$ and the vast majority of SG and RG stars. The former are of particular interest in, e.g., exoplanet studies (Lundkvist et al. 2016; Huber et al. 2019), and the latter for probing Galactic populations due to their visibility at great distances (Miglio et al. 2013; Mathur et al. 2016). In addition, with the automation in PBjam we also aim to make asteroseismology accessible to nonspecialists, thereby enabling its use in a wider range of contexts. The following refers to PBjam v1.0.0.<fnref n="15" />

## 2. Peakbagging with PBjam

The purpose of peakbagging is the measurement of the oscillation mode frequencies of a star. An example spectrum is shown in Figure 1, where the modes of the RG star KIC4448777¹⁶ are shown in terms of their S/N, which is the height of the mode peaks relative to the surrounding noise

floor. The modes are centered around a characteristic frequency, $v_{max}$, and their S/N drops off rapidly toward higher and lower frequencies. This is usually called the *p*-mode envelope.

The modes that are sufficiently excited to become visible vary depending on the physical properties of the star, as do the frequencies of the modes themselves. The number of visible radial orders depends on the S/N of the observations but typically range from just a few in the low S/N cases, to ≈10–15 in the best cases (Lund et al. 2017; Chaplin et al. 2020). The visibility of the modes also decreases with increasing *l* (see, e.g., Christensen-Dalsgaard & Gough 1982), and so usually only *l* = 0, 1 and 2 are observed and rarely *l* = 3 (Appourchaux et al. 2012). For each *l* there are 2*l* + 1 modes with an azimuthal order of $-l \le m \le l$, which may become split in terms of frequency due to rotation of the star (Gizon & Solanki 2003; Ballot et al. 2007). The spectrum model therefore typically has a large number of variables, and some kind of parameterization of the model is usually necessary.

Several methods of parameterization have been used previously to study stellar power spectra (see, e.g., Ballot et al. 2011b; Appourchaux et al. 2012; Handberg & Campante 2011; Lund et al. 2017; Corsaro 2019). With this version of PBjam the objective is to measure the mode frequencies of the *l* = 0, *m* = 0, and *l* = 2, *m* = 0 modes, which are highlighted in Figure 1. These modes contain information about the stellar properties for a large range of spectral types and evolutionary stages. The parameterization that is used in PBjam is therefore greatly simplified compared to other studies, as effects such as rotation and mode asymmetry are ignored (see, e.g., Davies et al. 2015; Benomar et al. 2018b, 2018a). Most significantly however, is that the *l* = 1 modes are not included in the PBjam peakbagging. These modes require special treatment to be fit accurately in an automated fashion across all evolutionary stages of a solar-like oscillator. This is particularly relevant for SG and RG stars, where *g*-modes in the deep stellar interior start to couple with the *l* = 1 *p*-modes, resulting in a complicated pattern of mixed mode frequencies, which makes even manual mode identification challenging and an automated approach more so (see, e.g., Appourchaux 2020). We aim to include the *l* = 1 treatment in future releases of PBjam, as they carry useful information about rotation and help to further constrain stellar ages (Metcalfe et al. 2010; Deheuvels & Michel 2011). However, the *l* = 0 cannot couple to any gravity mode, and strong coupling of the *l* = 2 *p*- and *g*-modes is only rarely observed (B. Mosser et al. 2021, in preparation), and so these mode pairs consistently appear with a near-regular pattern, at a rough interval called the large frequency separation, $\Delta\nu$. This makes them far easier to identify in an automated fashion, and only using these mode frequencies is sufficient to determine the stellar mass and radius to a precision of a few percent and to estimate the stellar age to within ~10%–20% (Davies & Miglio 2016; Buldgen et al. 2019; McKeever et al. 2019). Figure 1 shows the spectrum of the star KIC4448777 where the *l* = 0, 2 pairs are highlighted. The repeating pattern appears clearly by showing the mode frequencies modulo the large separation, where the *l* = 0 and *l* = 2 modes align as ridges. Only some of the *l* = 1 modes adhere to the same repeating pattern, while others are shifted due to the coupling with the internal gravity waves.

<footnote n="15">¹⁵ For the latest version see https://github.com/grd349/PBjam.</footnote>

<footnote n="16">¹⁶ Additional examples can be found in the PBjam GitHub repository.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left: S/N spectrum of the star KIC4448777 shown in light gray and the smoothed spectrum in dark gray. The spectrum consists of a regular series of radial ($l=0$) and quadrupole ($l=2$) modes shown by dashed lines. These are the modes that PBjam attempts to model. The overall power in the p-mode envelope can be approximated by a Gaussian centered on $\nu_{max}$. Right: the spectrum shown as an echelle diagram, where the frequency axis is shown modulo the large separation, $\Delta\nu = 16.97$ µHz, showing the repeating pattern of modes as distinct ridges. The mode frequencies shown as points in the left frame, with the same color code.</figcaption></figure>

The choice of spectrum model and the mode identification determines the set of fit parameters, $\theta$, and so to find the best-fit model we map the posterior probability

$$P(\theta|D) \propto P(D|\theta)P(\theta), \quad (1)$$

where $P(\theta)$ is the constraint on the fit parameters given any prior knowledge, and $P(D|\theta)$ is the probability of observing the data $D$ given the model. If we use a more informal syntax then $P(\theta)$ is what we know already about what the mode frequencies should be, and $P(D|\theta)$ is what the power spectrum or any other observations tell us about what the mode frequencies are. The posterior probability, $P(\theta|D)$, is what we ultimately care about: the combination of our existing knowledge and our new knowledge.

### 3. Core PBjam Classes

The main functionality in PBjam is the mode identification and model fitting, and it is currently comprised of three connected classes: KDE and Asy_peakbag, which perform the mode identification and encode the prior knowledge of the modes, and Peakbag which performs the final model fit but with minimal influence from the prior. Provided with a set of basic inputs shown in Table 1, these classes are run in sequence to produce an estimate of the mode frequencies in the spectrum. The three steps can be executed automatically by PBjam, but each may still be executed individually if needed. The automated process in PBjam is illustrated in Figure 2, with detailed descriptions in the following sections. The outputs of each step are summarized in Figure 4.

#### 3.1. KDE—Encoding Prior Knowledge

This section discusses the details of KDE. The two main functions of KDE are to construct a function that approximates our prior knowledge and to provide initial guesses for the fit parameters used in the following step in Asy_peakbag (see Section 3.2). The parameters, $\theta$, that KDE defines a prior for are therefore dependent on which inputs Asy_peakbag requires

<table><caption>**Table 1**
Fit Parameters in KDE and Asy_peakbag</caption><thead><tr><th>Parameter</th><th>Description</th></tr></thead><tbody><tr><td>$T_{\text{eff}}^{*\dagger}$</td><td>Effective surface temperature of the star.</td></tr><tr><td>$G_{\text{BP}} - G_{\text{RP}}^{*\dagger}$</td><td>Gaia photometric color index.</td></tr><tr><td>$\Delta\nu^{\dagger}$</td><td>Frequency difference of modes with the same $l$, but consecutive $n$.</td></tr><tr><td>$\nu_{\text{max}}^{\dagger}$</td><td>Frequency of maximum power of the p-mode envelope.</td></tr><tr><td>$\varepsilon$</td><td>Radial order phase term.</td></tr><tr><td>$\delta\nu_{02}$</td><td>Frequency difference of $l=0$ and $l=2$ modes.</td></tr><tr><td>$\alpha$</td><td>Scale (curvature) of radial order variation with frequency.</td></tr><tr><td>$H_{\text{max}}$</td><td>Height of the p-mode envelope.</td></tr><tr><td>$W_{\text{env}}$</td><td>Width in frequency of the p-mode envelope.</td></tr><tr><td>$w$</td><td>Width in frequency of the modes.</td></tr></tbody></table>

**Note.** Parameters marked with $\dagger$ are required input values which must be computed or estimated by the user in order to start KDE. Variables marked with $*$ are only used to remove degeneracies in the asymptotic relation, and are not part of the spectrum model itself.

(see Table 1). The methodology may, however, be applied to a variety of parameterizations of the spectrum model and we will therefore keep the notation in this section more general.

#### 3.1.1. Using a Kernel Density Estimate as a Prior

The first part of KDE is to construct a prior function that encodes our knowledge of how the mode frequencies behave at various stages of stellar evolution. This is done by computing a multivariate kernel density estimate (KDE) of $\theta$, based on fits to previously observed targets from Kepler.<fnref n="17" /> This sample is shown in Figure 3 for a subset of the parameters used in Asy_peakbag. Provided the sample of previous observations covers the physically meaningful range in parameter space, the resulting KDE captures the covariance of the different parameters in the fit. This yields a continuous function which approximates the prior, $P(\theta)$.

<footnote n="17">¹⁷ Future versions of PBjam will include CoRoT and TESS observations.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Diagram of the PBjam peakbagging process. The input parameters are shown in Table 1. The output oscillation mode parameters include the mode frequencies, heights, and widths along with their uncertainties.</figcaption></figure>

Constructing a KDE is fundamentally a data smoothing problem. Here we use a multivariate KDE with a bandwidth for each of the fit parameters, $\theta$, of the asymptotic relation and a few others (see Section 3.2). This gives full control over the degree of smoothing of the prior data. The KDE is constructed as

$$P(\theta) \propto \frac{1}{K} \sum_{i=1}^{K} Q_{H}(\theta - \theta_{i}), \quad (2)$$

where $\theta_i$ is the parameters of one of K previously fit stars used to construct the prior. $Q_H$ is the bandwidth matrix for which we use a multivariate Gaussian covariance matrix with all off-axis elements set to zero (sometimes called a D-type kernel):

$$Q_H(\theta) = (2\pi)^{-d/2} |\mathbf{H}|^{-1/2} \exp\left(-\frac{1}{2} \theta^T \mathbf{H}^{-1} \theta\right), \quad (3)$$

where d is the number of dimensions of the multivariate KDE. Here, $\mathbf{H}$ is a $d \times d$ diagonal matrix

$$\mathbf{H} = \mathrm{diag}(q_1^2, \dots, q_d^2), \quad (4)$$

where $q_i$ is a scalar and typically referred to as the bandwidth of one of the model parameters.

Computing the KDE is done using the statsmodels<fnref n="18" /> Python package (Seabold & Perktold 2010). This package also determines the optimal bandwidth for each fit parameter using a cross-validated maximum-likelihood approach. The cross validation is performed on a subset of the ~100 nearest neighbors of the target star in terms of $\nu_{max}$. This subset is chosen from a range of up to ±20σ, where σ is the uncertainty on the input $\nu_{max}$. This approach creates a KDE where the optimal bandwidth varies inversely with the local density of

prior data points, depending on the target star. The effect is to reduce the bandwidth in regions of the prior sample with many stars, while increasing the bandwidth in regions where we have few samples, i.e., little knowledge. The bandwidth is optimized every time a new fit is performed, and the prior is therefore expected to become more informative as we populate the different regions of the prior sample. We intend for this to be an ongoing process as new stars are continuously being observed and investigated. However, the bandwidth can be scaled by the user, which may be necessary where the prior sample is particularly sparse, for example at low $\nu_{max}$ (≈20 µHz).

An example of the resulting probability distributions of the $l=0$ modes determined by the KDE class are shown in the top frame of Figure 4.

### 3.1.2. Estimating the Most Probable Starting Point

The second function of KDE is to estimate the most probable region in parameter space for Asy_peakbag to start. This is done by sampling $P(\theta|D)$ shown in Equation (1). The prior, $P(\theta)$, is approximated using a KDE as shown above, which leaves the likelihood $P(D|\theta)$. At this stage we use the input parameters, shown in Table 1, as the observational data, $D$. These input parameters are: the frequency of the maximum power of the p-mode envelope, $\nu_{max}$; the frequency difference of the consecutive overtones, $\Delta\nu$; the effective temperature of the star, $T_{eff}$; and the Gaia photometric color index $G_{BP} - G_{RP}$ (Evans et al. 2018). The likelihood is the joint probability estimated by a series of normal distributions given by these input parameters and their uncertainties. The uncertainties on the inputs are thereby accounted for when estimating the initial starting location. Since $T_{eff}$ and $G_{BP} - G_{RP}$ contain much of the same information for characterizing a star, only one of these parameters needs to be provided. In the event that only one of these parameters is available, PBjam assumes a very wide, uninformative prior on the missing parameter.

Sampling the posterior, $P(\theta|D)$, is done using the affine-invariant Markov Chain Monte Carlo (MCMC) sampler from the emcee<fnref n="19" /> (Foreman-Mackey et al. 2013) package. The resulting percentile values of the marginalized posterior of each parameter are then passed to Asy_peakbag as the initial starting location in parameter space.

### 3.2. Asy_peakbag—Initial Guesses from the Asymptotic Relation

This section discusses the details of Asy_peakbag, which fits the asymptotic relation (see, e.g., Mosser et al. 2015) to the spectrum, using the output of the KDE class. This in turn provides the most probable frequency intervals for the final stage of peakbagging (see Section 3.3).

Asy_peakbag essentially performs a highly constrained peakbagging fit, which may in principle be used to estimate the mode frequencies on its own. However, the asymptotic relation used here is not a complete description of the mode frequencies, and so will not capture detailed variations related to, e.g., acoustic glitches (Mazumdar et al. 2014; Vrard & Cunha 2019) and coupling to g-modes in the deep interior of the star (Mosser et al. 2017). Furthermore, any error estimates on the asymptotic fit parameters will be highly correlated.

<footnote n="18">¹⁸ https://www.statsmodels.org/stable/index.html</footnote>

<footnote n="19">¹⁹ https://emcee.readthedocs.io/en/stable/</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. Corner plot of the parameters in the Asy_peakbag model. The diagonal shows the distribution of each model parameter for stars in the current sample used by the KDE class to generate a prior. The off-diagonal frames show the correlations between the different parameters. All parameters except ε are on logarithmic scales. While all 13,288 stars in our sample are shown, for clarity we only show a subset of the parameters. For a full list and description of the parameters see Table 1.</figcaption></figure>

Asy_peakbag is therefore best used to simply define a credible frequency range of the individual modes, and to identify the angular degree, for more detailed peakbagging later (for example using the Peakbag class).

In this case the spectrum model, M(θ, ν), consists of a sum of Lorentzian profiles, one for each visible mode in the spectrum (Anderson et al. 1990). For the purposes of PBjam we treat (n, l = 0) and (n − 1, l = 2) as a pair, where n represents the visible radial orders in the spectrum. The model for Asy_peakbag is then

$$
\begin{equation}
\begin{aligned}
M(\theta, \nu) = & b + \sum_{n=1}^{N} \frac{h_{n,0}}{1 + \frac{4}{w^2}(\nu - \nu_{n,0})^2} \\
& + \frac{h_{n-1,2}}{1 + \frac{4}{w^2}(\nu - \nu_{n-1,2})^2}.
\end{aligned}
\tag{5}
\end{equation}
$$

We fit Equation (5) to the S/N spectrum of a star and we therefore assume that the constant term b = 1. The S/N spectrum is the power spectrum divided by the background noise level, which is approximated by a running median of the

power spectrum.<fnref n="20" /> The precise treatment of the background noise level caused by granulation and long period variability may influence the measurement of the mode frequencies. However, the method applied here is sufficient for the purposes of setting credible frequency ranges for later detailed peakbagging.

Rather than fitting a series of free parameters for each mode in Equation (5), we use the asymptotic relation to encode our expectation of the mode pattern in solar-like oscillators. The mode frequencies are then given by

$$
\begin{align}
\nu_{n,0} &= \left( n + \varepsilon + \frac{\alpha}{2} (n - n_{\text{max}})^2 \right) \Delta\nu \nonumber \\
\nu_{n-1,2} &= \nu_{n,0} - \delta\nu_{02}, \tag{6}
\end{align}
$$

where ε is commonly referred to as a frequency offset, or phase term, and n_max = ν_max/Δν − ε (e.g., Kjeldsen et al. 2005). The parameter α is the scale of the second-order variation of

<footnote n="20">²⁰ See https://docs.lightkurve.org/api/lightkurve.periodogram.
Periodogram.html.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4. Outputs of the three main classes in PBjam: KDE, Asy_peakbag, and Peakbag. In all the frames the gray curve is the smoothed S/N spectrum of KIC4448777. Top: the output of KDE is a probability distribution showing the location of the radial ($l=0$) modes, based on what the prior sample predicts. This stage is independent of the observed spectrum. Middle: the asymptotic relation is fit to the spectrum, using the probability distribution from KDE as an input. A sample of 50 models drawn from the posterior distribution of the fit is shown in red. Bottom: the mode frequencies define the frequency ranges to be considered in the final peakbag step, where each $l=0, 2$ pair is then fit independently. For each pair 50 models drawn from the posterior distributions are shown in red.</figcaption></figure>

the modes frequencies between each radial order. The frequency of the quadrupole, $l=2$, modes are offset from the radial modes by $\delta\nu_{02}$, which is assumed to be constant for all the radial orders in the fit.

The relative mode heights, in terms of S/N, are approximated by a Gaussian envelope,

$$h_{n,0} = H_{\max} \exp(-0.5(\nu - \nu_{\max})^2 / W_{\text{env}}^2), \\ h_{n-1,2} = 0.7 h_{n,0}, \qquad (7)$$

where $W_{\text{env}}$ is the p-mode envelope width and $H_{\max}$ is the envelope height. The heights of the $l=2$ modes are scaled by a factor of 0.7, relative to those of the $l=0$ modes (Ballot et al. 2011a).

This approximates the reduced visibility due to geometric cancellation as the degree $l$ of the modes increases.

The widths of the Lorentzian profiles are known to be a function of $T_{\text{eff}}$, $\nu_{\max}$, and mode frequency (see, e.g., Appourchaux et al. 2014). However, precisely modeling these requires several additional parameters. This adds complexity to the model, which for the purposes of Asy_peakbag is unnecessary, and we therefore simply approximate the mode widths by a single constant value of $w$ for all modes.

With Asy_peakbag we also evaluate the posterior shown in equation Equation (1), however here we add the additional constraints from the power spectrum, i.e., the model fit. This means that the logarithm of the posterior probability can be written as

$$\ln P(\theta|D) \propto \ln \mathcal{L}(\theta) + \ln P(\theta), \qquad (8)$$

where $P(\theta)$ is given by Equation (2) and $\mathcal{L}(\theta)$ is the observational constraints given by

$$\ln \mathcal{L}(\theta) = \ln \mathcal{L}_{\text{S}}(\theta) + \ln \mathcal{L}_{\text{O}}(\theta). \qquad (9)$$

The log-likelihood $\ln \mathcal{L}(\theta)$ is the constraint from the power spectrum. The power in a frequency bin of the power spectrum is gamma distributed with a shape parameter $\alpha = 1$ and scale parameter $\beta = 1/M(\theta, \nu)$, and so the log-likelihood may be computed by (see, e.g., Woodard 1984; Duvall & Harvey 1986)

$$\ln \mathcal{L}_{\text{S}}(\theta) = -\sum_{j=1}^{J} \left( \ln M(\theta, \nu_j) + \frac{S_j}{M(\theta, \nu_j)} \right), \qquad (10)$$

where $S_j$ is the power in frequency bin $j$, and $J$ is the total number of frequency bins in the spectrum.

Here, the log-likelihood $\ln \mathcal{L}(\theta)$ is from the additional observational parameters where, similarly to KDE, Asy_peakbag also uses $T_{\text{eff}}$ and $G_{\text{BP}} - G_{\text{RP}}$. However, $\nu_{\max}$ and $\Delta\nu$ are now only used as fit parameters, with the only constraint on these parameters coming from the prior and the spectrum, and not from the user input. This prevents double counting the information from the input $\nu_{\max}$ and $\Delta\nu$ and the spectrum itself, as the former is often derived from the same spectrum that is being fit. The constraints from $T_{\text{eff}}$ and $G_{\text{BP}} - G_{\text{RP}}$ are again simply a sum of normal distributions with mean values and standard deviations equal to the inputs provided by the user.

The KDE that represents the prior, $P(\theta)$, is constructed from a sample of 13,288 Kepler stars ranging from RG stars with a $\nu_{\max} \approx 30$ µHz to MS stars at $\nu_{\max} \approx 4000$ µHz (see Figure 3)
that were previously fit using Equation (5). This sample<fnref n="21" /> is compiled from the White et al. (2011), Silva Aguirre et al. (2015), Serenelli et al. (2017), Lund et al. (2017), and Yu et al. (2018) catalogs, with the intent that more are added as they become available.

To sample the parameter space spanned by $\theta$, PBjam can use either the emcee or CPNest<fnref n="22" /> (Veitch et al. 2017) packages, either of which provide an estimate of the posterior probability distribution $P(\theta|D)$. The marginalized posterior distributions of the fit parameters are then passed to the next step, Peakbag.

An example of fitting the asymptotic relation to the spectrum of KIC4448777 can be seen in the middle frame of Figure 4.

### 3.3. Peakbag—Pair-by-pair Mode Fitting

The final part of the process is to relax the majority of the parameterization used in Asy_peakbag, thereby providing a relatively unconstrained estimate of the mode frequencies and their uncertainties. This allows PBjam to capture small variations in the mode frequencies that are not accounted for by the asymptotic relation.

We start our description of the Peakbag spectrum model by once again only considering the $l=0, 2$ pairs as is done in with Asy_peakbag, and we can therefore use Equation (5). However, in the Peakbag model we only consider the frequency range near the $l=0, 2$ pairs, and not the range in between the pairs that is typically occupied by the $l=1$ modes. This is done to speed up this final peakbagging step. In the following, the $l=0, 2$ pairs are therefore treated independently, and for simplicity we will simply denote each ($n, l=0$), ($n-1, l=2$) pair as ($n, 0$), ($n, 2$), so that $n$ is now the pair number and not the radial order. We then have a model, $M_n(\nu)$, for each mode pair and a small segment of the spectrum surrounding them (see bottom frame of Figure 4).

In Peakbag we use PyMC3 (Salvatier et al. 2016) to sample the posterior distribution. PyMC3 is a probabilistic programming language and so we will adopt a consistent notation here. This allows us to set up random variables to describe the model parameters. The mode frequencies of a pair $n$ are then given by

$$ \nu_{n,l} \sim \mathcal{N}(\nu_{n,l}^{\text{asy}}, (0.03\Delta\nu)^2) \quad \text{for } l=0, 2, \quad (11) $$

which denotes $\nu_{n,0}$ and $\nu_{n,2}$ as random variables that are distributed according to a normal distribution with mean mode frequencies, $\nu_{n,l}^{\text{asy}}$, from Asy_peakbag, and variances of $(0.03\Delta\nu)^2$.

If we were to turn Equation (11) into a PyMC3 model and draw samples from it, we would get posterior distributions that are normally distributed with mean values equivalent to the Asy_peakbag predictions and standard deviations that are 3% of $\Delta\nu$.

We use the same arrangement for mode heights and widths. However, these parameters can only take on positive values, and so to encode this knowledge we use log-normal distributions instead. As the prior mean values we again use the estimates from Asy_peakbag. However, for the mode heights we use a standard deviation of 0.4 dex, and 1 dex for the mode line widths. The large range for the latter is required since the line widths vary widely with frequency across the p-

mode envelope, and may be much smaller or larger than the approximate mean value estimated by Asy_peakbag.

The mode heights are then

$$ h_{n,l} \sim \log \mathcal{N}(\log h_{n,l}^{\text{asy}}, 0.16) \quad \text{for } l=0, 2. \quad (12) $$

Similarly for the line widths

$$ w_{n,l} \sim \log \mathcal{N}(\log w^{\text{asy}}, 1.0) \quad \text{for } l=0, 2. \quad (13) $$

We construct a prior for the background $b_n$ in a similar fashion. As above, we are working with the S/N spectrum and so the background should be unity. However, the flattening becomes an increasingly bad approximation as the S/N of the modes increases. This is typically a small effect, but to provide a more robust result we allow the background term for each pair of modes to vary independently, while still being consistent with the same log-normal distribution. We therefore set up the background prior as

$$ b_n \sim \log \mathcal{N}(0, 0.16). \quad (14) $$

The final step of building our PyMC3 model is adding the constraint from the spectrum itself. As with Asy_peakbag we assume that the power in each bin is drawn from a gamma distribution, however, here the scale parameter is $\beta = 1/M_n(\nu)$, which defines a model for each mode pair $n$. In PyMC3 terminology the spectrum constraint is then imposed by

$$ S_n \sim \gamma(1, 1/M_n(\nu)). \quad (15) $$

This completes our PyMC3 model, and we use the No U-turn Sampler (NUTS) to draw samples from the posterior probability distribution. An example of the models sampled by Peakbag is shown in the bottom frame of Figure 4. The output is in the form of summary statistics of each parameter posterior distribution. As a guide to the user, one of the output metrics is the ratio of the prior (Equation (11)) and the posterior widths of the mode frequencies. For modes where this ratio is ~1 the inference is influenced mainly by the prior, and likewise when the ratio is >1 the spectrum provides the most information, with the former often being the case for low S/N modes, and the latter for high S/N modes. While PBjam does not make any selection based on this metric, we suggest that modes with prior to posterior width ratios greater than ≈2 can be used to select modes where the spectrum dominates the inference. However, even when the inference is informed predominantly by the prior, the resulting mode frequency may still be used as a constraint on stellar models. It then simply reflects the estimate gained from the asymptotic relation shown in Equation (6), but with larger uncertainties than if the spectrum also contributes to the inference.

## 4. Conclusion

The current version of PBjam is suitable for measuring the $l=0$ and $l=2$ mode frequencies from MS stars up to RG stars. This is done using the expected pattern of the mode frequencies, provided by the asymptotic relation and the wealth of asteroseismic data provided by space missions such as Kepler. Apart from providing a few input parameters, the mode frequencies are peakbagged in a completely automated fashion. This allows nonexpert users to obtain seismic constraints for a wide range of problems that involve solar-like oscillators.

This sample of prior targets used by PBjam contains 13,288 stars observed by Kepler. As shown in Figure 3 this sample currently contains predominantly RG stars, and only a few MS stars. This is due to telemetry restrictions on the Kepler

<footnote n="21">²¹ Available in machine readable format https://github.com/grd349/PBjam/blob/master/pbjam/data/prior_data.csv.</footnote>

<footnote n="22">²² https://johnveitch.github.io/cpnest/</footnote>
spacecraft, which limited the number of targets that could be observed in high-cadence mode. MS stars typically have oscillation frequencies ≳1000 µHz, and so oscillations were only detected in a few hundred of these stars. It is our goal to keep updating the sample of observed targets as they are observed and peakbagged. The TESS mission is promising to contribute substantially in this respect, as the observing strategy covers almost the entire sky, and includes many bright SG and MS stars.

The main limitation of the current version of PBjam is that it does not consider the *l* = 1 modes. While the *l* = 0, 2 mode pairs alone are sufficient to constrain stellar models to a precision of a few percent in mass and radius and 10%–20% in age, the *l* = 1 mode frequencies will provide much tighter constraints. The coupling of the *l* = 1 *p*-modes to the internal gravity-dominated oscillations and how these change as the star evolves is the main issue with fitting these modes. Previous work by Kallinger (2019), Corsaro et al. (2020), and Appourchaux (2020) among others, have sought to address this issue in an automated fashion. These methods tend to use mode-by-mode significance tests of the peaks in the spectrum, followed by a comparison with the asymptotic relation to identify both the *l* = 0, 2 pairs and subsequently the *l* = 1 modes.

The aim of PBjam is to create a generative model that leverages the information from the whole *p*-mode envelope at once, and which scales evenly from the MS up to the RG branch. This has largely been achieved for the *l* = 0, 2 pairs in the current version, by first establishing the prior frequency range from the asymptotic relation, and then subsequently releasing those constraints. However, applying this approach for fitting the *l* = 1 modes requires that the knowledge of the expected mode pattern is accurate. While this is comparatively straightforward for MS stars, and to an extent RG stars (Vrard et al. 2016), it is less so for SGs. For these stars even small errors in the *l* = 1 peakbagging can have large effects on the resulting stellar parameters (see, e.g., Li et al. 2020). Extending the application of prior information, as is done in PBjam, to potentially incorporate the approaches mentioned above is one possible solution. We leave this investigation to the next release of PBjam.

M.B.N., W.H.B., and W.J.C. acknowledge support from the UK Space Agency. G.R.D., O.J.H., and W.J.C. acknowledge the support of the UK Science and Technology Facilities Council (STFC). E.C. acknowledges support from PLATO ASI-INAF agreement No. 2015-019-R.1-2018 and by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 664931. P.G. acknowledges funding from the German Aerospace Center (Deutsches Zentrum für Luft- und Raumfahrt) under PLATO Data Center grant 500O1501. J.O. acknowledges support from TESS GI grant G022092. R.A.G. acknowledges support from the GOLF and PLATO/CNES grants. This paper has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (CartographY GA. 804752). The authors acknowledge use of the BlueBEAR HPC service at the University of Birmingham. Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (grant agreement No. DNRF106). This paper includes data collected by the Kepler mission and

obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

**Software:** astropy (Astropy Collaboration et al. 2018), corner (Foreman-Mackey et al. 2016), CPNest (Veitch et al. 2017) emcee (Foreman-Mackey et al. 2013), matplotlib (Hunter 2007), numpy (Oliphant 2006), lightkurve (Barentsen et al. 2019), pandas (Reback et al. 2020), Python (Van Rossum & Drake 1995) pymc3 (Salvatier et al. 2016), scipy (Virtanen et al. 2020), sklearn (Pedregosa et al. 2011), statsmodels (Seabold & Perktold 2010).

## ORCID iDs

M. B. Nielsen https://orcid.org/0000-0001-9169-2599

G. R. Davies https://orcid.org/0000-0002-4290-7351

W. H. Ball https://orcid.org/0000-0002-4773-1017

A. J. Lyttle https://orcid.org/0000-0001-8355-8082

T. Li (李坦达) https://orcid.org/0000-0001-6396-2563

O. J. Hall https://orcid.org/0000-0002-0468-4775

W. J. Chaplin https://orcid.org/0000-0002-5714-8618

P. Gaulme https://orcid.org/0000-0001-8330-5464

L. Carboneau https://orcid.org/0000-0003-1001-5137

J. M. J. Ong

(王加冕) https://orcid.org/0000-0001-7664-648X

R. A. García https://orcid.org/0000-0002-8854-3776

B. Mosser https://orcid.org/0000-0002-7547-1208

I. W. Roxburgh https://orcid.org/0000-0002-7403-2764

E. Corsaro https://orcid.org/0000-0001-8835-2075

O. Benomar https://orcid.org/0000-0001-9405-5552

A. Moya https://orcid.org/0000-0003-1665-5389

M. N. Lund https://orcid.org/0000-0001-9214-5642

## References

Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Dordrecht: Springer)

Anderson, E. R., Duvall, T. L., Jr., & Jefferies, S. M. 1990, ApJ, 364, 699

Angelou, G. C., Bellinger, E. P., Hekker, S., & Basu, S. 2017, ApJ, 839, 116

Appourchaux, T. 2003, Ap&SS, 284, 109

Appourchaux, T. 2020, A&A, 642, A226

Appourchaux, T., Antia, H. M., Benomar, O., et al. 2014, A&A, 566, A20

Appourchaux, T., Chaplin, W. J., García, R. A., et al. 2012, A&A, 543, A54

Appourchaux, T., Michel, E., Auvergne, M., et al. 2008, A&A, 488, 705

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Baglin, A., Miglio, A., & Michel, E. 2009, in AIP Conf. Ser. 1170, Stellar Pulsation: Challenges for Theory and Observation, ed. J. A. Guzik & P. A. Bradley (Melville, NY: AIP), 310

Ballot, J., Barban, C., & van't Veer-Menneret, C. 2011a, A&A, 531, A124

Ballot, J., Brun, A. S., & Turck-Chièze, S. 2007, ApJ, 669, 1190

Ballot, J., Gizon, L., Samadi, R., et al. 2011b, A&A, 530, A97

Barentsen, G., Hedges, C., Vinčius, Z., et al. 2019, KeplerGO/lightkurve: Lightkurve, v1.5.0, Zenodo, doi:10.5281/zenodo.3549023

Benomar, O., Baudin, F., Campante, T. L., et al. 2009, A&A, 507, L13

Benomar, O., Bazot, M., Nielsen, M. B., et al. 2018a, Sci, 361, 1231

Benomar, O., Goupil, M., Belkacem, K., et al. 2018b, ApJ, 857, 119

Borucki, W. J., Koch, D., Basri, G., et al. 2010, Sci, 327, 977
Brown, T. M., Christensen-Dalsgaard, J., Weibel-Mihalas, B., & Gilliland, R. L. 1994, ApJ, 427, 1013

Buldgen, G., Rendle, B., Sonoi, T., et al. 2019, MNRAS, 482, 2305

Chaplin, W. J., & Miglio, A. 2013, ARA&A, 51, 353

Chaplin, W. J., Serenelli, A. M., Miglio, A., et al. 2020, NatAs, 4, 382

Christensen-Dalsgaard, J., & Gough, D. O. 1982, MNRAS, 198, 141

Corsaro, E. 2019, FrASS, 6, 21

Corsaro, E., & De Ridder, J. 2014, A&A, 571, A71

Corsaro, E., McKeever, J. M., & Kuszlewicz, J. S. 2020, A&A, 640, A130

Davies, G. R., Broomhall, A. M., Chaplin, W. J., Elsworth, Y., & Hale, S. J. 2014, MNRAS, 439, 2025

Davies, G. R., Chaplin, W. J., Farr, W. M., et al. 2015, MNRAS, 446, 2959

Davies, G. R., & Miglio, A. 2016, AN, 337, 774

Davies, G. R., Silva Aguirre, V., Bedding, T. R., et al. 2016, MNRAS, 456, 2183

Deheuvels, S., & Michel, E. 2011, A&A, 535, A91

Domingo, V., Fleck, B., & Poland, A. I. 1995, SoPh, 162, 1

Duvall, T. L., Jr., & Harvey, J. W. 1986, in Proc. NATO Advanced Research Workshop, Seismology of the Sun and the Distant Stars, ed. D. O. Gough (Dordrecht: D. Reidel), 105

Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Vousden, W., Price-Whelan, A., et al. 2016, corner.py: corner.py v2.0.0, Zenodo, doi:10.5281/zenodo.53155

García, R. A., & Ballot, J. 2019, LRSP, 16, 4

Gizon, L., & Solanki, S. K. 2003, ApJ, 589, 1009

Grundahl, F., Fredslund Andersen, M., Christensen-Dalsgaard, J., et al. 2017, ApJ, 836, 142

Hale, S. J., Howe, R., Chaplin, W. J., Davies, G. R., & Elsworth, Y. P. 2016, SoPh, 291, 1

Handberg, R., & Campante, T. L. 2011, A&A, 527, A56

Handler, G. 2013, in Planets, Stars and Stellar Systems (Vol. 4) ed. T. Oswalt & M. A. Barstow, Vol. 4 (Dordrecht: Springer), 207

Harvey, J. W., Hill, F., Hubbard, R. P., et al. 1996, Sci, 272, 1284

Huber, D., Chaplin, W. J., Chontos, A., et al. 2019, AJ, 157, 245

Hunter, J. D. 2007, CSE, 9, 90

Kallinger, T. 2019, arXiv:1906.09428

Kjeldsen, H., Bedding, T. R., Butler, R. P., et al. 2005, ApJ, 635, 1281

Lebreton, Y., & Goupil, M. J. 2014, A&A, 569, A21

Li, T., Bedding, T. R., Christensen-Dalsgaard, J., et al. 2020, MNRAS, 495, 3431

Lund, M. N., Silva Aguirre, V., Davies, G. R., et al. 2017, ApJ, 835, 172

Lundkvist, M. S., Kjeldsen, H., Albrecht, S., et al. 2016, NatCo, 7, 11201

Mathur, S., García, R. A., Huber, D., et al. 2016, ApJ, 827, 50

Mazumdar, A., Monteiro, M. J. P. F. G., Ballot, J., et al. 2014, ApJ, 782, 18

McKeever, J. M., Basu, S., & Corsaro, E. 2019, ApJ, 874, 180

Metcalfe, T. S., Chaplin, W. J., Appourchaux, T., et al. 2012, ApJL, 748, L10

Metcalfe, T. S., Monteiro, M. J. P. F. G., Thompson, M. J., et al. 2010, ApJ, 723, 1583

Miglio, A., Chiappini, C., Morel, T., et al. 2013, MNRAS, 429, 423

Mosser, B., Elsworth, Y., Hekker, S., et al. 2012, A&A, 537, A30

Mosser, B., Pinçon, C., Belkacem, K., Takata, M., & Vrard, M. 2017, A&A, 600, A1

Mosser, B., Vrard, M., Belkacem, K., Deheuvels, S., & Goupil, M. J. 2015, A&A, 584, A50

Oliphant, T. E. 2006, A Guide to NumPy, Vol. 1 (Spanish Fork, UT: Trelgol Publishing)

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, JMLR, 12, 2825

Rauer, H., Catala, C., Aerts, C., et al. 2014, ExA, 38, 249

Reback, J., McKinney, W., jbrockmendel, et al. 2020, pandas-dev/pandas: Pandas 1.0.0, v1.0.0, Zenodo, doi:10.5281/zenodo.3630805

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Proc. SPIE, 9143, 914320

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, e55

Seabold, S., & Perktold, J. 2010, in 9th Python in Science Conf., ed. S. van der Walt & J. Millman (Austin, TX: SciPy), 92

Serenelli, A., Johnson, J., Huber, D., et al. 2017, ApJS, 233, 23

Silva Aguirre, V., Davies, G. R., Basu, S., et al. 2015, MNRAS, 452, 2127

Silva Aguirre, V., Lund, M. N., Antia, H. M., et al. 2017, ApJ, 835, 173

Tassoul, M. 1980, ApJS, 43, 469

Van Rossum, G., & Drake, F. L., Jr 1995, Python Tutorial (Amsterdam: Centrum voor Wiskunde en Informatica Amsterdam)

Veitch, J., Del Pozzo, W., Cody, Pitkin, M. & ed1d1a8d 2017, johnveitch/cpnest: Minor Optimisation, version 0.1.4, Zenodo, doi:10.5281/zenodo.835874

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vrard, M., & Cunha, M. S. 2019, Influence of Structural Discontinuities Present in the Core of Red-giant Stars on the Observed Mixed-mode Pattern and Characterization of Their Properties, v1, Zenodo, doi:10.5281/zenodo.2528750

Vrard, M., Mosser, B., & Samadi, R. 2016, A&A, 588, A87

White, T. R., Bedding, T. R., Gruberbauer, M., et al. 2012, ApJL, 751, L36

White, T. R., Bedding, T. R., Stello, D., et al. 2011, ApJ, 743, 161

Woodard, M. F. 1984, PhD thesis, University of California, San Diego.

Yu, J., Huber, D., Bedding, T. R., et al. 2018, ApJS, 236, 42