EFFECT OF UNCERTAINTIES IN STELLAR MODEL PARAMETERS ON ESTIMATED MASSES AND RADII OF SINGLE STARS

It is extremely difficult to estimate masses and radii of stars, in particular masses of single stars. However, in this era of exoplanet studies the masses and radii of exoplanet hosts play an important role in establishing the properties of the planets.

The usual way to determine stellar properties is through spectroscopic or photometric analyses of the star combined with fitting to grids of stellar models (see, e.g., Takeda et al. 2007). The recent, dramatic advances in observational asteroseismology provide us with a relatively new way to determine masses and radii of solar-type stars. These advances have come in large part from new satellite observations, for example from the French-led CoRoT satellite (e.g., Baglin et al. 2002; Michel et al. 2008; Appourchaux et al. 2008), and in particular the NASA Kepler Mission (Chaplin et al. 2010; Gilliland et al. 2010). Kepler has provided detections of solar-like oscillations for hundreds of solar-type stars (Chaplin et al. 2011) when only around 20 such stars had detections prior to the Mission.

The Kepler Mission (Borucki et al. 2010) is currently the largest source of data on extra-solar planet candidates and their host stars. Properties of stars in the Kepler field of view are available in the Kepler Input Catalog (KIC; Brown et al. 2011). The catalog provides estimates of T_eff, log g, [Fe/H], and E(B − V) obtained from the analysis of multi-color photometry and also derived properties such as radius, mass, and luminosity. Verner et al. (2011) compared the seismically derived radii of a sample of solar-type stars with those in the KIC and found that although there was general agreement, there was an underestimation bias of up to 50% for stars with R_KIC < 2R_☉, showing that estimates of stellar properties are prone to uncertainties.

Stellar masses and radii obtained by modeling observed stellar parameters such as the effective temperature T_eff, metallicity [Fe/H], and surface gravity log g are subject to systematic errors caused by uncertainties in the inputs to stellar models, as well as errors that arise because of uncertainties in the measured stellar parameters. In this paper, we determine the effect of uncertainties in model parameters on derived stellar masses, radii, and average stellar densities. The paper is organized as follows: in Section 2, we discuss some possible sources of error; in Section 3, we describe how we estimate stellar masses and radii; in Section 4, we describe the proxy stars that we use to determine the effects of model uncertainties on estimated stellar parameters. Our results are presented and discussed in Section 5 and we present the main conclusions in Section 6.

Stellar properties estimated from observations of temperature, metallicity, and gravity depend on stellar models. As a result, all uncertainties in stellar models will affect the estimated mass and radius of the star under consideration. We consider three main sources of uncertainty: (1) the mixing length parameter, (2) the metallicity scale, and (3) the atmospheric model. Additionally, we examine what happens if there are systematic errors in the temperature scale.

One of the largest sources of uncertainty in stellar models is the treatment of convection. Given that it is not possible to follow both convective and evolutionary timescales in a stellar model, convection in stars is usually modeled using the so-called mixing length theory (MLT) or one of its variants. Unfortunately, there is a free parameter, the so-called mixing length parameter α, that must be specified before one can begin modeling. The usual practice is to assume that the solar mixing length applies to all stars. The solar mixing length (and the initial solar helium abundance) is determined by forcing a 1 M_☉ model to have the solar radius and solar luminosity at the current solar age. While this is the most commonly used approach, efforts to model stellar pulsation data show that stars do not necessarily have the solar value of α (see, e.g., Metcalfe et al. 2010; Benomar et al. 2010; Deheuvels & Michel 2011). Thus, assuming that the stars being analyzed have the solar value of α will necessarily affect our results, particularly since α controls the radii of stellar models.

Another potential source of error is the metallicity scale. Most stellar metallicity estimates are given with respect to solar metallicity. Converting to Z hence requires a value of the solar metallicity. Asplund et al. (2005; henceforth AGS) had suggested that Z/X for the Sun is 0.0165; this value is much lower than the usually accepted value of 0.023 (Grevesse & Sauval 1998) or 0.0245 (Grevesse & Noels 1993). Although Asplund et al. (2009) and Grevesse et al. (2010) revised the AGS results upward to Z/X = 0.0181, this is still lower than the Grevesse & Sauval (1998) value. This uncertainty in the solar metallicity scale adds uncertainties in evolutionary tracks and hence needs examining.

The atmospheric T–τ relations used in stellar models also play a role in defining the fundamental properties of models. It is common to use the Eddington T–τ relation to define the outer boundary of the star, but often other, semi-analytic T–τ relations, such as those of Vernazza et al. (1981) or Krishna Swamy (1966), are used. Since the T–τ relation effectively determines the layer at which T = T_eff, it plays a role in determining the radius of a model of a given mass and temperature, and can add to uncertainties to the estimated stellar parameters.

While the factors mentioned above increase uncertainties in stellar models, one observational issue that can affect stellar parameter determination is T_eff. Determinations of T_eff of the same star by different groups often result in vastly different temperature estimates (see, e.g., Metcalfe et al. 2010 for the case of one of the Kepler stars). There are often scale errors between one set of photometric temperature differences and another with the differences often being as large as 100 K (see, e.g., Casagrande et al. 2010; Pinsonneault et al. 2012). This will affect estimates of stellar parameters, and hence needs to be accounted for in the error budget.

Basu et al. (2010) described a way (the Yale–Birmingham or YB code) in which stellar radii can be determined using seismic parameters. The YB pipeline was modified by Gai et al. (2011) to obtain other stellar properties, such as mass and log g. Here, we use the same method for determining masses and radii of stars even in the absence of seismic data.

The YB method is based on finding the maximum likelihood of the set of input parameter data calculated with respect to the grid models. The estimate of the stellar property is obtained by taking an average of the properties on the grid that have the highest likelihood. We average all values with likelihoods over 95% of the maximum value of the respective likelihood function. For a given observational (central) input parameter set, the first key step in the method is to generate 10,000 input parameter sets by adding different random realizations of Gaussian noise to the actual (central) observational input parameter set. The distribution of any property, say radius, is then obtained from the central parameter set and the 10,000 perturbed parameter sets, which forms the distribution function. The final estimate of the property is the median of this distribution. We use 1σ limits from the median as a measure of the uncertainties.

The likelihood function is formally defined as

$$\begin{equation} \mathcal {L}=\prod ^{n}_{i=1}\left(\frac{1}{\sqrt{2\pi }\sigma _{i}}\times \exp (-\chi ^{2}/2)\right), \end{equation} \tag{ 1 }$$

where

$$\begin{equation} \chi ^{2}=\sum ^{n}_{i=1}\left(\frac{q^{{\rm obs}}_{i}-q^{{\rm model}}_{i}}{\sigma ^{i}}\right)^{2}, \end{equation} \tag{ 2 }$$

with q ≡ T_eff, [Fe/H], and log g in the non-seismic case, and T_eff, [Fe/H], the average large frequency separation, Δν, and the frequency of maximum power, ν_max, in the seismic case.

The two seismic observables, Δν and ν_max,  depend on global properties of a star. The large separation scales as the mean density of a star (e.g., Ulrich 1986; Christensen-Dalsgaard 1993), so that

$$\begin{equation} \frac{\Delta \nu }{\Delta \nu _{\odot }}\simeq \sqrt{\frac{M/M_{\odot }}{(R/R_{\odot })^{3}}}. \end{equation} \tag{ 3 }$$

Thus, assuming we know Δν for the Sun, we can find the mean density of other stars. Stello et al. (2009) have shown that this scaling holds over most of the HR diagram and errors are probably below 1%. The frequency of maximum power in the oscillations power spectrum, ν_max, is related to the acoustic cutoff frequency of a star (see, e.g., Kjeldsen & Bedding 1995; Bedding & Kjeldsen 2003; Chaplin et al. 2008), which in turn scales as M R⁻² T^−1/2_eff. Thus, if we know the solar values of ν_max, we can calculate ν_max for any star as

$$\begin{equation} {{\nu _{\rm max}}\over {\nu _{\rm max,\odot }}}\simeq {{M/M_{\odot }}\over {(R/R_{\odot })^2\sqrt{(T_{\rm eff}/T_{{\rm eff},\odot })}}}. \end{equation} \tag{ 4 }$$

If Δν, ν_max, and T_eff are known, Equations (3) and (4) represent two equations in two unknowns, M and R, which hence can be solved to obtain the mass and radius. This also allows us to calculate log g.

Although Equations (3) and (4) can be used to determine mass and radius directly, the errors in the derived quantities can be unphysically large, in the sense that the equations assume that all values of T_eff are possible for a star of a given mass and radius. However, the equations of stellar structure and evolution tell us otherwise—we know that for a given mass and radius, only a narrow range of temperatures are allowed. This is what led to the development of so-called grid-based methods, such as the YB pipeline, where the seismic observations are used in conjunction with a grid of models. This does however come at a cost, since it introduces some model dependence in the seismic estimates of the stellar properties.

We use four grids of models to determine stellar parameters. These "calibration" grids are (1) a grid of models that constitute the YY isochrones (Demarque et al. 2004); (2) a grid of models constructed with the Dartmouth Stellar Evolution code (DSEP; Dotter et al. 2008), as described by Dotter et al. (2008; these models can be downloaded from the DSEP Web site⁴); (3) the model grid of Marigo et al. (2008) constructed with the Padova stellar evolution code (Marigo et al. 2008; Girardi et al. 2000; these models can be downloaded from the Padova CMD Web site⁵); and (4) a grid of models constructed with the Yale Rotation and Evolution Code (YREC; Demarque et al. 2008) in its non-rotating configuration (these models are described in Gai et al. 2011). All models have been constructed with the solar-calibrated value of the mixing length for the particular code and physics used. Although the input physics in these grids is different, they are not sufficiently different from each other for us to use models from these grids as test models. Consequently we refer to models in these grids as "Normal models."

As can be seen from Equation (1), knowing the uncertainty in the observations is key. We use three sets of errors in this work. Set 1 corresponds to the uncertainties believed to be present in KIC for T_eff, [Fe/H], and log g. Set 2 corresponds to typical errors expected in (good) photometric estimates of T_eff, [Fe/H], and log g. Set 3 is what is expected from spectroscopic determinations. The uncertainties used are listed in Table 1. As can be seen in the table, for all three sets we assume identical errors in the seismic properties.

Gai et al. (2011) have already performed an in-depth analysis of asteroseismic radius and mass determinations. However, they did not consider all the sources of error that we shall consider here; consequently we also examine how the above sources of error affect seismic estimates of the stellar properties. For the seismic determinations we will use the same inputs as in Gai et al. (2011), i.e., T_eff, [Fe/H], the large separation Δν, and the frequency of maximum power, ν_max.

Basu et al. (2010) and Gai et al. (2011) include extensive discussions on how well the method works when seismic data are available. We find that the method works even for the combination of non-seismic data {T_eff, [Fe/H], log g} if we use synthetic stars drawn from the same family of models as the grid. In the absence of uncertainties in any of the parameters, the half-width at half-maximum (HWHM) for the distribution of differences between input and output masses and radii is less than 0.25%, while the semi-quartile distance (i.e., half the distance between the 1st and 3rd quartile points) is about 2%. If we use exact, error-free, data but assume that there could be errors as in Error Set 3, then in the case of mass the HWHM is still less than 0.25%, but the semi-quartile distance is somewhat larger, at 2.8%. For radius, the HWHM is 1%, while the semi-quartile distance is 7%. The probability distribution functions are shown in Figure 1. The results show that we can use the same technique for determining radii and masses of stars regardless of whether or not we use seismic data.

Zoom In Zoom Out Reset image size

In order to investigate the effects of the input parameters we have used several sets of simulated stars drawn from different grids of stellar models to determine the errors made in determining the masses, radii, and average densities of these stars.

In the first set of tests—which we call the "Identical" set—the proxy stars come from the same grid of models that is used to estimate their properties, an important implication being that both have the same physics. In the second set—the so-called Normal set—a mixture of proxy stars is drawn from all four grids, with the name reflecting the fact that these grids were constructed with usual stellar parameters. However, there are differences in detail between the different grids (i.e., differences in the physics inputs used).

We also use four other sets of proxy stars derived from grids constructed especially for this work. The first of these, "AGS," assumes that [Fe/H] = 0 corresponds to Z/X = 0.0165, the solar metallicity claimed by Asplund et al. (2005). This set should simulate the effects of getting the metallicity scale wrong (since the grids used to estimate the stellar properties assume the higher solar Z/X = 0.023 or similar). The AGS proxy stars have solar-calibrated values of the mixing length parameter. The next set, "MLT," comprises models constructed with non-solar values of the mixing length; in particular we use models with α = 1.5, 1.6, 1.7, 1.9, 2.0. 2.1, 2.2, 2.3, and 2.4. The solar α, for the same physics, is 1.826. A third set of so-called KS models was constructed with the Krishna Swamy (1966) T–τ relation in the atmosphere. The last set, which we call "T-shift," was derived from the four grids, just like the "Normal" models; however, the effective temperatures were increased by 100 K to simulate the effect of getting the temperature scale wrong. The different test cases are summarized in Table 2.

In all cases, we added to the seismic and non-seismic input parameters of the proxy stars' random errors consistent with the uncertainties listed in Table 1, and then used these perturbed data as the inputs to our grid pipeline.

Our aim is to see how well we recover the properties of the stars in each of the test cases. To this end, we look at the fractional deviations between the actual properties and the estimated properties of the proxy stars. If there are no systematic effects, the distribution of the fractional deviations will be symmetric and centered at zero. Our results for all test cases analyzed, and for all three sets of errors, are presented in Table 3. In each case we list the median of the distribution as well as the spread (standard deviation of the distribution).

These two values do not completely characterize the distributions and many of them are highly skewed, as we shall see below, and hence we also list the HWHM in both the negative error and positive error directions. The spread as well as the two HWHMs give information on the non-Gaussianity as well as the asymmetry of the distributions. The data in the table allow us to draw some important conclusions: (1) the errors in the results, as expected, decrease when the uncertainties in the simulated input data are decreased; (2) for the same set of errors, adding seismic inputs improves the results considerably; and (3) differences in physics (i.e., between the proxy stars and the model grids used to estimate the properties) lead to noticeable biases in the stellar properties estimated from non-seismic inputs alone, while the effects are much reduced when seismic input data are available.

The effect of errors is shown in Figure 2, where we plot the distribution of the errors for the Normal set results (all three error cases). In the non-seismic case, the best-determined property is the mass. This is not surprising since stellar model properties depend predominantly on (or are fixed predominantly by) mass. As can be seen from the figure, decreasing the input parameter uncertainties makes a dramatic difference in the non-seismic results: unless we have very precise estimates of T_eff, [Fe/H], and log g, mass errors can be as much as 50% in some cases. Radius determinations are much more uncertain, with density faring the worst (since the mass and radius errors propagate through). In our best-case scenario, i.e., Error Set 3, the HWHMs for mass, radius, and density are 8%, 12%, and 35%, respectively.

Zoom In Zoom Out Reset image size

In contrast, the best-determined quantity in the seismic case is density, which is not surprising since the input Δν scales as $\sqrt{\rho }$. The errors in the estimated densities are almost completely independent of the errors in the other input parameters, the error distributions being consistent with a Gaussian with σ = 2% (i.e., entirely consistent with the 1% error adopted for Δν). The next-best determined property is radius. Here, errors in T_eff and metallicity play only a minor role, as reported by Gai et al. (2011). In the seismic case, mass determinations fare the worst and reducing uncertainties in T_eff and [Fe/H] is essential. However, the mass results are still better than the non-seismic case with HWHM of 8%, 7%, and 5% for Error Sets 1, 2, and 3, respectively.

Figure 3 shows the error distribution when the other sets of proxy stars are used. We only show the distributions for Error Set 3 and have also plotted the "Normal" set results as a reference. What is clear from the figure is that there is a systematic shift in the results when the properties of the other proxy stars are deduced using normal stellar grids. The median values of the distributions listed in Table 3 attest to the shifts. We have not shown the results when the proxy stars are identical to the grid models; the results for identical stars are however listed in Table 3, and the spread of the results is actually a measure of the precision of the results. The spread of the other cases is, in some sense, a convolution of the precision and the accuracy of the method. While the results are equally precise in all cases, they are not as accurate, which gives rise to a shift of the peak of the distributions.

Zoom In Zoom Out Reset image size

The effect of uncertainties in the physics is most noticeable in the mass estimates, for both non-seismic and seismic cases. The peak of the distribution shifts away from zero and in many cases the distribution becomes wider. The effects on radius are smaller in a relative sense, and the effects on density are smaller again, relative to what we had seen in the Normal case. Given that mass is believed to be the most fundamental property of a star, and one that controls the evolution, this is somewhat disappointing although not very difficult to explain. In the absence of seismic data, the only observable that depends explicitly on mass is log g. Although the observed T_eff is determined by mass and metallicity, its explicit dependence is on radius. In the seismic case too, Δν and ν_max parameters depend more on radius than on mass. To determine masses of single stars better, we need to model the full oscillating spectrum of a star, not merely Δν and ν_max (see, e.g., Metcalfe et al. 2010; Deheuvels & Michel 2011).

Of all cases, an artificial shift in the T_eff scale produces the smallest effect. While the peak of the distribution shifts slightly, the width does not change for the non-seismic estimates, while the HWHM increases modestly, to only 6% (from 5%), for the seismic estimates. The relative insensitivity to the temperature scale is reassuring since it is not completely clear how the temperature scale can be improved in the short term. Of course, the result would be different if we shifted the temperatures by a much larger amount.

Differences in other physics inputs cause large changes, especially in the absence of seismic data. We find that masses of the low-Z proxy stars are generally overestimated by our calibration grids. One simple way to understand this is the following: masses of stars for which [Fe/H] = 0 is pegged to the low value of solar Z/X are overestimated with grids of Normal models. A low (Z/X) for a given helium abundance results in a star that is hotter than a higher Z/X star. For a given log g, these models have a higher T_eff, mimicking a higher-mass star in the grid of models that has an [Fe/H] scale pegged to the higher Z/X. Since the latest sets of abundances published by Asplund et al. (2009) and Grevesse et al. (2010) are higher than those of AGS (although still lower than the GS98 value), the errors we get for the "AGS" set of proxy stars could be considered to be an upper limit to the errors caused by the metallicity scale.

Unsurprisingly, proxy stars with non-solar mixing length also show systematic errors. Basu et al. (2010) had shown that even when seismic data are available, uncertainties in the mixing length can result in systematic errors in radius measurements. We find that in the absence of seismic data, the errors are much larger, particularly in the case of mass estimates. For a given mass, a change in α leads to a change in radius, with sub-solar α models showing a larger radius (and hence lower T_eff and log g) at a given luminosity than solar α models with the same metallicity. Super-solar α models show a smaller radius (and hence larger T_eff and log g). In most parts of the HR diagram this change in log g (at a given T_eff) seems to be interpreted more like a change in mass than the change in radius that it actually is, giving rise to errors in mass estimates that can be quite large. There are errors in the radius estimation too, but the resulting errors, relative to the errors in the case of Normal models, are not large.

While estimates of mass and radius that are made when seismic data are available also suffer from errors due to uncertainties in inputs to the stellar models, the effects are much smaller. The main reason for this is that the seismic input parameters depend directly on both radius and (albeit to a lesser extent) on mass.

The results for non-solar mixing length stars can be improved by using a grid that consists of several complete grids, each with a different mixing length. To test this we constructed YREC grids with five values of α, α = 1.4, 1.6, 1.826 (solar), 2.0, and 2.2, and re-derived the properties of the MLT stars using this grid. The results are indeed improved compared with the single α YREC grid, as can be seen from Figure 4, where we show results for the Error Set 3 case. As can be seen from the figure, all distributions are peaks at zero difference. Note that the single-α results seem worse that those in Figure 3; this is because the results in Figure 3 are a result of using four different grids, each with a slightly different value of α (α = 1.938 for Dotter et al., α = 1.68 for Marigo et al., α = 1.74 for YY, and α = 1.826 for YREC). Figure 4 demonstrates the need for making grids with a large range of values of the mixing length parameter, α, to obtain accurate properties of stars.

Zoom In Zoom Out Reset image size

The results for "KS" proxy stars, made with the Krishna Swamy atmospheres, can be understood to some extent in terms of the mixing length. A solar-calibrated model with a Krishna Swamy atmosphere needs a mixing length parameter of α = 2.1. The proxy stars we used were constructed with α = 1.826, which is the value of α we need to construct a solar model when the Eddington T–τ relation is used. Thus, the "KS" models in effect have sub-solar mixing length. This means that for the same mass and metallicity, the "KS" proxy stars have a higher log g at a given temperature. As in the case of the "MLT" proxy stars, this appears to be interpreted as a change in mass. The "KS" results also improve if we use a multi-α grid, as demonstrated in Figure 5, confirming that the effect of having a different model of stellar atmospheres is essentially the same as changing α.

Zoom In Zoom Out Reset image size

One might question whether it would be better in the seismic case to use Equations (3) and (4) directly to estimate the stellar properties, thereby removing the sensitivity of the results to the stellar models (and the uncertainties which arise from choices made concerning the input physics). While this will lead to more accurate results, such a strategy comes at the cost of giving much larger uncertainties (i.e., it reduces the precision in the estimated properties), for the reasons mentioned in Section 3 and shown in Gai et al. (2011). Our exercise does however indicate that unless we use a grid constructed with a wide range of values of the mixing length parameter, the formal uncertainties in the results only give a lower bound to the error in the result. Even with a multi-α grid, we will need to make an allowance because of uncertainty in the metallicity scale.

If we consider the complete ensemble of all proxy stars (but do not include the "T-shift" stars), the HWHMs for mass, radius, and density for the non-seismic estimates (Error Set 3) are 12.5%, 17.8%, and 36%, respectively. For the analagous seismic cases, they are 7.6%, 3.2%, and 2.3%, respectively. Again, these errors can be decreased if we use grids made to cover a large range in α.

Estimates of density are robust in all cases (in the sense that the deviation from the true density remains unchanged). Radius estimates change only slightly in the seismic case (HWHM rising to 3.2% from 2.5%), but rise somewhat more in the non-seismic case. It should be noted that the errors we find are in some sense a lower limit to the possible errors, particularly since modeling convection in stars remains uncertain. Additionally, we would expect all the sources of error to be present together, making the situation more complicated.

Uncertainties in inputs to stellar models add uncertainties in properties of single field stars estimated from observational input parameters T_eff, log g, and [Fe/H]. In the absence of seismic data, extremely precise estimates of these parameters are needed to make reasonable estimates of the masses and radii of the stars. With T_eff errors of only 50 K and errors of 0.1 dex in [Fe/H] and log g, we cannot determine masses to better than 8% and radii to better than 14%, even if we assume that we know the physics of stars perfectly. Note that we are talking about the accuracy of the results; the precision is often better. Uncertainties in stellar model parameters—in particular uncertainties in the mixing length parameter—can increase those errors by more than a factor of two. Uncertainty surrounding the solar metallicity (which establishes the [Fe/H] scale) also increases the errors. Small changes in the effective temperature scale do not affect the results substantially. Although results obtained by having seismic input parameter data available are also affected by uncertainties in stellar physics inputs, the effects are much smaller. Mass estimates are affected the most in the seismic case, with errors increasing to about ±7% for the cases considered. In the absence of seismic input data, the average stellar density cannot be determined very well; however with seismic data the average stellar density becomes the most well-determined property, with errors almost completely determined by errors in the seismic parameter Δν. We find that errors in mass and radius estimates can be decreased if we use underlying grids constructed with multiple values of the mixing length parameter.

This work is partially supported by NSF grant AST-1105930 to S.B. W.J.C., Y.E., and G.A.V. acknowledge financial support from the UK Science and Technology Facilities Council (STFC).