ISOCHRONES FOR OLD (>5 GYR) STARS AND STELLAR POPULATIONS. I. MODELS FOR −2.4 ⩽ [Fe/H] ⩽+0.6, 0.25 ⩽ Y ⩽ 0.33, AND −0.4 ⩽ [α/Fe] ⩽+0.4

To model the lowest masses (each grid has a minimum mass of $0.12 \,{\cal M_{\odot }}$), MARCS model atmospheres (Gustafsson et al. 2008) at an optical depth τ = 100 were attached to the interior structures following the procedures described by VandenBerg et al. (2008). In the case of 0.3–$0.4 \,{\cal M_{\odot }}$ models (or somewhat higher masses in the metal-rich or super-metal-rich regimes), the stellar photosphere was taken to be the outer boundary and the pressure at T = T_eff was determined by integrating the hydrostatic equation from very small optical depths to the photospheric value, assuming the semi-empirical Holweger & Müller (1974, hereafter HM74) T–τ structure (specifically, the fit to the latter given by VandenBerg & Poll 1989). (When using the Sun to calibrate such quantities as the convective mixing-length parameter, it is obviously important to assume the solar atmospheric structure instead of, say, a gray atmosphere (see, e.g., Morel et al. 1994). Encouragingly, the temperature stratification predicted by recent three-dimensional (3D) model atmospheres for solar parameters appears to satisfy observational constraints even better than the HM74 model, which is preferable to current one-dimensional (1D) model atmospheres in representing the surface layers of the Sun; see the discussion by Pereira et al. (2013).

The ramifications of different treatments of the atmospheric layers for low-mass, [Fe/H] =0.0 stellar models are shown in Figure 1. (For an instructive example of similar work carried more than 15 yr ago, see Brocato et al. 1998.) Using the properties of MARCS model atmospheres at τ = 100 to derive the outer BCs of stellar models clearly results in increasingly cooler temperatures and reduced luminosities with decreasing mass than attaching the same atmospheres to the interior structures at T = T_eff or deriving the boundary pressure using the scaled HM74 T–τ relationship (compare the locations of the open circles along the solid, dashed, and dot–dashed curves, respectively). For these computations (and, indeed, for all of the lower-MS models in which the MARCS atmospheres were attached at depth), the advanced equation-of-state (EOS) developed by A. Irwin was used in its most efficient "EOS4" mode.⁵ Because Irwin's EOS, even in the EOS4 mode, is slower by a factor of three to four than the EOS which is normally employed by the Victoria code (see VandenBerg et al. 2000), we have opted to use the latter for higher mass models (where the tracks are essentially independent of this choice, see V12).

Zoom In Zoom Out Reset image size

Since the lowest mass models are based a different EOS and surface BCs than those for higher masses, it is necessary to make a smooth transition between the two regimes. To accomplish this, evolutionary tracks for masses in the range 0.3–$0.8 \,{\cal M_{\odot }}$ were computed (for each combination of the chemical abundance parameters) to just past the zero-age main-sequence (ZAMS) location on the H-R diagram, assuming the EOS and atmospheric BCs employed in the lower-main-sequence (LMS) grids, on the one hand, and those used in the higher mass tracks, on the other. The mass for which the differences in log  L/L_☉ and log  T_eff of the respective ZAMS models were the smallest was taken to be the transition mass, and the mean luminosity and temperature differences, which were usually <0.003 in both log  L/L_☉ and log  T_eff, were calculated. These offsets were applied to a complete track for the transition mass (assuming the physics that has been employed for higher masses), which became the adopted track for that mass. No adjustments of any kind were made to the evolutionary sequences for lower or higher masses.

As shown in Figure 2 for a subset of the model grids, the resultant H-R diagrams (upper panels) and mass–luminosity relations (lower panels) are very smooth in the LMS region where this join has been made (and elsewhere). Only in magnified versions of this and similar plots are some very slight irregularities evident; e.g., the spacing between the fourth and fifth loci close to the $0.4 \,{\cal M_{\odot }}$ ZAMS models in the upper panel for Y = 0.33 is a bit larger than those between the third and fourth or fifth and sixth loci. However, they have no obvious impact on the predicted mass–luminosity relations. Note that some kind of a transition is unavoidable because proper model atmospheres, attached at depth, must be used as BCs for very low mass models in order to obtain the most realistic T_eff scale, whereas scaled HM74 atmospheric structures are the preferred choice for the Sun and solar-type stars (as discussed below).

Zoom In Zoom Out Reset image size

Figure 1 also plots the lower-MS portion of a 4 Gyr isochrone for solar abundances (the dotted curve) that was kindly provided to us by G. Feiden (2011, private communication). The BCs for these models (see Feiden et al. 2011) were derived from PHOENIX model atmospheres (Hauschildt et al. 1999): the latter were attached to the interior structures at τ = 100 if ${\cal M} \le 0.2 \,{\cal M_{\odot }}$, or at the photosphere, in the case of higher masses. Remarkably, the lowest mass models overlay those represented by the solid curve nearly perfectly, while at ${\cal M} \gtrsim 0.2 \,{\cal M_{\odot }}$, they are cooler than our atmosphere-interior models by ≲ 50 K. This is really quite good agreement between completely independent predictions of the properties of very low mass stars.

At lower metallicities, the differences between the present Victoria–Regina and published Dartmouth models (Dotter et al. 2007b) are even less, as shown in the left-hand plot in Figure 3. Although they employ different model atmospheres as BCs, the latter are apparently sufficiently similar to MARCS atmospheres that they yield nearly the same T_eff scale. Moreover, it is apparent from the near coincidence of the filled and open circles that the predicted mass–luminosity relations (especially at masses ${>} 0.15 \,{\cal M_{\odot }}$) are in excellent agreement as well. However, perhaps the most compelling demonstration of the reliability of our computations for LMS stars is provided in the right-hand plot, which illustrates how well our models reproduce the mass–radius (MR) relation that describes the lowest mass (0.213 and $0.241 \,{\cal M_{\odot }}$) components of the triple system, KOI-126 (Carter et al. 2011). A similar (equally successful) comparison, but using Dartmouth models, was reported by Feiden et al. (2011), who note that some low-mass binaries (notably CM Draconis) continue to be problematic, possibly due to the neglect of magnetic fields and activity effects (also see Feiden & Chaboyer 2013). Regardless, Figure 3 and additional plots presented later in this paper (see Section 5) provide encouraging support for the quality of our MS and LMS models. (Although discrepancies between predicted and observed CMDs are generally found at the faintest absolute magnitudes (e.g., Richer et al. 2008; Casagrande & VandenBerg 2014), they are likely due mostly to deficiencies in current color–T_eff relations, given the great difficulty of accounting for all of the sources of blanketing in model atmospheres and synthetic spectra for cool stars.)

Zoom In Zoom Out Reset image size

One question that warrants some discussion is the following: why has the scaled solar HM74 T–τ relation been used to determine the boundary pressures for the majority of our computed models instead of MARCS model atmospheres? As already mentioned, the former is the preferred choice for the calculation of a Standard Solar Model, and presumably for models which are relevant to stars with properties similar to that of the Sun. In fact, VandenBerg et al. (2008) showed that surface BCs derived in this way agree rather well with those obtained from scaled, differentially corrected (SDC) MARCS models over wide ranges in T_eff, gravity, and metallicity. (Note that the latter were constructed in order that the resultant solar model atmosphere reproduces the temperature structure derived by HM74.) While it is not necessarily the case that the SDC atmospheres provide a better representation of those applicable to, e.g., metal-deficient stars than standard MARCS models, the implied T_eff scale agrees quite well with that derived by Casagrande et al. (2010) (via the infrared flux method, IRFM) for field subdwarfs that have −2.0 ≲ [Fe/H] ≲ −0.6 (see VandenBerg et al. 2010, and Section 5.3 later in this paper).

As regards the VandenBerg et al. (2008) examination of the use of MARCS model atmospheres as BCs for interior structures, it is pertinent to note that an identical treatment of convection and the same chemical abundances were assumed in the Victoria and MARCS codes. Indeed, plots were included in the paper by VandenBerg et al. to illustrate the close agreement of, among other things, the predicted variations with depth (in the atmosphere) of the pressure, the adiabatic temperature gradient, the opacity, and the convective flux. It is not possible to obtain similar consistency with the large grids of model atmospheres published by Gustafsson et al. (2008) because, for one thing, they assumed a value of 1.5 for the mixing-length parameter α_MLT, whereas the evolutionary models presented in this paper required α_MLT = 2.007 to satisfy the solar constraint. There are also minor differences in the adopted heavy-element abundances, which could have some effect on the opacities at low temperature. In addition, the MARCS atmospheres were computed for a fixed helium abundance (log N(He) =10.93, on the scale log N(H) =12.0, or Y ≈ 0.26). In stellar models that take gravitational settling into account, the surface value of Y can fall to quite small values when the envelope convection zones become very thin (see V12). Such variations will affect the mean molecular weight in the atmospheric layers and presumably have some repercussions for the temperature structure.

Even though suitable model atmospheres for use as BCs in the computation of diffusive stellar models are not currently available, we did some limited explorations of the effects on the predicted T_eff scale of using the current MARCS grids in this way. The results of those experiments are shown in Figure 4, which plots evolutionary tracks for the indicated masses and chemical abundances, on the assumption of different treatments of the atmosphere. The dotted track for [Fe/H] =0.0 differs only slightly from one that passes through log T_eff = 3.7617 (T_eff = 5777 K) and M_bol = 4.75 (the solar temperature and bolometric luminosity) at the solar age (4.57 Gyr) insofar as it assumed Y = 0.25 whereas a Standard Solar Model requires Y = 0.2553. The solid and dashed loci for the same metallicity assume α_MLT = 2.007, as in the case of the dotted curve (and indeed, all other evolutionary sequences that have been computed), but instead of employing BCs based on the HM74 T–τ relation, interpolations in the tabulated properties of the MARCS model atmospheres at T = T_eff or at τ = 100, respectively, were carried out to fit the atmospheres to the interior structures.

Zoom In Zoom Out Reset image size

In contrast with the findings of VandenBerg et al. (2008), who found fairly small differences between the tracks (for similar masses and chemical compositions) when MARCS atmospheres were fitted to non-diffusive stellar models at the photosphere or at depth, the solid and dashed tracks are appreciably offset from each other. Because they run roughly parallel to one another, it seems more likely that differences in the low-T opacities or of the assumed value of α_MLT is responsible for this separation since, at the ZAMS location, diffusion has not had enough time to significantly alter the surface abundances. Regardless, one could apply suitable adjustments to the pressures predicted by the MARCS atmospheres, either at the photosphere or at τ = 100, to force consistency of these cases with the solar constraint. Doing so results in the long-dashed and dot–dashed curves, which are not very different from the dotted track. The solar calibration thus compensates for most of the differences in the assumed physics.

However, the computations for [Fe/H] = −2.40 show that there are systematic variations in the tracks as a function of metallicity. When MARCS atmospheres are used as BCs, with or without ad hoc adjustments to the pressures at the photosphere or at τ = 100, the resultant tracks are all hotter than the dotted track (by as much as 200 K, see Figure 4). Whether or not this is telling us that the isochrones presented in this study for low [Fe/H] values are too cool is hard to say (though it is tempting to conclude that this is probably the case). On the one hand, the use of the same T–τ structure for all stellar models, regardless of mass, metallicity, and evolutionary state, can hardly be realistic. On the other hand, our isochrones appear to be able to reproduce the properties of local subdwarfs with Hipparcos parallaxes quite well (see Section 5.3, and VandenBerg et al. 2010). It is always possible, for instance, that errors associated with the surface BCs are compensating for those arising from the treatment of convection or from other physics ingredients. The observed T_eff scale is simply not yet precise enough to provide good constraints on the temperatures predicted by stellar models.

The final point worth making here is that, as shown by V12, the Victoria and recent MESA (Paxton et al. 2011) evolutionary codes produce nearly identical tracks when very close to the same physics is assumed. Indeed, predictions of such quantities as the age, luminosity, and helium core mass at the RGB tip are also in excellent agreement. The same can be said of the tracks produced by the Dartmouth code (Dotter et al. 2007b), as those computations also appear to be nearly indistinguishable from ours when the same mass and chemical abundances, and very similar input physics, are adopted (see Brogaard et al. 2012, their Figure 3). Judging from, e.g., the plots provided by Bressan et al. (2012), such good consistency between the results of the various evolutionary codes currently in use is not always obtained, though the extent to which differences in the physics are responsible for this is not clear. Efforts should be made to understand the origin of any such discrepancies that are found.

M67 is younger than the age range focused on by the present investigation, but it is still worthwhile to fit isochrones to its CMD in order to see how well solar abundance models reproduce its MS and RGB. High-resolution spectroscopy has revealed that this system has nearly the same [Fe/H] value as the Sun and very close to the solar mix of the metals (Randich et al. 2006; Önehag et al. 2011).⁸ Moreover, as both the earlier dust maps by Schlegel et al. (1998) and the recent recalibration of them by Schlafy & Finkbeiner (2011, hereafter SF11) yield a reddening that is consistent with E(B − V) = 0.030 ± 0.003, the Sun may be used to provide quite an accurate estimate of the cluster distance modulus via the MS-fitting technique because of the similarity in their ages: according to most estimates, M67 is only ∼0.4–0.7 Gyr younger than the Sun (see, e.g., Richer et al. 1998; Michaud et al. 2004; Önehag et al. 2011).

As predictions of synthetic magnitudes for longer wavelength filters are likely to be less problematic than those for blue or ultraviolet bandpasses, we have opted to fit isochrones to the VK_S photometry for M67 that was analyzed by Brasseur et al. (2010). The latter combined a recent reduction of Johnson–Cousins V-band data for M67 (by one of the coauthors; namely, P. B. Stetson, see Stetson 2005) with Two Micron Sky Survey (2MASS) K_S photometry (Skrutskie et al. 2006) to produce a CMD that is well defined down to V ∼ 19. After the colors of the stars have been dereddened, assuming E(V − K_S) = 2.76 E(B − V) (Casagrande & VandenBerg 2014), and their apparent magnitudes have been decreased by 9.69 mag (the adopted apparent distance modulus), one obtains the CMD that is shown in Figure 9. The vertical offset was chosen so that the cluster MS lies just slightly fainter than the Sun (which has M_V = 4.82) at the solar color (V − K_S = 1.560; Casagrande et al. 2012), to be consistent with the aforementioned age difference. The isochrone that provides the best fit to the cluster subgiants has an age of 4.3 Gyr. (Interestingly, the age corresponding to a given magnitude difference between the SGB and the MS, at a fixed color, depends quite sensitively on the adopted value of Z. As shown by VandenBerg et al. (2007, their Figure 2) and Michaud et al. (2004), ages ≲ 4 Gyr are obtained for M67 if Z ≳ 0.017, and vice versa.) Encouragingly, the MARCS transformations yield the same V − K_S color for the Sun as that derived by Casagrande et al. to within 0.006 mag.

Zoom In Zoom Out Reset image size

In their study of M67, Brasseur et al. (2010) employed the 3.7 Gyr isochrone that was derived by Michaud et al. (2004) from models that took gravitational settling and radiative accelerations into account (and they made somewhat different assumptions about the solar normalization and initial abundances). Unlike our isochrone, it provides an excellent fit to the morphology of the TO, including the luminosity of the gap. In order for our computations to have similar success, it will be necessary to consider diffusion of the metals and to calibrate the extent of convective core overshooting. However, our interest here is not to derive the age of M67, but rather to ascertain how well our solar abundance models are able to reproduce those parts of a CMD (namely, the MS and RGB) that have no more than a weak dependence on age, and in addition, are virtually independent of whether or not diffusive processes and core overshooting are treated. Indeed, in these respects, our models fare quite well as they provide a very good match to both the cluster giant branch and the lower MS (down to M_V ∼ 8.5). This suggests that the predicted T_eff scale for solar metallicity stars is accurate over a wide range in luminosity. These results also provide support for stellar models (at least those for [Fe/H] ≈0.0) that use the scaled HM74 T–τ relation to derive the pressure at T = T_eff.

NGC 6791 is one of the oldest and most metal-rich open clusters known, and for these reasons it has been the subject of numerous investigations over the years (e.g., Garnavich et al. 1994; Origlia et al. 2006; Brogaard et al. 2011). Because its CMD is characterized by tight, well-defined photometric sequences (e.g., Stetson et al. 2003), and because two of its eclipsing binaries have been subjected to careful analyses (see Brogaard et al. 2012, 2011), NGC 6791 provides an especially powerful probe of the properties of old, super-metal-rich stars. In particular, the MR diagram for the observed binaries provides an important constraint on the helium content of NGC 6791 if the abundances of the metals are obtained from spectroscopic data. Moreover, an estimate of the binary, and hence cluster, distance that is completely independent of CMD considerations may be derived from the luminosities which are implied by their radii and, say, spectroscopically determined values of T_eff. Brogaard et al. (2012) concluded that NGC 6791 has an age of ≈8.3 Gyr if it has [Fe/H] =+0.35 (with the metals in the proportions given by Grevesse & Sauval 1998), Y = 0.30, E(B − V) = 0.14, and (m − M)_V = 13.51.

According to K. Brogaard (2013, private communication), a slightly lower [Fe/H] value (≲ 0.30) seems to be favored by the latest spectroscopic results (in agreement with the earlier findings of Boesgaard et al. 2009). If [Fe/H] =0.30 is adopted, it is a straightforward and relatively quick exercise to iterate between the fits of isochrones for different Y and age to the MR diagram for the binaries known as V18 and V20 (for numerical values of their properties, see Brogaard et al. 2012, their Table 1) and the cluster CMD to obtain the best possible consistency between them. We have opted to use the VJ photometry for NGC 6791 compiled by Brasseur et al. (2010): they collected new J observations, which were calibrated to the 2MASS system (Skrutskie et al. 2006) and then combined with V-band data for the same stars (e.g., Stetson 2005).

Initially, we found that an 8.5 Gyr isochrone for [Fe/H] =0.30 and Y = 0.28 provided quite a good fit to the photometry if E(B − V) = 0.14, which agrees well with the value of 0.133 from the SF11 dust maps, and (m − M)_V = 13.55. However, in a separate (concurrent) study, Casagrande & VandenBerg (2014) were able to obtain a consistent fit of the same isochrone to most of the CMDs that can be constructed for NGC 6791 from publicly available BVIJ and Sloan ugriz photometry if E(B − V) = 0.16 and the equivalent true distance modulus, (m − M)₀ = 13.05, are assumed.⁹ Curiously, the most problematic CMD turned out to be the same [(V − J), V] diagram that we have considered here. If E(B − V) = 0.16 is adopted, as implied by most of the observations considered by Casagrande & VandenBerg, the observed colors and/or the transformations to V − J must be corrected by a combined total of 0.04 mag (which is equivalent to a change of 0.018 mag in E(B − V) since E(V − J) ≈ 2.23 E(B − V); see Casagrande & VandenBerg 2014), in order to obtain a good fit of the isochrone to the Brasseur et al. CMD. (Further work will be needed to identify the cause of the color offset if, indeed, the foreground reddening in the direction of NGC 6791 truly is E(B − V) = 0.16.)

Figure 10 shows that the observed RGB is redder (at V ≲ 15.5) than the isochrone which otherwise does a fine job of reproducing the fainter photometry. This suggests that either the model temperatures along the upper giant branch are too hot, or the adopted color–T_eff relations for low gravities yield V − J colors that are too blue, or both. As Brasseur et al. (2010) did not obtain VJ data for fainter MS stars than those plotted in Figure 10, we are unable to comment on how well the isochrone fits near-IR data for LMS stars. However, the plots provided by Casagrande & VandenBerg (2014, see their Figures 12 and 13) indicate that the same isochrone tends to deviate to the blue of the cluster MS at 2–3 mag below the TO (depending on the selected color index), while providing a comparable fit to the upper MS, TO, and SGB stars as that shown in Figure 10. Presumably, the discrepancies at faint magnitudes are also indicative of errors in the model T_eff scale and/or the color transformations for cool, super-metal-rich stars.

Zoom In Zoom Out Reset image size

The dashed curve in this figure represents an isochrone for 8.0 Gyr, Y = 0.30, and [Fe/H] =+0.35 that has been computed assuming the m/Fe number abundance ratios given by Grevesse & Sauval (1998) instead of those determined by Asplund et al. (2009). As mentioned in the first paragraph of this section, Brogaard et al. (2012) derived a slightly higher age (8.3 Gyr) on the assumption of exactly the same chemical abundances. This is, however, the expected consequence of their adoption of a smaller value of (m − M)_V by 0.04 mag, which is easily within the 1σ uncertainty associated with the cluster distance modulus (see the Brogaard et al. paper for a discussion of this issue). Figure 10 shows that, with just a small difference in age, and minor changes to the adopted values of [Fe/H] and Y, isochrones based on either of the Grevesse & Sauval (1998) or Asplund et al. (2009) metals mixtures provide equally good fits to the CMD of NGC 6791 (as well as its binaries, see below). This reinforces the conclusions reached by Brogaard et al. from a similar analysis that such comparisons between theory and observations are not able to provide a clear preference for either solar abundance mixture, due in part to the compensating effects of the respective solar calibration.

The masses and radii that were determined for the components of the binaries V18 and V20 in NGC 6791 by Brogaard et al. (2012) are shown in Figure 11, along with the predicted MR relations from four different isochrones that assume scaled Asplund et al. (2009) metal abundances and one isochrone which assumes the same m/Fe ratios that were derived for the Sun by Grevesse & Sauval (1998). The solid and dashed curves represent the same isochrones that were plotted in the previous figure, and both provide reasonably good fits to the data. Indeed, very similar plots are given by Brogaard et al., who showed that it is only when 3σ error boxes are plotted that the observations can be intersected by a single isochrone. At this stage, it is not known whether the apparent discrepancies are due more to deficiencies in the models that have been compared with the observed masses and radii or to errors in the derived properties of the binaries. It would certainly be worthwhile to collect and analyze more observations of them and to add to the sample of completely eclipsing binaries that have been discovered to date in NGC 6791.

Zoom In Zoom Out Reset image size

Be that as it may, Figure 11 shows how the MR relation that is represented by the solid curve would be altered by, in turn, a 0.5 Gyr increase in age (the dot–dashed curve), a 0.05 dex reduction in [Fe/H] (the dotted curve), and a change in Y by +0.01 (the long-dashed curve). According to these results, we would have obtained a closer match of the dashed locus to the solid curve, with an equally good fit to the cluster CMD, if the former assumed [Fe/H] =+0.34 (instead of +0.35) or a larger helium abundance by δ Y ≈ 0.002. (The effects of variations in [α/Fe] have not been considered because, to within the uncertainties, the abundances of the α elements appear to be consistent with scaled solar values; see Brogaard et al. 2011.)

To conclude: aside from a possible zero-point error in the near-infrared photometry that we have used (Brasseur et al. 2010) or in the Casagrande & VandenBerg (2014) transformations to the J-band, our isochrones for [Fe/H] ∼0.3 are able to reproduce the observed [(V − J), V)] diagram of NGC 6791 quite well in a systematic sense—at least at V − J ≲ 2.2, which corresponds to T_eff ≳ 4600 K. (Especially encouraging comparisons between the same isochrones and many other CMDs of this open cluster, derived from available BVI and ugriz photometry, are provided by Casagrande & VandenBerg.) We have also demonstrated that it is easy to use the models presented in this study to iterate on the age and chemical abundance parameters until a consistent fit is found to both an observed CMD and the MR relation that can be obtained from observations of detached, eclipsing binaries that belong to the same cluster.

VandenBerg et al. (2010) have already shown that current Victoria–Regina stellar models satisfy the constraints provided by subdwarfs in the solar neighborhood that have well-determined M_V values from Hipparcos. In fact, good consistency between theory and observations is obtained on several different color–magnitude planes, particularly those involving red or near-infrared colors (also see Brasseur et al. 2010), or on the (log T_eff, M_V)-diagram if the temperatures of the Population II dwarfs are obtained from Casagrande et al. (2010). The T_eff scale derived by the latter is ∼150 K hotter than the one by Alonso et al. (1999, 1996), which was widely adopted during the last decade, though it agrees well with the hot temperature scale first proposed by King (1993), and three years later by Gratton et al. (1996). It may be recalled that the [Fe/H] values determined for GCs by Carretta & Gratton (1997) are based, in part, on a hot T_eff scale.

It is of some interest to revisit the work by R. G. Gratton and collaborators in the late 1990s, as their T_eff and [Fe/H] estimates for local subdwarfs are in remarkable agreement with the predictions of present-day isochrones. This is shown in Figure 12, which plots (in the bottom panel) the absolute visual magnitudes for the 10 subdwarfs with [Fe/H] ≲ −1.0 that have the smallest values of σ_π/π (van Leeuwen 2007), where π represents the trigonometric parallax, as a function of their effective temperatures. The sources of the T_eff (and [Fe/H]) determinations are Gratton et al. (1996), Gratton et al. (1997); Clementini et al. (1999), and R. G. Gratton (2001, private communication, as reported by Bergbusch & VandenBerg 2001, see their Section 4.1). If the models provided a perfect match to the observed stars, each of the subdwarfs would sit on the isochrone for its measured [Fe/H] value and the temperature implied by that isochrone would be identical to the spectroscopic estimate of T_eff. (Of course, even the best metallicity and T_eff determinations are uncertain by ∼ ± 0.2 dex and ∼ ± 70 K, respectively. Furthermore, the ages of the subdwarfs could well be higher or lower than 12 Gyr—though most of them are sufficiently faint that the effect of the age uncertainty will have negligible consequences for our comparisons with the observations.)

Zoom In Zoom Out Reset image size

The middle panel plots, as a function of the temperatures of the subdwarfs, the differences between the [Fe/H] values that were determined spectroscopically and those inferred from the isochrones that match the subdwarf locations in the (log T_eff, M_V) diagram. For the sample of 10 stars, the mean offset is only 0.04 dex, in the sense that the observed iron abundances are just slightly less than the values deduced from the isochrones, with a standard deviation of 0.29 dex. Interestingly, the differences between the observed ("Obs") and isochrone ("Iso") metallicities tend to be ≳ 0.0 for stars that have [Fe/H] ⩾−1.5 (those represented by open circles) whereas more metal-deficient stars, which are plotted as filled circles, all have "Obs − Iso" values <0.0. However, there is no obvious variation of δ [Fe/H] with temperature for either group of stars, which suggests that the models predict the correct lower-MS slopes.

One may alternatively interpolate in the isochrones to determine how much of an adjustment to the temperature of each subdwarf, at its observed M_V, would be required to locate it on the isochrone that has the same [Fe/H] as the subdwarf. The differences in T_eff so derived are plotted in the upper panel of Figure 12. Not surprisingly (because the abundance implied by a given line strength depends directly on the adopted temperature), stars with [Fe/H] ⩾−1.5 tend to have "Obs − Iso" values of δ T_eff ≳ 0.0, while the opposite is found for the most metal-poor stars. As in the middle panel, the level of agreement is surprisingly good: the mean offset and standard deviation are only −11 K and 65 K, respectively. Though the sample of stars is small, the models appear to fit the observations equally well over the entire temperature range encompassed by the stars. (Considering just the 10 subdwarfs in our sample, the temperatures and [Fe/H] values determined by R. G. Gratton and collaborators are, in the mean, 17 K and 0.09 dex higher, respectively, than the values tabulated by Casagrande et al. 2010.)

The same isochrones, when plotted on the [(V − I)₀, M_V] and [(V − K)₀, M_V] diagrams, provide equally satisfactory fits to the same subdwarfs—as shown in Figure 13. The differences in the predicted colors, which are based on the MARCS color–T_eff relations (Casagrande & VandenBerg 2014), are obviously in excellent agreement with those observed (from Casagrande et al. 2010) since, on both color planes, <δ (color) > = 0.00, with relatively little scatter about the horizontal dashed line. Moreover, the models appear to fit the brighter, bluer stars just as well as the reddest, faintest ones. It is important to appreciate that relatively high temperatures must be assumed for the subdwarfs in order for the MARCS transformations to yield the observed colors. Most broad-band colors (especially V − I and V − K) are much more dependent on T_eff than on [Fe/H] (or on gravity, which will, in any case, be close to log  g = 4.5 for the Population II dwarfs). That is, our isochrones are able to provide good fits to the observations only because they predict the particular T_eff scale that yields the observed colors when derived from the MARCS color–T_eff relations. Similar success would not have been obtained had the latter predicted much redder or bluer colors at the same T_eff or if we had adopted a significantly cooler empirical T_eff scale (e.g., Alonso et al. 1999, 1996).

Zoom In Zoom Out Reset image size

Despite the indications from the spectroscopic results described above, the recent calibration of the IRFM by Casagrande et al. (2010), the predictions of our stellar evolutionary models, and the color–temperature relations implied by the latest MARCS model atmospheres in support of a hot T_eff scale, it is important to remember that current 1D model atmospheres play a central role in each of these avenues of research. Indeed, as discussed by Magic et al. (2013), the very different temperature structures produced by 3D model atmospheres, particularly at low Z, are bound to impact T_eff and [m/H] determinations, as well as color transformations and the BCs employed by stellar models. Indeed, the importance of advancing our understanding of model atmospheres, which provide the interface between stellar interior models and observed stars and stellar populations, can hardly be understated.

In their extensive survey of GC ages, VandenBerg et al. (2013) found that isochrones generally provided reasonably good fits to the Hubble Space Telescope Advanced Camera for Surveys (ACS) photometry obtained by Sarajedini et al. (2007) when the cluster distances were determined from fits of zero-age horizontal-branch (ZAHB) loci to the observed horizontal branch (HB) stars. All of the models used in that investigation assumed the solar metal abundances given by Grevesse & Sauval (1998), with suitable enhancements to the abundances of α-elements and then scaled to the [Fe/H] values derived by Carretta et al. (2009a). As we have not yet computed ZAHBs for the chemical mixtures assumed in this study, we are unable to follow exactly the same procedure here in order to ascertain, in particular, how the inferred distance moduli will differ from those found by VandenBerg et al. However, the predicted ZAHB luminosities, at the same [Fe/H], are likely to be quite similar because the main difference in the solar mixtures given by Grevesse & Sauval and Asplund et al. (2009) are the abundances of CNO, which mainly affect the color of the HB.¹⁰

Because TO luminosity versus age relations depend sensitively on the absolute abundance of oxygen (see V12), and because our models assume a smaller value of [O/Fe] (by 0.1 dex) as well as a lower solar abundance of oxygen, it can be expected that we will obtain higher ages for metal-poor clusters than those derived by VandenBerg et al. (2013; if all other variables are kept constant). However, this is a moot point for the present discussion. Our main motivation for examining the CMDs of a few GCs is to check how well our models are able to reproduce the observed MS and RGB morphologies. To partially compensate for the expected effects of the different abundances of oxygen noted above (and of other metals), we have arbitrarily assumed slightly larger distance moduli (by ≲ 0.05 mag) and ages (by 0.25 Gyr) than the values derived by VandenBerg et al., and then matched the predicted and observed TOs. To accomplish this, it was necessary to apply a small blueward shift to the isochrones (by ≲ 0.02 mag) after the observed colors had been dereddened.¹¹ The result of this exercise is shown in Figure 14 for the GCs 47 Tuc, M3, M5, and M92. As in the VandenBerg et al. study, the [Fe/H] values derived by Carretta et al. (2009a) have been assumed.

Zoom In Zoom Out Reset image size

Except at M_F606W ≳ 8, where the solid curves deviate to the blue side of the observed lower-MS stars, the isochrones do quite a good job of matching the MSs of the GCs over the entire range in [Fe/H] sampled by them. The biggest differences between theory and observations occur along the lower RGB, where the models are too red. However, the tendency of photometric scatter due to blending to be preferentially blueward on the giant branch may explain some fraction of such offsets (see Bergbusch & Stetson 2009). Curiously, the discrepancies resemble the effect on the location of the RGB of varying the helium content: as shown in Figure 7, increasing Y causes a larger temperature shift at the base of the giant branch than near the tip. On the other hand, it is possible that our treatment of the atmospheric BC is responsible for the apparent difficulties (recall Figure 4), or perhaps they signal some problems with the color–T_eff relations that we have used or our treatment of convection. As noted in the introductory remarks given at the beginning of Section 5, predicted temperatures and colors are subject to many uncertainties, and it should not be a surprise to find some discrepancies between isochrones and observed CMDs.

To corroborate this point, we have generated isochrones for the case represented in Figure 4 by the dot–dashed curve. That is, a full set of evolutionary tracks has been computed for [Fe/H] =−2.40, Y = 0.25, and [α/Fe] =0.4 in which the surface BCs have been derived from the properties of MARCS model atmospheres at τ = 100, with the small increase in the pressure at that point implied by the corresponding Standard Solar Model. As shown in Figure 15, the 13 Gyr isochrone derived from these tracks, unlike the one plotted in the bottom, right-hand panel of the previous figure, provides a good fit to the lower RGB stars of M92, but not those at higher luminosities. (Granted, the predicted TO is slightly too blue, but an improved fit to the TO could be obtained, without affecting the location of the lower RGB, simply by assuming a somewhat higher oxygen abundance.) In this example, the discrepancies along the upper giant branch could be telling us that, e.g., our treatment of convection or the adopted color–T_eff relations are inadequate. We could force the models to provide an essentially perfect match to the data (by, for instance, suitable adjustments of the color transformations or the atmospheric BCs), but the assumed distance and chemical abundances of M92 may not be correct. Although isochrones may need to be "calibrated" for some investigations, not doing so enables one to retain the predictive power of stellar models. In fact, it is remarkable that current stellar models perform as well as they (appear to) do.

Zoom In Zoom Out Reset image size

As presented to the user, the evolutionary sequences are contained in EEP files which have been processed in such a way that track points with the same model number are equivalent in every track in every grid. Two caveats apply to this prescription. First, in grids that extend to masses lower than $0.4 \,{\cal M_{\odot }}$, tracks with masses ${\le } 0.4 \,{\cal M_{\odot }}$ are listed for equally logarithmically spaced ages from the ZAMS point up to a maximum age of ≈30 Gyr. Second, for those grids that contain tracks in which core contraction manifests itself after core hydrogen exhaustion, the MSTO point EEP (MSTO) becomes degenerate with the blue hook EEP (BLHK) for those lower mass tracks in which the blue hook is not present. We do not explicitly list these degenerate points: their presence (discussed below) is indicated by listing the primary EEPS in the header lines for each track.

The canonical EEP file names (with the extension *.eep) provide all the abundance information for the tracks contained within them: each one begins with a five character prefix that terminates with the underscore character followed by three abundance specifications, e.g., a0zz_p4y29m18.eep. In this example, a0 indicates the solar metals mixture (Asplund et al., 2009) and zz specifies the entire group of α-elements: decoding the rest of the name from left to right, p4 implies [α/Fe]= +0.4, y29 implies Y = 0.29, and m18 implies [Fe/H]= −1.8. A grid of tracks interpolated to [α/Fe]= +0.2, Y = 0.273, and [Fe/H]= −0.75 would have a0zz_p2y273m075.xeep as its file name, where the extension .xeep distinguishes it from the canonical grids. Had the interpolation been to [α/Fe] = −0.2, Y = 0.273, and [Fe/H] = +0.25, the file name would have been a0zz_m2y273p025.xeep—that is, the signs of the α-element and iron abundances are denoted by either p (+ve) or m (−ve).

For future reference with grids in which individual elements may be enhanced differently with respect to some basic [α/Fe] abundance ratio, such grids will have file names like a4xO_p1y25p02. The prefix a4xO decodes as "a basic [α/Fe] = +0.4 mixture with an extra degree of enhancement of the element oxygen." Decoding the rest of the name, p1 means that oxygen has been incrementally enhanced by +0.1 dex (above the amount in the basic [α/Fe] mixture) so that [O/Fe] = +0.5, and y25p02 means Y = 0.25 and [Fe/H] = +0.2. When the symbol for an element consists of a single letter (like C, N, or O), that letter appears just before the underscore, and x is used as a place-holder; otherwise, the third and fourth characters of the file name give the two-letter symbol of the metal (e.g., Ne, Mg) in question.

The contents of a *.eep file are illustrated in Figure 16. The header lines at the beginning of the file are reasonably straightforward to interpret (note, in particular, that the assumed [m/Fe] values are explicitly given for the main metals of interest), but the header lines for individual tracks require some explanation. The columns labeled Match, D(age), and D(log Teff) are redundant: they list information about how the models for the MS and SGB phases, which were obtained by solving the Lagrangian form of the stellar structure equations, were matched (at the base of the RGB) to the models for the subsequent evolution to the tip of the giant branch. As described in detail by VandenBerg (1992), a non-Lagrangian technique like the one developed by Eggleton (1971) was used to follow RGB evolution very efficiently. The indicated age and T_eff offsets were applied to the original track files for the RGB phase in order to obtain continuity with the Lagrangian models.

Zoom In Zoom Out Reset image size

The column labeled Zage lists the ZAMS age as we defined it in V12. Six evolutionary phase points are listed under Primary EEPs. The default setting for each EEP-point is 0, which should be interpreted to mean (except for the ZAMS) that that particular evolutionary phase does not occur, or has not been reached, in that particular track. The ZAMS EEP-point is listed as model 1 when the track is included in the full isochrone interpolation scheme. When listed as 0, it signals that an interpolation to the isochrone age is to be made directly within that track. (Generally this applies to tracks with masses ${\le } 0.4 \,{\cal M_{\odot }}$. Points on the isochrone between those corresponding to these track masses are obtained by spline interpolation.) If the track evolves sufficiently, the MSTO is listed at model 801, and if a blue hook occurs it is listed at model 921. In the absence of a BLHK EEP, as is the case for the $0.9 \,{\cal M_{\odot }}$ track in the grid shown, the base of the red-giant branch (BRGB) occurs at model 1421, the evolutionary pause on the giant branch (GBPS) at model 1621, and the tip of the giant branch (GBTP) at model 1921. When the BLHK is non-zero, as is the case for the tracks with masses ${>}0.9 \,{\cal M_{\odot }}$, the BRGB occurs at model 1541, the GBPS at 1741, and the RGBT at 2041.

In this example, tracks with masses ${>}0.9 \,{\cal M_{\odot }}$ have BLHK EEPs while those with masses ${\le } 0.9 \,{\cal M_{\odot }}$ do not. (A "nascent" BLHK EEP may be identified in some tracks to relax the interpolation scheme through the transition from lower mass tracks with radiative cores at central H exhaustion to those higher mass tracks with fully developed convective cores at the end of the MS phase; see Bergbusch & VandenBerg 2001.) Consequently, 120 degenerate EEP points would be inserted between the MSTO and BRGB for the lower mass tracks that evolve at least as far as the BRGB so that, for example, the BRGB in the $0.9 \,{\cal M_{\odot }}$ track would change from model 1421 to model 1541. The same approach works when interpolating grids of tracks to some set of target abundances.

The first four columns in the listing for each track give the model number, log L/L_☉, log T_eff, and the age in Gyr. Columns 5 and 6 are reserved for the surface helium and [Fe/H] abundances (not implemented for the example shown in Figure 16, while Columns 7 and 8 list the derivatives of luminosity and of effective temperature with respect to time. The latter are needed for the calculation of LFs and IPFs.

The interpolation of isochrones, LFs, and/or IPFs is made easy by the (FORTRAN) programs PBISO and PBIPF, which are improved versions of the MKISO and MKIPF codes, respectively, that were presented in previous publications (see V12 and references therein). We have now developed a new program, PBMIX, that can be used to interpolate within the canonical grids to produce a new grid of tracks at some set of the (up to) three abundance parameters—either ([α/Fe], Y, [Fe/H]) or ([m_i/Fe], Y, [Fe/H]) that is signaled by the file name prefix. PBMIX makes use of a parameter file, PBMIX.PAR, that resides in the working directory; it contains a listing of the track masses employed in constructing the grid, the abundance ranges spanned by the grids, and a listing of all the EEP files encompassed by the abundance ranges specified. If PBMIX is run in a directory that doesn't contain a PBMIX.PAR file, it will prompt you for all the requisite information and construct one automatically.

A sample file is shown in Figure 17. The first line lists the prefix for the file names associated with the grids contained in the working directory. (In this example, the grids have a basic [α/Fe] = +0.4 mixture with additional enhancements to [O/Fe].) The next line begins with the number of tracks that may be associated with each grid, followed by the mass values—since the masses are read in via a list-directed read statement, they need only be separated by blanks and can be spread over several lines (two lines in this example). The (fourth) line following the list of masses gives the number of and canonical values for the abundance ratio [O/Fe]. (If the file prefix were a0zz, this line would list the [α/Fe] values instead.) The next (fifth) line does the same thing for the helium abundances, as does the (sixth) line for [Fe/H]. Subsequently, the 24 canonical file names derived from the file prefix and the tabulated abundances—each with a unique set of [O/Fe], Y, and [Fe/H]—are listed.

Zoom In Zoom Out Reset image size

Running PBMIX is very simple. When PBMIX.PAR already exists, the program simply prompts the user for the target [α/Fe] value or the Δ[m_i/Fe] increment, and the target Y and [Fe/H] values, then it creates the default output file name based on those targets and writes the interpolated tracks to that file. As in the case of the canonical grids, isochrones may be generated for the interpolated grids for any age in the range from ∼5 Gyr to ∼15 Gyr by executing PBISO. The auxiliary code PBIPF may then be used to provide magnitudes and colors in the Johnson–Cousins, 2MASS, Sloan, or Hubble Space Telescope ACS or WFC3 photometric systems. (The color transformation tables that are needed to accomplish this must be generated using the computer programs and data provided by Casagrande & VandenBerg 2014.) The same code provides the option of producing LFs or IPFs. A brief, but quite thorough, description of how to use each code is provided in a manual that can also be retrieved from the CANFAR Web site.

To demonstrate that the interpolation errors are small, even though linear interpolation is employed for all three chemical abundance parameters, we have computed a set of evolutionary tracks for values of [Fe/H] =−0.1, [α/Fe] =+0.2, and Y = 0.27 that are midway between their respective grid values. These particular abundances were chosen because they lie along the [α/Fe] versus [Fe/H] relationship that has been derived for stars in the Galactic bulge (e.g., Ryde et al. 2010, their Figure 3), where the "knee" in that relation, which is believed to represent the point at which Type Ia supernovae began to contribute to the chemical evolution of the Bulge, occurs at an unusually high [Fe/H] value (≈ − 0.3). (To produce this grid, the necessary opacity data were generated in the usual way; see Section 3.) It can be expected that this case will provide quite a severe test of the interpolation scheme because (1) the difference between a line that connects any two of the points that are fitted by a parabolic curve, and that parabola, will be maximal at approximately the mid-point, and (2) the effects of opacity on mass–luminosity and age–mass relations are strongest at the highest metal abundances.

However, as shown in Figure 18, isochrones derived from this set of models differ only slightly from those for the same age that are obtained by interpolation in the grids which are being provided for general use via this paper. At a value of log  T_eff that is 0.01 dex cooler than the TO temperature, the SGBs of the interpolated isochrones are ≲ 0.04 mag brighter than the SGBs of isochrones that have been derived from evolutionary tracks computed for the same values of Y, [α/Fe], and [Fe/H]. This is not negligible, but as noted above, the interpolation errors are expected to be larger for this case than for most other choices of the abundance parameters that are representative of the observed abundances in stars. For instance, the transition from [α/Fe] ≈+0.4 to 0.0 typically begins at [Fe/H] <−1.0 in dwarf galaxies (Venn et al. 2004), and the second test case that we have considered (for Y = 0.26, [α/Fe] =+0.2, and [Fe/H] =−0.9) shows that the computed and interpolated isochrones which are applicable to such systems will be indistinguishable (see Figure 18).

Zoom In Zoom Out Reset image size

We would have preferred to compute grids of evolutionary sequences for δ[α/Fe] =0.2, but the model atmospheres that are needed to provide the surface BCs for the lowest mass models are not currently available for such a fine spacing of this quantity. It turns out, in fact, that most of the interpolation errors at high metallicities occur because a grid spacing of 0.4 dex range in [α/Fe] is too large. We reached this conclusion after comparing computed and interpolated isochrones for a third case; specifically, [α/Fe] =0.0 (one of the grid values), Y = 0.27, and [Fe/H] =+0.30, where the latter quantities are at the mid-points of the grid values, Y = 0.25, 0.29 and [Fe/H] =+0.2, +0.4, respectively. To avoid making Figure 8 too complicated, we elected not to plot these two isochrones, but we found that they superimpose each other so well that the differences between them are barely discernible.

V12 (see their Figures 4–6) had previously demonstrated that the errors associated with Y, [Fe/H] interpolations are negligible if the grid spacings of these variables are δY = 0.04 and δ[Fe/H] =0.2 dex. Although they used three-point interpolations in the tests that they conducted, we have found that linear interpolations work nearly as well (even at high metal abundances). Thus, there appears to be only a small regime of parameter space ([α/Fe] roughly halfway between the grid values, but only at relatively high [Fe/H] values) where minor interpolation errors can be expected. Elsewhere, such errors are of no consequence whatsoever.