MESA ISOCHRONES AND STELLAR TRACKS (MIST). I. SOLAR-SCALED MODELS

The accurate determination of solar abundances has been an ongoing effort for decades. Within the past decade, there has been a systematic downward revision of the solar metallicity from ${Z}_{\odot }\sim 0.02$ (e.g., Anders & Grevesse 1989) to ${Z}_{\odot }\lesssim 0.015$ (e.g., Asplund et al. 2009; Lodders et al. 2009; Caffau et al. 2011). Recently, the abundances of C, N, O, and Ne have experienced dramatic revisions as a collective result of improved atomic and molecular linelists and the introduction of 3D non-local thermodynamic equilibrium (NLTE) hydrodynamical modeling techniques (see Asplund et al. 2009 for a review). In this paper, we adopt the protosolar abundances recommended by Asplund et al. (2009) as the reference scale for all metallicities, unless noted otherwise. In other words, [Z/H] is computed with respect to $Z\;=\;{Z}_{\odot ,\mathrm{protosolar}}\;=\;0.0142$, not $Z\;=\;{Z}_{\odot ,\mathrm{photosphere}}\;=\;0.0134$. The difference between the two is a consequence of diffusion of heavy elements out of the photosphere over time. We emphasize that this is not an attempt to redefine ${Z}_{\odot ,\mathrm{photosphere}}$—the current photospheric abundances are determined by spectroscopy—but rather to clarify what a "solar metallicity model" entails.

To compute X and Y for an arbitrary Z, we use the following formulae:

$$\begin{eqnarray}&&{Y}_{{\rm{p}}}=0.249\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&Y={Y}_{{\rm{p}}}+\left(\displaystyle \frac{{Y}_{\odot ,\mathrm{protosolar}}-{Y}_{{\rm{p}}}}{{Z}_{\odot ,\mathrm{protosolar}}}\right)Z\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&X=1-Y-Z.\end{eqnarray} \tag{ 3 }$$

This approach adopts a primordial helium abundance ${Y}_{{\rm{p}}}\;=\;0.249$ (Planck Collaboration et al. 2015) determined by combining the Planck power spectra, Planck lensing, and a number of "external data" such as baryonic acoustic oscillations. In the above equations, we assume a linear enrichment law to the protosolar helium abundance, ${Y}_{\odot ,\mathrm{protosolar}}\;=\;0.2703$ (Asplund et al. 2009), such that ${\rm{\Delta }}Y/{\rm{\Delta }}Z\;=\;1.5$. Once Y is computed for a desired value of Z, X is trivial to compute. We assume the isotopic ratios ${}^{2}{\rm{H}}{/}^{1}{\rm{H}}\;=\;2\times {10}^{-5}$ and ${}^{3}\mathrm{He}{/}^{4}\mathrm{He}\;=\;1.66\times {10}^{-4}$ from Asplund et al. (2009).

3.2.1. Equation of State (EOS)

3.2.2. Opacities

3.2.3. Nuclear Network

The pressure and temperature in the outermost cell of a stellar model calculation must be specified as a set of boundary conditions in addition to the trivial boundary conditions at the center of the star. There are a multitude of options that range from simple analytic approximations to tables based on full atmospheric structure models.

The simplest choice, simple_photosphere, uses the Eddington $T(\tau )$ relation to obtain T:

$$\begin{eqnarray}&&{T}^{4}(\tau )\;=\;\displaystyle \frac{3}{4}{T}_{\mathrm{eff}}^{4}\left(\tau +\displaystyle \frac{2}{3}\right),\end{eqnarray} \tag{ 4 }$$

where ${T}_{\mathrm{eff}}$ is calculated directly from the MESA interior model. Similarly, P is computed as follows:

$$\begin{eqnarray}&&P\;=\;\displaystyle \frac{\tau g}{\kappa }\left[1+{P}_{0}\displaystyle \frac{\kappa }{\tau }\displaystyle \frac{L}{M}\displaystyle \frac{1}{6\pi {cG}}\right].\end{eqnarray} \tag{ 5 }$$

The second term in the square brackets accounts for the nonzero radiation pressure (e.g., Cox 1968), which can be significant in high-mass stars. P₀ is a dimensionless factor of order unity used to scale up the radiation pressure in order to help convergence in massive stars and post-asymptotic giant branch (post-AGB) stars radiating close to or at super-Eddington luminosities. We adopt ${P}_{0}\;=\;2$ for this work.

In most cases, the simple_photosphere option is a poor choice as there is no guarantee that κ and P from the interior model are consistent according to Equation (5), assumed to be the correct relation at the stellar surface. For this work, we adopt realistic model atmospheres as the outer boundary conditions for most locations in the Hertzsprung–Russell (HR) diagram. We have computed a new grid of 1D plane-parallel atmosphere models specifically for this project using the ATLAS12 code (Kurucz 1970, 1993). The atmospheres are computed to a Rosseland optical depth of 10³ with ${\alpha }_{\mathrm{MLT}}\;=\;1.25$ following the implementation of convection as outlined in Castelli et al. (1997).¹¹ We employed the latest atomic and molecular line lists from R. Kurucz, including molecules important for cool stars such as TiO and H₂O. Individual models are calculated for $\mathrm{log}(g)\;=\;0$ to 5 and ${T}_{\mathrm{eff}}\;=\;2500$ K to 50,000 K for $[{\rm{Z}}/{\rm{H}}]\;=\;-7.0$ to $+0.5$ on the Asplund et al. (2009) abundance scale. Beyond these limits the tables have been smoothly extrapolated to encompass all possible locations of the stellar tracks. This is a satisfactory solution since the few phases that fall into these extrapolated regimes (e.g., post-AGB) are typically very short-lived.

With model atmosphere tables in hand, one is left to choose where (in terms of Rosseland depth) to use the tables as boundary conditions for the models. The standard convention, which we adopt for most stars, is the photosphere, i.e., where $T\;=\;{T}_{\mathrm{eff}}$. However, for cooler dwarfs a more sensible choice is to set the boundary condition deeper in the atmosphere, i.e., $\tau \;=\;100$. This option will result in more realistic models for cool low-mass stars whose atmospheres are heavily influenced by molecules that are not included in the MESA interior model calculations. This issue is less critical for the cool giants because the structure of these stars is overall less sensitive to the boundary condition (i.e., the pressure at the photosphere for giants is much closer to zero than for dwarfs). We refer the reader to Section 6.2 for additional discussion on this topic.

For our grid of models, the tau_100_tables option is used for $0.1\mbox{--}0.3\;{M}_{\odot }$, photosphere_tables is used for $0.6\mbox{--}10\;{M}_{\odot }$, and simple_photosphere is used for $16\mbox{--}300\;{M}_{\odot }$. To facilitate a smooth transition between different regimes, we run both tau_100_tables and photosphere_tables for $0.3\mbox{--}0.6\;{M}_{\odot }$ and photosphere_tables and simple_photosphere for $10\mbox{--}16\;{M}_{\odot }$. The tracks in this transition region are then blended (see Section 5.2 for more details). At the highest masses, simple_photosphere is a sufficient approximation due to the flattening of opacity as a function of temperature for ${T}_{\mathrm{eff}}\gtrsim {10}^{4}$ K.

Microscopic diffusion and gravitational settling of elements are essential ingredients in stellar evolution models of low-mass stars, leading to a modification to the surface abundances and main-sequence (MS) lifetimes, as well as a shift in the evolutionary tracks toward lower luminosities and temperatures in the H-R diagram (e.g., Michaud et al. 1984; Morel & Baglin 1999; Salaris et al. 2000; Chaboyer et al. 2001; Bressan et al. 2012). Calculations of diffusion and gravitational settling are implemented in MESA following Thoul et al. (1994). All species are categorized into one of five "classes" according to their atomic masses, each of which has a representative member whose properties are used to estimate the diffusion velocities. MESA's default set of representative members for the five classes are ¹H, ³He, ⁴He, ¹⁶O, and ⁵⁶Fe. Atomic diffusion coefficients are calculated according to Paquette et al. (1986): the representative ionic charge for each class is estimated as a function of T, ρ, and free electrons per nucleon, while the diffusion velocity of the representative member is adopted. The diffusion equation is then solved using the total mass fraction within each class.

However, the inclusion of microscopic diffusion alone cannot reproduce observations of surface abundances in the Hyades open cluster and OB associations including the Orion association (e.g., Cunha & Lambert 1994; Varenne & Monier 1999; Daflon et al. 2001). Models with diffusion predict an overdepletion of helium and metals in the outer envelopes of stars with ${M}_{{\rm{i}}}\gt 1.4\;{M}_{\odot }$, a problem that appears to worsen with increasing mass due to a disappearing outer convection zone and a concomitant steepening of the temperature and pressure gradients (Morel & Thévenin 2002). The solution to this problem requires additional forces to counteract gravity. Radiative levitation (Vauclair 1983; Hu et al. 2011) can help to reduce the gravitational settling of highly charged elements, such as iron, via radiation pressure. However, it is thought to be mostly important in hot, luminous massive MS stars or helium-burning stars (e.g., hot subdwarf stars; Fontaine et al. 2008), and to have only a modest effect for solar-type stars (e.g., Alecian et al. 1993; Turcotte et al. 1998). We employ radiation turbulence (Morel & Thévenin 2002) to reduce the efficiency of diffusion in hot stars, though there exist other explanations for the observed surface abundances, including turbulent mixing due to differential rotation (e.g., Richard et al. 1996). We adopt the radiative diffusivity parameter ${D}_{\mathrm{RT}}\;=\;1$, which relates the strength of radiative diffusivity (the deposition of photon momentum into the fluid, resulting in radiative mixing) to the kinematic radiative viscosity.

Since the effects of elemental diffusion are most significant in the absence of more efficient mixing processes, such as convection, diffusion is neglected for fully convective MS stars in the MIST models. Additionally, diffusion is expected to play a reduced role in massive stars and during some post-MS phases (which are also associated with large convective envelopes) for which the evolutionary timescales are comparable to or much shorter than the timescale for diffusion and gravitational settling (e.g., Turcotte et al. 1998). Thus, the effects of microscopic diffusion are considered only for the MS evolution. However, it may have a notable impact on both the atmospheres and interiors of cooling WDs by modifying the surface abundances and lengthening cooling times through the release of gravitational energy. The Thoul et al. (1994) formalism, which assumes isolated interactions between two particles at a time, breaks down in the regime of strongly coupled plasmas. This is a particularly relevant issue for the interiors of cooling WDs, and there is ongoing effort in MESA to update the diffusion implementation to account for this. The inclusion of diffusion during the post-MS evolution, especially the WD cooling phase, is one of the priorities for Paper II.

To summarize, the implementation of diffusion is limited to MS stars above the fully convective limit for which it is most effective in terms of both the relevant timescales and the relative importance compared to other mixing processes.

The effect of rotation on stellar models has been studied for decades (e.g., Strittmatter 1969; Fricke & Kippenhahn 1972; Tassoul 1978; Zahn 1983; Pinsonneault 1997; Heger et al. 2000; Maeder & Meynet 2000; Palacios et al. 2003; Talon & Charbonnel 2005; Denissenkov & Pinsonneault 2007; Hunter et al. 2007; Chieffi & Limongi 2013; Cantiello et al. 2014), but it remains as one of the most challenging and uncertain problems in stellar astrophysics. Rotation is particularly important for massive stars, as rotationally induced instabilities, combined with the non-negligible effects of radiation pressure, can significantly alter their evolution (e.g., Heger et al. 2000; Meynet & Maeder 2000; Hirschi et al. 2004; Woosley & Heger 2006; Yoon et al. 2006; de Mink et al. 2010; Georgy et al. 2012; Langer 2012; Köhler et al. 2015). Although the overall importance of rotation in models—lifetimes, surface abundances, evolutionary fates, to name a few—has been explored, the details are not yet fully understood.

Rotation is inherently a 3D process, but the so-called "shellular approximation" allows stellar structure equations to be solved in 1D (Kippenhahn & Thomas 1970; Endal & Sofia 1976; Meynet & Maeder 1997; Heger et al. 2000; Paxton et al. 2013). This approximation is valid in the regime where strong anisotropic turbulence arises from differential rotation and smears out both chemical composition and velocity gradients along isobars (Zahn 1992; Meynet & Maeder 1997). As a result, the standard stellar structure equations are simply modified by centrifugal acceleration terms in the presence of rotation. More details on the implementation of rotation in MESA can be found in Paxton et al. (2013).

Our models are available in two varieties, with and without rotation. All rotating models are initialized with solid-body rotation on the zero-age main sequence (ZAMS), which is the standard choice in stellar evolution codes (Pinsonneault et al. 1989; Heger et al. 2000; Eggenberger et al. 2008). As discussed extensively in Heger et al. (2000), pre-main-sequence (PMS) stars achieve rigid rotation due to convection, and once they settle onto ZAMS, they establish close-to-rigid rotation mainly via Eddington–Sweet (ES) circulation and Goldreich–Schubert–Fricke (GSF) instability. However, it is worth noting that there are important exceptions. First, this rigid rotation approximation fails in the solar convection zone as inferred from helioseismology observations (e.g., Brown et al. 1989). Second, current detailed evolutionary models (e.g., Bouvier 2008; Gallet & Bouvier 2013) suggest that low-mass stars ($\lesssim 1.2\;{M}_{\odot }$), particularly those with slow and moderate rotation rates, have strong differential rotation profiles at ZAMS.

Currently, surface magnetic fields are not included in MESA calculations, which can couple to mass loss and give rise to magnetic braking (e.g., Weber & Davis 1967; Mestel 1968; ud-Doula & Owocki 2002; Meynet et al. 2011), a mechanism for winding down surface rotation over time in stars with appreciable convective envelopes (Kraft 1967; see also Section 3.6.5).

Since magnetic braking is not currently modeled in MIST, we do not include rotation for stars with ${M}_{{\rm{i}}}\leqslant 1.2\;{M}_{\odot }$ in order to reproduce the slow rotation rate observed in the Sun and in other low-mass stars. Over the mass range 1.2–$1.8\;{M}_{\odot }$, the rotation rate is gradually ramped up from 0 to the maximum value of ${v}_{\mathrm{ZAMS}}/{v}_{\mathrm{crit}}\;=\;{{\rm{\Omega }}}_{\mathrm{ZAMS}}/{{\rm{\Omega }}}_{\mathrm{crit}}\;=\;0.4$, where ${v}_{\mathrm{crit}}$ and ${{\rm{\Omega }}}_{\mathrm{crit}}$ are critical surface linear and angular velocities, respectively (See Equation (22)). This rotation rate, also adopted in the Geneva models (Ekström et al. 2012),¹² is motivated by both recent observations of young B stars (Huang et al. 2010) and theoretical work on rotation rates in massive stars (Rosen et al. 2012). A comparison with observed rotation rates of both MS and evolved stars in the mass range $1.2\mbox{--}1.5\;{M}_{\odot }$ (Wolff & Simon 1997; Canto Martins et al. 2011) reveals that our ramping scheme produces velocities that are reasonable ($\sim 10\mbox{--}25\;\mathrm{km}\;{{\rm{s}}}^{-1}$ during the main sequence for $1.3\mbox{--}1.35\;{M}_{\odot }$) even in the absence of magnetic braking.

Chemical mixing and angular momentum transport due to rotationally induced instabilities are discussed in Section 3.6.4, and rotationally enhanced mass loss is discussed in Section 3.7.3.

3.6.1. Convection

Mixing length theory (MLT), whose modern implementation in stellar evolution codes was pioneered by Böhm-Vitense (1958), describes the convective transport of energy in the stellar interior. There is a vexing yet crucial free parameter of order unity, ${\alpha }_{\mathrm{MLT}},$ that determines how far a fluid parcel travels before it dissolves into the background, ${l}_{\mathrm{MLT}}$, in units of the local pressure scale height, ${H}_{{\rm{P}}}$ (${l}_{\mathrm{MLT}}\;=\;{\alpha }_{\mathrm{MLT}}{H}_{{\rm{P}}}$). In other words, it parameterizes how efficient convection is, because a large ${\alpha }_{\mathrm{MLT}}$ means that the parcel travels a large distance before it deposits its energy into the ambient medium.

Convective mixing of elements is treated as a time-dependent diffusive process with a diffusion coefficient computed within the MLT formalism, which may later be modified by overshoot mixing across convection boundaries (see Section 3.6.2). Convective heat flux is computed by solving the MLT cubic equations to obtain the temperature gradients (e.g., Equation (42) in Henyey et al. 1965). We adopt the modified version of MLT from Henyey et al. (1965) instead of the standard MLT prescription (Cox 1968), as the latter assumes no radiative losses from fluid elements and is therefore applicable only at high optical depth.¹³ In addition to ${\alpha }_{\mathrm{MLT}},$ there are two free parameters, ν and y, which are multiplicative factors to the mixing length velocity and the temperature gradient in the convective element. The latter two parameters are set to 8 and 1/3, respectively (Henyey et al. 1965). This particular framework allows for convective efficiency to vary with the opacity of the convective element, an important effect to take into account in the layers near the stellar surface. The empirical calibration of ${\alpha }_{\mathrm{MLT}}$ is discussed in Section 4.1.

Classically, the location of the convective region is determined using the Schwarzschild criterion, which implies that a region is convectively stable if

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{T}\lt {{\rm{\nabla }}}_{\mathrm{ad}},\end{eqnarray} \tag{ 6 }$$

where ${{\rm{\nabla }}}_{T}$ is the local background temperature gradient (in practice, it is set to the radiative temperature gradient ${{\rm{\nabla }}}_{\mathrm{rad}}$) and ${{\rm{\nabla }}}_{\mathrm{ad}}$ is the adiabatic temperature gradient. Alternatively, the Schwarzschild criterion can be replaced by the Ledoux criterion, which also takes into account the composition gradient, ${{\rm{\nabla }}}_{\mu }$. In this case, a region is convectively stable if

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{T}\lt {{\rm{\nabla }}}_{{\rm{L}}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{{\rm{L}}}\equiv {{\rm{\nabla }}}_{\mathrm{ad}}-\displaystyle \frac{{\chi }_{\mu }}{{\chi }_{T}}{{\rm{\nabla }}}_{\mu }\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\chi }_{\mu }\equiv {\left[\displaystyle \frac{\partial \mathrm{ln}(P)}{\partial \mathrm{ln}(\mu )}\right]}_{\rho ,T}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\chi }_{T}\equiv {\left[\displaystyle \frac{\partial \mathrm{ln}(P)}{\partial \mathrm{ln}(T)}\right]}_{\rho ,\mu },\end{eqnarray} \tag{ 10 }$$

where the thermodynamic derivatives ${\chi }_{\mu }$ and ${\chi }_{T}$ are equal to −1 and 1 for an ideal gas, respectively. We adopt the Ledoux criterion for convection in our models to account for the composition effects in the stellar interiors.

3.6.2. Overshoot Mixing

Rather unsurprisingly, the MLT framework, which relies on a 1D diffusive model in place of a full 3D hydrodynamical treatment, offers an incomplete description of convection. To model the mixing occurring at convective boundaries, also known as overshoot mixing, one must turn to yet another set of parameterizations. Typically, a convective region is extended beyond the fiducial boundary determined by either the Schwarzschild or Ledoux criterion in order to account for the nonzero momentum of the fluid element approaching the edge of the convective zone, as well as its subsequent penetration into the non-convective region (e.g., Böhm 1963; Shaviv & Salpeter 1973; Maeder 1975; Roxburgh 1978; Bressan et al. 1981). This overshoot action leads to enhanced mixing and can account for both the observed properties of AGB and post-AGB stars (Herwig 2000), the observed MS width (e.g., Schaller et al. 1992), and the main-sequence turnoff (MSTO) morphology in clusters such as M67 (e.g., Magic et al. 2010).

There are two prescriptions for convective overshoot available in MESA. The first method, which we call step overshoot, is to simply extend the fiducial convective boundary by a fraction, usually ∼0.2, of the local pressure scale height. This instantaneous treatment is often calibrated to fit the observed MSTO of stellar clusters and associations and is a commonly adopted scheme in many stellar evolution codes (e.g., Demarque et al. 2004; Pietrinferni et al. 2004; Dotter et al. 2008; Brott et al. 2011a; Bressan et al. 2012; Ekström et al. 2012).

The second method, adopted in the present work, was motivated by the plume-like nature of convective elements seen in 2D and 3D radiation hydrodynamic simulations, where coherent downward and upward flows were observed rather than a hierarchy of blob-like eddies (e.g., Freytag et al. 1996). This led to a picture in which the turbulent velocity field decays exponentially away from the fiducial convective boundary and the convective element eventually disintegrates in the overshoot region through a diffusion process. Following the parametrization discussed in Herwig (2000), the resulting diffusion coefficient in the overshoot region is given by

$$\begin{eqnarray}&&{D}_{\mathrm{OV}}\;=\;{D}_{0}\mathrm{exp}\left(\displaystyle \frac{-2z}{{H}_{{\rm{v}}}}\right)\;;\;{H}_{{\rm{v}}}\;=\;{f}_{\mathrm{ov}}{H}_{{\rm{P}}},\end{eqnarray} \tag{ 11 }$$

where ${H}_{{\rm{v}}}$ is the velocity scale height of the overshooting convective elements at the convective boundary, ${f}_{\mathrm{ov}}$ is a free parameter that essentially determines the efficiency of overshoot mixing, ${H}_{{\rm{P}}}$ is the local pressure scale height, and D₀ is the diffusion coefficient in the unstable region "near" the convective boundary (more specifically at a depth of ${f}_{0,\mathrm{ov}}{H}_{{\rm{P}}}$ from the convective boundary). For simplicity, we adopt two sets of (${f}_{\mathrm{ov}}$, ${f}_{0,\mathrm{ov}}$) values, one for the core and another for shell/envelope, irrespective of the type of burning taking place in the overshoot region. For further simplicity, ${f}_{0,\mathrm{ov}}$ is set to $0.5{f}_{\mathrm{ov}}$. The temperature gradient in the overshoot region is assumed to be equal to the radiative gradient as in the step overshoot approach.

Other models in the literature make use of an additional parameterization to both avoid a physically unrealistic size of the overshoot region outside a small convective core and account for the possibility that convective overshoot efficiency is smaller in lower-mass stars (Demarque et al. 2004; Pietrinferni et al. 2004; Dotter et al. 2008; Bressan et al. 2012). In the step overshoot formalism where the size of the convective core is extended by a fraction ${f}_{\mathrm{ov},\mathrm{step}}$ of ${H}_{{\rm{P}}}$, one could end up with a physically unrealistic situation where the size of the overshoot region exceeds the size of the convective core itself. This can occur when the convective core is very small, e.g., for the critical mass around $1.1\mbox{--}1.2\;{M}_{\odot }$ when the CNO cycle begins to dominate over the pp-chain and the hydrogen-burning core becomes convective instead of radiative. Since the convective core boundary is not far from the center and ${H}_{{\rm{P}}}$ formally diverges as $r\to 0$, the size of the overshoot region, ${f}_{\mathrm{ov},\mathrm{step}}{H}_{{\rm{P}}}$, also diverges. Thus, to avoid the excessive growth of the convective core for low-mass MS stars, the common solution is to gradually ramp up the overshoot efficiency from ${M}_{1}\sim 1\;{M}_{\odot }$ to ${M}_{2}\sim 1.7\;{M}_{\odot }$, with no convective overshoot below M₁ and maximum efficiency above M₂. These boundary masses vary with metallicity due to opacity effects (Demarque et al. 2004; Bressan et al. 2012). In the exponential overshoot formalism adopted in this work, we bypass this issue and thus do not require a secondary parameterization involving M₁ and M₂.

We adopt a modest overshoot efficiency ${f}_{\mathrm{ov},\mathrm{core}}\;=\;0.016$ for the core that is roughly equivalent to ${f}_{\mathrm{ov},\mathrm{step}}\;=\;0.2$ in the step overshoot scheme (Magic et al. 2010). This value is calibrated to reproduce the shape of the MSTO in the open cluster M67 (Section 8.3.1). However, it is worth noting that the strength of core convective overshoot depends on numerous other factors as well. For instance, Stothers & Chin (1991) explored the role of opacities in models with core overshoot and found that the overall increase in radiative opacities from the OPAL group (Iglesias & Rogers 1991) compared to the older values from the Los Alamos groups (e.g., Cox & Stewart 1965, 1970) reduced the overshoot efficiency required to reproduce observations of intermediate- and high-mass stars. Magic et al. (2010) studied how variations in the solar abundances (Asplund et al. 2005 versus Grevesse & Sauval 1998), element diffusion, overshooting, and nuclear reaction rates influence the MSTO morphology in M67. The authors concluded that the appearance of the characteristic MSTO hook (also known as the "Henyey hook") depends sensitively on the choice of the input solar abundances and that the effects of uncertain input physics on the Henyey hook morphology are degenerate.

We emphasize that the strength of convective overshoot is calibrated purely empirically: the overshoot efficiency in the core is constrained by matching the MSTO in M67, and the overshoot efficiency in the envelope, ${f}_{\mathrm{ov},\mathrm{env}}$, is chosen along with ${\alpha }_{\mathrm{MLT}}$ during solar calibration (Section 4.1). As noted in Bressan et al. (2012), envelope overshoot has a negligible effect on the evolution, e.g., the MS lifetime, though it is believed to affect the surface abundances of light elements, the location of the red giant branch (RGB) bump, and the extension of the blue loops. We also remind the reader that the overshoot efficiency in shells, e.g., hydrogen-burning shells during the RGB, is set to ${f}_{\mathrm{ov},\mathrm{env}}$ for simplicity.

3.6.3. Semiconvection and Thermohaline Mixing

As noted in Section 3.6.1, we adopt the Ledoux criterion for convection in our models. Due to the additional composition gradient term, a region that is formally convectively unstable to Schwarzschild criterion may be stable according to the Ledoux criterion (i.e., a thermally unstable medium with a stabilizing, positive composition gradient), which leads to a type of mixing called semiconvection. The importance of semiconvection on the evolution of massive stars has been studied for many decades, e.g., during the core helium burning phase (CHeB) (Stothers & Chin 1975; Langer et al. 1985; Grossman & Taam 1996). The fraction of a star's core helium burning lifetime spent on the Hayashi line relative to that in the blue part of the H-R diagram, in other words, the ratio of red supergiants to blue supergiants, is found to depend sensitively on the inclusion of semiconvection in the model. Additionally, the resulting core mass has significant implications for the supernova progenitors, from their ability to undergo a successful explosion (e.g., Sukhbold & Woosley 2014) to the actual nucleosynthetic yields (e.g., Langer et al. 1989; Heger & Woosley 2002; Rauscher et al. 2002). Semiconvection also operates in lower-mass stars with convective cores on the MS, and it can have an important effect on the actual size and appearance of the core (e.g., Faulkner & Cannon 1973; Silva Aguirre et al. 2011; Paxton et al. 2013).

Alternately, a thermally stable medium may have a negative, destabilizing composition gradient, which triggers a different type of instability called thermohaline mixing. An inverted chemical composition gradient is rare in stars, since fusion usually occurs inside out and synthesizes lighter elements into heavier products. This phenomenon can occur due to mass transfer in binaries (Stothers & Simon 1969; Stancliffe et al. 2007), off-center ignition in semi-degenerate cores (Siess 2009), or the ³He(³He,2p)⁴He reaction taking place just beyond the hydrogen-burning shell during the RGB, horizontal branch (HB), and AGB (Eggleton et al. 2006; Charbonnel & Zahn 2007; Cantiello & Langer 2010; Charbonnel & Lagarde 2010; Stancliffe 2010). Thermohaline mixing is thought to be responsible for the modification of surface abundances of RGB stars near the luminosity bump that otherwise cannot be explained using standard models. However, more recent work suggests that thermohaline mixing alone cannot account for the observed surface abundance anomalies (Denissenkov 2010; Traxler et al. 2011; Wachlin et al. 2011, 2014).

In MESA, semiconvection and thermohaline mixing are both implemented as time-dependent diffusive processes. The diffusion coefficient for semiconvection is computed following Langer et al. (1983):

$$\begin{eqnarray}&&{D}_{\mathrm{sc}}\;=\;{\alpha }_{\mathrm{sc}}\left(\displaystyle \frac{K}{6{C}_{P}\rho }\right)\displaystyle \frac{{{\rm{\nabla }}}_{T}-{{\rm{\nabla }}}_{\mathrm{ad}}}{{{\rm{\nabla }}}_{{\rm{L}}}-{{\rm{\nabla }}}_{T}},\end{eqnarray} \tag{ 12 }$$

where K is the radiative conductivity, C_P is the specific heat at constant pressure, and ${\alpha }_{\mathrm{sc}}$ is a dimensionless free parameter. Similarly, the diffusion coefficient for thermohaline mixing is computed following Ulrich (1972) and Kippenhahn et al. (1980):

$$\begin{eqnarray}&&{D}_{\mathrm{th}}\;=\;{\alpha }_{\mathrm{th}}\left(\displaystyle \frac{3\;K}{2{C}_{P}\rho }\right)\displaystyle \frac{-\tfrac{{\chi }_{\mu }}{{\chi }_{T}}{{\rm{\nabla }}}_{\mu }}{{{\rm{\nabla }}}_{T}-{{\rm{\nabla }}}_{\mathrm{ad}}},\end{eqnarray} \tag{ 13 }$$

where ${\alpha }_{\mathrm{th}}$ is a dimensionless number that describes the aspect ratio of the mixing blobs or "fingers" (a large ${\alpha }_{\mathrm{th}}$ corresponds to slender fingers).

As summarized in Paxton et al. (2013), the range of ${\alpha }_{\mathrm{sc}}$ and ${\alpha }_{\mathrm{th}}$ adopted by various authors spans several orders of magnitude, partially due to differences in the implementation in various codes. We adopt ${\alpha }_{\mathrm{sc}}\;=\;0.1$, though values as small as 0.001 or as large as 1 can be found in the literature (Langer 1991; Yoon et al. 2006). For thermohaline mixing we adopt ${\alpha }_{\mathrm{th}}\;=\;666$, as this value has been shown to reproduce the surface abundances anomalies in RGB stars past the luminosity bump (Charbonnel & Zahn 2007; Cantiello & Langer 2010). Note, however, that in the literature ${\alpha }_{\mathrm{th}}$ spans the range 1–1000 (Kippenhahn et al. 1980; Charbonnel & Zahn 2007; Cantiello & Langer 2010; Stancliffe 2010; Wachlin et al. 2011). There are ongoing theoretical efforts aimed at eliminating these free parameters with full 3D simulations (e.g., Traxler et al. 2011; Brown et al. 2013; Spruit 2013; Wood et al. 2013).

3.6.4. Rotationally Induced Instabilities

3.6.5. Magnetic Effects

The implementation of mass loss in stellar evolution calculations is based on a number of observationally and theoretically motivated prescriptions. It is frequently cited as one of the most uncertain ingredients in stellar evolution and is thought to play a crucial role in the advanced stages of stellar evolution for low-mass stars and in all phases of evolution for massive stars (see Smith 2014 for a recent review on this topic). In this section we review our treatment of mass loss across the H-R diagram. We note that the total mass-loss rate is capped at ${10}^{-3}\;{M}_{\odot }\quad {\mathrm{yr}}^{-1}$ in all models to prevent convergence problems.

3.7.1. Low-mass Stars

Mass loss for stars with masses below $10\;{M}_{\odot }$ is treated via a combination of the Reimers (1975) prescription for the RGB and Blöcker (1995) for the AGB. Both mass-loss schemes are based on global stellar properties such as the bolometric luminosity, radius, and mass:

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{R}}}=4\times {10}^{-13}{\eta }_{{\rm{R}}}\displaystyle \frac{(L/{L}_{\odot })(R/{R}_{\odot })}{(M/{M}_{\odot })}\;{M}_{\odot }\quad {\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{B}}}=4.83\times {10}^{-9}{\eta }_{{\rm{B}}}\displaystyle \frac{{(L/{L}_{\odot })}^{2.7}}{{(M/{M}_{\odot })}^{2.1}}\displaystyle \frac{{\dot{M}}_{{\rm{R}}}}{{\eta }_{{\rm{R}}}}\;{M}_{\odot }\quad {\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 15 }$$

where ${\eta }_{{\rm{R}}}$ and ${\eta }_{{\rm{B}}}$ are factors of order unity. These free parameters have been tuned to match numerous observational constraints, including the initial–final mass relation (IFMR) (Section 8.2; see also Kalirai et al. 2009), AGB luminosity function (Section 8.5.1; see also Rosenfield et al. 2014), and asteroseismic constraints from open cluster members in the Kepler fields (Miglio et al. 2012). The Blöcker (1995) mass-loss scheme was proposed as an alternative to the classic Reimers (1975) prescription to account for the onset of the superwind phase found in dynamical simulations of atmospheres of Mira-like variables. However, we emphasize that these are still empirically motivated recipes and therefore remain agnostic on the subject of the actual mechanism driving the winds (see a review by Willson 2000 for a discussion on, e.g., dust-driven winds).

For simplicity, we turn on Reimers mass loss at the beginning of the evolution, but a negligible amount of mass loss occurs throughout the MS ($\sim {10}^{-13}\;{M}_{\odot }\quad {\mathrm{yr}}^{-1}$ for a solar-metallicity $1\;{M}_{\odot }$ star). Once core helium is depleted, the mass-loss rate is chosen to be max[${\dot{M}}_{{\rm{R}}}$, ${\dot{M}}_{{\rm{B}}}$] at any given time. We adopt ${\eta }_{{\rm{R}}}\;=\;0.1$ and ${\eta }_{{\rm{B}}}\;=\;0.2$ in order to reproduce the IFMR (Section 8.2) and the AGB luminosity functions in the Magellanic Clouds (Section 8.5.1).

3.7.2. High-mass Stars

3.7.3. Rotationally Enhanced Mass Loss

Observations of O- and B-type stars have long argued for rotationally enhanced mass-loss rates (Gathier et al. 1981; Vardya 1985; Nieuwenhuijzen & de Jager 1988). It is now a standard ingredient in modern stellar evolution models with rotation (e.g., Heger et al. 2000; Brott et al. 2011a; Potter et al. 2012). In MESA, mass-loss rates are enhanced in models as a function of surface angular velocity Ω as follows:

$$\begin{eqnarray}&&\dot{M}({\rm{\Omega }})=\dot{M}(0){\left(\displaystyle \frac{1}{1-{\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit}}}\right)}^{\xi },\end{eqnarray} \tag{ 21 }$$

where $\dot{M}(0)$ is the standard mass-loss rate (Reimers, Blöcker, or "Dutch"), ξ is assumed to be 0.43 (Friend & Abbott 1986; Bjorkman & Cassinelli 1993; Langer 1998), and ${{\rm{\Omega }}}_{\mathrm{crit}}$ is the critical angular velocity at the surface:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{crit}}^{2}\;=\;\left(1-\displaystyle \frac{L}{{L}_{\mathrm{Edd}}}\right)\displaystyle \frac{{GM}}{{R}^{3}}.\end{eqnarray} \tag{ 22 }$$

The Eddington luminosity, ${L}_{\mathrm{Edd}}$, is a mass-weighted average over the optical depth τ between 1 and 100. For stars close to the Eddington limit, ${{\rm{\Omega }}}_{\mathrm{crit}}$ approaches 0 and therefore even a small Ω will result in a dramatic boost according to Equation (21). To prevent the mass loss from becoming too catastrophic, we cap the rotational boost to 10⁴.

As mentioned in Section 3.6.1, it is customary to calibrate the mixing length parameter, ${\alpha }_{\mathrm{MLT}},$ using helioseismic data and surface properties of the Sun. We utilize the MESA test suite solar_calibration, which conducts an extensive parameter search using the simplex method to obtain a set of input parameters that reproduces the solar observations. We vary the initial composition of the Sun, ${\alpha }_{\mathrm{MLT}},$ and convective overshoot in the envelope. For each iteration, a new set of parameters is drawn and the star is evolved from the PMS to 4.57 Gyr.¹⁷ A global ${\chi }^{2}$ value is computed by summing over the $\mathrm{log}L$, $\mathrm{log}R$, surface composition, ${R}_{\mathrm{cz}}$ (the location of the base of the convection zone), and $\delta {c}_{{\rm{s}}}$ (model−observed sound speed) terms with non-uniform user-defined weights. This process is repeated until ${\chi }^{2}$ ceases to change considerably and the tolerance parameters are met. For simplicity and consistency, the target solar values we adopt are the same nominal values recommended by the IAU.

Table 2 summarizes the solar calibration results, and Figure 1 shows $\delta {c}_{{\rm{s}}}/{c}_{{\rm{s}}}$, the fractional error in sound speed compared to that from helioseismic observations (Rhodes et al. 1997), as a function of radius. The best-fit MIST solar-calibrated model is shown in black, and two Serenelli et al. (2009) models adopting Grevesse & Sauval (1998) (GS98) and Asplund et al. (2009) (AGSS09) protosolar abundances are shown in red and blue, respectively.

Zoom In Zoom Out Reset image size

Table 2.  Solar Calibration Results

Parameter	Target	Model Value	Fractional Error (%)
${L}_{\odot }({10}^{33}\;\mathrm{erg}\quad {{\rm{s}}}^{-1})$	3.828a	3.828	$4.1\times {10}^{-3}$
${R}_{\odot }({10}^{10}\;\mathrm{cm})$	6.957a	6.957	$1.2\times {10}^{-3}$
${T}_{\mathrm{eff},\odot }$ (K)	5772a	5772	$2.4\times {10}^{-3}$
${X}_{\mathrm{surf}}$	0.7381b	0.7514	1.8
${Y}_{\mathrm{surf}}$	0.2485c	0.2351	5.4
${Z}_{\mathrm{surf}}$	0.0134b	0.0134	0.3
${R}_{\mathrm{cz}}$	0.7133d	0.7321	2.6
${\alpha }_{\mathrm{MLT}}$	...	1.82	...
${f}_{\mathrm{ov},\mathrm{env}}$	...	0.0174	...
X${}_{\mathrm{initial}}$	...	0.7238	...
Y${}_{\mathrm{initial}}$	...	0.2612	...
Z${}_{\mathrm{initial}}$	...	0.0150	...

Notes.

^aXXIXth IAU resolutions B2 and B3. ^bAsplund et al. (2009). ^cBasu & Antia (2004). ^dBasu & Antia (2004).

Download table as:  ASCIITypeset image

Although the model $\mathrm{log}L$ and $\mathrm{log}{T}_{\mathrm{eff}}$ are in excellent agreement with the observed values, there are noticeable discrepancies in the surface helium abundance, the location of the base of the convection zone, and the sound speed profile. Although many initial guesses and different weighting schemes were explored, we were unable to obtain a solar model that satisfies all available observational constraints. This is a well-known problem for solar models that adopt the AGSS09 abundances (e.g., Asplund et al. 2009; Serenelli et al. 2009).

The helium surface abundance at 4.57 Gyr in the best solar model is much lower compared to the helioseismologically inferred value of 0.2485 ± 0.0034 (Basu & Antia 2004). While most elemental abundances are determined through 3D spectroscopic modeling, the helium abundance is inferred indirectly from helioseismology, relying on the change in the adiabatic index in the He ii ionization zone near the surface (Asplund et al. 2009). The tension between the helioseismic value and the inferred abundance from interior modeling was noted in Asplund et al. (2009) and remains an unsolved problem to this day.

The sound speed profile comparison shows that the GS98 model is in good agreement while the models adopting the AGSS09 abundances show a large deviation at $\sim 0.6\;{R}_{\odot }$. This is partly due to the discrepancy between the predicted and observed locations of the convective boundary. In particular, the lower oxygen and neon abundances in AGSS09 relative to the older models like GS98 (or an even newer model like Caffau et al. 2011) imply a smaller mean opacity below the convective zone, which shifts the inner convective boundary farther out in radius. The abundance of oxygen, one of the most abundant and important elements, has undergone a striking overall downward revision over the past few decades. Still, there are likely lingering uncertainties in the surface abundance determinations due to the challenges associated with spectroscopic modeling, such as non-LTE effects, line blending, and uncertainties in the atomic and molecular data (Asplund et al. 2009).

The Serenelli et al. (2009) AGSS09 sound speed profile is in better agreement with the observed profile, likely because their model matches the location of the base of the convection zone more closely than the MIST model predicts. However, their model prefers an even lower surface helium abundance compared to the best-fit helium abundance in our model, and the present-day luminosity, radius, and effective temperature are not included as part of their fit.

Several explanations have been offered to reconcile the mismatch between the standard AGSS09 solar model and helioseismology results. One resolution invokes increased opacity, an idea that has been quantitatively explored by several authors (e.g., Bahcall et al. 2005a; Christensen-Dalsgaard et al. 2009; Serenelli et al. 2009). These authors concluded that a $\sim 10\%$ increase is required to match the observations, but that current atomic physics calculations do not leave room for such substantial change in opacities. However, there was a recent upward revision of iron opacities, based on new experimental data, that might account for roughly half the increase in the total mean opacity required to resolve this problem (Bailey et al. 2015). Other possible resolutions include more efficient diffusion processes in the radiative zone (e.g., Asplund et al. 2004), increased neon abundance to compensate for decreased oxygen abundance (e.g., Antia & Basu 2005; Bahcall et al. 2005b; Drake & Testa 2005), the introduction of new physics currently missing from stellar evolution calculations (e.g., convectively induced mixing in radiative zone; Young & Arnett 2005), and improved implementations of the current input physics (e.g., a replacement for MLT; Arnett et al. 2015).

Although the widely adopted practice is to fix the solar-calibrated ${\alpha }_{\mathrm{MLT}}$ across all masses, evolutionary phases, and abundances, there has been a recent effort to map out ${\alpha }_{\mathrm{MLT}}$ as a function of $\mathrm{log}g$ and ${T}_{\mathrm{eff}}$ (Trampedach 2007; Trampedach et al. 2014) as well as metallicity (Magic et al. 2015) from full 3D radiative hydrodynamic calculations of convection in stellar atmospheres. Recently, Salaris & Cassisi (2015) included variable ${\alpha }_{\mathrm{MLT}}$ and $T(\tau )$ boundary condition from Trampedach et al. (2014) in their stellar evolution calculations and found that varying ${\alpha }_{\mathrm{MLT}}$ has a small effect on the evolution and surface properties. We adopt a constant value of ${\alpha }_{\mathrm{MLT}}\;=\;1.82$ for the present work, but discuss the implications of this assumption in more detail in Section 8.3.1.

In summary, we adopt solar-calibrated ${\alpha }_{\mathrm{MLT}}$ and convective overshoot efficiency in the envelope ${f}_{\mathrm{ov},\mathrm{env}}$ (${f}_{0,\mathrm{ov},\mathrm{env}}$ is fixed to $0.5{f}_{\mathrm{ov},\mathrm{env}}$). As noted earlier in the Introduction, we adopt the Asplund et al. (2009) abundance mixture, and "solar metallicity" in this paper refers to the initial bulk $Z\;=\;{Z}_{\odot ,\mathrm{protosolar}}\;=\;0.0142$.

Lithium is a very fragile element that is burned via proton capture at temperatures as low as $2.5\times {10}^{6}$ K. Mixing processes such as convection that transport lithium from the outer layers to the interior, where the temperature is sufficiently high, lead to the depletion of surface lithium on very short timescales. The time evolution of surface lithium abundance therefore depends very sensitively on the initial stellar mass (proxy for temperature) and mixing physics.

Standard stellar evolution models so far have not been able to successfully reproduce the solar surface lithium abundance, indicating the need to include extra physical mechanisms. One way to account for missing physics in the models is to vary ${\alpha }_{\mathrm{MLT}}$ (see, e.g., Lyon models; Baraffe et al. 1998, 2003, 2015). Alternatively, models that incorporate the effects of rotation and internal gravity waves (e.g., Charbonnel & Talon 2005) are able to reproduce both the solar interior rotation profile and surface lithium abundances for the Sun and other galactic cluster stars. Relatedly, Somers & Pinsonneault (2014) found that radius dispersion on the PMS (correlated with rotation and chromospheric activity) can explain the spread in lithium abundances in young clusters such as the Pleiades.

In Figure 2 we show the evolution of surface lithium abundance relative to hydrogen, ${\rm{A}}{(}^{7}\mathrm{Li})\;=\;\mathrm{log}({N}_{{}^{7}\mathrm{Li}}/{N}_{{}^{1}{\rm{H}}})+12$, as a function of time for several $1\;{M}_{\odot }$ models. The purple square and circle are surface lithium abundance for the present-day Sun (Asplund et al. 2009) and the typical surface lithium abundance for nearby solar-metallicity, young clusters (e.g., Jeffries & Oliveira 2005; Sestito & Randich 2005; Juarez et al. 2014). Each pair of numbers in parentheses corresponds to the two envelope overshoot efficiency parameters ${f}_{\mathrm{ov},\mathrm{env}}$ and ${f}_{0,\mathrm{ov},\mathrm{env}}$, respectively. The surface lithium abundance decreases over time in all of the displayed models due to the inclusion of diffusion. The solid black line represents the fiducial model that adopts the solar-calibrated envelope overshoot parameters (0.0174, 0.0087) and ${\alpha }_{\mathrm{MLT}}\;=\;1.82$. The fiducial model burns too much lithium early on and then does not deplete lithium efficiently on the MS. The two dashed maroon and blue lines are models with less and more efficient overshoot, resulting in reduced and enhanced lithium depletion, respectively. The three dot-dashed lines show additional variations in input physics: the red lines correspond to ${\alpha }_{\mathrm{MLT}}\;=\;1.7$, while the orange and green lines are models that include PMS rotation with $v/{v}_{\mathrm{crit}}\;=\;0.01$ and 0.10, respectively. As expected, a lower ${\alpha }_{\mathrm{MLT}}$ produces a puffier, cooler star, resulting in less lithium depletion. The inclusion of rotation during the PMS produces very different, potentially more promising behavior, though the current absence of magnetic braking (see Sections 3.6.5 and 3.5) in MESA results in an unrealistically high rotation speed ($v/{v}_{\mathrm{crit}}$ between 0.1 and 1) on the MS. The models presented here fail to simultaneously match both the young and present-day solar surface lithium abundances. The inclusion of rotation on the PMS, the implementation of magnetic braking, and the exploration of a variable ${\alpha }_{\mathrm{MLT}}$ to mimic the effects of non-standard physics are planned for follow-up investigations in the near future.

Zoom In Zoom Out Reset image size

One of the main goals of the MIST project is to produce extensive grids of stellar evolutionary tracks and isochrones that cover a wide range in stellar masses, ages, evolutionary phases, and metallicities. The stellar mass of evolutionary tracks ranges from $0.1$ to $300\;{M}_{\odot }$ for a total of $\gtrsim 100$ models, and the ages of isochrones cover $\mathrm{logAge}\;=\;5$ to $\;10.3$ in 0.05 dex steps. For ${M}_{{\rm{i}}}\leqslant 0.7\;{M}_{\odot }$, the models are terminated at TAMS, i.e., central ¹H abundance drops to 10⁻⁴. For a $0.7\;{M}_{\odot }$ star at Z_⊙, this limit is typically reached at an age >35 Gyr. For ${M}_{{\rm{i}}}\gt 0.7\;{M}_{\odot }$, the models are evolved through either the WD cooling phase ("low-mass" type) or the end of carbon burning ("high-mass" type), depending on which criterion is satisfied first. We adopt this flexible approach to take into account the blurry boundary—further complicated by its metallicity dependence—between low- and intermediate-mass stars that end their lives as WDs and high-mass stars that continue to advanced stages of burning. In particular, the "low-mass" type models are terminated when Γ, the central plasma interaction parameter, also known as the Coulomb coupling parameter, exceeds 20. Γ is defined to be ${\bar{Z}}^{2}{e}^{2}/{a}_{i}{k}_{b}T$, where $\bar{Z}$ is the average ion charge, e is the electron charge, a_i is the mean ion spacing, k_b is the Boltzmann constant, and T is the temperature. A large Γ corresponds to a departure from the ideal gas limit toward solidification (crystallization of a pure oxygen WD occurs at ${\rm{\Gamma }}\approx 175;$ Paxton et al. 2011). The "high-mass" type models—stars that are sufficiently massive to burn carbon—are stopped when the central ¹²C abundance drops to 10⁻². The metallicity ranges from $[\mathrm{Fe}/{\rm{H}}]\;=\;-2.0$ to $+0.5$, with 0.25 dex spacing. We also provide an additional set of models evolved from the PMS to the end of core helium burning for $-4.0\leqslant [{\rm{Z}}/{\rm{H}}]\lt -2.0$ for modeling ancient, metal-poor populations. We provide a limited set at low metallicities at this time due to computational difficulties we have encountered. In particular, mixing between convective boundaries during the thermally pulsating AGB phase results in the ingestion of protons into a burning region, resulting in dramatically higher nuclear burning luminosities (e.g., Lau et al. 2009; Stancliffe et al. 2011; Woodward et al. 2015). Non-solar-scaled abundance grids will be presented in Paper II.

We note that in any grid, there is a subset of models that does not run to completion due to convergence issues. This is not generally problematic because the mass sampling is sufficiently fine such that there are enough models to smoothly interpolate a new EEP track and/or construct isochrones. We also note that there are interesting features in the tracks and isochrones that may appear to be numerical issues at first glance, but in fact a number of them are real phenomena captured in the MESA calculations. We refer the reader to the Appendix for a discussion of these features. Of these, an example feature that may be a numerical artifact is an extremely short-lived glitch that appears during the early post-AGB phase for a subset of the models. Since it has zero bearing on the evolution of the star, we post-process this feature out of the final evolutionary tracks in order to facilitate the construction of smooth isochrones.

We briefly describe the isochrone construction process and defer a detailed discussion of this topic to Paper 0 (Dotter 2016). In each evolutionary track, we identify a set of the so-called primary equivalent evolutionary points (primary EEPs). These correspond to specific stages of evolution defined by a set of physical conditions, such as the terminal-age main sequence (TAMS; central hydrogen exhaustion) and the tip of the RGB (RGBTip; a combination of lower limits on the central helium abundance and luminosity). Next, each segment between adjacent primary EEPs is further divided into so-called "secondary-EEPs" according to a distance metric that evenly samples the tracks in certain relevant variables, such as ${T}_{\mathrm{eff}}$ and L. Put another way, primary EEPs serve as physically meaningful reference locations along the evolutionary track, and secondary EEPs finely sample the track between primary EEPs. This method maps a set of evolutionary tracks from ordinary time coordinates onto uniform EEP coordinates. The primary EEPs and corresponding evolutionary phases are listed in Table 3. In Figures 3 and 4, we show example $1\;{M}_{\odot }$ and $30\;{M}_{\odot }$ evolutionary tracks in the EEP format, with colored dots marking the locations of the primary EEPs. Note that we require both IAMS and TAMS points in order to properly resolve the MSTO for stars that burn hydrogen convectively in their cores during the MS, i.e., the Henyey hook. Also note that although "RGBTip" does not have the same morphological significance in high-mass stars as in low-mass stars since high-mass stars ignite helium under non-degenerate conditions, we retain this terminology for consistency reasons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As described in Section 3.3, we use three different boundary conditions depending on the initial mass of the model ($\tau \;=\;100$ tables for $0.1\mbox{--}0.3\;{M}_{\odot }$, photosphere tables for $0.6\mbox{--}10\;{M}_{\odot }$, and simple photosphere for $16\mbox{--}300\;{M}_{\odot }$). To facilitate a smooth transition from one regime to another, we run both $\tau \;=\;100$ and photosphere tables for $0.3\mbox{--}0.6\;{M}_{\odot }$ and photosphere tables and simple photosphere for $10\mbox{--}16\;{M}_{\odot }$. For every mass in the transition regime, the two EEP tracks are blended with a smooth weighting function to create a hybrid EEP track:

$$\begin{eqnarray}&&w\;=\;0.5\left[1-\mathrm{cos}\left(\pi \displaystyle \frac{{M}_{\ast }-{M}_{1}}{{M}_{2}-{M}_{1}}\right)\right],\end{eqnarray} \tag{ 23 }$$

where M₁ and M₂ are the transition masses (e.g., ${M}_{1}\;=\;0.3\;{M}_{\odot }$ and ${M}_{2}\;=\;0.6\;{M}_{\odot }$ for the transition from $\tau \;=\;100$ to photosphere tables).

To generate an isochrone at age t₀ with all of the EEP tracks now in hand, we first cycle through all masses and construct a piecewise monotonic function between ${M}_{{\rm{i}}}$ and t for each EEP point.¹⁸ Next, we interpolate to obtain ${M}_{{\rm{i}}}({t}_{0})$. Once we have ${M}_{{\rm{i}}}({t}_{0})$ for every EEP, we can now construct an isochrone for any parameter, e.g., L, by interpolating that parameter as a function of ${M}_{{\rm{i}}}$, e.g., $L({M}_{{\rm{i}}}({t}_{0}))$, at every EEP. The EEP framework is superior to a direct interpolation scheme in time coordinates as it can properly treat evolutionary phases with short timescales (e.g., post-AGB) or complex trajectories (e.g., thermally pulsating AGB, TPAGB).

A sensible approach is to construct a monotonic relationship between mass and age assuming that "increasing mass = decreasing phase lifetime" is always true. However, interesting non-monotonic behaviors begin to appear in certain evolutionary phases over a narrow age (or equivalently, mass) interval if the mass resolution is sufficiently high (see Figure 13 and also Girardi et al. 2013). Put another way, two stars of different initial masses are at the same evolutionary stage in terms of their EEPs over a special narrow age interval. This effect will be explored in future work, but for this work, we enforce monotonicity in the age–mass relationship.

The published tracks and isochrones include a wealth of information, ranging from basic parameters such as $\mathrm{log}(L)$, $\mathrm{log}({T}_{\mathrm{eff}})$, $\mathrm{log}(g)$, and surface abundances of 19 elements to more detailed quantities such as $\beta \equiv {P}_{\mathrm{gas}}/{P}_{\mathrm{total}}$ and asteroseismic parameters (the full list of available parameters is available on the project Web site). To highlight this fact, we show examples of isochrones in several different projections in Figure 5. From left to right, the top three panels feature isochrones in the $\mathrm{log}(g)-\mathrm{log}({T}_{\mathrm{eff}})$, ${\rm{\Delta }}\nu -\mathrm{log}({T}_{\mathrm{eff}})$, and ${\nu }_{\mathrm{max}}-\mathrm{log}({T}_{\mathrm{eff}})$ planes. Δν and ${\nu }_{\mathrm{max}}$ are asteroseismic quantities that correspond to the large frequency separation for p-modes and the frequency of maximum power, respectively, which can be readily obtained from power spectra of, e.g., Kepler light curves. We clarify that Δν and ${\nu }_{\mathrm{max}}$ in the MIST models are computed from simple scaling relations (e.g., Ulrich 1986; Brown et al. 1991; Kjeldsen & Bedding 1995) and not from full pulsation analysis, though the pulsation code GYRE (Townsend & Teitler 2013) is integrated into MESA. The bottom left panel shows isochrones in the $\mathrm{log}({T}_{c})$–$\mathrm{log}({\rho }_{c})$ plane, where the dashed, dot-dashed, and solid lines show thresholds for hydrogen, helium, and carbon ignition. The 10 Myr isochrone contains massive stars that are able to ignite carbon, whereas at the older ages, the isochrones cannot reach sufficiently high central densities and temperatures. At 10 Gyr, only the low-mass stars remain and helium core flash at RGBTip shows up as a discontinuous sharp feature. The bottom middle and right panels show $\mathrm{log}(L/{L}_{\mathrm{Edd}})$, the ratio of total luminosity to Eddington luminosity, and $\dot{M}$, the mass-loss rate, as a function of initial mass. As expected, both quantities generally increase as the initial mass increases. A very prominent increase immediately followed by a sharp decrease at intermediate ages (0.1 and 1 Gyr) in both panels is due to the TPAGB phase, followed by the post-AGB and WD cooling phases.

Zoom In Zoom Out Reset image size

Bolometric corrections are necessary to transform theoretical isochrones into magnitudes that allow for direct comparisons with observations. The bolometric corrections are largely based on a new grid of stellar atmosphere and synthetic spectra created with the ATLAS12 and SYNTHE codes (C. Conroy et al. 2016, in preparation). These same models are used for the surface boundary conditions discussed in Section 3.3. They include the latest atomic line list from R. Kurucz (including both laboratory and predicted lines) and many molecules including CH, CN, TiO, H₂, H₂O, SiO, C₂, SiH, MgH, CrH, CaH, FeH, CO, NH, VO, and OH. We have also computed model atmosphere and spectra for carbon stars with ${\rm{C}}/{\rm{O}}\;=\;1.05$ over the range $2400\quad {\rm{K}}\lt {T}_{\mathrm{eff}}\lt 4700$ K and $-1.0\lt \mathrm{log}g\;[{\rm{g}}\;{\mathrm{cm}}^{-3}]\lt 0.5$. Our carbon star spectra agree well with the models of Aringer et al. (2009). The primary differences arise at $\gt 2\;\mu {\rm{m}}$ as our models do not currently include the important molecules C₃, HCN, and C₂H₂. We chose to create our own carbon star spectra in order to have models covering the full wavelength range and at the same resolution as our main spectral library. The ATLAS12/SYNTHE spectra are combined with synthetic spectra for H-rich WDs with ${\rm{6000}}\quad {\rm{K}}\leqslant {T}_{\mathrm{eff}}\leqslant {\rm{50,000}}$ K from Koester (2010). These are supplemented by a set of blackbody spectra with ${\rm{200,000}}\quad {\rm{K}}\leqslant {T}_{\mathrm{eff}}\leqslant 1$ million K. The coverage of the different synthetic libraries is shown in Figure 6. Example isochrones at $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;7.0$, 7.5, 8.5, 9.0, and 10.0 are overplotted for reference.

Zoom In Zoom Out Reset image size

Bolometric corrections are computed from the synthetic spectra following Equation (1) of Girardi et al. (2008). As noted in the Introduction, we adopt as a zero point ${L}_{\circ }\;=\;3.0128\times {10}^{35}$ erg s⁻¹ to define ${M}_{\mathrm{bol}}\;=\;0$. This is equivalent to adopting solar values of ${M}_{\mathrm{bol},\odot }\;=\;4.74$ and ${L}_{\odot }\;=\;3.828\times {10}^{33}$ erg s⁻¹. The bolometric corrections include a range of extinction values, as characterized by both A_V and R_V, following the extinction curve of Cardelli et al. (1989). We provide A_V = 0 to 6 with R_V = 3.1, though other R_V values can be made upon request. We emphasize that Z of the bolometric correction is matched to the surface Z (and surface ${\rm{C}}/{\rm{O}}$ ratio where relevant) for each star along the isochrone.

The photometric systems included in the initial MIST release are summarized in Table 4. This is only an initial set and will expand over time. Photometric systems define their magnitude scales according to a flux standard.

Table 5.  Convergence Test Results

$1\;{M}_{\odot }$
C: end of H burn	4	2	1	1/2	1/4	C: end of He burn	4	2	1	1/2	1/4
${N}_{\mathrm{cells}}$	347	486	821	1593	3174	${N}_{\mathrm{cells}}$	586	1139	2208	3931	6931
${N}_{\mathrm{timesteps}}$	335	500	818	1432	2726	${N}_{\mathrm{timesteps}}$	9607	12185	17826	32928	66044
$M/{M}_{\odot }$	1.0	1.0	1.0	1.0	1.0	$M/{M}_{\odot }$	0.95	0.951	0.951	0.951	0.951
$\xi \;(\times {10}^{-3})$	0.006	0.009	0.004	0.001	0.0	$\xi \;(\times {10}^{-3})$	−0.615	−0.148	0.063	0.046	0.0
$\mathrm{log}(L/{L}_{\odot })$	0.261	0.256	0.256	0.256	0.256	$\mathrm{log}(L/{L}_{\odot })$	2.031	2.024	2.008	2.001	2.001
$\xi \;(\times {10}^{-3})$	10.64	1.151	0.586	1.054	0.0	$\xi \;(\times {10}^{-3})$	70.953	53.714	14.967	−0.799	0.0
$\mathrm{log}({T}_{\mathrm{eff}})\quad [{\rm{K}}]$	3.764	3.766	3.767	3.767	3.767	$\mathrm{log}({T}_{\mathrm{eff}}\;[{\rm{K}}])$	3.63	3.634	3.636	3.637	3.637
$\xi \;(\times {10}^{-3})$	−5.652	−1.65	−0.493	−0.159	0.0	$\xi \;(\times {10}^{-3})$	−14.812	−6.501	−1.779	−0.06	0.0
$\mathrm{log}(t)$	9.94	9.939	9.941	9.942	9.943	$\mathrm{log}(t)$	10.056	10.057	10.059	10.06	10.06
$\xi \;(\times {10}^{-3})$	−5.617	−8.577	−3.747	−0.836	0.0	$\xi \;(\times {10}^{-3})$	−10.478	−7.686	−3.576	−1.376	0.0
$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	7.295	7.294	7.294	7.294	7.294	$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	8.282	8.28	8.275	8.273	8.273
$\xi \;(\times {10}^{-3})$	1.165	−0.103	0.037	0.202	0.0	$\xi \;(\times {10}^{-3})$	20.975	16.362	4.178	0.649	0.0
$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	2.805	2.788	2.786	2.786	2.784	$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	4.754	4.742	4.733	4.727	4.726
$\xi \;(\times {10}^{-3})$	48.7	10.331	4.975	4.094	0.0	$\xi \;(\times {10}^{-3})$	66.899	37.961	15.172	2.266	0.0
$5\;{M}_{\odot }$
C: end of H burn	4	2	1	1/2	1/4	C: end of He burn	4	2	1	1/2	1/4
${N}_{\mathrm{cells}}$	345	547	1015	2003	3992	${N}_{\mathrm{cells}}$	603	1248	2148	4582	7283
${N}_{\mathrm{timesteps}}$	383	584	1019	1919	3699	${N}_{\mathrm{timesteps}}$	644	1117	2123	4262	8454
$M/{M}_{\odot }$	4.997	4.997	4.997	4.997	4.997	$M/{M}_{\odot }$	4.984	4.985	4.985	4.985	4.986
$\xi \;(\times {10}^{-3})$	−0.057	−0.051	−0.02	−0.019	0.0	$\xi \;(\times {10}^{-3})$	−0.355	−0.124	−0.033	−0.07	0.0
$\mathrm{log}(L/{L}_{\odot })$	3.168	3.154	3.139	3.137	3.136	$\mathrm{log}(L/{L}_{\odot })$	3.294	3.23	3.173	3.175	3.17
$\xi \;(\times {10}^{-3})$	74.552	40.912	5.28	0.766	0.0	$\xi \;(\times {10}^{-3})$	330.852	147.317	7.776	12.006	0.0
$\mathrm{log}({T}_{\mathrm{eff}})\quad [{\rm{K}}]$	4.159	4.16	4.16	4.16	4.164	$\mathrm{log}({T}_{\mathrm{eff}}\;[{\rm{K}}])$	3.614	3.624	3.63	3.63	3.631
$\xi \;(\times {10}^{-3})$	−12.433	−10.03	−7.975	−8.023	0.0	$\xi \;(\times {10}^{-3})$	−37.609	−15.963	−1.504	−1.559	0.0
$\mathrm{log}(t)$	8.025	8.021	8.01	8.009	8.008	$\mathrm{log}(t)$	8.077	8.076	8.071	8.072	8.071
$\xi \;(\times {10}^{-3})$	40.171	29.525	5.393	2.926	0.0	$\xi \;(\times {10}^{-3})$	13.82	12.067	−0.536	2.725	0.0
$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	7.663	7.658	7.657	7.657	7.657	$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	8.379	8.371	8.368	8.371	8.369
$\xi \;(\times {10}^{-3})$	13.611	2.911	0.973	0.874	0.0	$\xi \;(\times {10}^{-3})$	21.884	3.362	−2.952	4.452	0.0
$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	2.001	2.001	2.006	2.007	2.007	$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	4.243	4.245	4.246	4.247	4.243
$\xi \;(\times {10}^{-3})$	−13.275	−14.014	−0.887	0.314	0.0	$\xi \;(\times {10}^{-3})$	2.19	6.022	7.495	9.271	0.0
$20\;{M}_{\odot }$
C: end of H burn	4	2	1	1/2	1/4	C: end of He burn	4	2	1	1/2	1/4
${N}_{\mathrm{cells}}$	1467	2683	5024	9912	17655	${N}_{\mathrm{cells}}$	2325	4598	9133	18303	36666
${N}_{\mathrm{timesteps}}$	1803	2493	4460	8438	16409	${N}_{\mathrm{timesteps}}$	2951	3635	6218	11591	22501
$M/{M}_{\odot }$	19.279	19.202	19.184	19.188	19.212	$M/{M}_{\odot }$	12.156	12.783	13.115	14.655	14.856
$\xi \;(\times {10}^{-3})$	3.456	−0.532	−1.475	−1.287	0.0	$\xi \;(\times {10}^{-3})$	−181.746	−139.532	−117.17	−13.479	0.0
$\mathrm{log}(L/{L}_{\odot })$	5.083	5.03	5.02	5.017	5.015	$\mathrm{log}(L/{L}_{\odot })$	5.174	5.131	5.121	5.101	5.099
$\xi \;(\times {10}^{-3})$	169.241	35.182	11.427	5.929	0.0	$\xi \;(\times {10}^{-3})$	190.493	77.659	53.169	4.741	0.0
$\mathrm{log}({T}_{\mathrm{eff}})\quad [{\rm{K}}]$	4.468	4.442	4.437	4.434	4.435	$\mathrm{log}({T}_{\mathrm{eff}}\;[{\rm{K}}])$	3.568	3.535	3.532	3.537	3.537
$\xi \;(\times {10}^{-3})$	78.496	17.202	4.979	−1.602	0.0	$\xi \;(\times {10}^{-3})$	74.031	−4.799	−10.998	−0.389	0.0
$\mathrm{log}(t)$	6.998	6.961	6.953	6.949	6.947	$\mathrm{log}(t)$	7.034	7.002	6.995	6.993	6.991
$\xi \;(\times {10}^{-3)}$	123.993	31.461	12.596	5.123	0.0	$\xi \;(\times {10}^{-3})$	103.348	25.481	9.44	4.182	0.0
$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	7.819	7.815	7.815	7.815	7.814	$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	8.511	8.505	8.504	8.503	8.503
$\xi \;(\times {10}^{-3})$	10.992	2.764	1.188	0.834	0.0	$\xi \;(\times {10}^{-3})$	17.555	5.504	2.622	0.795	0.0
$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	1.375	1.393	1.396	1.397	1.397	$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	3.525	3.545	3.548	3.55	3.551
$\xi \;(\times {10}^{-3})$	−48.836	−10.201	−2.678	−0.329	0.0	$\xi \;(\times {10}^{-3})$	−58.531	−13.623	−6.043	−1.536	0.0
$100\;{M}_{\odot }$
C: end of H burn	4	2	1	1/2	1/4	C: end of He burn	4	2	1	1/2	1/4
${N}_{\mathrm{cells}}$	1468	2727	5448	10818	...	${N}_{\mathrm{cells}}$	1569	3152	6486	12465	...
${N}_{\mathrm{timesteps}}$	275	323	323	328	...	${N}_{\mathrm{timesteps}}$	696	572	559	584	...
$M/{M}_{\odot }$	65.916	25.712	24.583	23.752	...	$M/{M}_{\odot }$	33.813	11.914	11.664	11.454	...
$\xi \;(\times {10}^{-3})$	1775.107	82.495	34.981	0.0	...	$\xi \;(\times {10}^{-3})$	1952.055	40.126	18.353	0.0	...
$\mathrm{log}(L/{L}_{\odot })$	6.274	5.805	5.774	5.749	...	$\mathrm{log}(L/{L}_{\odot })$	6.162	5.504	5.489	5.475	...
$\xi \;(\times {10}^{-3})$	2352.624	139.574	58.767	0.0	...	$\xi \;(\times {10}^{-3})$	3867.197	67.71	31.926	0.0	...
$\mathrm{log}({T}_{\mathrm{eff}})\quad [{\rm{K}}]$	4.603	4.882	4.877	4.874	...	$\mathrm{log}({T}_{\mathrm{eff}}\;[{\rm{K}}])$	5.228	5.2	5.199	5.198	...
$\xi \;(\times {10}^{-3})$	−463.387	18.996	8.502	0.0	...	$\xi \;(\times {10}^{-3})$	71.311	4.274	3.703	0.0	...
$\mathrm{log}(t)$	6.476	6.546	6.553	6.56	...	$\mathrm{log}(t)$	6.52	6.595	6.601	6.608	...
$\xi \;(\times {10}^{-3})$	−174.392	−29.716	−15.581	0.0	...	$\xi \;(\times {10}^{-3})$	−183.629	−28.855	−14.912	0.0	...
$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	7.894	7.87	7.868	7.867	...	$\mathrm{log}({T}_{{\rm{c}}})\quad [{\rm{K}}]$	8.56	8.529	8.53	8.528	...
$\xi \;(\times {10}^{-3})$	65.67	8.314	3.609	0.0	...	$\xi \;(\times {10}^{-3})$	78.495	3.91	4.801	0.0	...
$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	1.03	1.208	1.216	1.223	...	$\mathrm{log}({\rho }_{{\rm{c}}})\quad [{\rm{g}}\;{\mathrm{cm}}^{-3}]$	3.151	3.383	3.392	3.392	...
$\xi \;(\times {10}^{-3})$	−359.052	−35.601	−15.776	0.0	...	$\xi \;(\times {10}^{-3})$	−426.172	−19.84	−0.399	0.0	...

Download table as:  ASCIITypeset images: 1 2

As described in detail in Section 5.1, the MIST models cover a wide range in stellar masses, ages, metallicities, and evolutionary phases. The stellar mass of evolutionary tracks ranges from $0.1$ to $300\;{M}_{\odot }$, the ages of isochrones cover $\mathrm{logAge}\;=\;5$ to $\;10.3$, and the metallicity ranges from $[\mathrm{Fe}/{\rm{H}}]\;=\;-4.0$ to $+0.5$. The evolution is continuously computed from the PMS phase to the end of hydrogen burning, WD cooling phase, or the end of carbon burning, depending on the initial stellar mass and metallicity.

Figure 7 illustrates the range of stellar masses, ages, and evolutionary phases, showing the evolutionary tracks and isochrones in the left and right panels, respectively. As noted in Section 3, the models include rotation with $v/{v}_{\mathrm{crit}}\;=\;0.4$ by default. Figure 8 shows the effect of rotation on both the evolutionary tracks and isochrones, where the rotating and non-rotating models are shown in solid and dashed lines, respectively. For display purposes, we omit the post-AGB phase where relevant. Rotation makes stars more luminous during the MS because rotational mixing (see Section 3.6.4) promotes core growth. It makes the star appear hotter or cooler depending on the efficiency of rotational mixing in the envelope: if rotational mixing introduces a sufficient amount of helium into the envelope and increases the mean molecular weight, the star becomes more compact and hotter. However, in the absence of efficient rotational mixing, the centrifugal effect dominates, making the star appear cooler and more extended. This is because in a rotating system, ordinary gravitational acceleration g is replaced by ${g}_{\mathrm{eff}}$ that includes both the gravitational and centrifugal terms. Since ${T}_{\mathrm{eff}}\propto {g}_{\mathrm{eff}}^{1/4}$ for a star with a radiative envelope and $| {g}_{\mathrm{eff}}| \lt | g|$, ${T}_{\mathrm{eff}}$ is thus lower in a rotating system. Generally speaking, rotation increases ${T}_{\mathrm{eff}}$ in massive stars where rotational mixing operates efficiently, and it decreases ${T}_{\mathrm{eff}}$ in low-mass stars, as well as at all ZAMS locations where the centrifugal effect dominates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The top two panels of Figure 9 show a series of $2.2\;{M}_{\odot }$ (left) and $6.4\;{M}_{\odot }$ (right) evolutionary tracks for three metallicities. A portion of the PMS phase in both panels and the evolution following the CHeB phase in the right panel are omitted for display purposes. The tracks become hotter and more luminous with decreasing metal content due to a lower Rosseland mean opacity. In the bottom panels, we zoom in on the TPAGB and CHeB to further highlight the effects of metallicity on these evolutionary phases. The blue loops become hotter and more prominent with decreasing metallicity, and the entire TPAGB phase also shifts to hotter temperatures as metallicity decreases. To scrutinize the effect of metallicity on the TPAGB phase in more detail, we plot the number of TPs executed by each model as a function of mass at a number of metallicities in Figure 10. There are two notable features: the maximum number occurs at around $2\;{M}_{\odot }$ independent of metallicity; and the number of thermal pulses increases with [Z/H] from −1.0 to 0.0, and there is a hint that the trend reverses for $[{\rm{Z}}/{\rm{H}}]\gt 0$. We leave a more detailed discussion of the number of TPs as a function of uncertain physical parameters (e.g., mixing, mass loss) and comparisons to other databases for future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to illustrate the range of metallicities available in the current MIST models, we show 10 Gyr isochrones at $[{\rm{Z}}/{\rm{H}}]\;=\;-4$, −3, −2, −1, 0, and 0.5 in Figure 11. For clarity, we only show phases up to the RGBTip. As the metallicity decreases, the MSTO becomes hotter and more luminous (and the MSTO mass decreases), and the RGBTip becomes fainter due to the helium ignition occurring at lower core masses. Note that the isochrone changes more subtly with metallicity in the very metal-poor regime, i.e., $[{\rm{Z}}/{\rm{H}}]\lesssim -2$.

Zoom In Zoom Out Reset image size

In the top panel of Figure 12 we show phase lifetimes as a function of initial mass for $[{\rm{Z}}/{\rm{H}}]\;=\;-0.25$, 0.0, and $+0.25$ in solid, dot-dashed, and dotted lines, respectively. The bottom panel shows the ratio of phase lifetimes to the MS lifetime. Note that the "RGB" label refers to the phase between TAMS and helium ignition, which includes the short subgiant branch (SGB) evolution, and "post-AGB" includes the white dwarf cooling phase up to ${\rm{\Gamma }}\;=\;20$. The post-AGB timescales in the MIST models (adopting the definition from Miller Bertolami 2016) are consistent with those reported by Miller Bertolami (2016) and Weiss & Ferguson (2009), which are a factor of $3-10$ shorter compared to the older post-AGB stellar evolution models (Vassiliadis & Wood 1994; Blöcker 1995). High-mass stars are not included because they do not go through the same set of evolutionary phases featured here. The TPAGB and post-AGB phases are not shown for a subset of the models that do not completely evolve through those evolutionary stages. Unsurprisingly, the lifetimes generally decrease with increasing mass, though there are some notable exceptions, including the peak in CHeB and AGB lifetimes at $\sim 2\;{M}_{\odot }$ (see the discussion below).

Zoom In Zoom Out Reset image size

The left panel of Figure 13 is a slight variation of the previous plot, where we now show the cumulative age as a function of mass for $[{\rm{Z}}/{\rm{H}}]\;=\;0.0$. In the right panel, we zoom in on a particularly interesting mass range around $2\;{M}_{\odot }$, where there is a noticeable increase in the CHeB lifetime. This effect, explored in detail in Girardi et al. (2013), is due to the transition from the explosive ignition of helium in degenerate cores of low-mass stars, i.e., helium core flash at RGBTip, to quiescent ignition of helium in more massive stars. This is because more massive stars start burning helium at lower core masses and therefore have smaller initial helium-burning luminosities compared to those that undergo the helium core flash (see also Girardi 1999). Note that this non-monotonicity in the CHeB lifetime is present even in models computed without RGB mass loss. Thus, the simple "increasing mass = decreasing MS lifetime" rule of thumb is violated over a very narrow mass range. In other words, there is a special age at which two stars of different masses are at the same evolutionary stage in terms of their EEPs. This notion of "double EEPs" is discussed in more detail in Paper 0 (Dotter 2016) and will be explored in a future paper.

Zoom In Zoom Out Reset image size

To explore the effects of uncertain input physics and choices of free parameters on the resulting tracks and isochrones, we ran several sets of models varying one ingredient at a time. In Figure 14 we present four such variations. The black solid lines represent solar-metallicity models with fiducial parameters listed in Table 1. We note that several of these parameters affect the lifetimes of different evolutionary phases, though this is not explicitly shown in the figures.

Zoom In Zoom Out Reset image size

In the top left panel, we show two sets of isochrones, which are identical except for their surface boundary conditions. The three sets of ages correspond to $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;9.0,9.6$, and 10.2. The solid lines are isochrones with our default implementation of boundary conditions: $\tau \;=\;100$ and photosphere tables from the ATLAS12/SYNTHE atmosphere models plus the simple Eddington gray $T(\tau )$ approximation for the hottest stars. The dashed lines are models with the Eddington gray $T(\tau )$ relation applied across the entire mass range. For clarity, we only plot the isochrones up to the RGBTip. As expected, the choice of boundary conditions has the largest effect on the lower MS populated by cool, compact dwarfs. The shift in ${T}_{\mathrm{eff}}$ on the RGB is smaller but important, amounting to $\sim 50\mbox{--}60$ K. Differences both on the SGB and RGB and on the lower MS are larger at non-solar metallicities. In particular, the isochrones become more discrepant on the SGB and RGB at $[{\rm{Z}}/{\rm{H}}]\;=\;-1$ and on the lower MS at $[{\rm{Z}}/{\rm{H}}]\;=\;+0.3$.

In the top right panel, we show three sets of isochrones at $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;6.7,7.4$, and 9.0 to illustrate the impact mass loss has on various stages of evolution. The solid lines correspond to isochrones with ${\eta }_{{\rm{R}}}\;=\;0.1,\;{\eta }_{{\rm{B}}}\;=\;0.2$, and ${\eta }_{\mathrm{Dutch}}\;=\;1.0$, whereas the dashed lines and dotted lines represent mass-loss rates that are twice and half as efficient, respectively. For display purposes, we omit the PMS and post-AGB phases. The temperature evolution is especially sensitive to the mass-loss rates at the youngest ages because massive star evolution is strongly affected by the choice of input physics. The morphology of the CHeB is also directly influenced by mass-loss rates: lower and higher $\dot{M}$ values yield hotter and cooler CHeB, respectively. The morphology and lifetime of the notorious TPAGB phase are also affected mass-loss rates, with more efficient winds resulting in fainter and fewer TPs as expected. The mass-loss efficiency, parameterized by the η parameter, is calibrated empirically to match various observational constraints, including the number ratios of different stellar types for the high-mass stars (Section 9.3) and the AGB luminosity functions for the low-mass stars (Section 8.5.1).

In the bottom left panel, we highlight the importance of core convective overshoot. The three sets of ages shown are $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;7.0,7.5$, and 8.0. The solid lines are isochrones with the default core overshoot parameter ${f}_{\mathrm{ov},\mathrm{core}}\;=\;0.016$ in the exponential diffusive overshoot formalism. The dashed and dotted lines represent decreased and increased (${f}_{\mathrm{ov},\mathrm{core}}\;=\;0.012,\;0.020$) overshoot efficiency, respectively. For clarity, we plot the evolution up through helium burning only. Since convective overshoot enhances the mixing of fresh fuel into the core, a higher overshoot efficiency results in longer MS lifetimes and systematically higher MSTO masses and luminosities. Likewise, SGB and CHeB luminosities are higher for more efficient overshoot due to the larger resulting core masses. Note that the overshoot efficiency in the core is constrained by matching the MSTO in M67, and the overshoot efficiency in the envelope is determined during solar calibration.

In the bottom right panel, we show isochrones with different values of the mixing length parameter ${\alpha }_{\mathrm{MLT}}$. The four ages shown are $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;7.0,8.0,8.7$, and 10.0. The solid lines are isochrones with the solar-calibrated value ${\alpha }_{\mathrm{MLT}}\;=\;1.82$. The dashed and dotted lines correspond to isochrones with ${\alpha }_{\mathrm{MLT}}\;=\;1.6$ and ${\alpha }_{\mathrm{MLT}}\;=\;2.0$, respectively. For display purposes, we only show the evolution up through the RGBTip. The physical interpretation of a small value of ${\alpha }_{\mathrm{MLT}}$ is a fluid parcel traveling a short radial distance (in units of ${H}_{{\rm{P}}}$) before it deposits its internal energy and blends into the surrounding medium. The net effect is reduced convective efficiency, thus cooler temperatures and more inflated radii. For this reason, ${\alpha }_{\mathrm{MLT}}$ is used to mimic the effects of physical ingredients, e.g., the inhibition of convection by a magnetic field, that are missing from the majority of current stellar evolution models (but see Feiden & Chaboyer 2013). For example, ${\alpha }_{\mathrm{MLT}}$ that is smaller compared to the solar-calibrated value is commonly used to bring models into agreement with observations of inflated radii in low-mass stars (see Section 8.1).

We first compare isochrones at Z_⊙ adopted by each model. In Figure 15, we compare MIST (black; this work), PARSEC v1.2S (red; Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014), Y² (orange; Demarque et al. 2004), DSEP (green; Dotter et al. 2008), BaSTI "non-canonical" with $\eta \;=\;0.4$ (turquoise; Pietrinferni et al. 2004), and Lyon (navy; Baraffe et al. 1998, 2003, 2015) for $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;7.5$, 8.0, 9.0, and 10.0. The DSEP and Lyon models are not included in the top two panels because they are not available at young ages. Note that of the models featured here, only the MIST models include the effects of rotation.¹⁹ We choose MIST models with rotation instead of those without to make the comparison because the fiducial models include rotation and all calibrations are performed on this set. Overall, the MIST isochrones are in broad agreement with other isochrones. Although the absolute metal contents differ by as much as 30% between various models due to differences in the preferred definition of Z_⊙, the isochrones are less discrepant than one might imagine because they have been calibrated to match the properties of the Sun. There is more noticeable discrepancy at the young ages due to the complex and uncertain physics—such as core convective overshoot—governing the evolution of massive stars. In particular, the CHeB phase (e.g., the development of the blue loop) is notoriously sensitive to the details of input physics (e.g., McQuinn et al. 2011) though there are ongoing efforts to address this issue (e.g., Tang et al. 2016). Moreover, the PARSEC isochrones depart notably from the rest of the models on the lower MS due to their recent implementation of $T\mbox{--}\tau$ boundary conditions that have been empirically calibrated to match the observed mass–radius relations for cool dwarfs (Chen et al. 2014).

Zoom In Zoom Out Reset image size

In Figure 16, we now compare isochrones at fixed Z. Note that although Z is the same, there are still element-to-element variations due to the different solar abundance scales adopted by each group. We plot $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;10.0$ isochrones at Z = 0.0001 and Z = 0.03 for PARSEC, Y², DSEP, BaSTI, and MIST. Note that only the MIST models follow the evolution continuously from the helium ignition in the degenerate core (RGBTip) to the CHeB through a series of helium flashes (see also Appendix A). The models are broadly in agreement, though there are some differences in the lower MS and the extent of the CHeB. The former is likely mostly due to differences in the adopted boundary conditions in the models, and the latter is possibly due to differences in the adopted Reimers mass-loss efficiency (BaSTI, PARSEC, and MIST adopt ${\eta }_{{\rm{R}}}\;=\;0.4$, 0.2, and 0.1, respectively, while Y² does not include mass loss).

Zoom In Zoom Out Reset image size

Figure 17 presents a comparison between the MIST, PARSEC, and Lyon models for a $0.3\;{M}_{\odot }$ evolutionary track. As noted above, the PARSEC models now adopt a modified $T\mbox{--}\tau$ relation for low-mass stars, which likely explains the relatively large difference between that model and the others. The Lyon and MIST models largely agree, although several sharp features are noticeable in the Lyon models during the PMS phase that do not appear in the MIST models. The origin of these small differences is unclear to us.

Zoom In Zoom Out Reset image size

In Figure 18 we show the evolution of integrated colors of a simple stellar population for MIST (black), PARSEC/COLIBRI (red; Bressan et al. 2012; Marigo et al. 2013; Rosenfield et al. 2014), and BaSTI (blue; Pietrinferni et al. 2004) isochrones at solar metallicity. The colors are computed by integrating along the isochrone at a given age with weights provided by the Kroupa IMF (Kroupa 2001). They are calculated using the Python bindings²⁰ to the Flexible Stellar Population Synthesis code (FSPS, v2.6; Conroy et al. 2009; Conroy & Gunn 2010). We turn the AGB circumstellar dust option off to enable a more direct comparison between the three models. There are no predictions from the BaSTI models at $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\lt 7.5$ because they only go up to $10\;{M}_{\odot }$ in mass. Note that we used the same bolometric corrections for all three cases so any variation in color is purely due to differences in the isochrones.

Zoom In Zoom Out Reset image size

Overall, the models are in good agreement with each other, especially in B − V, though there are some noticeable differences between the models in other colors. At $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\gtrsim 9$ in FUV − V, the MIST prediction turns over toward bluer colors while the PARSEC/COLIBRI and BaSTI predictions continue to get redder. This qualitative difference is due to the inclusion of the post-AGB and WD phases in the MIST models. In V − K and J − K, the large spikes at young ages ($\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\sim 7$) are due to the onset of the RSG phase from massive stars. This feature appears at a slightly later time in MIST compared to in PARSEC/COLIBRl, which points to the differences in the lifetimes of massive stars in the two databases. The inclusion of rotational mixing in the MIST models may explain the longer MS lifetimes. Finally, significant differences between the three models occur at intermediate ages in redder colors, where TPAGB stars are expected to contribute a significant fraction of the total luminosity. The MIST color predictions fall between the BaSTI and PARSEC predictions. We note that the full luminosity and temperature variations—the actual thermal pulses—are included in the MIST isochrones.

We compare MIST and Geneva (Ekström et al. 2012) evolutionary tracks for a wide range of masses in Figure 19. The Geneva models, shown in pink, also include rotation with $v/{v}_{\mathrm{crit}}\;=\;0.4$ and adopt a similarly low-metallicity solar abundance scale—${Z}_{\odot }\;=\;0.014$ to be exact—with the elemental mixture from Asplund et al. (2005) combined with the Ne abundance from Cunha et al. (2006). At fixed stellar mass, the Geneva models are hotter and more luminous at TAMS, which implies that rotational mixing is more efficient in their models compared to that in the MIST models. As discussed in Section 6.1, efficient rotational mixing gives rise to hotter temperatures and higher luminosities due to larger core sizes and increased μ in the envelope.

Zoom In Zoom Out Reset image size

Some MIST models, such as the $120\;{M}_{\odot }$ star in Figure 19, lose their hydrogen-rich envelope very promptly as they reach the so-called ${\rm{\Omega }}{\rm{\Gamma }}$-limit (Maeder & Meynet 2000). As evident from Equations (21) and (22), massive stars with ${\rm{\Gamma }}\;=\;L/{L}_{\mathrm{Edd}}\to 1$ only require the smallest amount of rotation Ω to receive a large boost in mass-loss rates. As the star evolves, its surface metallicity increases due to a combination of mixing processes and mass loss, and as a result, its surface Rosseland mean opacity increases. This in turn decreases ${L}_{\mathrm{Edd}}$, which makes it easier for a star to experience a large rotational boost. The star then may enter a positive feedback loop where mass loss leads to even more efficient mass loss until it removes all of its envelope and becomes a very compact star almost completely devoid of angular momentum. At the moment, it is not clear whether nature produces such stars, perhaps because of more complex behavior not included in the current 1D models.

Differences in the efficiency of rotational mixing between the MIST and Geneva models are further explored in the left panel of Figure 20, which shows the ratio of MS lifetimes for rotating and non-rotating models as a function of initial mass. This ratio is expected to be greater than unity since rotational mixing channels additional fuel into the core. The solid black, green, blue, and dashed pink lines correspond to the default model at solar metallicity with $v/{v}_{\mathrm{crit}}\;=\;0.4$ and ${f}_{\mu }\;=\;0.05$, solar metallicity model with $v/{v}_{\mathrm{crit}}\;=\;0.6$ and ${f}_{\mu }\;=\;0.05$, solar metallicity model with $v/{v}_{\mathrm{crit}}\;=\;0.4$ and ${f}_{\mu }\;=\;0.01$, and the Geneva model, respectively. Recall that f_μ is the parameter that captures the sensitivity of rotational mixing to the mean molecular weight gradient, such that a small f_μ corresponds to efficient mixing even in the presence of a stabilizing composition gradient.²¹ The default MIST model experiences only a modest enhancement in the MS lifetime. In contrast, the Geneva model experiences an overall $\sim 25\%$ increase in the MS lifetime for stars more massive than $2\;{M}_{\odot }$ (Ekström et al. 2012; Georgy et al. 2013).

Zoom In Zoom Out Reset image size

Although the MIST and Geneva models experience quantitively different amounts of MS lifetime boost, this is not entirely surprising given their different implementations of rotational mixing. Moreover, massive star evolution, regardless of the inclusion of rotation, is highly uncertain and very sensitive to small changes in the input physics. At fixed $v/{v}_{\mathrm{crit}}$, the efficiency of rotational mixing depends sensitively on f_μ. As expected, MS lifetime boost in the MIST models is increased for a higher rotational mixing efficiency via increased rotation velocity or decreased f_μ. The default values ${f}_{c}\;=\;1/30$ and ${f}_{\mu }\;=\;0.05$ are adopted from Heger et al. (2000). This combination, though not unique, is able to reproduce many of the observational constraints such as the high-mass star ratios (see Section 9.3) and observed surface nitrogen enrichment (see Section 9.4). The fact that both MIST and Geneva models broadly reproduce observational constraints in spite of the different lifetime enhancements implies that current observations are not uniquely constraining. For reference, we note that the MS lifetimes for the rotating models in Geneva and MIST agree to within 10%–15% at solar metallicity: for low-mass stars ($\lesssim 1.5\;{M}_{\odot }$), the MS lifetime is shorter in the Geneva models, whereas for higher-mass stars, the MIST models have MS lifetimes that fall between those of non-rotating and rotating Geneva models.

In the right panel, we compare the ratio of MS lifetimes among MIST models with different metallicities. Since the primary mass-loss mechanism for massive stars is strongly metallicity-dependent line-driven winds, rotational mixing becomes more important at low metallicities due to the lowered efficiency of angular momentum loss, as expected.

Relations between mass, radius, luminosity, and temperature provide powerful and fundamental tests of stellar evolution models. In the past two decades, there have been enormous improvements in measuring these quantities to high precision from a variety of techniques, including eclipsing binaries and interferometry (see Torres et al. 2010 for a recent review on this topic).

In Figure 21, we plot $\mathrm{log}(R)$, $\mathrm{log}(L)$, and $\mathrm{log}({T}_{\mathrm{eff}})$ as a function of stellar mass for the DEBCat²² sample, an online catalog of detached eclipsing binaries (DEBs) with well-measured parameters compiled from the literature (Southworth 2015), and a sample of DEBs selected from the literature that was homogeneously reanalyzed by Torres et al. (2010). Note that the Torres et al. (2010) sample appears to show smaller scatter, especially around $\sim 1\;{M}_{\odot }$. We applied a $\mathrm{log}(g)$ cut—$\mathrm{log}(g)\gt 4.1\quad \mathrm{cm}\;{{\rm{s}}}^{-2}$ and $3.4\;\mathrm{cm}\;{{\rm{s}}}^{-2}$ for ${M}_{{\rm{i}}}\gt 1.2\;{M}_{\odot }$ and $\lt 1.2\;{M}_{\odot }$, respectively, as estimated from our model isochrones—to remove evolved stars from the sample of likely MS stars. Furthermore, we removed from the final sample a few conspicuous outliers identified as PMS or RGB stars in the literature. The predicted ZAMS relations for solar metallicity are shown as black solid lines in each of the panels. Since the ages of the stars are unknown, we also show the full range of possible MS values as the gray shaded region. The vertical dashed, dotted, and dot-dashed lines demarcate the initial masses for which ${t}_{\mathrm{MS}}\sim {t}_{\mathrm{Hubble}}$ for $[{\rm{Z}}/{\rm{H}}]\;=\;-1.0$, 0.0, and $+0.5$. Below these masses, we do not expect stars to have evolved off the ZAMS relations. Overall, the observed points fall comfortably within the region bounded by the ZAMS and TAMS relations. However, the observed stars start to deviate from the predicted relations below ${M}_{{\rm{i}}}\lesssim 0.7\;{M}_{\odot }$. In the insets, we zoom in on the low-mass range to show that the models systematically underpredict radii by ∼0.03 dex and overpredict temperatures by ∼0.05 dex, for a total deficit of ∼0.2 dex for the predicted luminosity.

Zoom In Zoom Out Reset image size

There is a well-known discrepancy between observed and predicted effective temperature, radius, and luminosity relations for stars with appreciable convective layers, most notably M dwarfs (e.g., Casagrande et al. 2008; Torres et al. 2010; Kraus et al. 2011; Feiden & Chaboyer 2013; Spada et al. 2013; Torres 2013; Chen et al. 2014). At fixed stellar mass, models tend to predict stars that are 5%–10% hotter and 10%–20% smaller in radius compared to observations. This disagreement is present in both field stars and DEBs, suggesting that this is an effect intrinsic to dwarfs (Boyajian et al. 2012; Spada et al. 2013). However, there are systematic errors of a few percent expected from DEB light-curve analysis due to variations in the spot size and coverage (Morales et al. 2010). A proposed explanation for this mismatch invokes magnetic activity and rotation effects that are not currently modeled accurately (Spruit & Weiss 1986; Morales et al. 2008, 2010; Irwin et al. 2011; Kraus et al. 2011; Feiden & Chaboyer 2012; MacDonald & Mullan 2014; Jackson & Jeffries 2014). Large-scale magnetic fields are thought to both inhibit the upwelling of hot convective bubbles and generate more starspots on the surface (e.g., Feiden & Chaboyer 2012). In order to conserve flux, the stellar radius is inflated, causing a subsequent decrease in the surface temperature. Rotation may play a role since it is believed to generate a dynamo and has been linked to magnetic activity (see Section 3.6.5). Furthermore, the choice of surface boundary conditions in stellar models has a non-trivial effect on the mass–radius relation and the optical color–magnitude diagrams (CMDs) at the lowest masses (e.g., Baraffe et al. 1995; Chabrier et al. 1996; Baraffe et al. 1997; Spada et al. 2013; Chen et al. 2014). When it comes to modeling cool dwarfs, it is especially important to use accurate boundary conditions—such as those computed from atmosphere models, e.g., PHOENIX (Hauschildt et al. 1999) and ATLAS12/SYNTHE (Kurucz 1970, 1993)—in place of simple models that assume gray atmospheres (see Section 3.3).

Low- and intermediate-mass stars shed a nontrivial fraction of their mass via winds during the course of their lifetime, eventually terminating their lives as WDs. Total mass loss integrated over the lifetime directly connects the initial mass to the remnant mass through the IFMR (e.g., Reimers 1975; Weidemann 1977; Renzini & Fusi Pecci 1988; Weidemann 2000). It is an important diagnostic for the cumulative effect of mass loss occurring at various stages of evolution. The expectation is that stars with higher initial masses produce more massive WD remnants (e.g., Claver et al. 2001; Dobbie et al. 2004; Williams & Bolte 2007; Catalán et al. 2008a; Kalirai et al. 2008; Williams et al. 2009; Zhao et al. 2012).

We compare the predicted IFMR to a sample of eight young open clusters, three older open clusters with ages $\gt 1\;\mathrm{Gyr}$, and Sirius B (see Ferrario et al. 2005 and Kalirai et al. 2008 for references therein). It is useful to study clusters of a variety of ages because it allows us to probe a large range of initial masses. Observed initial and final masses for each WD in a cluster are inferred using the following method (see, e.g., Kalirai et al. 2008 for details). A combination of WD spectral analysis and modeling yields both the WD mass (final mass) and cooling age (age since the end of shell helium burning on the TPAGB). The WD progenitor age up to the end of the TPAGB is simply the difference between the total age of the system as estimated from the cluster turnoff and the WD cooling age from the previous step. Finally, stellar evolution models provide the progenitor mass (initial mass) corresponding to the WD progenitor age. Note that initial and final masses are not directly observed but instead are inferred from modeling: the final mass comes from the spectral analysis, while the initial mass depends on stellar evolution theory and CMD analysis.

In Figure 22 we plot the predicted IFMRs for three values of [Z/H] that altogether encompass the metallicities of the systems represented here. Individual measurements within a single cluster have been binned to represent a weighted mean. Overall, the models are in excellent agreement with the data, though there are some notable outliers like NGC 6819 and Sirius B. The low-mass plateau at ${M}_{{\rm{i}}}\lesssim 2\;{M}_{\odot }$ (${M}_{{\rm{f}}}\sim 0.6\;{M}_{\odot }$) is marginally consistent with the peak of the galactic disk WD mass function near $\sim 0.6\;{M}_{\odot }$ (Liebert et al. 2005; Kleinman et al. 2013; Kepler et al. 2015), although a full model that folds in the age and metallicity distributions as well as the IMF weights will be required for a robust comparison against the observed mass function (see also Catalán et al. 2008b).

Zoom In Zoom Out Reset image size

Furthermore, the steep slope predicted around ${M}_{{\rm{i}}}\sim 3\quad -\quad 4{M}_{\odot }$ ($\mathrm{slope}\sim 0.16$) is consistent with the empirical estimate from Cummings et al. (2015), who found a slope of ${M}_{{\rm{f}}}\;=\;(0.163\pm 0.022){M}_{{\rm{i}}}+(0.238\pm 0.071)\quad {M}_{\odot }$ from a combined sample of newly identified WDs in NGC 2099 and the WDs in the Hyades and Praesepe from Kalirai et al. (2014). Interestingly, in agreement with Romero et al. (2015) but in contrast with Marigo & Girardi (2007), our theoretical relations show a clear systematic trend with metallicity above $3\;{M}_{\odot }$. As discussed in Marigo & Girardi (2007), the core can grow/erode through shell burning/third dredge-up (TDU) during the TPAGB, while mass loss more or less determines when the TPAGB phase terminates. Metallicity is predicted to affect all these processes: at low metallicities, TDU and hot bottom burning efficiencies, as well as the core mass at the onset of the first thermal pulse, increase, but mass-loss efficiency is believed to decrease. The measurement of the IFMR as a function of metallicity therefore has great potential for constraining these uncertain evolutionary phases.

It is worth noting that calculations of ${M}_{{\rm{i}}}$ and ${M}_{{\rm{f}}}$ that rely on, e.g., WD spectral analysis and isochrone fitting and the quality of the data vary between cluster to cluster. As noted in Kalirai et al. (2008), the sample is likely affected by small systematic offsets and nonzero field contamination. In particular, precision measurements of ages with the MSTO method in young systems are particularly challenging due to the vertical placement of the MSTO in an optical CMD (note the larger error bars toward younger ages and higher ${M}_{{\rm{i}}}$). Moreover, the thickness of the hydrogen and helium layer assumed in the WD spectral model can also have a non-negligible effect on the inferred cooling age and thus the initial mass (Prada Moroni & Straniero 2002; Catalán et al. 2008b). Catalán et al. (2008b) found that reasonable variations in the envelope thickness can lead to differences as large as $1\;{M}_{\odot }$ for progenitors with ${M}_{{\rm{i}}}\gtrsim 5\;{M}_{\odot }$ but a more modest $\sim 0.1\;{M}_{\odot }$ difference for the lower masses. Other model uncertainties include the assumed core composition. A stringent test of the IFMR should be entirely self-consistent; a single set of isochrones should be used to estimate the "observed" masses ${M}_{{\rm{i}}}$ and ${M}_{{\rm{f}}}$.

We conclude this section with a comparison between the final mass and the core mass at the first thermal pulse in Figure 23. It is useful to consider the core mass before the star experiences significant core growth from the subsequent thermal pulses, because this comparison is devoid of large uncertainties in mass loss, TDU, and hot bottom burning that strongly influence the TPAGB phase (e.g., Wagenhuber & Groenewegen 1998; Weidemann 2000). In the left panel, we show the predicted IFMR from Figure 22 and the core mass at the beginning of the TPAGB phase as a function of initial mass in solid and dashed lines, respectively. The right panel shows the fractional growth in core mass during the TPAGB phase (${M}_{{\rm{f}}}-{M}_{1\mathrm{stTP}}/{M}_{{\rm{f}}}$) as a function of initial mass. Overall, there is considerable growth in core mass occurring during this phase, with a broad peak over $2\mbox{--}3\;{M}_{\odot }$. The location of maximum growth coincides with the peak in TPAGB lifetime as shown in Figure 13. This result is in good agreement with what Bird & Pinsonneault (2011) and Kalirai et al. (2014) found using the Pietrinferni et al. (2004) and PARSEC/COLIBRI (Bressan et al. 2012; Marigo et al. 2013) models, respectively.

Zoom In Zoom Out Reset image size

8.3.1. Star Clusters

In this section we present comparisons with observed CMDs of star clusters. Models shown here include a reddening correction according to the standard ${R}_{{\rm{V}}}\equiv {A}_{{\rm{V}}}/E(B-V)\;=\;3.1$ reddening law from Cardelli et al. (1989). The majority of the systems presented here are metal-rich and we will provide comparisons with metal-poor clusters with non-solar-scaled abundance patterns in Paper II. The CMD comparisons here are by-eye fits to check that the models yield a reasonable agreement. We plan to perform more robust CMD fitting with MATCH (Dolphin 2002) in future work.

M67 (NGC 2682), an intermediate-age (4 Gyr) solar-metallicity open cluster at a distance of ∼800 pc, is a benchmark system for stellar evolution models (Taylor 2007; Sarajedini et al. 2009). In particular, its well-developed Henyey hook on the MSTO is used to calibrate core convective overshoot in low- and intermediate-mass stars (e.g., Michaud et al. 2004; VandenBerg et al. 2006; Magic et al. 2010; Bressan et al. 2012).

Figure 24 shows optical and near-infrared (NIR) CMDs and $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;9.58$ (3.80 Gyr) isochrones with $[{\rm{Z}}/{\rm{H}}]\;=\;0.0$, A_V = 0.2, and $\mu \;=\;9.7$, where μ is the distance modulus. The Sandquist (2004) BV sample (blue points) was carefully selected using proper motion, radial velocity, variability, and CMD-location information to yield cluster members that are most likely to be single stars. The Sarajedini et al. (2009) Two Micron All Sky Survey (2MASS) sample (mauve points) reflects a membership probability cut of $\gt 20\%$.

Zoom In Zoom Out Reset image size

The MS, MSTO morphology, and the RC luminosity are well matched in all three colors. However, the isochrone begins to peel away from the MS ridge line at fainter than $V\sim 17$, which corresponds to ${M}_{{\rm{i}}}\sim 0.7\;{M}_{\odot }$. This is a well-known issue: models that successfully reproduce the NIR colors (e.g., Sarajedini et al. 2009) predict MS colors that are too blue in the optical (e.g., An et al. 2009; Chen et al. 2014). One of the goals of future work is to revisit and address this problem.

Praesepe (M44; NGC 2632) is a young (∼757 Myr) and moderately super solar cluster at a distance of 180 pc, making it one of the nearest open clusters to the Sun (Taylor 2006; Gáspár et al. 2009; Carrera & Pancino 2011). A combination of its rich cluster membership ($N\gtrsim 1000;$ Kraus & Hillenbrand 2007), large proper motion, and proximity makes Praesepe a favorable target for stellar population studies.

In Figure 25, we show 2MASS (mauve points) and UKIDSS Galactic Clusters Survey photometry data (blue points) from Wang et al. (2014) and Boudreault et al. (2012), respectively. The Wang et al. (2014) sample consists of proper-motion-selected cluster members that are likely to be single stars and includes stars with masses as low as $\sim 0.15\;{M}_{\odot }$. The Boudreault et al. (2012) sample was obtained using astrometric and five-band photometric selection criteria and consists mostly of low-mass stars with ${M}_{{\rm{i}}}\lesssim 0.8\;{M}_{\odot }$. We applied an additional cut to remove objects that were flagged as variable stars and/or had $\lt 50\%$ membership probability. We overplot $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;8.8$ (630 Myr) isochrones with $[{\rm{Z}}/{\rm{H}}]\;=\;+0.15$, A_V = 0.08, and $\mu \;=\;6.26$. Overall, the MIST isochrones provide good fits in all three colors.

Zoom In Zoom Out Reset image size

Pleiades (M45) is a young (∼100 Myr) solar-metallicity open cluster at a distance of only 133 pc (Soderblom et al. 2005, 2009; Melis et al. 2014). Like Praesepe, its richness ($N\sim 1400$) and proximity make it a popular choice for testing stellar evolution models. In fact, the tension between observed and predicted CMD locations of K and M dwarfs dates back to Herbig (1962), who proposed non-coeval evolution, i.e., age spread, as a solution to the discrepancy. Theoretical isochrones were systematically offset toward higher luminosities and cooler temperatures in B − V but predicted fainter and hotter stars in redder colors such as V − I (Stauffer et al. 2003; Kamai et al. 2014). A more recent proposed explanation attributes this discrepancy to magnetic activity (e.g., spots) and/or rotation (e.g., Stauffer et al. 2003).

In Figure 26, we show optical and NIR CMDs constructed from the Stauffer et al. (2007) compilation of 2MASS and BVI_C photometry from the literature (mauve points), a sample of newly identified 2MASS candidates by Stauffer et al. (2007) (green points), and the Kamai et al. (2014) sample of proper-motion members from the Stauffer catalog with updated BVI_C photometry (blue points). We overplot $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;8.0$ (100 Myr) isochrones with $[{\rm{Z}}/{\rm{H}}]\;=\;0.0$, A_V = 0.1, and $\mu \;=\;5.62$. Overall, the isochrones are well-matched to the observed CMDs in all filters. However, as seen in M67, the isochrones depart blueward from the MS ridge line in V − I_C and B − V at fainter than $V\sim 14$, corresponding to ${M}_{{\rm{i}}}\sim 0.6\;{M}_{\odot }$.

Zoom In Zoom Out Reset image size

NGC 6791 is one of the most well-known and well-studied open clusters in the Milky Way. Its unusually old age ($\sim 8\;\mathrm{Gyr}$) and high metallicity ($[\mathrm{Fe}/{\rm{H}}]\sim 0.3-0.5$), combined with its rich membership, make it a unique system for studying extreme stellar populations and their chemical evolution (e.g., Bedin et al. 2005; Carney et al. 2005; Origlia et al. 2006; Kalirai et al. 2007; Brogaard et al. 2012). NGC 6791 is particularly suitable for testing the present MIST models given its solar-scaled [α/Fe] abundances.

In Figure 27, we show optical and NIR CMDs in $B-V$ and $V-I$ (mauve points; Brogaard et al. 2012) and in $J-K$ (blue points; Carney et al. 2005). The Brogaard et al. (2012) sample consists of photometry from Stetson et al. (2003) that has been empirically corrected for differential reddening effects. We overplot $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;9.93$ (8.5 Gyr) isochrones with $[{\rm{Z}}/{\rm{H}}]\;=\;+0.47$, A_V = 0.32, and $\mu \;=\;13.1$.

Zoom In Zoom Out Reset image size

The agreement between data and MIST isochrones is generally good, in particular the MSTO and SGB morphologies and the CHeB (red clump) luminosity in all three colors. However, we see the same issue in B − V and V − I as in Figures 24 and 26: MIST isochrones predict colors that are too blue for stars fainter than $V\sim 21$. We return to this point at the end of this section.

Ruprecht 106 is a relatively low mass ($M\sim {10}^{4.8}\;{M}_{\odot }$) globular cluster with a metallicity of $[\mathrm{Fe}/{\rm{H}}]\sim -1.5$ (Kaluzny et al. 1995; Brown et al. 1997; Francois et al. 1997; Dotter et al. 2011; Villanova et al. 2013). Its peculiar properties include the lack of α-enhancement and the absence of abundance spread in light elements, e.g., Na-O anticorrelation, both of which are typical of globular clusters (Carretta et al. 2009). The spread in abundances found in globular clusters has been proposed to be a signature of self-enrichment through multiple generations of stellar populations (Kraft 1994; Gratton et al. 2004; D'Antona & Caloi 2008; Piotto et al. 2012; Gratton et al. 2012). Thus, the modest mass of Ruprecht 106 combined with its "stubby" HB morphology (consistent with the lack of helium variation via self-pollution) suggest that it is an archetypical single-population globular cluster (Caloi & D'Antona 2011; Villanova et al. 2013). It is an excellent choice for testing our low-metallicity models because we can bypass the issue of α-enhancement. We plan to perform many more tests against a large sample of globular clusters with our α-enhanced models in Paper II.

Figure 28 shows the optical CMD from the HST ACS observations in the ${\rm{F}}606{\rm{W}}$ and F814 broadband filters (Dotter et al. 2011). We plot a $\mathrm{log}(\mathrm{Age})\ [\mathrm{year}]\;=\;10.08$ (12.0 Gyr) isochrone with $\mu \;=\;16.7$, A_V = 0.55, and $[{\rm{Z}}/{\rm{H}}]\;=\;-1.50$ in black. The lower MS, SGB, and RGB are very well matched though the model fails to accurately predict the extreme blue extent of the HB.

Zoom In Zoom Out Reset image size

We conclude this section by addressing a recurring issue raised from the CMD comparisons. Our models predict V − I and B − V colors that are too blue for stars below $\lesssim 0.6-0.7\;{M}_{\odot }$. These discrepancies are due to missing and/or inaccurate atomic and molecular line opacity data used in our bolometric corrections. Efforts are ongoing to address these shortcomings.

8.3.2. The Quadruple System LkCa 3

LkCa 3, a quadruple system of PMS stars in the Taurus–Auriga star-forming region (Torres et al. 2013), offers an excellent opportunity to test PMS evolution models. Operating under the assumption that the system is coeval, we expect all four objects to fall on a single-age isochrone. However, there is evidence that early accretion episodes affect the location and evolution of PMS stars on the H-R diagram (e.g., Baraffe et al. 2009; Hosokawa et al. 2011), possibly complicating the comparison to PMS at these young ages. Torres et al. (2013) concluded that the predicted CMD locations of the four components in the Dartmouth models (Dotter et al. 2008) are better matched than those from the Lyon models (Baraffe et al. 2003), though their recently updated models (Baraffe et al. 2015) show good agreement with the observations as well. Their updates include a new solar abundance scale (a combination of Asplund et al. 2009 and Caffau et al. 2011), improved linelists, and recalibrated mixing length parameter for the treatment of convection.

In Figure 29, we show the observed LkCa 3 stars from Torres et al. (2013) with our 1, 1.4, and 3 Myr solar-metallicity isochrones in the V − H versus M_V plane. The observations as reported by Torres et al. (2013) were already corrected for interstellar extinction assuming ${A}_{{\rm{V}}}\;=\;0.31$. The best-fitting isochrone shown in solid black line indicates that the age of the LkCa 3 system is ∼1.4 Myr, consistent with previous estimates.

Zoom In Zoom Out Reset image size

Upper Scorpius (Upper Sco) is one of three subgroups (Upper Centaurus Lupus and Lower Centaurus Crux) in Scorpius-Centaurus (Sco-Cen), the nearest OB association from the Sun with $d\sim 145\;\mathrm{pc}$ (de Zeeuw et al. 1999; Preibisch et al. 2002). The three subgroups altogether constitute a rich environment to study the formation and evolution of massive stars, circumstellar disks, low-mass stars, and brown dwarfs (e.g., Preibisch et al. 2002; Chen et al. 2011; Lodieu 2013). There has been some recent tension in the literature over the ages of these subgroups (Song et al. 2012). In particular, the age of Upper Sco is quoted to be either ∼5 or ∼11 Myr depending on the analysis method and the spectral types of stars used to estimate the age (de Geus et al. 1989; Preibisch et al. 2002; Lodieu et al. 2008; Slesnick et al. 2008; Pecaut et al. 2012; Herczeg & Hillenbrand 2015). This discrepancy poses a problem for other studies that rely on accurate age measurements, e.g., inferred mass functions (Preibisch et al. 2002; Lafrenière et al. 2008; Lodieu 2013).

Pecaut et al. (2012) used isochrones to determine the age of Upper Sco from a kinematic sample of PMS F-type stars. The authors compared luminosities and temperatures derived from Hipparcos and 2MASS photometry to four different sets of models (D'Antona & Mazzitelli 1997; Siess et al. 2000; Demarque et al. 2004; Dotter et al. 2008). They inferred an age of ∼13 Myr—much older than the previous estimates of ∼5 Myr—regardless of the model used in the analysis. This intriguing result prompted the authors to repeat their analysis using a number of other data sets from the literature, namely, the MSTO B-type stars, M supergiant Antares (α Sco), MS A-type stars, and PMS G-type stars. The resulting ages were all consistently older than 5 Myr albeit with a large scatter, and the authors concluded that Upper Sco has a median age of $11\pm 1\pm 2$ (statistical, systematic) Myr.

Figure 30 compares 2, 5, 10, and 20 Myr MIST isochrones at solar metallicity with observations from Preibisch & Zinnecker (1999), Preibisch et al. (2002), and Pecaut et al. (2012), in blue, green, and red symbols, respectively. For spectral types later than G (with the exception of Antares, a RSG), the models yield an age of $\lesssim 5\;\mathrm{Myr}$, which is consistent with the earlier results (e.g., de Geus et al. 1989; Preibisch et al. 2002). For the hotter stars, however, the best-fit model has an older age of ∼10 Myr, consistent with the conclusion from Pecaut et al. (2012).

Zoom In Zoom Out Reset image size

Recent work by Kraus et al. (2015) on UScoCTIO 5, a low-mass spectroscopic binary consisting of nearly identical stars that was recently observed by Kepler as part of the K2 mission (Howell et al. 2014), offered yet another perspective on this issue. Thanks to direct mass and radius measurements, the authors were able to avoid making a comparison in ${T}_{\mathrm{eff}}$ and therefore exclude potential problems with temperature scales as the culprit of this discrepancy. They found that none of the considered models—Lyon (Baraffe et al. 2015), Dartmouth (Dotter et al. 2008), Pisa (Tognelli et al. 2011), Siess (Siess et al. 2000), and PARSEC (Chen et al. 2014)—predicted an age around ∼11 Myr given the system's stellar properties.²³ A simple exercise where they horizontally shifted the evolutionary tracks in the $\mathrm{log}L\mbox{--}\mathrm{log}{T}_{\mathrm{eff}}$ plane to bring the predicted ${T}_{\mathrm{eff}}$ into agreement with the observations yielded an age of ∼11 Myr, consistent with the age derived from hotter stars. As a result, the authors concluded that the low age predicted from low-mass stars is likely problematic and recommended ∼11 Myr instead of ∼5 Myr as the probable correct age for Upper Sco. This again emphasizes the imperative need for more robust and detailed modeling of the PMS and low MS stars, which we plan to explore in future work.

Low- and intermediate-mass stars ($1\;{M}_{\odot }\lesssim {M}_{{\rm{i}}}\lesssim 8\;{M}_{\odot }$, depending on metallicity) ascend the Hayashi track for the second time and enter the AGB phase after they exhaust their central helium supply. During the first part of the AGB phase (the Early AGB; EAGB), the star contains at its center a degenerate carbon and oxygen core surrounded by helium-burning and hydrogen-burning shells. As the star ascends the AGB, the helium-burning shell moves outward until it reaches the hydrogen-rich zone and shuts off. Meanwhile, a thin helium-rich shell starts to grow in mass due to the helium ash raining down from the hydrogen-burning shell, which now dominates the total energy output of the star. The TPAGB phase begins when the helium shell reaches a critical mass and ignites in a thermonuclear runaway as a consequence of thin shell instability (Schwarzschild & Härm 1965). The resulting expansion of the overlying material quenches the hydrogen-burning shell while the helium-burning shell settles into a period of quiescent burning. Next, the outer envelope begins to contract, causing the bottom of the hydrogen-rich shell to heat up and ignite. The helium-burning shell moves outward in mass until it eventually becomes extinguished. The entire cycle repeats as a series of thermal pulses until the star sheds its envelope and becomes a post-AGB star (observationally, a planetary nebula). The entire TPAGB phase lasts approximately 10⁵–10⁶ yr depending on mass and metallicity (see Figure 12), and the majority of that time is spent in the quiescent "interpulse" state.

The AGB phase plays a significant role in the chemical evolution of galaxies due to its rich nucleosynthetic processes coupled with its high typical mass-loss rate, which can be as large as $\sim {10}^{-3}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}$ during the "superwind" phase (Willson 2000; Herwig 2005). In particular, mixing through repeated dredge-up episodes can enrich the surfaces of AGB stars with heavy elements formed from the slow neutron capture process (s-process). In massive (${M}_{{\rm{i}}}\gtrsim 4-8\;{M}_{\odot }$) super-AGB stars, products of hot bottom burning through the CNO, NeNa, and MgAl cycles are transported up to the surfaces as well. Super-AGB stars are interesting in their own right because they occupy the blurry mass boundary within which ONe and ONeMg WDs can form as a consequence of advanced burning in the core (e.g., Doherty et al. 2015). From a stellar population synthesis perspective, the combination of high luminosities and relatively long lifetimes of AGB stars implies that stars in the AGB phase contribute a large fraction to the integrated light in intermediate-age (a few Gyr) galaxies (Frogel et al. 1990; Maraston 2005; Henriques et al. 2011; Melbourne et al. 2012; Conroy 2013; Noël et al. 2013).

Given its importance, it is therefore disconcerting that the AGB phase is still one of the most poorly understood stellar evolutionary stages. This is because a number of very uncertain and complex physical processes such as mixing and mass loss operate simultaneously and contribute significantly to the evolution (see Lattanzio 2007 for an excellent short overview of these issues). As succinctly summarized in Cassisi & Salaris (2013), "AGB stars are fascinating objects, where a complicated interplay between physical and chemical processes takes place; an occurrence that still makes computing reliable stellar models for this evolutionary phase a challenge." Nevertheless, there has been steady progress toward improving our understanding of AGB stars by calibrating the models against optical and NIR photometry and spectroscopy of the AGB population in nearby galaxies, including the Magnellanic Clouds, dwarf spheroidals, and spirals (e.g., Marigo & Girardi 2001, 2007; Girardi & Marigo 2007; Girardi et al. 2010; Boyer et al. 2011; Rosenfield et al. 2014). Since the duration of this phase strongly influences the time evolution of the rest-frame optical and NIR integrated light in intermediate-age stellar populations, it is important to carefully calibrate the models against observations in the local universe so that integrated light from unresolved systems, e.g., high-redshift galaxies, can be accurately interpreted.

8.5.1. AGB Luminosity Functions

We calibrate the AGB phase in our models via carbon star (C star) and total AGB luminosity functions (LFs) in the Magellanic Clouds. C stars, the most evolved subset of AGB stars, are formed when the atmosphere becomes carbon-rich (${\rm{C}}/{\rm{O}}\gt 1$) as a consequence of recurrent TDU episodes (Iben 1983). The observed C star LF is a popular diagnostic used to calibrate the TDU efficiency (Groenewegen & de Jong 1993; Marigo et al. 1999; Marigo & Girardi 2007). The faint cutoff contains information about the minimum initial mass that experiences TDU, the maximum near ${M}_{\mathrm{bol}}\sim -5.0$ is sensitive to the efficiency of mass loss as well as TDU, and the bright end reflects the decreasing numbers due to both mass loss leading to the termination of the TPAGB phase and hot bottom burning that converts carbon to nitrogen in intermediate-mass AGB stars (Cassisi & Salaris 2013).

The observed LFs are constructed using photometric catalogs of cool, evolved stars in the Magellanic Clouds from the Surveying the Agents of Galaxy Evolution (SAGE) SMC and LMC surveys (private communication, M. Boyer; see also Meixner et al. 2006; Gordon et al. 2011; Boyer et al. 2011 for more details). The very wide baseline in wavelength—optical UBVI from the Magellanic Clouds Photometric Survey (MCPS; Zaritsky et al. 2002) all the way out to $160\;\mu {\rm{m}}$ from Spitzer/MIPS—allows for accurate photometric classification of AGB subtypes and identification of contaminants such as unresolved background galaxies, compact H ii regions, and young stellar objects. From the entire catalog, we select only those located within the MCPS footprint to ensure consistency with the model prediction, which folds in star formation histories derived from the MCPS observations. Next, we select AGB (x-, O-, C-, aO-AGB according to the scheme introduced in Boyer et al. 2011) and C stars (x-, C-AGB) to construct cumulative LFs.

To create the predicted LFs, we utilize the FSPS (v2.6; Conroy et al. 2009; Conroy & Gunn 2010). We first convolve the isochrones with star formation and metallicity histories (SMC, Harris & Zaritsky 2004; LMC, Harris & Zaritsky 2009) to generate composite CMDs in various filters, including the effects of circumstellar dust around AGB stars that can strongly influence the flux at longer wavelengths (Villaume et al. 2015). Next, we select the AGB and C stars using the same CMD cuts that were applied to the SAGE observed sample and construct the LFs assuming a Kroupa IMF. As a consistency check, we ensured that the convolution of the adopted star formation histories with the integrated luminosities and stellar masses reproduces the observed integrated light in the NIR and total stellar mass quoted in Harris & Zaritsky (2004) and Harris & Zaritsky (2009) to within a factor of two.

In Figures 31 and 32, we plot the observed cumulative AGB and C star LFs in four different bands for the SMC and LMC, respectively, in thick black lines and the associated Poisson uncertainties in gray bands. We overplot in red the MIST prediction assuming the same CMD cuts as applied to the data. Overall, the MIST models do a reasonably good job, but the models slightly overpredict and underpredict the numbers for the SMC and LMC, respectively. Originally, our goal was to use these LF comparisons as a means of calibrating various input parameters, e,g., mass loss and TDU efficiencies, that influence AGB star lifetimes, luminosities, and the formation of C stars. Along the way, we encountered several factors that have rendered this task more challenging than initially expected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, there are uncertainties in the adopted star formation and metallicity histories from Harris & Zaritsky (2004) and Harris & Zaritsky (2009). These were derived using a different set of isochrones (Padova; Girardi et al. 2002, to be exact) compared to the MIST isochrones used for the AGB LF predictions. There is thus a fundamental inconsistency that can be addressed by reconstructing the LMC and SMC SFHs with MIST isochrones. Moreover, the recovered star formation histories are sensitive to the adopted dust attenuation model and to crowding, which affect the completeness in high-density areas like 30 Doradus. Also important is the recovered metallicity history, which is far more imprecise than the star formation history itself, since the predicted AGB colors and evolution are extremely metallicity sensitive.

Second, since predicted AGB and C stars are selected using the observed CMD cuts, a small mismatch in the locations of the isochrones, especially in color, may strongly influence the comparison with the data, e.g., separation of C stars from O stars (${\rm{C}}/{\rm{O}}\lt 1$). Undertaking the comparison in multiple bands makes this a particularly demanding task because obtaining good agreement across all wavelengths leaves little room for error in each component: star formation and metallicity histories, stellar evolutionary models, bolometric corrections, and AGB circumstellar dust models. A perhaps more straightforward test is to compare the predicted and observed CMDs directly, though the same uncertainties will make this a challenging task as well.

Finally, although Boyer et al. (2011) took great care to ensure a high-fidelity sample, there is most likely a nonzero amount of contamination in the final sample of AGB stars by, e.g., RSGs. By adopting the observed CMD cuts as the selection criteria, we can, to some degree, account for the possibility of contamination from RSGs in the AGB sample and O stars in the C star sample.

To illustrate the challenges of comparing samples based on CMD cuts, we have also computed LFs by identifying AGB and C stars directly in the isochrone according to their evolutionary stages. It is trivial to tag the phase of every star in the predicted CMD because we have all of the necessary evolutionary information, e.g., stellar mass and surface C/O abundance. The MIST model predictions are shown in blue in Figures 31 and 32. There are a few interesting and revealing differences. Interestingly, objects at the bright end of the J- and ${K}_{{\rm{s}}}$-band LFs are absent in the phase-selected MIST models, which suggests that RSG contamination may be important for the bright end of the AGB LFs. The C star LFs show a more dramatic difference. Possible reasons for the discrepancy include inaccurate modeling of surface abundance enrichment through TDU and winds, the current lack of C/O-variable molecular opacities in the envelope in the models, and any deficiencies in the AGB circumstellar dust models. The implementation of low-temperature molecular opacities in MESA, which have been shown to play an important role in AGB evolution (Marigo 2002; Marigo et al. 2003), is a high priority for the MIST project.

We conclude this section by comparing the predicted LFs from the MIST isochrones to those from other widely used isochrones. In Figures 33 and 34, we plot MIST, PARSEC/COLIBRI (Bressan et al. 2012; Marigo et al. 2013; Rosenfield et al. 2014), and BaSTI (Pietrinferni et al. 2004) predictions in red, dark blue, and green, respectively, for J and [3.6]. The LFs were computed with stars selected from CMD cuts. The $\dot{M}$ required for the computation of AGB circumstellar dust effects came directly from the isochrone files for MIST and PARSEC/COLIBRI, while for BaSTI $\dot{M}$ was computed using the Vassiliadis & Wood (1993) AGB mass-loss prescription. We emphasize that all the rest of the input ingredients for constructing the LFs, including the bolometric corrections in FSPS, are identical for the three isochrones showcased here in order to isolate the effects of varying the isochrones alone. Overall, the models are in broad agreement with each other and the observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Section 8.3.1 we used the MSTO morphology of M67 to calibrate the efficiency of core overshoot. Another popular calibration option is to match the observed width of the MS (e.g., Ekström et al. 2012). We check to see if the MSTO-calibrated overshoot efficiency predicts MS width that is consistent with the observed MS width reported by Castro et al. (2014). Following Langer & Kudritzki (2014), the authors performed their analysis on a so-called "spectroscopic H-R diagram." It differs from an ordinary H-R diagram in that it still has $\mathrm{log}({T}_{\mathrm{eff}})$ on the x-axis but a new quantity ${\mathcal{L}}\equiv {T}_{\mathrm{eff}}^{4}/g$ on the y-axis. The main advantage of a spectroscopic H-R diagram is that all of the relevant quantities can be obtained from spectroscopic analysis without having to worry about ambiguities in distance or extinction.

In Figure 35, we plot three lines demarcating the MS region (ZAMS, TAMS, UPPER) from Castro et al. (2014), which are empirical fits to the probability density distribution constructed from a sample of more than 600 stars. The top portion of the TAMS line is missing because there is no clean break from the MS to cooler temperatures in the observed distribution of stars for $\mathrm{log}{\mathcal{L}}/{{\mathcal{L}}}_{\odot }\gtrsim 4$. This continuous distribution could be explained by the inflation of stars approaching the Eddington limit or the presence of helium-burning stars. Castro et al. (2014) compared evolutionary tracks from Ekström et al. (2012) and Brott et al. (2011a) both with and without rotation and found that the ZAMS loci are generally well reproduced by all models except at the highest masses. However, these massive objects could be missing from the observed sample simply due to their rarity or obscuration by their birth clouds. The authors concluded that no model can reproduce the broad MS width at high ${\mathcal{L}}$ and suggested that core overshoot efficiency may be mass dependent.

Zoom In Zoom Out Reset image size

We overplot a series of solar-metallicity MIST evolutionary tracks with masses ranging from 10 to $80\;{M}_{\odot }$ with (solid red) and without (solid blue) rotation in Figure 35. The MIST models also correctly predict the ZAMS line but fail to reproduce the broad MS width at the highest masses. These models suggest that rotation alone cannot explain the full extent of the MS width. A mass-dependent core overshoot efficiency is a topic we plan to explore in future work.

When a high-mass star runs out of hydrogen in its core, it may migrate toward the Hayashi track and become an RSG. For a long time, the observed RSGs were found to be too luminous and cool compared to the predicted RSGs (Massey & Olsen 2003), which was problematic as the region redward of the theoretical Hayashi track represents a "forbidden zone." At fixed metallicity, each point along the Hayashi line corresponds to the coolest possible model in hydrostatic equilibrium.

Levesque et al. (2005) resolved this problem for Galactic stars, demonstrating that the new effective temperature scale computed from improved MARCS atmospheric models yielded much better agreement between the Geneva evolutionary tracks and the observed H-R diagram locations. Shortly after, Levesque et al. (2006) utilized these new models to analyze a sample of RSGs in the Magellanic Clouds and confirmed previous results from Elias et al. (1985) that these RSGs belong to earlier spectral subtypes compared to their galactic counterparts. This is also consistent with the theoretical expectation that the Hayashi line should shift toward warmer temperatures at lower metallicities. Levesque et al. (2006) also found that RSGs in the SMC span a wide range of ${T}_{\mathrm{eff}}$ for a given ${M}_{\mathrm{bol}}$, possibly due to the increased importance of rotational mixing at lower metallicities (see Section 7.3 for more details).

In Figure 36, we compare the MIST evolutionary tracks in black and the sample of observed RSGs from Massey et al. (2009). Only for the top left panel (M31), we show additional tracks at $[{\rm{Z}}/{\rm{H}}]\;=\;0.0$ in gray since the metallicity of M31 is still under debate (Venn et al. 2000; Trundle et al. 2002; Sanders et al. 2012; Zurita & Bresolin 2012). The observed ${L}_{\mathrm{bol}}$ shown here is derived from M_K rather than M_V since the former is less sensitive to extinction. The typical measurement uncertainty in ${T}_{\mathrm{eff}}$ ranges from ∼100 K for the warmest stars to ∼20 K for the coolest stars, and the uncertainty in $\mathrm{log}(L)$ is negligible (∼0.05 dex). We expect the high density of observed stars to coincide with the location of the Hayashi line (Drout et al. 2012). For the LMC and the SMC, ${T}_{\mathrm{eff}}$ and maximum luminosity are both reproduced by the models. For M31 and the MW, the predicted slopes of the RGB tracks are too shallow compared to the observations, but it is still encouraging that no observed RSGs fall in the forbidden zone. We plan to investigate this further in the future.

Zoom In Zoom Out Reset image size

One of the most popular tests of massive-star models is to compare the observed and predicted ratios of stars belonging to different evolutionary stages as a function of metallicity. The ratio of the IMF-weighted sums of phase A and B lifetimes serves as a proxy for the observed number ratio of stars in phases A and B, ${N}_{\mathrm{obs},{\rm{A}}}/{N}_{\mathrm{obs},{\rm{B}}}$:

$$\begin{eqnarray}&&\displaystyle \frac{\int {t}_{{\rm{A}}}(M)\phi (M)\;{dM}}{\int {t}_{{\rm{B}}}(M)\;\phi (M)\;{dM}}\;=\;\displaystyle \frac{{N}_{\mathrm{obs},{\rm{A}}}}{{N}_{\mathrm{obs},{\rm{B}}}},\end{eqnarray} \tag{ 24 }$$

where t is the phase lifetime and ϕ is the IMF weight. This implicitly assumes that the star formation history is constant over the range of ages considered, which is likely a reasonable approximation for massive stars with ${M}_{{\rm{i}}}\gtrsim 10\;{M}_{\odot }$ and MS lifetimes $\lesssim 20\;\mathrm{Myr}$ (but see also Dohm-Palmer & Skillman 2002, where the authors examined the ratio of blue supergiants (BSGs) to RSGs as a function of age in Sextans A).

Here we present three such tests: the ratio of W-R subtypes WC to WN, the ratio of W-R to O-type stars, and the ratio of BSGs to RSGs. We convert metallicities reported in the literature—$\mathrm{log}({\rm{O}}/{\rm{H}})+12$, Z, and [Fe/H]—to a common scale in [Z/H] to enable comparison with the models. There is an estimated ∼0.1 dex uncertainty in our converted [Z/H] values since there are spatial metallicity gradients within the galaxies and variations in the solar abundances adopted by different groups. Note that these models do not include the effects of binary evolution (see Eldridge et al. 2008).

9.3.1. WC/WN

A W-R star is an evolved massive star with little to no hydrogen in its outer layers. It is formed once the star sheds its hydrogen-rich envelope through mass loss, revealing hydrogen-burning products such as helium and nitrogen (WN subtype) and helium-burning products such as carbon and oxygen (WC subtype). Since the predominant mass-loss mechanism in hot massive stars is likely radiative momentum transfer onto metal ions in the atmosphere, mass loss is predicted to increase with metallicity (Vink & de Koter 2005). As a result, the ratio of WC to WN subtypes is expected to increase with increasing metallicity of the environment (Maeder & Conti 1994), which makes this ratio a useful calibrator for metallicity-dependent mass loss in massive star evolutionary models.

There has been a long-standing mismatch between the predicted and observed WC/WN ratios, especially at high metallicities (Meynet & Maeder 2005; Neugent & Massey 2011; Neugent et al. 2012). However, it was unclear whether this disagreement was due to poor models or completeness issues with observations. On the observations front, Neugent et al. (2012) completely revised the WN/WC ratio in M31 by discovering more than 100 new W-R stars with an estimated completeness fraction of ∼0.95. A comparison between the observed WC/WN ratios in M31, M33, SMC, and LMC and the ratios predicted by non-rotating and rotating Geneva models (Meynet & Maeder 2005; Ekström et al. 2012) revealed only a marginal improvement from the new rotating models. Furthermore, they concluded that additional models at different metallicities (full grids of models were only available for two metallicities at that time) were required for a more informative comparison.

We identify W-R stars in our models following the classification scheme introduced in Georgy et al. (2012). We group the WNL and WNE stars (late and early subtypes of WN) as part of the WN stars, exclude the ambiguous WNC stars (WN to WC transition), and include the WO subtype with the WC stars. We emphasize that this theoretical classification scheme—based on the average surface abundances—is technically not equivalent to the classification scheme used by observers who rely on the spectroscopic detection of emission lines (see e.g., van der Hucht 2001).

First, we compute the phase lifetimes for each evolutionary track. For each phase, we sum up the lifetimes for all stellar masses with weights provided by the IMF (Kroupa 2001). The total lifetime for a given phase is a theoretical proxy for the observed number of stars in that phase, which means that we can now take the ratio of WC to WN lifetimes and directly compare to the observations.

In the left panel of Figure 37, we plot the predicted WC-to-WN ratio as a function of metallicity in a solid black line. The red circles are observations from Neugent & Massey (2011) and Neugent et al. (2012), and the blue diamond point is the observed ratio computed by Georgy et al. (2012), estimated in the 3 kpc radius volume around the solar neighborhood (see their Section 4.4 for references therein). Although the MIST model is currently unable to reproduce the observed trend with metallicity, the predicted ratios are in agreement with the observed values to within a factor of 2–3. We plan to improve this further in future work.

Zoom In Zoom Out Reset image size

9.3.2. W-R/O

9.3.3. BSG/RSG

Rotation has been proposed as a viable mechanism for enriching the surfaces of stars (Heger et al. 2000; Meynet & Maeder 2000; see also Sections 3.6.4 and 3.7.3) by inducing extra mixing and enhancing mass loss.²⁵ Models with rotation generally predict a surface enrichment of helium and nitrogen along with a concomitant depletion of carbon and oxygen during the MS as the products of the CNO cycle get dredged up to the surface through rotational mixing (Maeder & Meynet 2000; Yoon & Langer 2005). For this reason, observed surface abundances have been used to calibrate the efficiency of rotational mixing in the models (f_μ and f_c; Pinsonneault et al. 1989; Heger et al. 2000; Brott et al. 2011a).

Following Ekström et al. (2012), we check that our rotating models are able to match the range of observed surface nitrogen enrichment in galactic O- and B-type stars. Figure 38 shows the predicted surface nitrogen enrichment midway through the MS (${X}_{{\rm{c}}}\sim 0.5;$ red) and at TAMS (blue) as a function of initial mass. Models with and without rotation are plotted in solid and dashed lines, respectively. The red shaded box corresponds to the mean surface nitrogen excess observed in a sample of galactic MS B-type stars with masses below $20\;{M}_{\odot }$, and the blue shaded box on top corresponds to the maximum observed excess. These observed numbers come from Table 2 of Maeder & Meynet (2012), which is a compilation of data from Gies & Lambert (1992), Kilian (1992), Morel et al. (2008), and Hunter et al. (2009). Overall, our models are in excellent agreement with the observed range of nitrogen enrichment on the surfaces of B-type MS stars, and in marginal agreement with the maximum observed excess. They are also broadly in agreement with the predictions from the Geneva models (see Figure 11 of Ekström et al. 2012), though there is a noticeable difference for stars below $\sim 2\;{M}_{\odot }$. This is due to the inclusion of magnetic braking effects in the Geneva models; stars with masses below $\sim 1.7\;{M}_{\odot }$ experience extra surface nitrogen enrichment due to enhanced mixing induced by large shear in the outer layers, though the authors caution that this effect may be overestimated in their current implementation. In the MIST models, we turn off rotation for stars with ${M}_{{\rm{i}}}\lt 1.2\;{M}_{\odot }$ and gradually increase the rotation rate from ${v}_{\mathrm{ZAMS}}/{v}_{\mathrm{crit}}\;=\;0.0$ to 0.4 over the mass range ${M}_{{\rm{i}}}\;=\;1.2\mbox{--}1.8\;{M}_{\odot }$. As discussed in Section 3.5, the purpose of this ramping scheme is to compensate for the absence of magnetic braking in MESA, which is important for low-mass stars with appreciable convective envelopes. At low masses where the MS lifetimes are long and rotational mixing is inefficient or non-existent, the predicted surface nitrogen abundances actually decrease during the MS due to diffusion. At higher masses without the inclusion of rotation, the predicted nitrogen enrichment during the MS is zero. Additional nitrogen enhancement measurements in stars with ${M}_{{\rm{i}}}\lt 10\;{M}_{\odot }$ would provide a valuable constraint on the models.

Zoom In Zoom Out Reset image size

In the left panel of Figure 39, we show 7, 16, and $40\;{M}_{\odot }$ solar-metallicity evolutionary tracks with (solid) and without (dashed) rotation in the $\mathrm{log}(g)\mbox{--}\mathrm{log}({T}_{\mathrm{eff}})$ plane. In the right panel, we show the evolution of surface nitrogen abundance for the same three models. The symbols correspond to O- and B-type galactic dwarfs and supergiants (Takeda & Takada-Hidai 2000; Villamariz et al. 2002; Villamariz & Herrero 2005; Crowther et al. 2006; Morel et al. 2008; Searle et al. 2008; Przybilla et al. 2010). To limit the comparison to O- and B-type stars, we exclude the A, F, and Cepheid stars from the Takeda & Takada-Hidai (2000) sample. We adopt the symbol and color scheme from Figure 12 of Ekström et al. (2012) to aid comparison with their figure. As in Ekström et al. (2012), we do not compare surface rotation velocities because only $v\mathrm{sin}i$ is known for most observed stars. The left panel suggests that the observed stars have initial masses ranging roughly between 7 and $40\;{M}_{\odot }$ and that the observed $\mathrm{log}(g)$ values for the dwarfs (open symbols) are in good agreement with the MS location. The right panel demonstrates that models can broadly reproduce the range of observed surface nitrogen abundances. A fair number of observed points fall below the initial N/H ratio in the models, which Ekström et al. (2012) suggest is due to abundance variations in the birth cloud.

Zoom In Zoom Out Reset image size

We also note that for some samples consisting exclusively of slow rotators (e.g., the Morel et al. 2008 sample of B stars, whose $v\mathrm{sin}i$ values range from 10 to $60\;\mathrm{km}\;{{\rm{s}}}^{-1}$), our default rotating models may be inappropriate for a direct comparison. The observed rotational distribution function is quite broad (Huang et al. 2010), and also the observed stars were most likely born with a range of initial metallicities, which would increase scatter in the observed nitrogen abundances. Moreover, there are physical mechanisms that have not been taken into account in our models, such as mass loss and gain from binary mass transfer (see, e.g., Eldridge et al. 2008; de Mink et al. 2009, 2013), and typical measurement uncertainties are quite large: 1000 K for ${T}_{\mathrm{eff}}$, $0.1\mbox{--}0.2$ dex for $\mathrm{log}(g)$, and 0.2 dex for the N/H ratio. This is meant to be a coarse-grained comparison to assess whether or not the models can reasonably reproduce the range of observed surface nitrogen abundances. Finally, we note that this is still a controversial topic: there are known slow/fast rotators with/without surface nitrogen enrichment, and some authors even find that de-projected velocities show no statistically significant correlation with surface nitrogen abundances (e.g., Hunter et al. 2008; Brott et al. 2011b; Rivero González et al. 2012; Bouret et al. 2013; Aerts et al. 2014). We defer a more detailed comparison and discussion of these issues to future work.