BAFFLES: Bayesian Ages for Field Lower-mass Stars

Calcium emission strength, as given by the index ${R}_{{HK}}^{{\prime} }$, is connected to the strength of a star's magnetic field through the stellar dynamo (Noyes et al. 1984). The rotation of the star and convection within induces a magnetic field whose strength is proportional to the rate of rotation (Skumanich 1971; Noyes et al. 1984). Over time, the star's rotation inevitably slows as it ejects ionized particles in its stellar wind, which carry away angular momentum (Kraft 1967; Weber & Davis 1967). As a result, the magnetic field strength—and thus, calcium emission strength—generally decrease with age.

The index ${R}_{{HK}}^{{\prime} }$ is a measure of the flux in the narrow emission line in the core of the Calcium II H and K absorption lines at ∼3968 Å and ∼3934 Å, respectively (Noyes et al. 1984; Wright et al. 2004). Index ${R}_{{HK}}^{{\prime} }$ is derived from an intermediate index, the S index, which represents the ratio of the narrow emission flux to the background continuum flux. The S index provides a relative comparison of emission strength, yet includes both chromospheric and photospheric contributions and is dependent on B − V (as well as age). Therefore, to remove the dependencies on B − V, the S index is transformed by two empirically determined polynomials in B − V, resulting in ${R}_{{HK}}^{{\prime} }$ (Noyes et al. 1984; Wright et al. 2004), where the polynomials have been calibrated over a B − V range of 0.45–0.90, corresponding to an approximate spectral type range of F6 to K2. In addition to the long-term decline in activity over time, the S value for a single star also varies by ∼10% over that star's activity cycle (Wright et al. 2004).

The strength of the lithium absorption line traces the amount of lithium present in the photosphere of a star. When stars initially form, their primordial lithium abundances are similar to the abundance from Big Bang nucleosynthesis, with number densities of ∼10⁻⁹ that of hydrogen (Sestito & Randich 2005). Over time, stars deplete their primordial lithium via nuclear burning in the core and convective mixing, such that measurements of remaining surface lithium correlate with stellar age (Skumanich 1971; Soderblom 2010). For stars cooler than ∼7000 K, lithium abundance can be measured based on the EW of the absorption of the lithium doublet at 6708 Å (Soderblom 2010); hotter stars (OBA spectral types) have ionized their lithium and have negligible 6708 Å absorption even with no lithium burning.

Lithium's two isotopes, ⁶Li and the more abundant ⁷Li, burn at temperatures of 2.2 million K and 2.6 million K, respectively. Since stellar surface temperatures are much lower (∼2500 K for low-mass M stars and ∼46,000 K for high-mass O-stars), in order to burn, lithium must be brought into hotter layers via convection (Soderblom et al. 1990). As a result, the rate of lithium depletion largely depends on the depth of the convection zone, allowing lower-mass stars—which, while having lower surface temperatures, have much deeper convective layers—to deplete lithium faster than higher-mass stars (Soderblom et al. 1990). In addition to convection, it is thought that slow mixing induced by rotation and angular momentum loss may affect lithium depletion (Sestito & Randich 2005), such that lithium abundance is a function of age, spectral type, and the initial rotation rate and rotational evolution of the star.

We seek an expression that returns an age PDF for a single star given an ${R}_{{HK}}^{{\prime} }$ measurement, which is the posterior

$$\begin{eqnarray}&&p(t| \hat{r}),\end{eqnarray} \tag{ 1 }$$

where t is the age and $\hat{r}$ is the measured value of ${R}_{{HK}}^{{\prime} }$ for a single star, with measurement uncertainty of ${\sigma }_{\hat{r}}$. We evaluate this posterior using Bayes' rule

$$\begin{eqnarray}&&p(\theta | D)=\displaystyle \frac{1}{Z}p(D| \theta )p(\theta ),\end{eqnarray} \tag{ 2 }$$

where the four terms are posterior ($p(\theta | D)$), evidence (Z), likelihood ($p(D| \theta )$), and prior (p(θ)), functions of the data (D) and parameters of the model (θ).

For calcium, the parameters of our model, θ, are the age t and the true value of ${R}_{{HK}}^{{\prime} }$ for the star, r, while the data, D, are our measured value of ${R}_{{HK}}^{{\prime} }$, $\hat{r}$. We also assume the evidence, Z, is a constant. With these terms, Bayes' rule becomes

$$\begin{eqnarray}&&p(r,t| \hat{r})\propto p(\hat{r}| r,t)p(r,t).\end{eqnarray} \tag{ 3 }$$

 Our knowledge of the true value of ${R}_{{HK}}^{{\prime} }$ for the star, r, comes from a measurement with an associated measurement error: $\hat{r}$ and ${\sigma }_{\hat{r}}$. In the case of ${R}_{{HK}}^{{\prime} }$, the astrophysical scatter among stars in a single cluster is generally much larger than the measurement uncertainty for any one star. Thus, our model should incorporate both the overall trend that, for clusters of different ages, average r (μ_r) decreases with increasing age, and that there is a scatter about this mean at a single age (σ_r). We expect both these terms to evolve with time, and express them as functions μ_r = f(t) and σ_r = g(t). If the scatter is fit by a Gaussian, our prior on r then becomes

$$\begin{eqnarray}&&p(r| t)={ \mathcal N }(r| f(t),g(t)),\end{eqnarray} \tag{ 4 }$$

while the prior on t, p(t), is flat for a uniform star formation rate, uniform in linear age between 1 Myr and 13 Gyr. Although the star formation rate increases at ages older than ∼8 Gyr, this prior is a reasonable approximation for ages <5 Gyr (Snaith et al. 2015), which also corresponds to the oldest benchmark clusters we utilize. Higher-mass stars have main-sequence lifetimes shorter than the full range of our prior. A stellar lifetime prior is a complicated function of B − V, especially because stars of a given mass evolve in color over time. Rather than commit to a particular set of isochrones, we choose to keep BAFFLES as empirically driven as possible. An isochrone-based age prior can be applied to a BAFFLES posterior once generated, and we advise caution when considering an age posterior with significant probability at very large ages for higher-mass stars. Together, these define a joint prior for our problem

$$\begin{eqnarray}&&p(r,t)=p(r| t)p(t)={ \mathcal N }(r| f(t),g(t))p(t).\end{eqnarray} \tag{ 5 }$$

 In the general case of measurements with Gaussian error bars, likelihood would be given by a normal distribution,

$$\begin{eqnarray}&&{ \mathcal L }(\hat{r}| r)={ \mathcal N }(\hat{r}| r,{\sigma }_{\hat{r}}).\end{eqnarray} \tag{ 6 }$$

However, for ${R}_{{HK}}^{{\prime} }$, we assume that the uncertainty is negligible, especially given the larger astrophysical scatter, σ_r. Therefore, we instead take the likelihood to be a delta function,

$$\begin{eqnarray}&&p(\hat{r}| r)={ \mathcal L }(\hat{r}| r)=\delta (r-\hat{r}).\end{eqnarray} \tag{ 7 }$$

We have no direct data on the age, t, but it is a parameter of our model, so we rewrite the likelihood as

$$\begin{eqnarray}&&p(\hat{r}| r)=p(\hat{r}| r,t)=\delta (r-\hat{r}).\end{eqnarray} \tag{ 8 }$$

We can now rewrite Equation (8), the joint posterior over r and t, as

$$\begin{eqnarray}&&p(r,t| \hat{r})\propto \delta (r-\hat{r}){ \mathcal N }(r| f(t),g(t))p(t),\end{eqnarray} \tag{ 9 }$$

and after marginalizing over r and taking p(t) to be a constant, we solve for $p(t| \hat{r})$,

$$\begin{eqnarray}&&p(t| \hat{r})\propto \int p(r,t| \hat{r}){dr}={ \mathcal N }(r| f(t),g(t)).\end{eqnarray} \tag{ 10 }$$

 If the astrophysical scatter is Gaussian, then by determining functional forms for f(t) and g(t) from our cluster data, we can evaluate the likelihood and produce a posterior for any star given a measurement $\hat{r}$. In Section 3.1.4, however, we present evidence that the scatter is not well-modeled by a Gaussian, and introduce a new numerical function to describe the prior on r.

The derived quantity ${R}_{{HK}}^{{\prime} }$ is formulated to be independent of B − V color, which is accomplished by using two polynomials in B − V to convert the raw S_HK value into ${R}_{{HK}}^{{\prime} }$. To determine the extent to which ${R}_{{HK}}^{{\prime} }$ is in fact independent of color, we initially considered using a two-parameter linear fit to the cluster ${R}_{{HK}}^{{\prime} }$ as a function of B − V, similar to that in Mamajek & Hillenbrand (2008), since the slopes seemed non-negligible. However, since our data set included many clusters with only a handful of calcium measurements, the fit slopes were poorly determined and the fits crossed frequently, the latter being a nonphysical outcome. As in the right panel of Figure 1, linear fits to the clusters resulted in nonmonotonic changes in ${R}_{{HK}}^{{\prime} }$ over time, especially in the reddest and bluest regions of our B − V range. Although Mamajek & Hillenbrand (2008) used linear fits for each cluster, they interpolated cluster means for solar B − V (∼0.65) only. For solar B − V, the cluster means are still monotonic, something not true for other B − V values that were included in our study.

Zoom In Zoom Out Reset image size

There is a significant improvement in ${\chi }^{2}$ from the linear fit to the constant fit, dropping from 554 (constant) to 423 (linear), assuming a constant measurement error for each star of 0.1 dex, as estimated by Mamajek & Hillenbrand (2008). Based on the Bayesian information criterion, this presents very strong evidence in favor of the linear model (ΔBIC = 83.3). Nevertheless, we find the behavior of the linear fits in the right panel of Figure 1 to be unphysical: at the reddest and bluest ends, the evolution in ${R}_{{HK}}^{{\prime} }$ is nonmonotonic and implies wild swings in calcium activity as a function of age, based on a handful of data points in each cluster, as well as poor sampling across the entire B − V range. As a result, to avoid overfitting sparse data, we adopted a constant fit for ${R}_{{HK}}^{{\prime} }$, where each cluster is represented by the median value of ${R}_{{HK}}^{{\prime} }$, with no B − V dependence. A constant fit has the advantage of capturing the monotonic decrease in ${R}_{{HK}}^{{\prime} }$ while remaining the simplest fit. Mamajek & Hillenbrand (2008) advocate determining age from ${R}_{{HK}}^{{\prime} }$ through an age–activity–rotation relation, the effect of which is a significant B − V dependence on ${R}_{{HK}}^{{\prime} }$ for objects of similar ages (see their Figure 11), which varies by ∼0.15 dex across B − V. As there are limited ${R}_{{HK}}^{{\prime} }$ measurements in benchmark clusters, it is currently difficult to confirm this behavior of ${R}_{{HK}}^{{\prime} }$ as a function of color. In fact, more direct solutions to a B − V dependence of ${R}_{{HK}}^{{\prime} }$ would be to either redetermine the polynomial parameters or to fit directly in S_HK, and either would likely require a larger data set than that presented here.

From the fits above, we have nine cluster ages and their respective mean log(${R}_{{HK}}^{{\prime} }$) values, which we use to find the mean log(${R}_{{HK}}^{{\prime} }$) at all ages covered by our prior, μ_r = f(t). We fit log(${R}_{{HK}}^{{\prime} }$) as a function of age with a second-order polynomial, constrained to be monotonically decreasing, and where each cluster in the fit is weighted by the number of stars it contains. Figure 2 shows this fit against the median value of each cluster, with plotted error bars indicating the standard deviation in each cluster. Our fit is consistent with polynomial fits from previous authors. Mamajek & Hillenbrand (2008) use linear fits for finding each cluster's mean ${R}_{{HK}}^{{\prime} }$ as a function of B − V, and then fit a third-order polynomial to age, based on the value of each cluster's linear fit evaluated at solar B − V of 0.65. The largest discrepancies between the two fits are, unsurprisingly, at ages lower than that of the youngest benchmark cluster (Upper Sco) and larger than that of the oldest (M67). Soderblom et al. (1991) experimented with several different second-order polynomials, correcting for disk heating and a uniform star formation rate.

Zoom In Zoom Out Reset image size

We next examine the astrophysical scatter of ${R}_{{HK}}^{{\prime} }$ about the mean, σ_r = g(t). We begin by computing the residuals of ${R}_{{HK}}^{{\prime} }$ for every star in a cluster to the median value for all stars in the cluster. The standard deviations of these residuals are plotted in Figure 3, where uncertainty in the standard deviation (σ_m) of the mth cluster with N_m stars is given by the equation appropriate for Gaussian scatter, $\pm \tfrac{{\sigma }_{m}}{\sqrt{2{N}_{m}-2}}$. There is some evidence that the scatter between 20 and 200 Myr is larger than the scatter for younger or older stars. This is reminiscent of Figure 1 of Gallet & Bouvier (2013), where solar-type stars spin up between ∼10 and 50 Myr as they contract when approaching the main sequence, and the dispersion in rotation rate between the fast rotators and slow rotators in a single cluster increases, compared to stars younger than 20 Myr or older than 200 Myr. We investigated using a Gaussian or inverted parabola to fit the data, but we do not have strong evidence that such a fit is justified. The clusters with the most measurements, Pleiades, Hyades, and M67, show the strongest evidence for a change in the standard deviation with time. However, the clusters Upper Sco, UCL+LCC, β Pic, Tuc/Hor, α Per, each have only a few stars (8, 8, 6, 6, and 12 respectively), meaning the standard deviations are not well-determined. As a result, we treat g(t) as a constant and note that larger data sets with more stars in these young clusters would be needed in order to precisely measure time-dependent behavior.

Zoom In Zoom Out Reset image size

 We instead compute the residuals for each star with respect to the fit f(t), and observe that these residuals, while somewhat symmetric, are not well-fit by a Gaussian (the dashed gray curve in Figure 4). In particular, the best-fit Gaussian underestimates the peak and is significantly wider at ∼1σ compared to the data. We thus return to Equation (4) and replace the normal distribution in our prior with a new numerical function,  

$$\begin{eqnarray}&&p(r| t)=\ { \mathcal H }(r| f(t)),\end{eqnarray} \tag{ 11 }$$

which has a mean f(t). The amplitude of the scatter is encoded by ${ \mathcal H }$ itself, and since we take g(t) to be a constant, the shape of ${ \mathcal H }$ does not change over time, only its mean.

Zoom In Zoom Out Reset image size

To evaluate ${ \mathcal H }$, we fit a function to the smoothed CDF of the residuals, as in the right panel of Figure 4, capturing the non-Gaussian shape of the astrophysical scatter. The tails of the PDF are constrained to decrease exponentially out to ∼4 standard deviations and then are fixed at zero. We perform the smoothing with a Savitzky–Golay filter, which fits successive windows with a third-order polynomial so as to remove jumps in the function from star to star without significantly increasing the width of the distribution. The final function is a good fit to the data (see left panel of Figure 4). We normalize the final ${ \mathcal H }$ distribution so that it has an integral of unity. Then, since we find ${ \mathcal H }$ to be slightly asymmetric with a median of 0.0026 dex, we shift the distribution so that its median has a value of zero. 

The numerical fit ${ \mathcal H }$ is one of many possible implementations of the prior function. Other choices with wider tails (such as the Student's-t distribution or the Lorentzian distribution) also partially capture the non-Gaussian behavior. We found the best fit with the Student's-t distribution, which came closest to matching the residual distribution. We also found no significant difference between it and the empirical function ${ \mathcal H }$ on our final age posteriors.

Rewriting Equation (8), we then have an expression for our posterior given by  

$$\begin{eqnarray}&&p(r,t| \hat{r})\propto \delta (r-\hat{r}){ \mathcal H }(r| f(t))p(t),\end{eqnarray} \tag{ 12 }$$

where f  and ${ \mathcal H }$ are determined from our cluster data sets above, and p(t) is constant for a uniform star formation rate. We can then rewrite our calcium age posterior in Equation (10) using our function ${ \mathcal H }$: 

$$\begin{eqnarray}&&p(t| \hat{r})\propto { \mathcal H }(\hat{r}| f(t)).\end{eqnarray} \tag{ 13 }$$

We implement this method with an array of 1000 elements uniformly sampled in log age from 1–13,000 Myr, and we evaluate Equation (13) at each point in the array for the $\hat{r}$ of a single star. These probabilities are then normalized to integrate to unity (accounting for uneven bin sizes) and provide the age posterior for that star.

As with calcium, lithium EW decreases with time, with an astrophysical scatter about this trend for objects of the same age. Following the framework we developed for calcium, we define functions for the mean EW as a function of time (i), the standard deviation about that mean (j), and the shape of the distribution function about the mean (${ \mathcal K }$). These three functions are the lithium equivalents of f, g, and ${ \mathcal H }$ used above for calcium. The mean i(t, b) is decidedly a function of both age (t) and B − V color (b). However, when we consider the log of the EW (l), the scatter about this mean appears to be independent of color, so we define j(t) as a function of time only. The parameters of our model are l, b, and t, requiring a joint prior in all three for Bayes' equation, p(l, b, t),  

$$\begin{eqnarray}&&p(l,b,t)=p(l| b,t)p(b)p(t)={ \mathcal K }(l| i(t,b),j(t))p(b)p(t),\end{eqnarray} \tag{ 14 }$$

where p(b) is the prior on B − V color, which we take to be flat, since we will generally have a precise measurement of color for a given star, and p(t) is the prior on age, again flat for a constant star formation rate.

We assume a Gaussian likelihood for both l and b, given measurements of $\hat{l}$ and $\hat{b}$, and measurement errors of ${\sigma }_{\hat{l}}$ and ${\sigma }_{\hat{b}}$,

$$\begin{eqnarray}&&{ \mathcal L }(\hat{l},\hat{b})={ \mathcal N }({10}^{\hat{l}}| {10}^{l},{\sigma }_{\hat{l}}){ \mathcal N }(\hat{b}| b,{\sigma }_{\hat{b}}),\end{eqnarray} \tag{ 15 }$$

where 10 is raised to the power of l and $\hat{l}$ because, while l is a log quantity, measurement errors are typically quoted in linear units (e.g., mÅ). Combining likelihood and prior, and again assuming the evidence to be constant, we obtain an expression for the posterior  

$$\begin{eqnarray}&&p(l,b,t| \hat{l},\hat{b})\propto { \mathcal N }({10}^{\hat{l}}| {10}^{l},{\sigma }_{\hat{l}}){ \mathcal N }(\hat{b}| b,{\sigma }_{\hat{b}}){ \mathcal K }(l| i(t,b),j(t)),\end{eqnarray} \tag{ 16 }$$

which, when marginalized over l and b, gives the marginalized posterior on age,

$$\begin{eqnarray}&&p(t| \hat{l},\hat{b})=\int \int p(l,b,t| \hat{l},\hat{b}){dldb}.\end{eqnarray} \tag{ 17 }$$

As with calcium, all that remains is to define the functions i(t, b), j(t), and ${ \mathcal K }(l| i(t,b),j(t))$ from our cluster data.

For a single cluster, the log of the EW, l, appears as a Gaussian or parabola as a function of B − V, as shown in Figure 5. The reddest and bluest stars in the cluster tend to have the smallest values for lithium EW, while intermediate B − V stars (G stars) have the largest lithium EW. This behavior is the result of two primary processes. First, redder, lower-mass stars have deeper convective envelopes, so they more quickly convect lithium to deeper, hotter layers of the star, where it is fused, resulting in faster depletion of lithium. Meanwhile, blue stars have hotter photospheres, so there are fewer lithium atoms in the ground state to absorb 6708 Å light. Stars are expected to have uniform lithium abundance (N(Li)) at formation, but this translates to a range of Li EW values as a function of color, given the different photospheric temperatures across this range. In addition to the Gaussian shape, the "lithium dip" is observed for stars between B − V of ∼0.36 and ∼0.42 (6900 K and 6600 K) for stars that are ≳500 Myr, where there is a decrease in lithium abundance in this narrow range compared to stars on either side of the dip (Boesgaard & Tripicco 1986; Balachandran 1995). A suggested explanation for the lithium dip is that, at the hot end of the dip, magnetic field strength is increasing with decreasing stellar mass, spinning down the outer layers of the star and creating turbulent mixing from internal shear between these layers and the faster-rotating core. Moving to the cooler end of the dip thus corresponds to the rise of internal gravity waves, which more efficiently spin down the core, such that there is less turbulent mixing (Talon & Charbonnel 2010). Under this model, surface lithium is preferentially destroyed in the narrow region of the lithium dip, while it is preserved on either side.

Zoom In Zoom Out Reset image size

For each cluster, we simultaneously fit both the mean  and the standard deviation of l at a single value of t = t_m, where t_m is the age of the given cluster m. We take these fits, $i^{\prime} ({t}_{m},b)$ and $j^{\prime} ({t}_{m})$, as preliminary values for the mean, i(t, b), and standard deviation, j(t), evaluated at the age of the cluster. We assume a functional form of a second-order polynomial for $i^{\prime} ({t}_{m},b)$, while at a single age the standard deviation, $j^{\prime} ({t}_{m})$, is a constant that does not depend on color. To fit these parameters, we assumed a Gaussian likelihood, which for lithium detections takes the form  

$$\begin{eqnarray}&&p(\hat{l}| {t}_{m},b)={ \mathcal N }(\hat{l}| l,{\sigma }_{\hat{l}})=\displaystyle \frac{1}{\sqrt{2\pi }j^{\prime} ({t}_{m})}{e}^{-\tfrac{{\left(\hat{l}-i^{\prime} ({t}_{m},b)\right)}^{2}}{2j^{\prime} {\left({t}_{m}\right)}^{2}}}.\end{eqnarray} \tag{ 18 }$$

For lithium upper limits ($\hat{u}$), we represent the likelihood as the integral of the Gaussian function from $-\infty$ to the upper limit $\hat{u}$,  

$$\begin{eqnarray}&&p(\hat{u}| {t}_{m},b)={\int }_{-\infty }^{\hat{u}}\displaystyle \frac{1}{\sqrt{2\pi }j^{\prime} ({t}_{m})}{e}^{-\tfrac{{\left(l-i^{\prime} ({t}_{m},b)\right)}^{2}}{2j^{\prime} {\left({t}_{m}\right)}^{2}}}{dl},\end{eqnarray} \tag{ 19 }$$

and then fit these four parameters—three for the polynomial in b that defines $i^{\prime} ({t}_{m},b)$, and one for the standard deviation, $j^{\prime} ({t}_{m})$. The fit itself is performed by assigning one of these likelihoods to each star, based on whether there is a lithium measurement or upper limit, then maximizing the product of likelihoods over all cluster stars.

The lithium dip is clear in the ∼700 Myr Hyades data set, so we fit an inverted Gaussian to the dip (0.39 < B − V < 0.52) and a second-order polynomial to the stars outside the dip. There is no clear evidence for a lithium dip at younger ages in M34 (∼200 Myr). By ∼4 Gyr, stars have evolved off the main sequence, leaving no stars bluer than B − V ≈ 0.5 in M67. Fits to each cluster are shown in Figures 5 and 6.

Zoom In Zoom Out Reset image size

Unlike ${R}_{{HK}}^{{\prime} }$, Li EW varies substantially across both age and B − V. As a result, the decline in lithium EW as a function of age must be fit across multiple slices of B − V. The polynomials fit to each cluster (Figure 5) define our preliminary fits $i^{\prime} ({t}_{m},b)$ as a function of B − V at the age for each individual cluster (t_m). We examine 64 B − V slices uniformly spaced between B − V of 0.35 and 1.9. At each slice n, then, we have 10 values of $i^{\prime} ({t}_{m},b)$ corresponding to our 10 cluster data sets to which we add two additional points from primordial Li EW and Blue Lithium Depletion Boundary (BLDB, described below). From these, we determine the 64 fits, $i(t,{b}_{n})$, as a function of age.

To help constrain the young end of the fits of Li EW and age, we approximate primordial Li EW using the MIST model isochrones (Choi et al. 2016) in conjunction with our NGC 2264 fit. In particular, we seek to extend the fit to this ∼5 Myr cluster to the first age point in our grid, 1 Myr. At every B − V value, we determine the corresponding effective temperature using the conversions from Pecaut & Mamajek (2013); we then find the Li abundance, N(Li), and initial stellar mass at 5 Myr from the MIST isochrones. Next, we find the Li abundance from the same initial mass star using the 1 Myr isochrones. We convert T_eff and the Li abundance to Li EW using the curve of growth in Soderblom et al. (1993) for T_eff > 4000 K, and that in Zapatero Osorio et al. (2002) for T_eff ≤ 4000 K. The difference in Li EW between 1 and 5 Myr is added to the fit to NGC 2264 in order to determine the primordial Li EW at every B − V value (Figure 6). The change found in Li EW between 1 and 5 Myr is only significant between 0.8 ≲ B − V ≲ 1.4, and is negligible elsewhere.

We define the BLDB as the B − V color for which stars redder than this boundary have no detectable lithium absorption, which we use to help constrain the older and redder range of fits to Li EW against age. Since the redder stars have deeper convective envelopes, they burn lithium faster than the bluer stars in the cluster. As a result, the nested polynomials of Figure 6 generally get narrower and move blueward over time, and thus the BLDB point moves blueward with increasing age. The BLDB is distinct from the classical LDB, which moves redder over time as a cluster's high-mass brown dwarfs deplete their lithium, while at the same time, all brown dwarfs evolve to redder colors as they cool over time, outside the brief deuterium-burning phase. We have defined the BLDB in order to add an additional data point to our fits for B − V > 0.7, as these are most important for constraining the ages of stars with B − V > 1.4, for which there are fewer literature measurements, especially at older ages. For each B − V slice redward of 0.7, the fit to BLDB points gives an approximation of maximum age associated with log(Li EW) = 0.5, or Li EW = 3.2 mÅ (Figure 7).

Zoom In Zoom Out Reset image size

For each value of B − V, we use the mean value of l, $i^{\prime} ({t}_{m},{b}_{n})$, from each cluster, in addition to the primordial lithium point and BLDB point, to fit the intermediate ages between the cluster ages and complete our grid. Unlike calcium, where the fits to individual clusters were independent of B − V, for lithium there is a strong B − V dependence, and for redder regions, the fit to the mean EW reaches unphysically small values. When cluster means drop below ∼3 mÅ (0.5 on the log scale), we do not expect any detections, and clusters with $i^{\prime} ({t}_{m},b)$ below this value are not included in the fitting process.

As in Figure 8, we fit a 2–4 segment piecewise function to the cluster means, primordial Li EW, and BLDB point. The first segment is always between the Primordial Li EW value and NGC 2264, and the fit is constrained to decrease monotonically with age. Additionally, for B − V slices inside the lithium dip (0.41 ≤ B − V ≤ 0.51), the final piecewise segment is constrained to go through the Hyades point. The locations of the segment breaks (except for the first break at NGC 2264) were free parameters. Weights for the cluster means were determined based on the relative proportion of stars the cluster had at a given B − V slice in relation to the total number of stars. The BLDB point is given an uncertainty of about 0.15 dex, compared to 1 dex for clusters in poorly constrained regions. Although different functional forms were good fits to the decrease in lithium over time for some B − V ranges, only the piecewise function was flexible enough to capture the shape more generally.

Zoom In Zoom Out Reset image size

With our grid of 64 B − V slices and 1000 age slices for mean EW of lithium, i(t, b), we next empirically determine the distribution of the residuals, ${ \mathcal K }(l| i(t,b),j(t))$, as we did with calcium (Figure 9). Residuals are with respect to the value of i(t, b) evaluated at the age of each cluster and the B − V value of the star, and upper limits are not considered in this step. As with calcium, we smooth the CDF of the residuals with Savitzky–Golay filters and take the derivative to convert to a PDF. We next fit exponential functions to the two tails, which we connect with the smoothed PDF, then normalize to have integral unity, defining ${ \mathcal K }(l| i(t,b),j(t))$. We also center the distribution at zero by subtracting off the residual median value of 0.033, which ensures that ${ \mathcal K }$ does not introduce a systematic bias toward older ages. Unlike calcium, we find no evidence for even a weak dependence on time of the standard deviation of the residuals (j(t)). As a result, the shape of ${ \mathcal K }$ is not a function of age or color, while the mean value is.

Zoom In Zoom Out Reset image size

Since we see no evidence for an age dependence in the astrophysical scatter, we take j(t) to be a constant, and slightly rewrite our posterior from Equation (16),  

$$\begin{eqnarray}&&p(l,b,t| \hat{l},\hat{b})\propto { \mathcal N }({10}^{\hat{l}}| {10}^{l},{\sigma }_{\hat{l}}){ \mathcal N }(\hat{b}| b,{\sigma }_{\hat{b}}){ \mathcal K }(l| i(t,b)).\end{eqnarray} \tag{ 20 }$$

To determine the age posterior for a single star, we construct a dense grid covering B − V from 0.35 to 1.9 and age from 1 to 13,000 Myr. We use a grid of 64,000 elements (64 B − V × 1000 age, logarithmically spaced in age), with the mean lithium abundance i(t, b) calculated at each gridpoint. At each combination of (t, b), we first marginalize over l by multiplying ${ \mathcal K }(l| i(t,b))$ (our prior) by ${ \mathcal N }({10}^{\hat{l}}| {10}^{l},{\sigma }_{\hat{l}})$, a Gaussian representing the measurement of lithium EW and the associated measurement error, over an array of 1000 elements, logarithmically spaced in l between 0.5 and 1585 mÅ. To do this multiplication, however, we first convert ${ \mathcal K }(l| i(t,b))$ to a function (similar to a log-normal) in linear space because the likelihood for $\hat{l}$, ${ \mathcal N }({10}^{\hat{l}}| {10}^{l},{\sigma }_{\hat{l}})$, is defined in linear space. If no measurement uncertainty is given, we use a default error for ${\sigma }_{\hat{l}}$ of 15 mÅ, which is noted as a typical error by Soderblom et al. (1993). The products of these functions evaluated at all 1000 points are then summed, which gives the probability at that specific (t, b) gridpoint. We then marginalize over B − V color by weighting each (t, b) gridpoint by the Gaussian likelihood for b, ${ \mathcal N }(\hat{b}| b,{\sigma }_{\hat{b}})$, (assuming ${\sigma }_{\hat{b}}$ = 0.01 mag if no error is given) and summing over the product. To minimize computation time, instead of computing this probability at all (t, b) locations, we only evaluate gridpoints at 15 sampled values of b within $4\cdot {\sigma }_{\hat{b}}$ of $\hat{b}$, with all others set to 0. Having marginalized over both l and b, we are left with a marginalized posterior over only age, $p(t| \hat{l},\hat{b})$.

If the measurement $\hat{l}$ is an upper limit $\hat{u}$, we instead integrate ${ \mathcal K }(l| i(t,b))$ from $-\infty$ to $\hat{u}$ to find the probability at each (t, b) gridpoint. Thus, upper limits result in a plateau of probability at old ages, with a rapid drop-off toward younger ages.