The BAyesian STellar algorithm (BASTA): a fitting tool for stellar studies, asteroseismology, exoplanets, and Galactic archaeology

ABSTRACT

1 INTRODUCTION

Obtaining reliable properties of field stars is of paramount importance for many fields in astrophysics. An accurate characterization of exoplanets requires precise knowledge of the parent star radius and mass, the ultimate fate and evolutionary remnant of a star can only be determined if we know its initial mass, and the study of the formation and evolution of our Galaxy is incomplete without the distribution of stellar ages across the Milky Way. These are just examples of how determining fundamental properties of stars has become the cornerstone of research for a variety of subjects.

Empirical methods to determine physical characteristics of stars are restricted to targets where e.g. years of monitoring are feasible (eclipsing binaries), or their bright apparent magnitude allows the measurement of the angular diameter (interferometry). This severely limits the number of stars where these techniques are applicable, and demands the development of methods in which some measured quantities (e.g. stellar effective temperature, surface composition, and luminosity) are compared to model predictions to infer stellar properties (such as age).

The advent of large-scale stellar surveys providing a myriad of data for thousands of stars across the Galaxy has led to the further development of algorithms that can combine different data sets to determine stellar properties. These algorithms vary in the approach used to extract the final parameters of stars (e.g. machine-learning, neural networks, Bayesian inference), the method to determine uncertainties (e.g. confidence intervals, Gaussian errors, Monte Carlo sampling), and the set of stellar tracks or isochrones considered in the analysis. Moreover, the possible combinations of input data are different for these codes: while some rely on spectroscopic, photometric, and astrometric data (e.g. StarHorse (Queiroz et al. 2018), MADE (Das & Sanders 2019), others have the capability of including asteroseismic information [e.g. PARAM (Rodrigues et al. 2017), Isoclassify (Huber et al. 2017), and AIMS (Rendle et al. 2019)]. The latter point is of key importance: due to their dependence on the internal stellar structure, reproducing the observed pulsation properties allows for a determination of the stellar radius, mass, and age of solar-type stars and red giants to a level of precision that cannot be achieved when fitting only atmospheric properties. Asteroseismology has therefore become an invaluable tool for a large variety of studies thanks to the rapidly increasing amount of data available since the launch of the space-missions CoRoT, Kepler, and TESS, which will continue to ramp up as future missions such as PLATO 2.0 (Rauer et al. 2014) begin to acquire data.

The irruption of asteroseismology in the scene of stellar properties determination poses tremendous challenges to fitting algorithms due to the large variety of oscillation quantities that one can try to reproduce. In stars whose driving mechanism is stochastic excitation from their outer convective envelopes (called solar-like oscillators, see e.g. Chaplin & Miglio 2013, for a review), the so-called global asteroseismic parameters can almost always be determined if oscillations are detected. If the data are of sufficient quality, individual frequencies of oscillation (or combinations of them) can be reproduced in main-sequence stars, as well as modes of mixed character that dominate the information content in subgiants. More evolved red giants present a much richer spectrum of pulsation, where state-of-the-art fitting algorithms reproduce only a subset of the observed frequencies while also considering a probe of the stellar core in the form of a gravity-mode period spacing. It is clear that fully exploiting the richness of asteroseismic data from solar-like oscillators requires algorithms capable of fitting a large variety of oscillation properties, which become relevant at different evolutionary stages and are strongly dependent on the data quality available. Moreover, these data should be supplemented with knowledge of the stellar effective temperature, a measurement of chemical composition, and ideally a determination of luminosity or absolute magnitude from astrometric and photometric data.

With this in mind we have developed the BAyesian STellar Algorithm (BASTA), originally introduced in Silva Aguirre et al. (2015). BASTAis a fitting tool written in Python (Van Rossum & Drake 2009) designed to take advantage of the large variety of data obtained by large-scale ground-based surveys and space missions to precisely characterize stars. It has been built in a flexible way that allows the user to choose any combination of input data to be fit and, to the best of our knowledge, it is the fitting code that includes the largest number of global asteroseismic quantities and individual frequency diagnostics for solar-type stars, subgiants, and red giants. Other codes do not include information from individual oscillation frequencies (e.g. Isoclassify and PARAM), or glitch properties (e.g. AIMS), or simply do not take any asteroseismic input (e.g. StarHorse and MADE). Moreover, BASTAallows to simultaneously reproduce spectroscopic, photometric, astrometric, and asteroseismic data in a self-consistent manner, and is the only code where parallaxes can be fitted directly in addition to e.g. individual oscillation frequencies without the need of transforming the astrometric information into a luminosity estimate. BASTA can run using publicly available compilations of stellar isochrones or tailor-made sets of evolutionary tracks with a wide combination of input physics.

BASTAhas been extensively used to determine stellar properties of both asteroseismic and non-seismic exoplanet host stars discovered by Kepler (Silva Aguirre et al. 2015; Lundkvist et al. 2016; Bonomo et al. 2019), K2 (Johnson et al. 2018; Persson et al. 2018; Van Eylen et al. 2018a; Hjorth et al. 2019a ; Korth et al. 2019; Lund et al. 2019), tess(Gandolfi et al. 2018; Huber et al. 2019), and MASCARA (Talens et al. 2018; Hjorth et al. 2019b). The sample of precise asteroseismic parameters originally derived in Silva Aguirre et al. (2015) has enabled detailed studies of e.g. exoplanet eccentricities (Van Eylen & Albrecht 2015) and the radius gap (Van Eylen et al. 2018b; see also Fulton et al. 2017). Since then, the applications of BASTAhave been extended to a large variety of studies across fields of astrophysics such as characterization of asteroseismic targets (e.g. Serenelli et al. 2017; Silva Aguirre et al. 2017, 2020b; Stokholm et al. 2019), Galactic archaeology (e.g. Casagrande et al. 2016; Silva Aguirre et al. 2018; Spitoni et al. 2019, 2020; Nissen et al. 2020), open clusters (e.g. Lund et al. 2016; Stello et al. 2016; Arentoft et al. 2019), and the study of physical processes in stars such as rotation, convective overshoot, and magnetic activity (e.g. van Saders et al. 2016; Booth et al. 2017; Hjørringgaard et al. 2017). It has also been shown to be one of the most accurate pipelines available in tests using artificial data of main-sequence stars (Reese et al. 2016).

In this paper we introduce the public version of BASTA. We describe the Bayesian approach followed when determining stellar properties, the compilations of stellar track and isochrones available, the main features and capabilities included in the code, and present validation results of the fitting algorithm using artificial data.

2 THE BAYESIAN FRAMEWORK

The user can specify an informative prior on stellar mass as given by the initial mass function (IMF). There are several versions of IMF included in BASTA(Salpeter 1955; Miller & Scalo 1979; Kennicutt, Tamblyn & Congdon 1994; Scalo 1998; Kroupa 2001; Baldry & Glazebrook 2003; Chabrier 2003). To decrease the computation time, the user can pre-select a region of the grid for which the likelihoods are computed. This selection can be made on any available grid properties with user-defined tolerances. Technically, this is equivalent to assuming specific non-informative priors on certain stellar parameters.

3 GRIDS OF STELLAR MODELS

The functionalities of BASTArely on the use of collections of stellar evolution tracks or isochrones to extract the properties of stars by means of Bayesian inference. In its current version BASTA runs over publicly available sets of isochrones and tracks as well as custom-computed ones, which we process and store in Hierarchical Data Format version 5 (HDF5) and make available upon request. The functionalities and adopted input physics for each case are described in the following subsections.

3.1 BaSTI isochrones and tracks

Stellar properties can be determined with BASTA making use of the recently updated library ‘a Bag of Stellar Tracks and Isochrones’ (BaSTI¹). This compilation consists of evolutionary tracks and isochrones that are available for four distinct science cases defined by the inclusion of different physical processes as given in Table 1. We give a brief description of its main features in this section, and refer the reader to Hidalgo et al. (2018) and Pietrinferni et al. (2021) for additional details.

The prescription for convective core overshooting consists of instantaneous mixing beyond the region formally defined by the Schwarzschild criterion, keeping the radiative temperature gradient in this region. In the case of main-sequence models with convective core, the overshoot region is defined by the distance λ_ov × H_p, where H_p is the local pressure scale height and λ_ov is a free parameter. It has been set equal to 0.2, decreasing to zero when the mass decreases below a certain value. The approach used for decreasing λ_ov from its maximum value to zero depends on both the chemical composition and stellar mass (see section 2.2 of Hidalgo et al. 2018, for details). During the core-helium burning stage, regardless of the considered science case, core mixing is modelled by accounting for semiconvection and suppression of breathing pulses.

The science case including atomic diffusion follow the prescription of Thoul, Bahcall & Loeb (1994), while mass-loss is taken into account in the formulation of Reimers (1975) with an efficiency of the free parameter η set to 0.3. The temperature stratification in the outer stellar layers is given by the Vernazza, Avrett & Loeser (1981) formulation, but for the case of very-low mass stellar models for which outer boundary conditions based on accurate model atmospheres have been adopted (see Hidalgo et al. 2018; Pietrinferni et al. 2021, for a detailed discussion on this topic). For the adopted physical framework and solar heavy element distribution (see Section 4.2.1 below), the calibration of a Standard Solar Model (SSM) sets the value of the mixing length parameter α_MLT, and of the initial solar metallicity and He abundance. At the solar age, the BaSTI SSM² matches the solar luminosity and radius as well as the present (Z/X)_⊙ abundance ratio with a value for the mixing length parameter α_MLT = 2.006 (we refer to Hidalgo et al. 2018, for details about the properties of the BaSTI SSM).

The BaSTI stellar models are available both for the Caffau et al. (2011) solar heavy element mixture, and for an α −element enhanced mixture³ (⁠|$\rm [\alpha /Fe]=+0.4$|⁠). For each selected metallicity, the corresponding initial helium abundance has been obtained by adopting an He− enrichment ratio equal to ΔY/ΔZ = 1.31, and a primordial He abundance equal to Y_P = 0.247. From the calibration of an SSM the resulting initial solar abundances are Z_ini = 0.01721 and Y_ini = 0.2695. From the compilation of 21 metallicities initially available with the release of the BaSTI models, we have increased the resolution in chemical composition and age by interpolating metallicity points for all science cases using the available online routine provided by the BaSTI team.

The current release of the BaSTI library also contains evolutionary tracks spanning masses from 0.1 |$\rm M_{\odot }$| to 15 |$\rm M_{\odot }$|⁠. For a subset of these (ranging from |$\sim 0.7~\rm M_{\odot }$| to 4.0 |$\rm M_{\odot }$| with the lower limit being metallicity dependent) the full interior structure is also provided throughout the evolution. The number of structures stored⁴ at each mass varies with metallicity, ranging from ∼300 individual models at the high-mass end to more than 9000 structures at the low-mass end. These are stored in the standardized fgong format described in the website⁵ of the Aarhus Red Giants Challenge (Silva Aguirre et al. 2020a). The availability of interior structures allows us to compute oscillation frequencies and determine asteroseismic observables as described in Section 4.1.

Our grids of BaSTI isochrones and tracks reach ages up to 16 Gyr, which allow us to properly construct the posterior distributions of old stars without risking the appearance of an edge effect (see e.g. Valle et al. 2014, 2015).

3.2 Custom-computed evolutionary tracks

While the use of publicly available compilations of tracks and isochrones (such as BaSTI) makes it easy to compare the BASTA results with those from other fitting codes, it poses a limit on our flexibility to explore different combinations of input physics and numerical implementations. For this reason, we developed algorithms that calculate and process user-defined grids of evolutionary tracks with any combination of input physics and make them readily usable with BASTA. The limitations on the input physics are given by the features available in an evolutionary code, and currently we support the usage of GARSTEC (Weiss & Schlattl 2008) and the latest publicly available version of MESA (Paxton et al. 2011, 2013, 2015).

We have two different approaches for grid sampling. In the first, which is the conventional approach, we compute tracks on a predefined mesh of stellar parameters. The mesh points are typically equally spaced along each parameter. The stellar model grids calculated in this manner are known as Cartesian grids. In the second approach, we sample the parameter space uniformly using a quasi-random number generator described in Sobol (1967). Note that quasi-random number generators perform better than pseudo-random number generators, specifically in high-dimensional spaces, as the latter provides more clumpy distributions. We refer to grids computed in this way as Sobol grids.

We can generate Cartesian and Sobol grids over a number of stellar parameters including mass (M_ini), initial helium abundance (Y_ini), initial metallicity ([Fe/H]_ini), alpha enhancement (⁠|$\rm [\alpha /Fe]$|⁠), mixing-length (α_MLT), and overshoot and mass-loss (η). The parameters used to generate the grid define the dimension of the weight |$w_{\boldsymbol {\Theta }}$| included in our marginalized posterior distribution (see equation 5). Note that we construct grids over [Fe/H]_ini (instead of initial metal mass fraction), because this quantity is well constrained by the observed metallicity, and hence allows convenient choice of the parameter space over which one needs to calculate a grid to model the observed star. We can either treat Y_ini as a free parameter similar to other stellar parameters, or determine it from the [Fe/H]_ini assuming values for the primordial helium abundance (default Y_p = 0.248; Fields et al. 2020), and the helium-to-metal enrichment ratio (default ΔY/ΔZ = 1.4). The default values can be changed by the user. In the case of GARSTEC, we have implemented alpha enhanced and depleted mixtures of the Asplund et al. (2009) solar abundances, ranging from |$\rm [\alpha /Fe]=-0.2$| to |$\rm [\alpha /Fe]=0.6$| in steps of 0.1. We adopt consistent opacity tables from OPAL for high-temperatures (Iglesias & Rogers 1996) and the Ferguson et al. (2005) opacities in the low-temperature regime.

4 AVAILABLE FITTING PARAMETERS

One of the core features of BASTAis its ability to handle inhomogeneities in the available data as long as they are contained in the grid model parameters. BASTA can easily deal with heterogeneous input and provide a robust set of stellar properties based on the likelihood of the models in its grids. In the following sections, we describe the quantities available for fitting asteroseismic, photometric, spectroscopic, and astrometric data. A complete and up-to-date list of code parameters can be found in the code documentation (available on GitHub).

4.1 Asteroseismology

The wealth of data from main-sequence and red giant stars provided by asteroseismic space missions has driven the development of BASTAtowards the study of solar-like oscillators. These are stars whose oscillations are excited by the same mechanism as in the Sun, and comprise the vast majority of targets where asteroseismic quantities are available. For solar-like oscillators, there are several observables that the user can select to be reproduced by the models and fit with BASTA. In all grids of stellar models currently supported (see Section 3), theoretical oscillation frequencies have been computed using the Aarhus adiabatic oscillation package (adipls; Christensen-Dalsgaard 2008).

4.1.1 Global asteroseismic quantities

There is extensive literature devoted to the testing and validating the scaling relations by using independent constraints in stellar masses and radii from e.g. binary stars, clusters, interferometry, and parallaxes (see e.g. Huber et al. 2012, 2017; Miglio et al. 2012; Silva Aguirre et al. 2012; White et al. 2013; Gaulme et al. 2016; Brogaard et al. 2018; Sahlholdt & Silva Aguirre 2018; Sahlholdt et al. 2018, to name a few). There is general agreement that the scaling relations are accurate to within a few percent, but as it is clear from equations (8) and (9) these extrapolations from the solar values do not take into account variations with chemical composition nor the evolutionary stage of the star.

To address some of these issues, several prescriptions have been proposed to correct the 〈Δν〉 scaling relation, while the equation for ν_max is purely empirical and no corrections across the range of interest in T_eff and |$\rm [Fe/H]$| have been derived yet (but see Belkacem et al. 2011, for an initial explanation on the theory behind equation 9). We have implemented in BASTAfour additional determinations of the average large frequency separation that aim at further decreasing the level of systematic deviation in the scaling relation equation (8). A detailed comparison between the performances of each prescription can be found in Viani et al. (2019).

The first two determinations consist of a linear fit as a function of radial order to the model individual frequencies of ℓ = 0 weighted by a Gaussian centred at ν_max, the only difference being the adopted Full Width at Half Maximum of 0.25ν_max (White et al. 2011) or |$0.66\nu _\mathrm{max}^{0.88}$| (Mosser et al. 2012a). The average large frequency separation 〈Δν〉 is determined from the slope of this fit and is meant to mimic as close as possible the manner in which this quantity is derived from the observations. This average large frequency separation is available for the grids of stellar models where we compute the individual oscillation frequencies (see Sections 3.1 and 3.2). As a by-product of this calculation we obtain the dimensionless offset ϵ that can be used to correct for differences between the observed and model radial order (see Section 4.1.2 below).

The other two determinations are those of Sharma et al. (2016) and Serenelli et al. (2017) who, based on the same principle of the White et al. (2011) approach, have computed a correction factor across the Hertzprung–Russell Diagram for the value of 〈Δν〉 obtained from equation (8) that depends on the mass, metallicity, effective temperature, and evolutionary state (hydrogen-core or-shell burning or core-helium burning) of the star. These corrections can be computed for any grid of stellar tracks or isochrones, as they are independent of the availability of theoretical oscillation frequencies.

As a final remark, we note that when defining the solar reference values 〈Δν_⊙〉 and ν_{max, ⊙} in equations (8) and (9) there is an implicit assumption that fitting a target with those values of 〈Δν〉 and ν_max (and solar temperature and metallicity) will result in a star of 1 M_⊙ and 1 R_⊙ (but not necessarily solar age, as this depends on the input physics used to construct the models). To ensure this level of consistency, all theoretical values of 〈Δν〉 in our grids of models are re-scaled by a factor given by the fraction 〈Δν_⊙〉/〈Δν_⊙〉_grid, where 〈Δν_⊙〉_grid is the average large frequency separation computed from the individual ℓ = 0 modes of a solar model computed with the same input physics of the corresponding grid using the full width at half maximum (FWHM) of White et al. (2011).

4.1.2 Individual oscillation frequencies

If the time-series of observations is of sufficient signal-to-noise ratio and resolution, it is possible to extract individual oscillation frequencies characterized by radial order n and angular degree ℓ (see e.g. Davies et al. 2016; Lund et al. 2017). For stellar disc-integrated observations, as it is the case of the space missions, geometrical cancellation suppresses the signal from all modes except those with low degree ℓ ≤ 3. If the asymptotic theory can be applied to describe the oscillation (e.g. Tassoul 1980), modes of odd and even degree are separated by ∼〈Δν〉/2. As a result, the observed oscillation spectrum contains one mode of each degree ℓ = 1, 2, 3 between two consecutive ℓ = 0 modes. Departures from this asymptotic description occur e.g. in the presence of mixed modes (see Section 4.1.5).

Most current stellar evolutionary models use a rudimentary description of the outer-most layers of stars with convective envelopes (e.g. the mixing-length theory of Böhm-Vitense 1958), leading to systematic frequency shifts in the oscillation modes when compared to observations (Christensen-Dalsgaard et al. 1996). In order to correct for this effect, a so-called surface correction needs to be applied to the model frequencies. There are three prescriptions for the surface corrections implemented in BASTA and all include power-law dependence on frequency. They aim at obtaining the corrected model frequencies |$\nu _{n,\ell }^{\mbox{cor}}$| from the original model frequencies |$\nu _{n,\ell }^{\mbox{model}}$| by determining the corresponding fitting coefficients of the power-law correction to make them as close as possible to the observed frequencies |$\nu _{n,\ell }^{\mbox{obs}}$|⁠.

When fitting individual oscillation frequencies, it is necessary to correctly match each observed mode with its corresponding theoretical counterpart of identical radial order n as the surface correction and thus the computation of the likelihood of the model depend on the difference in frequency between the observed and modelled frequency of the same radial order and angular degree. Identifying the radial order in the observations of main-sequence stars is relatively straight forward using the dimensionless offset ϵ (White et al. 2012), and BASTA has an option for the user to apply a suitable correction to the radial order of the observed frequencies, if desired. An example of the best-fitting model found by BASTA when fitting the individual oscillation frequencies of the Kepler target 16 Cyg A is shown in Fig. 1.

Figure 1.

Échelle diagram of the main-sequence star 16 Cyg A and the subgiant star HR 7322. The coloured circles show the observed oscillation modes (red: ℓ = 0, green ℓ = 1, blue: ℓ = 2) from Lund et al. (2017) and Stokholm et al. (2019), respectively, while the coloured symbols with a black outline show the modes predicted from model with the highest likelihood (same colour coding), corrected using the combined surface correction from Ball & Gizon (2014) and scaled using the observed Δν of the star. The size of the symbols from the model is scaled inversely with their normalized mode inertias: the larger the symbol, the greater the probability of the mode being observed. The lighter coloured symbols with no outline are not matched to any observation, but still predicted by the model.

Open in new tabDownload slide

4.1.3 Frequency combinations

There are strong correlations between combinations including five individual frequencies, and Deheuvels et al. (2016) showed that the set r₀₁₀ results in almost singular covariance matrices with large condition numbers that can lead to overfitting the data as recently suggested by Roxburgh (2018). Instead, the latter study recommends the usage of r₀₁ or r₁₀ in combination with r₀₂ to form the series r₀₁₂ or r₁₀₂, respectively. Here, the underlying problem is in the use of conventional numerical methods to estimate the inverse of covariance matrices, which turns out to be highly inaccurate if these are ill-conditioned. We overcome this issue by using the Moore–Penrose pseudo-inverse (see e.g. Strang 2006) of the covariance matrices in the likelihood. The pseudo-inverse is calculated using singular value decomposition and sets singular values to zero when these are below a threshold that is defined relative to the largest singular value.

The approach devised in BASTAto handle frequency combinations attempts to give the user as much flexibility as possible while keeping the statistical approach robust, and allows fitting any of these quantities (d₀₁, d₁₀, d₀₂, r₀₁, r₁₀, r₀₂, r₀₁₀, r₀₁₂, r₁₀₂) as desired. In all cases the user only needs to provide the individual oscillation frequencies and BASTA will calculate the needed combinations. It is possible to supply the code with the necessary correlations across terms, and if these are not given then BASTA calculates them by doing 10 000 Monte Carlo realizations drawn from random Gaussian distributions of the individual frequencies. An example of the fit to the frequency ratios obtained for 16 Cyg A is shown in the top panel of Fig. 2.

Figure 2.

Ratios and glitch signatures for 16 Cyg A as a function of oscillation frequencies. In the top panel, the black and green circles represent the observed ratios r₀₂ and r₀₁, respectively, while the black and green diamonds show the corresponding quantities for the best-fitting model (see the legend). In the bottom panel, the circles represent the sum of the observed glitch signatures from the helium ionization zone and the base of convection zone, while the diamonds show the same for the best-fitting model.

Open in new tabDownload slide

4.1.4 Acoustic glitches

There are regions inside solar-like oscillating stars where the sound speed changes at length scales substantially shorter than the local wavelengths of the acoustic waves. These regions are known as acoustic glitches in the stellar structure, and the two most prominent are the helium ionization zone and the base of the envelope convection zone. Glitches in the acoustic structure of stars leave tiny signatures in the observed oscillation frequencies currently detectable from space borne missions (see e.g. Miglio et al. 2010; Mazumdar et al. 2012) as well as from ground-based facilities (see e.g. Bedding et al. 2010; Grundahl et al. 2017). The bottom panel of Fig. 2 shows an example of glitch signatures as seen in the oscillation frequencies of 16 Cyg A observed by the Kepler satellite.

In the current implementation, we do not use the parameters associated with the base of convection zone glitch in the stellar properties determination as it is typically difficult to reliably extract them from the contemporary frequency precision for two reasons: (1) the small amplitude of the base of convection zone signature, and (2) the issue caused by aliasing (Mazumdar & Antia 2001). Having said that, it is straightforward to modify the current implementation to use the glitch parameters related to the base of convection zone.

4.1.5 Mixed modes and the frequency matching routine

As stars evolve beyond the core-hydrogen burning phase their oscillation pattern develop irregularities: when the helium core contracts as a product of stellar evolution, the frequency of the g modes increases and interactions between the p-mode behaviour near the surface and the g-mode behaviour near the core take place. Modes can experience mixed properties and exchange nature, causing these so-called mixed modes to deviate from the regular oscillation pattern of the p modes and thus to be visible in an échelle diagram as avoided crossings.

Mixed modes have a substantial diagnostic potential as they are sensitive to properties of the stellar core. However, mixed modes complicate the analysis of individual frequencies as this departure from simple asymptotic theory results in the presence of more than one non-radial mode of a given angular degree between two consecutive ℓ = 0 modes, and some of these theoretically predicted mixed modes do not reach observable amplitudes. This complicates the matching of model modes to observed modes, which can lead to incorrect computations of the surface correction and of the likelihood evaluation of the given model.

To address this issue we note that the amplitude of a mode can be roughly estimated from the mode inertia (see equation 13). The frequency fitting procedure in BASTAuses a mode matching algorithm, where the observed modes are matched to their most-likely counterpart in the model using their frequencies as well as their inertias as a proxy of the likelihood of observability. We describe the matching procedure and assumptions in the following paragraphs.

The frequency matching routine is based on mode counting. Even though the extracted frequencies can be uncertain at times as they depend on the power spectrum background and systematic effects in the pipeline, these effects do not change the relative ordering of the modes. An ℓ = 1 mode and an ℓ = 0 mode will not exchange order due to effects such as the surface effect. Due to the physical nature of the radial modes, avoided crossings do not occur in the pattern of the ℓ = 0 modes. We therefore use the observed radial modes to define a number of frequency intervals. If higher angular degree modes are observed outside the frequency range encompassed by the radial modes, the frequency binning is extended to lower and higher frequencies in steps of Δν.

BASTA counts how many modes of a given angular degree ℓ ≠ 0 are present in each frequency bin between two consecutive radial modes. When there is an equal number of model and observed modes, the modes are matched based only on frequency. We note in passing that this is the typical case for main-sequence stars that do not present avoided crossings. If there are more observed modes than modelled modes in a given range, the model is rejected as it fails to accurately describe the observed pattern.

If, between two ℓ = 0 modes, j modelled modes exist with inertias I₁ < … < I_j and k observed modes of same angular degree with j > k, we need to select k modelled modes to match one-to-one to the observed modes. This is usually the case when avoided crossing takes place in the model, but not all mixed modes have high enough amplitudes to be observed. Intuitively, one might think to just pick the k modelled modes with the lowest inertia as they should have the highest likelihood of being observed. However, small differences in inertia might cause this to result in a miss-matching. Instead, BASTAselects two inertia thresholds a and b (with a < b) and subdivides the modelled modes into a set A with inertias less than a, a set B with inertia between a and b, and a set C for inertias greater than b. The set A thus contains modes that are likely to be detected, while set C contains modes that are unlikely to be observed, and set B contains modes that are somewhere in between.

By picking all modes in A and a subset of B such that k modes are selected in total, the modes can be matched one-to-one in frequency to the observed modes. Specifically, the thresholds a and b are chosen based on the k’th smallest modelled inertia I_k: a is I_k/10 and b is 10I_k. This ensures that |A| ≤ k ≤ |A| + |B|, where |A| is the number of modes in the set. Each possible match is evaluated based on the total L₁ distance in frequency space and the subset of B that minimizes this metric is chosen. This match between the observed and modelled modes will then be used in the following surface effect and model likelihood computation.

Fig. 1 shows the échelle diagrams for two examples of asteroseismic fitting to individual frequencies using this matching algorithm. The code can nicely reproduce the observed oscillation pattern of the main-sequence star 16 Cyg A, and also follow the rapid evolution of a distinct dipole mixed mode in the bright F6 subgiant star HR 7322. Further examples of matching the mixed-mode pattern in subgiant stars are given in Section 6 and Appendix  B below.

4.1.6 Period spacing

The observed power spectra of red-giant stars exhibit a complex pattern due to the presence of mixed modes. As described above, mixed modes behave as gravity modes in the inner regions of the star and as acoustic modes in the outer layers, and hence their observations provide a unique opportunity to probe the conditions deep in the stellar core (see e.g. Beck et al. 2011; Bedding et al. 2011; Mosser et al. 2012b). The detection of the mixed modes makes it possible to measure the gravity mode period spacing, and currently measurements of the asymptotic period spacing for the dipole modes are available for several thousands of Kepler red-giant stars (see e.g. Stello et al. 2013; Mosser et al. 2014; Datta et al. 2015; Vrard, Mosser & Samadi 2016).

4.2 Atmospheric properties

4.2.1 Surface chemical composition

Under these assumptions, all custom-constructed grids used in BASTA are mapped from the input |$\rm [Fe/H]_\mathrm{ini}$| and |$\rm [\alpha /Fe]$| into |$\rm [M/H]$| using equation (28), and then the ratio (Z/X)_model is determined using equation (27). The BaSTI models on the other hand follow the inverse procedure where we use their initial mass fractions of X and Z as given in their data base to determine |$\rm [M/H]$| following equation (27), and |$\rm [Fe/H]$| is then calculated for the corresponding case of |$\rm [\alpha /Fe]$| using equation (28) and the solar mixture of Caffau et al. (2011). This procedure ensures that all surface abundance ratios are computed in a consistent manner, but it is up to the user to ensure that the combination of input values of |$\rm [Fe/H]_\mathrm{ini}$| and |$\rm [\alpha /Fe]$| are reasonable for observed stars. We note in passing that the surface abundance by mass of elements such as ³He, ¹²C, ¹³C, ¹⁴N, and ¹⁶O are also available in our custom-computed grids.

4.2.2 Synthetic photometry, parallaxes, and distances

All grids of stellar models available in BASTA(cf. Section 3) can be mapped from the theoretical H–R diagram to various photometric systems using bolometric corrections (BCs) tables provided by Hidalgo et al. (2018). This allows us to determine synthetic magnitudes in more than 15 photometric systems including those that are of relevance for asteroseismic and exoplanet studies (e.g. Kepler and TESS passbands), compilations of bright stars (Tycho-2 and Hipparcos), large photometric surveys (e.g. 2MASS, Skymapper, sloan, and VISTA), and naturally all Gaia data releases. A list of the photometric systems currently available in BASTAis provided in Table A1.

BASTAsupports the inclusion of parallaxes as a fitting parameter together with at least one apparent magnitude m_ζ. The evaluation is then done by comparing the grid-model absolute magnitude M_ζ, computed from the bolometric luminosity and effective temperature using bolometric corrections, to the measured absolute magnitude computed from the apparent magnitude m_ζ, an estimate of extinction, and the distance modulus from the observed parallax μ = 5log (d) − 5. The computation of the observed absolute magnitude undergoes multiple transformations and thus the assumption of this value being normally distributed like the other fitting parameters is weak. Instead this parameter is evaluated by constructing the likelihood distribution of M_ζ and including it in the calculation of the posterior using equation (2).

The procedure when parallax is included as a fitting parameter is as follows. In the first step, BASTA constructs prior distributions in distance and apparent magnitude using their observed values. The distance distribution is analytically calculated from the measured parallax and its associated uncertainty using the exponentially decreasing space density prior with a scale length of 1.35 kpc described in Bailer-Jones (2015) and Astraatmadja & Bailer-Jones (2016). For the apparent magnitude m_ζ, we assume a normal distribution with the observed values of m_ζ and its uncertainty as the mean and standard deviation, respectively.

In the second step, BASTAsamples over the distance and apparent magnitude distributions to calculate the reddening. When multiple distributions are considered it is difficult to properly sample the tails of all distributions if one just draws samples following the distributions and let the number density of the samples determine the probability density function. If we draw samples from two normal distributions, then the odds of drawing values at e.g. 3 standard deviations away from the mean in both distributions is less than 1 in 100 000. This means that the tails of the resulting distribution would be artificially steep. To overcome this, BASTA draws the samples systematically over a large range of values for each given parameter.

For a number N of samples in distance, BASTAdraws half of them linearly across a range of d ∈ [10^{(− 0.4)}, 10^(4.4)] pc, and the other half as quantiles of the normal distribution around the mode of the distance distribution derived from the observed parallax as described above. Similarly, we produce K apparent magnitude samples distributing half of them linearly across the limiting magnitudes published by each relevant survey (or else assume m_ζ ∈ [−10, 25]), and the rest from a normal distribution around the observed apparent magnitude.

For each of the N × K pairs of distance and apparent magnitude we determine the reddening E(B − V) using the latest version of the Bayestar dust map (currently Bayestar19, see Green et al. 2019), and if the target falls outside of its coverage we apply the simpler map of Schlegel, Finkbeiner & Davis (1998). Since Bayestar provides multiple estimates of the colour excess at each distance, we determine the individual reddening values using the median and standard deviation of those samples. To transform the reddening values into absorption A_ζ estimates at a given filter ζ, we use extinction coefficients from table 6 in Schlafly & Finkbeiner (2011) assuming the Cardelli, Clayton & Mathis (1989) relation A_ζ = 3.1 × E(B − V). For filters not contained in that compilation, we default to the temperature and metallicity dependent extinction coefficients from Casagrande & Vandenberg (2014).

We compute the absolute magnitudes M_ζ from our N × K groups of distance, apparent magnitude, and absorption using the distance modulus. This sample is converted into a probability distribution function by weighting each obtained absolute magnitude by the underlying observed distance and apparent magnitude probability distributions. This distribution is then included in the computation of the posterior (see equation 2). We note that in cases where the parallax uncertainty is smaller than 5 per cent and the absolute magnitude prior is symmetric (i.e. if the distance between the median and each quantile is within a predefined threshold), we fit a Gaussian function to the probability distribution of absolute magnitudes and use the analytical expression in the computation of the likelihood.

Fig. 3 shows an example of the distance sampling procedure. The fitting includes effective temperature, metallicity, and parallax using the 2MASS filters. For the purpose of this example we have modified the apparent J-magnitude to be in disagreement with the other two magnitudes, as it can be seen in the Kiel diagram. The right-hand panel of Fig. 3 shows the probability density function of |$\log \, (g)$| predicted by each of the filters. If BASTA did not sample the tails of the magnitude distributions far from their median values and standard deviations, the full likelihood as defined in equation (2) will be zero. The designed sampling scheme avoids this singularity and provides a robust statistical solution.

Figure 3.

Left: Kiel diagram depicting BaSTI isochrones colour-coded according to their correspondence to each fitted parameter. Right: Resulting |$\log \, (g)$| probability density function as predicted in each filter, as well as the joint distribution obtained from equation (2; dashed line). See the text for details.

Open in new tabDownload slide

In addition to fitting parallax directly, distances can be determined in BASTAindependent of any parallax input as long as one photometric apparent magnitude is provided to the code. In this case, we solve for the distance and absorption iteratively using the magnitudes and the dust map until convergence is reached (normally within three iterations). If multiple apparent magnitudes are given as input we derive a distance in each passband and determine a joint distance by multiplying the individual probability distribution functions.

5 INTERPOLATION

If the resolution of the grid used with BASTAis lower than what is desired, the grid can be interpolated to match a user-defined resolution. BASTA includes this option and can interpolate along or across stellar evolution tracks or isochrones (or both along and across, which we refer to as the combined option), and every parameter of the grid can be included in this interpolation. The individual procedures for each type of interpolation are described in the following sections.

To minimize computation time, the interpolation of tracks or isochrones can be limited to only a section of the original grid given by limits in any of the grid parameters. When fitting multiple stars simultaneously, these limits can be applied on a star-by-star basis (and thus producing one interpolated subgrid per star), or as a single set of limitations for the full collection of stars (producing only one new subgrid of tracks or isochrones).

5.1 Interpolation along tracks or isochrones

For interpolation along the tracks or isochrones the user must define two relevant quantities. The first is the desired resolution in a stellar property between consecutive points in the track, normally an observed quantity that will be fitted (e.g. large frequency separation, or individual oscillation frequencies). The second relevant quantity is a (smoothly varying) grid parameter to be used as the independent variable in the interpolation. We refer to this parameter as the ‘base parameter’, where typical examples are age, central hydrogen content, or central density for stellar evolution tracks, or the initial mass for the case of isochrones. The user can choose which base parameter will be used, or BASTAwill consider age or initial mass as default for stellar tracks or isochrones.

Once these two quantities are defined, the number of points along each interpolated track or isochrone is estimated as the number of points needed to satisfy the desired resolution assuming an equal spacing in the base parameter. The interpolation is performed as a 1-dimensional function using either a linear or cubic method via scipy.interpolate.interp1d (from the scipy package; Virtanen et al. 2020) on a track-by-track basis. The new base along with the interpolated parameters are stored and replace the original track. Therefore, one should consider these as completely new evolutionary tracks, instead of simple refinements of the original.

5.2 Interpolation across tracks or isochrones

Interpolation across tracks or isochrones can be applied to grids of identical microphysics. The base for this interpolation is formed by the parameters used to construct the original grid (see Section 3), and their spacing gives the original grid resolution. For a Cartesian interpolated grid a desired resolution can be set for each of these parameters, and the code determines the minimal amount of tracks needed to satisfy this resolution one parameter at the time. For a Sobol-sampled interpolated grid the user defines a multiplicative increase in the number of tracks of the selected section of the grid. BASTAthen automatically determines a homogeneous distribution of new tracks in the base parameters that meets the required increase in the number of tracks. To retain this homogeneity, the resulting grid consists solely of interpolated tracks and the original ones are excluded. An example of these tracks assignments can be seen in Fig. 4. We note that a Sobol-sampled interpolated grid can be constructed starting from a Cartesian grid, as shown in the right-hand panel of the figure.

Figure 4.

Example of how new tracks are assigned for interpolation across tracks for the case of a base consisting of M_ini and |$\rm [Fe/H]_\mathrm{ini}$|⁠. The simplices show the tessellation of the original grid, where each new track is interpolated to by using the tracks forming the simplex it is contained within. The left-hand panel shows the assignment of tracks for a Cartesian input grid with a starting resolution of |$0.05 \, \mathrm{dex}$| in initial metallicity and |$0.01\, \rm M_{\odot }$| in initial mass, and a desired interpolated resolution of |$0.01\, \mathrm{dex}$| and |$0.005\, \rm M_{\odot }$|⁠, respectively. The middle panel shows the assignment of tracks for a Sobol input grid with an increase in the number of tracks by a factor of 2. The right-hand panel shows the assignment of tracks for a Cartesian input grid but with Sobol assignment of tracks with an increase in resolution by a factor of 2. See the text for details.

Open in new tabDownload slide

The enveloping tracks used to interpolate each new track are determined from a tessellation of the base parameters using scipy.spatial.Delaunay (based on the qhull package; Barber, Dobkin & Huhdanpaa 1996). The number of tracks considered to envelop a new track corresponds to the number of parameters in the base plus one. In the example shown in Fig. 4 this requires three enveloping tracks for each new track (there are two parameters in the base, M_ini and [Fe/H]_ini). The tessellation divides the original grid into triangles, as it can be seen in the example.

In addition to the enveloping tracks, the user must define a single additional smoothly varying quantity that runs along the evolution of the track to perform the interpolation (e.g. age, central hydrogen content, central density). We note that, when applying the combined interpolation method, this quantity does not need to be the same as the one selected for interpolation along the track. For example, the user can choose central hydrogen content when interpolating across tracks, and the subsequent interpolation along tracks can be performed with age as the independent variable.

To avoid extrapolations, the range of the smoothly varying parameter in the new track is limited to an interval that is contained by all enveloping tracks, and its spacing is determined as the mean of the spacing in the enveloping tracks. Using this basis, each parameter in the new track is interpolated separately using scipy.interpolate.LinearNDInterpolater (that relies on a new tessellation of the points along each enveloping track). Individual oscillation frequencies are treated as independent parameters, and are therefore interpolated individually.

As a test of our interpolation procedure, Fig. 5 shows the results of reconstructing an evolutionary track extracted from the Sobol grid presented in the middle panel of Fig. 4. The top panel depicts the evolution of the large frequency separation Δν as a function of age for the three enveloping tracks determined from the tessellation of the base parameters, as well as the track reconstructed with our combined interpolation approach. The bottom panel presents the outcome of the combined interpolation procedure as the fractional difference in Δν between the interpolated track and the original track. Deviations are at the 10⁻⁴ level, an order of magnitude smaller than the uncertainty in Δν measured for the best Keplertargets (see fig. 6 in Yu et al. 2018).

Figure 5.

Example of the performance of our interpolation scheme. Top: large frequency separation as a function of age for the interpolated track (mint green) and the enveloping tracks used for the interpolation. The latter have been colour-coded according to the base quantity used along the track (the central density ρ_cen). Bottom: fractional difference in Δν between the original and interpolated track. See the text for details.

Open in new tabDownload slide

Figure 6.

Posterior distributions of stellar parameters for the Kepler target 16 Cyg A obtained fitting the set (T_eff, |$\rm [Fe/H]$|⁠, Δν, ν_max) with a Sobol grid of different resolutions. Left: original grid. Right: interpolated grid. See the text for details.

Open in new tabDownload slide

The impact of our interpolation procedure can be seen in Fig. 6, where we compare the obtained posterior PDFs for mass, radius, age, and density when fitting data of the Keplertarget 16 Cyg A. For this example we considered as input parameter the effective temperature and metallicity from Ramírez, Meléndez & Asplund (2009) and the large frequency separation and frequency of maximum power derived by Lund et al. (2017). These data were then fitted to a Sobol grid in its original resolution, and an interpolated grid constructed with the combined method and an increase in resolution by a factor of 5 in the number of tracks (across tracks) and spacing of 0.1μHz between the lowest observed ℓ = 0 frequency (along tracks). The derived quantities are in very good agreement (and certainly within their respective uncertainties), but the resulting distributions are significantly smoother after the interpolation procedure.

Before closing this section we note that when the combined method is chosen, the interpolation across tracks or isochrones is performed before interpolation along the tracks or isochrones to increase computation efficiency. We have tested the inverse case (along before across), and confirmed the differences in the recovery procedure presented in Fig. 5 are of the same magnitude regardless of the order of the interpolation. This inverse case (along before across) is still available for usage in BASTA at a much larger computational cost.

6 VALIDATION WITH ARTIFICIAL DATA

We performed a thorough end-to-end validation of BASTA to quantify the robustness of our pipeline in retrieving stellar properties. For this purpose we produced artificial data from models extracted from a grid and used the same grid to fit the observables, which is equivalent to assuming that the underlying stellar models are true representations of the observations. By quantifying the deviations in our derived parameters from the true values we can estimate the level of accuracy of our pipeline, as these will depend exclusively in the reliability of BASTA. In addition, our obtained uncertainties in stellar properties provide an estimate of the typical statistical precision, which in turn depends on the assumed observational errors and the combination of input quantities fitted. We emphasize that this exercise allows us to test the accuracy and statistical precision obtained for a given set of observables and their measured uncertainty, but we cannot account for deviations between the physical and observed properties of a real star and those predicted by a grid of stellar models. This additional systematic uncertainty undoubtedly exists and remains to be quantified, but it depends on our many shortcomings in theory of stellar evolution and goes well beyond the scope of this paper.

For this particular exercise we constructed a Sobol grid of stellar models using GARSTEC comprising 3000 evolutionary tracks, a mass range between 0.8 and 1.5 M_⊙, initial metallicity −0.5 <[Fe/H]_ini < +0.5, the Asplund et al. (2009) solar mixture, and an enrichment law ΔY/ΔZ = 1.4. The tracks were evolved from the pre-main sequence to a value of the large frequency separation Δν = 10 μHz, roughly corresponding to the lower RGB. We determined individual oscillation frequencies of angular degree ℓ = 0, 1, 2, 3 for stars down to a value of Δν = 30 μHz, and only radial modes for more evolved stars. The coverage of the grid was selected to ensure that we can test our fitting procedures in various types of stars (main sequence, subgiants, and red giants) which have different observed quantities to be fitted (individual frequencies, mixed-modes, period spacing).

We assume the following observational uncertainties (taken from Serenelli et al. 2017 and Lund et al. 2017 for main-sequence and subgiant targets, and Yu et al. 2018 for red giant stars): 70 K in T_eff, 0.1 dex in |$\rm [Fe/H]$|⁠, 0.1 dex in |$\log \, (g)$|⁠, 0.5 per cent of the observed value in Δν, and 2 per cent in ν_max. These are typical uncertainties in the global asteroseismic properties for Kepler stars observed for more than 50 d (see Serenelli et al. 2017). The asymptotic period spacing uncertainty is assumed to be 1.5 per cent, which corresponds to the average uncertainty measured by Vrard et al. (2016) for RGB stars in the range 10 μHz ≤Δν ≤ 30 μHz. For the individual frequencies of oscillation we adopt a two-step process where we first identify which modes would be detected given an assumed observation length, and then determine their uncertainties following a recipe derived from Kepler targets.

For a given theoretical mode calculation we estimate which modes would likely be measurable based on results from the main-sequence stars included in the LEGACY study (Lund et al. 2017). From the LEGACY data we first estimate the typical minimum and maximum frequencies of measured ℓ = 0 modes in units of ν_max, and derive simple linear relations of these frequencies as a function of ν_max. The typical frequency intervals of measurable ℓ = 0 modes are found to range from |$\pm 0.4\, \nu _\mathrm{max}$| at |$\nu _\mathrm{max}=\rm 1000\, \mu Hz$| to approximately |$\pm 0.2\, \nu _\mathrm{max}$| at |$\nu _\mathrm{max}=\rm 4000\, \mu Hz$|⁠. We then estimate the expected relative amplitude of ℓ = 0 modes at these frequency limits, using a relation for the envelope width of the assumed Gaussian modulation of mode amplitudes around ν_max. This relation is determined from fits to LEGACY data, and includes both a dependence on ν_max and T_eff (Lund et al., in preparation). By including mode visibilities from Lund et al. (2017), we can then estimate the corresponding relative amplitudes of non-radial modes and assess which of these exceeds the limit for detectability set by the ℓ = 0 modes. Following Ball et al. (2018), we include information on the mode inertia for non-radial modes, where we divide amplitudes by the square-root of the Q-factor (ratio of the mode inertia relative to the ℓ = 0 inertia at the corresponding frequency, see also equation 10).

Concerning frequency uncertainties for the modes deemed measurable from the above procedure we use a polynomial relation between the frequency uncertainties in units of ν_max of ℓ = 0 modes from the LEGACY data and the corresponding relative frequency away from ν_max. At ν_max the typical minimal uncertainty is found to be of the order |${\sim }7.6\times 10^{-5}\, \nu _\mathrm{max}$|⁠, which for a |$\nu _\mathrm{max}=\rm 2000\, \mu Hz$| star corresponds to |$\rm \sigma _{\nu }\sim 0.15\, \mu Hz$|⁠. Away from ν_max the typical uncertainties increase by factors of 5 (⁠|$\nu \sim -0.4\, \nu _\mathrm{max}$|⁠) to 10 (⁠|$\nu \sim +0.4\, \nu _\mathrm{max}$|⁠).

With the observational uncertainties in all relevant observational quantities defined, the validation procedure was designed as follows. We selected 443 models from the grid to be used as artificial targets mimicking the observed distribution in large frequency separation of the Keplermain-sequence and subgiant sample of Lund et al. (2017; see Fig. 7). This sample reaches values of Δν ≃ 20μHz, which we extended in same proportion as the last bin to Δν ≃ 10μHz to encompass the base of the RGB. We considered various sets of input quantities (see Table 2) and generated the synthetic data by drawing a sample from a normal Gaussian distribution with the model value of the relevant parameter as the mean and a standard deviation given by the observational uncertainty described above. In that manner, each quantity in the input set does not exactly correspond to the model value but it has been perturbed by typical observational uncertainties. This is the closest we can simulate observations of a real star, under the assumption that the underlying grid of models is a true representation of the physics at play in stellar evolution.

Figure 7.

Left: distribution in large frequency separation of the synthetic data used for the validation procedure. Right: Kiel diagram depicting the distribution of our synthetic targets across the grid, colour-coded by the value of the large frequency separation.

Open in new tabDownload slide

The distribution of the recovered fractional values in density, radius, mass, and age for the main-sequence artificial targets are shown in Fig. 8 (equivalent figures for the subgiant and RGB targets are given in Appendix  B). As we can see, the inferred properties are generally in good agreement with the corresponding underlying input parameters. The inclusion of asteroseismic data dramatically improves the quality of the recovered properties compared to the cases where only atmospheric parameters are used.

Figure 8.

Distribution of normalized differences between derived stellar properties with BASTAand the underlying stellar model for the main-sequence targets. Values are computed as (Solution-Model) or Model. See the text for details.

Open in new tabDownload slide

As is also visible in Fig. 8, the inferred properties do not match exactly with the underlying model values. To ensure that this is a result of finite precision of the data, we reduced the assumed standard uncertainties by a factor of 10⁻⁴ and repeated the above validation. BASTA recovered the original model in 440 out of 443 cases, and in 3 cases we pick the model just before or after which is expected due to the perturbation to the model parameters that defines the input observables. These results confirm that the differences are due to finite data precision and are not an artefact of combined effects of priors and weights.

Since we have determined stellar properties assuming typical observed uncertainties in various sets of input, we can use the BASTA results to quantify the statistical uncertainty obtained across evolutionary phases for various sets of input. Table 2 summarizes the results, where we have listed for each stellar property the average of the 16th and 84th percentiles normalized by the median derived by BASTA. It is clear from the table how the accuracy increases as asteroseismic quantities are included in comparison to just the spectroscopic input, as well as the additional gain from including individual frequencies. This exercise also shows that, if our stellar models are a good representation of the observed stars, statistical uncertainties below the 10 per cent level in age are within reach thanks to the advent of asteroseismology.

Another interesting point arising from the validation is the extremely high precision in the recovered stellar properties for subgiant stars when individual oscillation frequencies are fitted. As mixed-modes evolve rapidly and are very sensitive to the conditions in the stellar core, such a small statistical uncertainty is not surprising as long as the code is capable of identifying the correct underlying model. We remind the reader that in subgiant stars the number of model non-radial oscillation frequencies between two radial modes exceeds the number of observed frequencies as their visibility depends in the inertia. The procedure we described in Section 4.1.5 can successfully identify which mixed-modes should be visible and thus avoids an artificial increase in the statistical uncertainties due to missidentified modes. We include a few representative Échelle diagrams from the validation procedure in Appendix  B, where the accuracy of the algorithm is demonstrated.

7 NOTES ABOUT CODE DEVELOPMENT AND CONTRIBUTING

The core developers of BASTA have customized the code following the needs of our own research fields. As the code begins to be used by scientists around the world, we expect that the inclusion of additional features will become highly desirable. Naturally the amount of proposed improvements to the code will increase proportionally to the number of users, and it will go beyond the available time of the core developers to implement them and maintain a stable version of BASTA.

For these reasons we rely on the git version-control system and the GitHub repository⁶ to encourage users to develop and share their contributions with the rest of the users. New features, additions, or extensions should be in separate branches or forks so that they do not need the direct involvement of the core developing team. Once the new additions are completed, the code developers will be happy to handle the pull request and make it part of the stable version of BASTA after proper testing has been completed. A more detailed contribution guide can be found on the GitHub repository along with the installation instructions and a tracker for issues/bugs/suggestions. We encourage the users to follow the project on GitHub (please ‘watch’ and select ‘custom or releases’ as a minimum) to get notifications and updates on new developments. We welcome any participation including pull requests.

8 CONCLUDING REMARKS

We have introduced the public version of the BASTA, an open-source code that determines stellar properties using a set of observables and a grid of stellar models or isochrones. It is flexible in its input and can combine a large number of spectroscopic, photometric, astrometric, and asteroseismic input to extract properties of stars under a statistically robust Bayesian scheme. The large number of functionalities included in BASTA, combined with the various sets of publicly available or custom-computed grids of stellar models or isochrones, make it the most versatile pipeline for stellar analysis currently available.

We have thoroughly described the type of input that can be given to retrieve stellar properties, and discussed the assumptions made when predicting these quantities from the underlying grids of evolutionary models. We have performed an extensive validation test that confirms the reliability of BASTAin determining accurate stellar properties, and use these results to quantify the typical statistical uncertainties obtained for various combinations of fitting parameters. Our results show that asteroseismic ages with statistical uncertainties below the 10 per cent level are achievable for data sets of the quality obtained by the Kepler satellite, and the only limiting factor is the reliability of our stellar evolution and pulsation models.

We hope to have provided the community with a useful analysis tool for stellar properties, which is specifically designed to meet the challenges of large inhomogeneous data sets available in the era of large-scale stellar surveys. Its capability of combining a wide range of input data makes it an invaluable tool to fully exploit the potential of current stellar catalogues (by simultaneously fitting e.g. TESS, APOGEE, 2MASS, and Gaia eDR3 data), and its future developments will make it ready for the challenge of the next data deluge from surveys such as PLATO 2.0 (Rauer et al. 2014), WEAVE (Dalton et al. 2014), 4MOST (de Jong et al. 2019), and the Legacy Survey of Space and Time at the Vera Rubin Observatory (LSST; Ivezić et al. 2019).

ACKNOWLEDGEMENTS

Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant agreement No. DNRF106). VAB-K acknowledges support from VILLUM FONDEN (research grant 10118), the Independent Research Fund Denmark (Research grant 7027-00096B), and the Carlsberg foundation (grant agreement CF19-0649). SC acknowledges support from Premiale INAF MITiC (PI: B. Garilli), from Istituto Nazionale di Fisica Nucleare (INFN) (Iniziativa specifica TAsP), from Progetto Mainstream INAF (PI: S. Cassisi) and grant AYA2013- 42781P from the Ministry of Economy and Competitiveness of Spain. AMS is partially supported by grant PID2019-108709GB-I00 of the MICINN of Spain. AS acknowledge support from the European Research Council Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, G.A. n. 772293, http://www.asterochronometry.eu).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

The code is available via https://github.com/BASTAcode/BASTA.

The following softwares were used during this research: astropy (Astropy Collaboration 2013), corner (Foreman-Mackey 2016), dustmaps (Green 2018), emcee (Foreman-Mackey et al. 2013), h5py (Collette 2013), healpy (Górski et al. 2005), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scikit-learn (Pedregosa et al. 2011), scipy (Virtanen et al. 2020), and tqdm (da Costa-Luis et al. 2021).

Footnotes

REFERENCES