This is the first of a series of papers presenting the Modules for Experiments in Stellar Astrophysics MESA Isochrones and Stellar Tracks MIST project a new comprehensive set of stellar evolutionary tracks and isochrones computed using MESA a stateoftheart opensource 1D stellar evolution package. In this work we present models with solarscaled abundance ratios covering a wide range of ages 5slant logAge yearslant 10.3 masses 0.1slant MMSUB SUBslant 300 and metallicities 2.0slant ZHslant 0.5. The models are selfconsistently and continuously evolved from the premain sequence PMS to the end of hydrogen burning the white dwarf cooling sequence or the end of carbon burning depending on the initial mass. We also provide a grid of models evolved from the PMS to the end of core helium burning for 4.0slant ZHlt 2.0. We showcase extensive comparisons with observational constraints as well as with some of the most widely used existing models in the literature. The evolutionary tracks and isochrones can be downloaded from the project website at A hrefhttpwaps.cfa.harvard.eduMISThttpwaps.cfa.harvard.eduMISTA.
['stars: evolution', 'stars: general', 'stars: interiors', 'Astrophysics - Solar and Stellar Astrophysics']
Mesa Isochrones and Stellar Tracks MIST. I. Solarscaled Models
{'MESA ISOCHRONES AND STELLAR TRACKS (MIST). I: SOLAR-SCALED MODELS': 'JIEUN CHOI 1 , AARON DOTTER 2 , CHARLIE CONROY 1 , MATTEO CANTIELLO 3 , BILL PAXTON 3 , AND BENJAMIN D. JOHNSON 1 Draft version May 2, 2016', 'ABSTRACT': 'This is the first of a series of papers presenting the Modules for Experiments in Stellar Astrophysics (MESA) Isochrones and Stellar Tracks (MIST) project, a new comprehensive set of stellar evolutionary tracks and isochrones computed using MESA, a state-of-the-art open-source 1D stellar evolution package. In this work, we present models with solar-scaled abundance ratios covering a wide range of ages (5 ≤ log(Age) [yr] ≤ 10 . 3), masses (0 . 1 ≤ M / M glyph[circledot] ≤ 300), and metallicities ( -2 . 0 ≤ [Z / H] ≤ 0 . 5). The models are self-consistently and continuously evolved from the pre-main sequence (PMS) to the end of hydrogen burning, the white dwarf cooling sequence, or the end of carbon burning, depending on the initial mass. We also provide a grid of models evolved from the PMS to the end of core helium burning for -4 . 0 ≤ [Z / H] < -2 . 0. We showcase extensive comparisons with observational constraints as well as with some of the most widely used existing models in the literature. The evolutionary tracks and isochrones can be downloaded from the project website at http://waps.cfa.harvard.edu/MIST/.  \nKeywords: stars: evolution, stars: general', '1. INTRODUCTION': "Stars are ubiquitous objects. Individually, they are both hosts to rich exoplanet systems and progenitors of some of the most spectacular transients in the distant Universe. As an ensemble, they are both cosmic engines that transformed the state of the early Universe as well as fossils bearing clues of galaxy formation and evolution. The interpretation of these systems and phenomena hinges on our understanding of stellar physics and well-calibrated stellar evolution models across a wide range of masses, metallicities, and evolutionary states.  \nOn the observational front, large quantities of precise data from recent and ongoing space- and ground-based programs have initiated an explosive growth in fields such as asteroseismology (e.g., CoRoT ; Baglin et al. 2006, Kepler ; Gilliland et al. 2010), time-domain astronomy (e.g., Palomar Transient Factory ; Law et al. 2009, Pan-STARRS ; Kaiser et al. 2010), galactic archaeology (e.g., APOGEE ; Zasowski et al. 2013, GALAH ; De Silva et al. 2015), and resolved stellar populations (e.g., the Hubble Space Telescope (HST) program PHAT ; Dalcanton et al. 2012). Moreover, future surveys, e.g., LSST (Ivezic et al. 2008), and next-generation observatories such as JWST and Gaia will provide an unprecedented volume of high-quality data whose analysis demands uniform models covering all relevant phases of stellar evolution.  \nDecades since the golden era of stellar astrophysics (e.g., Burbidge et al. 1957; Böhm-Vitense 1958; Schwarzschild 1958; Henyey et al. 1964; Paczy'nski 1970), the field has enjoyed a renaissance in recent years, largely due to technological advances in both computing and observational astronomy. Improvements in computers and numerical algorithms have resulted in a tremendous speedup in solving the nonlinear, coupled differential equations of stellar structure and evolution. Another important factor was the availability of increasingly precise tabulated opacities, nuclear reaction rates, and equations of state. Accordingly, a large number of stellar evolution models have been published to tackle a wide variety of problems in astrophysics. Studies of old, low-mass stellar populations in globular clusters and quiescent galaxies have relied on models such as BaSTI (Pietrinferni et al. 2004), DSEP (Dartmouth; Dotter et al. 2008), GARSTEC (Weiss & Schlattl 2008), Lyon (Baraffe et al. 1998, 2003, 2015), Padova/PARSEC (Girardi et al. 2002; Marigo et al. 2008; Bressan et al. 2012), Y 2 (Yi et al. 2001; Kim et al. 2002; Yi et al. 2003; Demarque et al. 2004), Victoria-Regina (VandenBerg et al. 2006), and more. On the other hand, studies of young, massive stellar populations in clusters and starforming environments have made use of e.g., Geneva (Ekström et al. 2012; Georgy et al. 2013), STARS (Eggleton 1971; Pols et al. 1995; Eldridge & Tout 2004), and STERN (Brott et al. 2011a; Köhler et al. 2015) stellar evolution models.  \nIn order to fully exploit these new and upcoming datasets, we have set out to construct stellar evolution models within a single, self-consistent framework using Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015), a popular open-source 1D stellar evolution package. 4 MESA is well-suited for this purpose due to its flexible architecture and its capability to self-consistently model the evolution of different types of stars to advanced evolutionary stages in a single continuous calculation. Furthermore, its  \nthread-safe design enables parallel computing, which greatly reduces the computation time and makes the large-scale production of models feasible.  \nThe paper is organized as follows. Section 2 gives an overview of MESA, focusing on the high-level architecture of the code and its time step and spatial mesh controls. Section 3 reviews the treatment of the relevant physical processes as implemented in a 1D framework, all of which is summarized in Table 1. Calibration of the model against the properties of the Sun is discussed in Section 4, and a short overview of the model outputs, the method for constructing isochrones, and the treatment of bolometric corrections is presented in Section 5. Section 6 presents an overview of the properties of the models, and Section 7 features comparisons with some of the most widely used existing databases. Sections 8 and 9 separately present comparisons with data for low-mass ( glyph[lessorsimilar] 10 M glyph[circledot] ) and high-mass stars ( glyph[greaterorsimilar] 10 M glyph[circledot] ). Finally, Section 10 concludes the paper with a discussion of caveats and future work.  \nThis is the first of a series of papers presenting MESA Isochrones and Stellar Tracks (MIST), a new set of stellar evolutionary tracks and isochrones. In this paper, we present models with solar-scaled abundances covering a wide range of ages, masses, phases, and metallicities computed within a single framework. We will subsequently present models with non-scaled-solar abundances, including e.g., alpha-enhanced, carbon-enhanced metal-poor, and CNONa-enhanced for modeling globular clusters, in Paper II.  \nWe define some conventions and assumptions adopted in the paper. We use M i throughout to refer to the initial stellar mass of a model. Both Z and [Z / H] are used to refer to metallicities, where Z is the metal mass fraction. For the models presented in this paper, [Z / H] = [Fe / H] since we adopt solar-scaled abundances. The [Z / H] notation assumes the Asplund et al. (2009) protosolar birth cloud bulk metallicity, not the current photospheric metallicity, as the reference value (see 3.1 for more details). Magnitudes are quoted in the Vega system unless noted otherwise. Where necessary, we adopt a Kroupa initial mass function (IMF) (Kroupa 2001). Lastly, in accordance with the XXIX th IAU resolutions B2 and B3, 5 we adopt the following two conventions. First, we use the following nominal values to express stellar properties in solar units: M glyph[circledot] = 1 . 988 × 10 33 g, L glyph[circledot] = 3 . 828 × 10 33 erg s -1 , R glyph[circledot] = 6 . 957 × 10 10 cm, and T eff , glyph[circledot] = 5772 K. Formally, the IAU published these values in SI units, but we report them in cgs units here to be consistent with the convention adopted in this work. Second, the zero point of the absolute bolometric magnitude scale is set by enforcing that M bol = 0, which corresponds exactly to L · = 3 . 0128 × 10 35 erg s -1 . This zero point was chosen such that the nominal solar luminosity L glyph[circledot] has an M bol , glyph[circledot] = 4 . 74 mag, a value commonly adopted in the literature.", '2. THE MESA CODE': "MESA 6 is an open-source stellar evolution package that is undergoing active development with a large user base worldwide (Paxton et al. 2011, 2013, 2015). Its 1D stellar evolution module, MESAstar, has been thoroughly tested against existing stellar evolution codes and databases, including BaSTI/FRANEC (Pietrinferni et al. 2004), DSEP (Dartmouth; Dotter et al. 2008), EVOL (Blöcker 1995; Herwig  \n2004; Herwig & Austin 2004), GARSTEC (Weiss & Schlattl 2008), Lyon (Baraffe et al. 1998, 2003, 2015), KEPLER (Heger et al. 2005), and STERN (Petrovic et al. 2005; Yoon & Langer 2005; Brott et al. 2011a). Its highly modular and therefore flexible infrastructure combined with its robust numerical methods enable its application to a wide range of problems in computational stellar astrophysics, from asteroseismology to helium core flash in low-mass stars, as well as the evolution of giant planets, accreting white dwarfs (WD), and binary stars.  \nHere we provide an overview of some of the features in MESAstar, namely its time step controls, adaptive mesh refinement, and parallelization. We refer the reader to Paxton et al. (2011, 2013, 2015) for more detailed information.  \nMESAstar simultaneously solves the fully coupled Lagrangian structure and composition equations using a Newton-Raphson solver. The required numerics (e.g., matrix algebra, interpolation) and input physics (e.g., opacity, mass loss) are organized into individual thread-safe 'modules,' each of which is an independently functional unit that generates, tests, and exports a library to the general MESA libraries directory. This modular structure is unique among stellar evolution programs. One of its main advantages is that experimentation with different available physics or even implementation of new physics is easy and straightforward. The user input is given at runtime via the inlist file, which contains the user's choice for parameters for input physics, time step, mesh, and output options. The run\\_star\\_extras.f file, a Fortran module that is compiled at runtime, allows the user to introduce new routines that hook into the source codes in order to adapt MESA to the problem of interest. Examples include the introduction of new physics routines, modification of model outputs, and customization of time step adjustments.  \nChoosing an appropriate time step throughout various stages of stellar evolution is critical to the accurate evolution of a model. It must be both sufficiently small to allow convergence and sufficiently large to carry out evolution calculations in a reasonable amount of time. In MESA, a new time step is first proposed using a scheme based on digital control theory (Söderlind & Wang 2006). Next, the proposed time step undergoes a series of tests to check if the hypothetical changes to various properties of the model (e.g., nuclear burning rate, luminosity, central density) in a single time step are adequately small, as excessively large changes can lead to convergence or accuracy issues in subsequent evolution.  \nAt the beginning of each evolution step, MESAstar checks whether or not a spatial mesh adjustment is necessary. Similar to time step controls, there is a trade-off between having sufficiently small cells to properly resolve physical processes locally while avoiding an unnecessarily fine grid that is computationally expensive. Remeshing involves the splitting and merging of cells, and each remesh plan aims to minimize the number of cells affected in order to reduce numerical diffusion and improve convergence. At the same time, the remesh plan must meet the criteria for the tolerated cellto-cell changes in relevant variables, which are specifiable by the user in addition to the basic variables, e.g., mass, radius, and pressure. For instance, cells near a convective boundary might be split in order to better resolve its location, while cells well within the convective zone might be merged if they are sufficiently similar in, e.g., composition and temperature. We refer the reader to the appendix for temporal and spatial convergence tests.  \nMESA is optimized for parallelization and uses OpenMP 7 to carry out computations in parallel. Since Paxton et al. (2011), MESA's performance has been greatly improved mainly due to the implementation of a new linear algebra solver-the single-most computationally expensive component-that is compatible with multicore processing. As demonstrated in Figure 48 of Paxton et al. (2013), many key components of MESAstar, such as the linear algebra solver and the evaluation of the nuclear reaction network, closely obey the ideal scaling law.  \nFor this work, we use MESA version v7503 compiled with GNU Fortran version 4.9.3 installed as part of the MESA SDK version 245. 8", '3. ADOPTED PHYSICS': 'In this section we review the relevant physics adopted in the models and their implementation in MESA. Readers who are interested in the most salient points can skip to Table 1, which presents a summary of the adopted physics. For the effects of varying some key physics ingredients, see Section 6.2.', '3.1. Abundances': "The accurate determination of solar abundances has been an ongoing effort for decades. Within the last decade, there has been a systematic downward revision of the solar metallicity from Z glyph[circledot] ∼ 0 . 02 (e.g., Anders & Grevesse 1989) to Z glyph[circledot] glyph[lessorsimilar] 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 glyph[circledot] , protosolar = 0 . 0142, not Z = Z glyph[circledot] , 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 glyph[circledot] , photosphere-the current photospheric abundances are determined by spectroscopy-but rather to clarify what a 'solar metallicity model' entails.  \nTo compute X and Y for an arbitrary Z , we use the following formulae:  \nY p = 0 . 249 (1)  \nY = Y p + ( Y glyph[circledot] , protosolar -Y p Z glyph[circledot] , protosolar ) Z (2)  \nX =1 -Y -Z . (3)  \nThis approach adopts a primordial helium abundance Y 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 glyph[circledot] , protosolar = 0 . 2703 (Asplund et al. 2009), such that ∆ Y / ∆ Z = 1 . 5. Once Y is computed for a desired value of Z , X is trivial to compute. We assume the isotopic ratios 2 H / 1 H = 2 × 10 -5 and 3 He / 4 He = 1 . 66 × 10 -4 from Asplund et al. (2009).  \nTable 1 Summary of the Adopted Physics.", '3.2.1. Equation of State (EOS)': 'The EOS tables in MESA are based on the OPAL EOS tables (Rogers & Nayfonov 2002), which smoothly transition to the SCVH tables (Saumon et al. 1995) at lower temperatures and densities. The extended MESA EOS tables cover X = 0 . 0 , 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8, and 1.0, and Z = 0 . 0 , 0 . 02, and 0.04. At higher metallicities, MESA switches to the MacDonald EOS tables (MacDonald & Mullan 2012) for Z = 0 . 2 and 1.0, which, unlike the HELM EOS tables (Timmes & Swesty 2000) used in the earlier versions of MESA, allow for partially ionized species. The HELM and PC tables (Potekhin & Chabrier 2010) are used at temperatures and densities outside the range covered by the OPAL + SCVH + MacDonald tables, and assume full ionization. The EOS tables in MESA also cover the late stages of WD cooling, during which the ions in the core crystallize, although the current MIST models do not reach such conditions. The low-mass models are evolved until Γ , the central plasma interaction parameter or Coulomb coupling parameter, reaches 20 (see Section 5.1), and crystallization occurs at Γ ≈ 175 for pure oxygen.', '3.2.2. Opacities': 'MESA divides the radiative opacity tables into two temperature regimes, high (log T glyph[greaterorsimilar] 4) and low (log T glyph[lessorsimilar] 4), and treats them separately. This system allows for the user to choose, for the low temperature opacities, between either Ferguson et al. (2005) or Freedman et al. (2008) with updates to ammonia opacity from Yurchenko et al. (2011) and the pressureinduced opacity for molecular hydrogen from Frommhold et al. (2010). The high temperature opacity tables come from either OPAL (Iglesias & Rogers 1993, 1996) or OP (Seaton 2005). The OPAL tables are split into two types, Type I and Type II: Type I tables are used for 0 . 0 ≤ X ≤ 1 . 0 -Z and 0 . 0 ≤ Z ≤ 0 . 1 for a fixed abundance pattern; Type II tables are optionally available which allow for enhanced carbon and  \noxygen abundances in addition to those already accounted for in Z , covering 0 . 0 ≤ X ≤ 0 . 7 and 0 . 0 ≤ Z ≤ 0 . 1. Type II opacities are particularly important for helium burning and beyond. The electron conduction opacity tables are originally based on Cassisi et al. (2007), but they have been extended to cover temperatures up to 10 10 K and densities up to 10 11 . 5 g cm -3 (Paxton et al. 2013).  \nWe use the Ferguson et al. (2005) low temperature tables and the OPAL Type I tables, then gradually switch to the OPAL Type II opacities starting at the end of hydrogen burning, smoothly interpolating between X < 10 -3 and X < 10 -6 . Note that we use the Asplund et al. (2009) protosolar mixture where available to be consistent with our choice of the solar abundance scale, but the opacity tables implemented in MESA were computed for the Asplund et al. (2009) photospheric abundances.', '3.2.3. Nuclear Network': 'We import the nuclear reaction rates directly from the JINA REACLIB database, 9 a compilation of the latest reaction rates in the literature (Cyburt et al. 2010). For example, the 15 N(p, α ) 12 C reaction rate comes from Angulo et al. (1999), while the tripleα reaction rate comes from Fynbo et al. (2005). We use the JINA reaction rates for p-p chains, cold and hot CNO cycles, tripleα process, α -capture up to 32 S, Ne-Na and Mg-Al cycles, and C/O burning. We adopt the mesa\\_49.net nuclear network in MESA.  \nThe nuclear network tracks and solves for the abundances of the following 52 species: n, 1 H, 2 H, 3 He, 4 He, 7 Li, 7 Be, 9 Be, 10 Be, 8 B, 12 C, 13 C, 13 N, 14 N, 15 N, 14 O, 15 O, 16 O, 17 O, 18 O, 17 F, 18 F, 19 F, 18 Ne, 19 Ne, 20 Ne, 21 Ne, 22 Ne, 21 Na, 22 Na, 23 Na, 24 Na, 23 Mg, 24 Mg, 25 Mg, 26 Mg, 25 Al, 26 Al, 27 Al, 27 Si, 28 Si, 29 Si, 30 Si, 30 P, 31 P, 31 S, 32 S, 33 S, 34 S, 40 Ca, 48 Ti, 56 Fe. The three heaviest elements 40 Ca, 48 Ti, and 56 Fe are our modifications to the default mesa\\_49.net network and are inert  \nspecies-they do not participate in any nuclear reactionsthat are thus only affected by mixing and diffusion processes.  \nElectron screening is included for both the weak and strong regimes. We use the default option in MESA which computes the screening factors by extending the classic Graboske et al. (1973) method with that of Alastuey & Jancovici (1978), and adopting plasma parameters from Itoh et al. (1979) for strong screening. 10', '3.3. Boundary Conditions': '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 is a multitude of options that ranges from simple analytic approximations to tables based on full atmospheric structure models.  \nThe simplest choice, simple\\_photosphere , uses the Eddington T ( τ ) relation to obtain T :  \nT 4 ( τ ) = 3 4 T 4 eff ( τ + 2 3 ) , (4)  \nwhere T eff is calculated directly from the MESA interior model. Similarly, P is computed as follows:  \nP = τ g κ [ 1 + P 0 κ τ L M 1 6 π cG ] . (5)  \nThe 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 0 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 superEddington luminosities. We adopt P 0 = 2 for this work.  \nWith 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.,  \nIn 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 3 with α MLT = 1 . 25 following the implementation of convection as outlined in Castelli et al. (1997). 11 We employed the latest atomic and molecular line lists from R. Kurucz, including molecules important for cool stars such as TiO and H2O. Individual models are calculated for log( g ) = 0 to 5 and T eff = 2500 K to 50 , 000 K for [Z / 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.  \nwhere T = T eff. However, for cooler dwarfs a more sensible choice is to set the boundary condition deeper in the atmosphere, i.e., τ = 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.  \nFor our grid of models, the tau\\_100\\_tables option is used for 0 . 1-0 . 3 M glyph[circledot] , photosphere\\_tables is used for 0 . 6-10 M glyph[circledot] , and simple\\_photosphere is used for 16-300 M glyph[circledot] . To facilitate a smooth transition between different regimes, we run both tau\\_100\\_tables and photosphere\\_tables for 0 . 3-0 . 6 M glyph[circledot] and photosphere\\_tables and simple\\_photosphere for 10-16 M glyph[circledot] . 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 eff glyph[greaterorsimilar] 10 4 K.', '3.4. Diffusion': "Microscopic diffusion and gravitational settling of elements are essential ingredients in stellar evolution models of lowmass 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 HR 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 1 H, 3 He, 4 He, 16 O, and 56 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.  \nHowever, 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 over-depletion of helium and metals in the outer envelopes of stars with M i > 1 . 4 M glyph[circledot] , 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  \nradiation 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 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.  \nTo 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.  \nSince 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 postMS 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.", '3.5. Rotation': "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.  \nRotation 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).  \nOur 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 GoldreichSchubert-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 ( glyph[lessorsimilar] 1 . 2 M glyph[circledot] ), particularly those with slow and moderate rotation rates, have strong differential rotation profiles at ZAMS.  \nSince magnetic braking is not presently modeled in MIST, we do not include rotation for stars with M i ≤ 1 . 2 M glyph[circledot] 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 glyph[circledot] , the rotation rate is gradually ramped up from 0 to the maximum value of v ZAMS / v crit = Ω ZAMS / Ω crit = 0 . 4, where v crit and Ω 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), 12 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-1 . 5 M glyph[circledot] (Wolff & Simon 1997; Canto Martins et al. 2011) reveals that our ramping scheme produces velocities that are reasonable ( ∼ 10-25 km s -1 during the main sequence for 1 . 3-1 . 35 M glyph[circledot] ) even in the absence of magnetic braking.  \nCurrently, 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).  \nChemical 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. Mixing Processes': '3.6.1. Convection  \nMixing Length Theory (MLT), whose modern implementation in stellar evolution codes was pioneered by BöhmVitense (1958), describes the convective transport of energy in the stellar interior. There is a vexing yet crucial free parameter of order unity, α MLT, that determines how far a fluid parcel travels before it dissolves into the background, l MLT, in units of the local pressure scale height, H P ( l MLT = α MLT H P). In other words, it parametrizes how efficient convection is, because a large α MLT means that the parcel travels a large distance before it deposits its energy into the ambient medium.  \nConvective mixing of elements is treated as a timedependent 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. 13 In addition to α 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 α MLT is discussed in Section 4.1  \nClassically, the location of the convective region is determined using the Schwarzschild criterion, which implies that a region is convectively stable if  \n∇ T < ∇ ad , (6)  \n∇ T < ∇ L (7)  \nwhere ∇ T is the local background temperature gradient (in practice, it is set to the radiative temperature gradient ∇ rad) and ∇ 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, ∇ µ . In this case, a region is convectively stable if  \n∇ L ≡∇ ad -χ µ χ T ∇ µ (8)  \nχ µ ≡ [ ∂ ln( P ) ∂ ln( µ ) ] ρ , T (9)  \nχ T ≡ [ ∂ ln( P ) ∂ ln( T ) ] ρ , µ , (10)  \nwhere the thermodynamic derivatives χ µ and χ 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 it  \ncan account for both the observed properties of asymptotic giant branch (AGB) and post-AGB stars (Herwig 2000), the observed MS width (e.g, Schaller et al. 1992), and the main sequence turn off (MSTO) morphology in clusters such as M67 (e.g., Magic et al. 2010).  \nThe 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  \nThere 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).  \nD OV = D 0 exp ( -2 z H v ) ; H v = f ov H P , (11)  \nwhere H v is the velocity scale height of the overshooting convective elements at the convective boundary, f ov is a free parameter that essentially determines the efficiency of overshoot mixing, H P is the local pressure scale height, and D 0 is the diffusion coefficient in the unstable region 'near' the convective boundary (more specifically at a depth of f 0 , ov H P from the convective boundary). For simplicity, we adopt two sets of ( f ov, f 0 , 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 , ov is set to 0 . 5 f ov. The temperature gradient in the overshoot region is assumed to be equal to the radiative gradient as in the step overshoot approach.  \nOther 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 as well as to 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 ov , step of H 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-1 . 2 M glyph[circledot] 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 P formally diverges as r → 0, the size of the overshoot region, f ov , step H P, also diverges. Thus, to avoid the excessive growth of the convective core for lowmass MS stars, the common solution is to gradually ramp up the overshoot efficiency from M 1 ∼ 1 M glyph[circledot] to M 2 ∼ 1 . 7 M glyph[circledot] , with no convective overshoot below M 1 and maximum effi-  \nabove M 2. 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 1 and M 2.  \nWe 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 ov , env, is chosen along with α 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 ov , env for simplicity.  \nWe adopt a modest overshoot efficiency f ov , core = 0 . 016 for the core which is roughly equivalent to f ov , 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 vs. 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.", '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 HR 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).  \nAlternately, 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 3 He( 3 He,2p) 4 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).  \nIn 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):  \nD sc = α sc ( K 6 CP ρ ) ∇ T -∇ ad ∇ L -∇ T , (12)  \nD th = α th ( 3 K 2 CP ρ ) -χ µ χ T ∇ µ ∇ T -∇ ad , (13)  \nwhere K is the radiative conductivity, CP is the specific heat at constant pressure, and α sc is a dimensionless free-parameter. Similarly, the diffusion coefficient for thermohaline mixing is computed following Ulrich (1972) and Kippenhahn et al. (1980):  \nwhere α th is a dimensionless number that describes the aspect ratio of the mixing blobs or 'fingers' (a large α th corresponds to slender fingers).  \nAs summarized in Paxton et al. (2013), the range of α sc and α th adopted by various authors spans several orders of magnitude, partially due to differences in the implementation in various codes. We adopt α 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 α 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 α 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': 'In MESA, the transport of both chemicals and angular momentum arising from rotationally induced instabilities are treated in a diffusion approximation (Endal & Sofia 1978; Zahn 1983; Pinsonneault et al. 1989; Heger et al. 2000; Yoon &Langer 2005) in place of the alternative diffusion-advection approach (Maeder & Meynet 2000; Meynet & Maeder 2000; Eggenberger et al. 2008; Potter et al. 2012). The five rotationally induced instabilities included in our models are: dy-  \nnamical shear instability, secular shear instability, SolbergHøiland (SH) instability, ES circulation, and GSF instability (Heger et al. 2000; Paxton et al. 2013). Of these, ES circulation and shear instabilities have the largest impact on the evolution of a rotating star. We refer the reader to Heger et al. (2000) and Maeder & Meynet (2000) for excellent overviews of these phenomena.  \nOnce the diffusion coefficients for these rotational mixing processes are computed, they are combined with the diffusion coefficients for convection, semiconvection, and thermohaline. This grand sum enters the angular momentum and abundance diffusion equations solved at each time step. There are two free parameters in this implementation, first introduced by Pinsonneault et al. (1989) to model the Sun: fc , a number between 0 and 1 that represents the ratio of the diffusion coefficient to the turbulent viscosity, which scales the efficiency of composition mixing to that of angular momentum transport; and f µ , a factor that encodes the sensitivity of rotational mixing to the mean molecular weight gradient, ∇ µ . A small fc corresponds to a process that transports angular momentum more efficiently than it can mix material, and a small f µ means that rotational mixing is efficient even in the presence of a stabilizing ∇ µ . We adopt fc = 1 / 30 and f µ = 0 . 05 following Pinsonneault et al. (1989), Chaboyer & Zahn (1992), and Heger et al. (2000). As we demonstrate in Section 9.4, these parameters produce surface nitrogen enhancements that are in reasonable agreement with the observations.', '3.6.5. Magnetic Effects': 'There is a growing body of evidence that our understanding of internal angular momentum transport in stars is not complete. For example, the observed spin rates of WDs and neutron stars (Heger et al. 2005; Suijs et al. 2008) and the angular velocity profiles inferred from asteroseismic observations of red giants (Eggenberger et al. 2012; Cantiello et al. 2014) cannot be reproduced with models that only include rotational mixing from hydrodynamic instabilities and circulations.  \nMagnetic fields are also observed near the stellar surface, which are thought to be either of fossil origin (e.g., Braithwaite & Spruit 2004) or generated through dynamo operating in convective zones in the outer layers of low-mass stars (e.g., Brandenburg & Subramanian 2005). However, as discussed in Section 3.5, magnetic braking due to the coupling between winds and surface magnetic fields is not yet included in MESA.  \nSpruit-Tayler (ST) dynamo is a mechanism for the amplification of seed magnetic fields in radiative stellar interiors in the presence of differential rotation (Spruit 2002). Stellar models including torques from ST dynamo fields can reproduce the flat rotation profile in the solar interior (Mestel & Weiss 1987; Charbonneau & MacGregor 1993; Eggenberger et al. 2005) and the observed spin rates of WDs and neutron stars (Heger et al. 2005; Suijs et al. 2008), although they still cannot explain the slow rotation rates of cores in red giants (Cantiello et al. 2014). The chemical mixing and the transport of angular momentum due to internal magnetic fields are not included in our models, though this is implemented in MESA following KEPLER (Heger et al. 2005) and STERN (Petrovic et al. 2005). We note that the very existence of the ST-dynamo loop is still under debate (Braithwaite 2006; Zahn et al. 2007), and there are ongoing efforts to understand the role of angular momentum transport via magnetic fields in radiative stellar regions (e.g. Rüdiger et al. 2015; Wheeler et al. 2015).', '3.7. Mass Loss': '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 HR diagram. We note that the total mass loss rate is capped at 10 -3 M glyph[circledot] yr -1 in all models to prevent convergence problems.', '3.7.1. Low Mass Stars': 'Mass loss for stars with masses below 10 M glyph[circledot] 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:  \n˙ M R = 4 × 10 -13 η R ( L / L glyph[circledot] )( R / R glyph[circledot] ) ( M / M glyph[circledot] ) M glyph[circledot] yr -1 , (14)  \n˙ M B = 4 . 83 × 10 -9 η B ( L / L glyph[circledot] ) 2 . 7 ( M / M glyph[circledot] ) 2 . 1 ˙ M R η R M glyph[circledot] yr -1 , (15)  \nwhere η R and η 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).  \nFor simplicity, we turn on Reimers mass loss at the beginning of the evolution, but a negligible amount of mass loss occurs throughout the MS ( ∼ 10 -13 M glyph[circledot] yr -1 for a solar metallicity 1 M glyph[circledot] star). Once core helium is depleted, the mass loss rate is chosen to be max[ ˙ M R, ˙ M B] at any given time. We adopt η R = 0 . 1 and η 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': 'For hot and luminous massive stars, mass loss is thought to arise from the absorption of ultraviolet photons by metal ions in the atmosphere, resulting in a preferentially outward momentum transfer from the absorbed photons to the plasma (line-driven winds; Lucy & Solomon 1970; Castor et al. 1975). For our models, mass loss for stars above 10 M glyph[circledot] uses a combination of radiative wind prescriptions, collectively called the Dutch mass loss scheme in MESA, inspired by Glebbeek et al. (2009). There is an option for an overall scaling factor η Dutch analogous to η R and η B for the low-mass stars, but we adopt η Dutch = 1 . 0. For prescriptions that include metallicity-scaling, we retain the reference Z glyph[circledot] adopted by each author. We expect this difference to have a negligible  \neffect relative to the large overall uncertainties in mass loss prescriptions.  \nWe now describe our mass-loss scheme for high-mass stars in each region of the HR diagram:  \n- 1. For T eff > 1 . 1 × 10 4 Kand X surf (surface hydrogen mass fraction) > 0.4, the mass loss prescription from Vink et al. (2000, 2001) is used, which is appropriate for the early phases of the evolution prior to the stripping of the hydrogen-rich envelope. The authors computed mass loss rates using a Monte Carlo radiative transfer code, taking into account multiple scatterings and assuming that the loss of photon energy is coupled to the momentum gain of the wind. The Vink mass loss rate is specified by five parameters; M , L , T eff, v ∞ / v esc, and Z surf.  \nFor 2 . 75 × 10 4 < T eff < 5 × 10 4 K,  \n˙ M V, hot = 10 . L / L .  \n10 10 . 92[log( T eff / 4 × 10 K)] × ( Z surf / Z glyph[circledot] ) 0 . 85  \n-6 697 × ( 10 5 glyph[circledot] ) 2 194 × ( M / 30 M glyph[circledot] ) -1 . 313 × ( v ∞ / v esc 2 . 0 ) -1 . 226 × ( T eff / 4 × 10 4 K) 0 . 933 × -4 2 . (16)  \nThe ratio of terminal flow velocity to the escape velocity increases with metallicity following v ∞ / v esc = 2 . 6( Z surf / Z glyph[circledot] ) 0 . 13 .  \nFor 1 . 1 × 10 4 < T eff < 2 . 25 × 10 4 K, 14  \n( M / 30 M glyph[circledot] ) -1 . 339 × ( v ∞ / v esc 2 . 0 ) -1 . 601 × T / 4 1 . 07 Z / Z 0 . 85 , (17)  \n˙ M V, cool = 10 -6 . 688 × ( L / 10 5 L glyph[circledot] ) 2 . 210 ×  \n( eff 4 × 10 K) × ( surf glyph[circledot] )  \nwhere v ∞ / v esc = 1 . 3( Z surf / Z glyph[circledot] ) 0 . 13 .  \nFor 2 . 25 × 10 4 ≤ T eff ≤ 2 . 75 × 10 4 K, either ˙ M V, hot or ˙ M V, cool is adopted depending on the exact position of the so-called bi-stability jump, a phenomenon in which ˙ M increases with decreasing T eff due to the recombination of metal lines:  \nT eff, jump = 61 . 2 + 2 . 59log 〈 ρ 〉 , (18)  \nwhere 〈 ρ 〉 corresponds to the characteristic wind density at 50% of the terminal velocity of the wind. The successful predictions of the mass loss rates and the bistability jump near T eff ∼ 2 . 5 × 10 4 K for the Galactic and SMC O-type stars have made the Vink prescription a popular choice among massive star modelers (but see the discussion below).  \n- 2. Once the star reaches T eff > 10 4 K and X surf < 0 . 4, it is formally identified as a Wolf-Rayet (WR) star and we switch over to the Nugis & Lamers (2000) empirical mass loss prescription which depends strongly on luminosity and chemical composition:  \n˙ M NL = 10 -11 ( L / L glyph[circledot] ) 1 . 29 X 1 . 7 surf Z 0 . 5 surf M glyph[circledot] yr -1 . (19)  \n14 In Vink et al. (2001), mass loss rates were computed for T eff ≥ 1 . 25 × 10 4 K, but the prescription is extended down to T eff = 1 . 1 × 10 4 K in MESA.  \nThis formula has been shown to reproduce the properties of a large sample of Galactic WR stars (WN, WC, and WO subtypes) that have well-constrained stellar and wind parameters (Nugis & Lamers 2000). With Equation 19, we are able to reproduce the observed ratio of WC to WN subtypes as a function of metallicity (see Section 9.3).  \n- 3. For all stars with T eff < 10 4 K, 15 including stars in the red supergiant (RSG) phase, we utilize the de Jager et al. (1988) empirically derived wind prescription:  \n˙ M dJ = 10 -8 . 158 ( L / L glyph[circledot] ) 1 . 769 T -1 . 676 eff M glyph[circledot] yr -1 . (20)  \nAlthough a quantitative model of mass loss in RSGs does not exist yet, it is believed that the main mechanism is dust-driven outflows. The low temperatures and pulsations in the outer layers lead to the condensation of dust at large radii, which is then driven out due to radiation pressure on grains (Mauron & Josselin 2011). In a recent work comparing different wind prescriptions with mass loss rates estimated from a sample of RSGs in the Galaxy and in the Magellanic Clouds, Mauron & Josselin (2011) found that the de Jager rates agree to within a factor of 4 with most estimates derived from the 60 µ m flux. Furthermore, the authors concluded that the de Jager wind recipe performs better overall compared to more recent prescriptions, though they recommended an additional metallicity scaling ( Z surf / Z glyph[circledot] ) 0 . 7 . For simplicity, we adopt the original de Jager prescription available in MESA.  \nA recent review by Smith (2014) explores some of the shortcomings of these prescriptions, one of which is that they fail to account for the clumpiness and inhomogeneity in outflows. 16 For instance, mass loss inferred from H α emission or free-free continuum excess that assumes a homogeneous wind (e.g., de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990) is believed to overestimate the true rate by a factor of 2 to 3. However, the reduced mass loss rate corrected for clumpiness may be problematic for the formation of WR stars, which requires the removal of their hydrogen-rich envelopes. Eruptive episodic mass loss episodes and/or binaries are likely to play a role (e.g., Smith & Owocki 2006; Yoon & Cantiello 2010; Sana et al. 2012), but neither phenomenon can be realistically captured in a simple recipe for implementation in a 1D stellar evolution code. Other forms of enhanced mass loss rates include super-Eddington winds when the stellar luminosity exceeds Eddington luminosity (Gräfener et al. 2011; Vink et al. 2011).', '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  \n- 15 In MESA, mass loss rate for 10 4 < T eff < 1 . 1 × 10 4 is computed by smoothly transitioning between the low temperature prescription (de Jager) and high temperature prescription (Vink or Nugis & Lamers).  \nfunction of surface angular velocity Ω as follows:  \n˙ M ( Ω ) = ˙ M (0) ( 1 1 -Ω / Ω crit ) ξ , (21)  \nwhere ˙ 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 Ω crit is the critical angular velocity at the surface:  \nΩ 2 crit = ( 1 -L L Edd ) GM R 3 . (22)  \nThe Eddington luminosity, L Edd, is a mass-weighted average over the optical depth τ between 1 and 100. For stars close to the Eddington limit, Ω 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 4 .", '4.1. Solar Calibration': "As mentioned in Section 3.6.1, it is customary to calibrate the mixing length parameter, α 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, α 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. 17 A global χ 2 value is computed by summing over the log L , log R , surface composition, R cz (the location of the base of the convection zone), and δ c s (model -observed sound speed) terms with non-uniform user-defined weights. This process is repeated until χ 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.  \nAlthough the model log L and log T 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).  \nTable 2 summarizes the solar calibration results and Figure 1 shows δ c s / c 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 solarcalibrated 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.  \nThe 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).  \n/circledot  \nTable 2 Solar calibration results.  \na XXIX th IAU resolutions B2 and B3  \nb Asplund et al. (2009)  \nc Basu & Antia (2004)  \n- d Basu & Antia (2004)  \nFigure 1. The fractional error in sound speed compared to helioseismology observations from Rhodes et al. (1997) for the present model (black) and two Serenelli et al. (2009) models each with Grevesse & Sauval 1998 (GS98; red) and Asplund et al. 2009 (AGSS09; blue) meteoritic abundances.  \n<!-- image -->  \nWhile 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.  \nThe 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 ∼ 0 . 6 R glyph[circledot] . 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 further 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 deter-  \nFigure 2. The evolution of surface lithium abundance relative to hydrogen, A( 7 Li) = log( N 7 Li / N 1 H ) + 12, as a function of time for several 1 M glyph[circledot] models. Each pair of numbers in parentheses corresponds to the two envelope overshoot efficiency parameters f ov , env and f ov , env , 0, respectively. The solid black line represents the fiducial model that adopts the solar-calibrated envelope overshoot parameters (0.0174, 0.0087) and α MLT = 1 . 82. 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 dotted-dashed lines show additional variations in input physics: the red line correspond to α MLT = 1 . 7, while the orange and green lines are models that include PMS rotation with v / v crit = 0 . 01 and 0.10, respectively. 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). The current models cannot simultaneously reproduce both observed lithium abundances.  \n<!-- image -->  \nnations 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).  \nSeveral 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 ∼ 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).  \nThe 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.  \nAlthough the widely adopted practice is to fix the solarcalibrated α MLT across all masses, evolutionary phases, and abundances, there has been a recent effort to map out α MLT as a function of log g and T 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 α MLT and T ( τ ) boundary condition from Trampedach et al. (2014) in their stellar evolution calculations and found that varying α MLT has a small effect on the evolution and surface properties. We adopt a constant value of α MLT = 1 . 82 for the present work, but discuss the implications of this assumption in more detail in Section 8.3.1  \nIn summary, we adopt solar-calibrated α MLT and convective overshoot efficiency in the envelope f ov , env ( f 0 , ov , env is fixed to 0 . 5 f ov , 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 glyph[circledot] , protosolar = 0 . 0142.", '4.2. Lithium Depletion': 'Lithium is a very fragile element that is burned via proton capture at temperatures as low as 2 . 5 × 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.  \nIn Figure 2 we show the evolution of surface lithium abundance relative to hydrogen, A( 7 Li) = log( N 7 Li / N 1 H ) + 12, as a function of time for several 1 M glyph[circledot] 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 ov , env and f 0 , ov , 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 α 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 dotted-dashed lines show additional variations in input physics: the red line correspond to α MLT = 1 . 7, while the orange and green lines are models  \nStandard 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 α 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.  \nthat include PMS rotation with v / v crit = 0 . 01 and 0.10, respectively. As expected, a lower α 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 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, as well as the exploration of a variable α MLT to mimic the effects of non-standard physics are planned for follow-up investigations in the near future.', '5. MODEL OUTPUTS AND BOLOMETRIC CORRECTIONS': 'In this section, we provide an overview of the two principal model outputs: evolutionary tracks and isochrones, which can be downloaded from http://waps.cfa.harvard.edu/MIST/. There are both theoretical and observational isochrones available. We offer both packaged models for download and a web interpolator that generates models with user-specified parameters.', '5.1. Ages, Masses, Phases, Metallicities': "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 M glyph[circledot] to 300 M glyph[circledot] for a total of glyph[greaterorsimilar] 100 models, and the ages of isochrones cover logAge = 5 to logAge = 10 . 3 in 0.05 dex steps. For M i ≤ 0 . 7 M glyph[circledot] , the models are terminated at TAMS, i.e., central 1 H abundance drops to 10 -4 . For a 0 . 7 M glyph[circledot] star at Z glyph[circledot] , this limit is typically reached at an age > 35 Gyr. For M i > 0 . 7 M glyph[circledot] , the models are either evolved through the WD cooling phase ('lowmass' 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 boundaryfurther 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 ¯ Z 2 e 2 / aikbT , where ¯ Z is the average ion charge, e is the electron charge, ai is the mean ion spacing, kb 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 Γ ≈ 175; Paxton et al. 2011). The 'high-mass' type modelsstars that are sufficiently massive to burn carbon-are stopped when the central 12 C abundance drops to 10 -2 . The metallicity ranges from [Fe / 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 ≤ [Z / H] < -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.  \nWe 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 shortlived 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.", '5.2. Isochrone Construction': "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 socalled 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 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 glyph[circledot] and 30 M glyph[circledot] 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.  \nAs described in Section 3.3, we use three different boundary conditions depending on the initial mass of the model ( τ = 100 tables for 0 . 1-0 . 3 M glyph[circledot] , photosphere tables for 0 . 6-10 M glyph[circledot] , and simple photosphere for 16-300 M glyph[circledot] ). To facilitate a smooth transition from one regime to another, we run both τ = 100 and photosphere tables for 0 . 3-0 . 6 M glyph[circledot] and photosphere tables and simple photosphere for 10-16 M glyph[circledot] . For every mass in the transition regime, the two EEP tracks are blended with a smooth weighting function to create a hybrid EEP track:  \nw = 0 . 5 [ 1 -cos ( π M ∗ -M 1 M 2 -M 1 )] , (23)  \nwhere M 1 and M 2 are the transition masses (e.g., M 1 = 0 . 3 M glyph[circledot] and M 2 = 0 . 6 M glyph[circledot] for the transition from τ = 100 to photosphere tables).  \nTable 3 Primary EEPs and corresponding evolutionary phases.  \na i.e., zero age horizontal branch; ZAHB for low-mass stars.  \nb terminal age horizontal branch; TAHB.  \n<!-- image -->  \nFigure 3. Left: an example 1 M glyph[circledot] evolutionary track in the equivalent evolutionary point (EEP) format, with the locations of the primary EEP points marked by colored circles. The gray box marks the zoomed-in region shown in the right panel. Right: a zoomed-in view of the track.  \n<!-- image -->  \nFigure 4. Same as Figure 3 but now for a 30 M glyph[circledot] evolutionary track. Note the different primary EEP point following the core helium burning phase. 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.  \n<!-- image -->  \nTo generate an isochrone at age t 0 with all of the EEP tracks now in hand, we first cycle through all masses and construct a piecewise monotonic function between M i and t for each EEP point. 18 Next, we interpolate to obtain M i( t 0). Once we have M 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 i, e.g., L ( M 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).  \nA 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  \nFigure 5. Isochrones from log(Age) [yr] = 5-10 in various physical quantities. Top left: isochrones in the log( g ) -log( T eff) plane. Top middle: asteroseismic isochrones in the ∆ ν -log( T eff) plane, where ∆ ν corresponds to the large frequency separation for p-modes. Note that ∆ ν is computed from the scaling relations. Top right: asteroseismic isochrones in the ν max -log( T eff) plane, where ν max corresponds to the frequency of maximum power. Note that ν max is computed from the scaling relations. Bottom left: isochrones in the log( Tc )-log( ρ c ) plane with the dashed, dotted-dashed, and solid lines showing thresholds for hydrogen, helium, and carbon ignition. The discontinuous feature at log(Age) [yr] ≥ 10 is due to the ignition of helium in a degenerate core, i.e, helium core flash at the tip of the RGB. Bottom middle: isochrones in the log( L / L Edd)M i plane, where L Edd is the Eddington luminosity. Bottom right: isochrones in the ˙ M -M i plane, where ˙ M is the mass loss rate. In the last two figures, isochrones become steadily truncated at the high-mass end as more stars evolve and die.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \ndifferent 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.", '5.3. Available Model Outputs': 'The published tracks and isochrones include a wealth of information, ranging from basic parameters such as log( L ), log( T eff), log( g ), and surface abundances of 19 elements to more detailed quantities such as β ≡ P gas / P total and asteroseismic parameters (the full list of available parameters is available on the project website). 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 log( g ) -log( T eff), ∆ ν -log( T eff), and ν max -log( T eff) planes. ∆ ν and ν max are asteroseismic quantities that correspond to the large frequency separation for pmodes 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 ν 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 log( Tc )-log( ρ c ) plane, where the dashed, dotted-dashed, and solid lines show thresholds for hydrogen, helium, and carbon ignition. The  \n10 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 log( L / L Edd), the ratio of total luminosity to Eddington luminosity, and ˙ 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.', '5.4. Bolometric Corrections': '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 (Conroy et al. 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, H2, H2O, SiO, C2, SiH, MgH, CrH, CaH, FeH, CO, NH, VO, and OH. We have also computed model atmosphere and spectra for carbon stars  \nTable 4 Current list of photometric systems.  \nFigure 6. Coverage of the synthetic spectra grids used to derive bolometric corrections. Example isochrones at log(Age) [yr] = 7 . 0, 7.5, 8.5, 9.0, and 10.0 are overplotted in gray for reference. The coverage in log g and log T eff is the same for all metallicities.  \n<!-- image -->  \nwith C / O = 1 . 05 over the range 2400 < T eff < 4700 K and -1 . 0 < log g [g cm -3 ] < 0 . 5. Our carbon star spectra agree well with the models of Aringer et al. (2009). The primary differences arise at > 2 µ m as our models do not currently include the important molecules C3, HCN, and C2H2. 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 6 , 000 ≤ T eff ≤ 50 , 000 K from Koester (2010). These are supplemented by a set of blackbody spectra with 200 , 000 ≤ T eff ≤ 1 million K. The coverage of the different synthetic libraries is shown in Figure 6. Example isochrones at log(Age) [yr] = 7 . 0, 7.5, 8.5, 9.0, and 10.0 are overplotted for reference.  \nBolometric corrections are computed from the synthetic  \nspectra following Equation 1 of Girardi et al. (2008). As noted in the Introduction, we adopt as a zeropoint L · = 3 . 0128 × 10 35 erg s -1 to define M bol = 0. This is equivalent to adopting solar values of M bol , glyph[circledot] = 4 . 74 and L glyph[circledot] = 3 . 828 × 10 33 erg s -1 . The bolometric corrections include a range of extinction values, as characterized by both AV and RV , following the extinction curve of Cardelli et al. (1989). We provide AV = 0 to 6 with RV = 3 . 1, though other RV values can be made upon request. We emphasize that Z of the bolometric correction is matched to the surface Z (and surface C / O ratio where relevant) for each star along the isochrone.  \nThe 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.', '6.1. Tracks and Isochrones': "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 M glyph[circledot] to 300 M glyph[circledot] , the ages of isochrones cover logAge = 5 to logAge = 10 . 3, and the metallicity ranges from [Fe / 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.  \nFigure 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 crit = 0 . 4 by default. Figure 8 shows the effect of rotation on both the evolutionary tracks and isochrones, where the rotating and nonrotating 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  \n<!-- image -->  \nFigure 7. An example solar metallicity grid of stellar evolutionary tracks (left) and isochrones (right) covering a wide range of stellar masses, ages, and evolutionary phases.  \n<!-- image -->  \n<!-- image -->  \nFigure 8. Similar to Figure 7 but now showing the effect of rotation on both the evolutionary tracks (left) and isochrones (right). Models with and without rotation are shown in solid and dashed lines, respectively. For rotating models, solid-body rotation with v ZAMS / v crit = Ω ZAMS / Ω crit = 0 . 4 is initialized at ZAMS. The PMS phase is not shown in the left panel for display purposes.  \n<!-- image -->  \nweight, 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 eff that includes both the gravitational and centrifugal terms. Since T eff ∝ g 1 / 4 eff for a star with a radiative envelope and | g eff | < | g | , T eff is thus lower in a rotating system. Generally speaking, rotation increases T eff in massive stars where rotational mixing operates efficiently, and it decreases T eff in low-mass stars as well as at all ZAMS locations where the centrifugal effect dominates.  \nThe top two panels of Figure 9 show a series of 2 . 2 M glyph[circledot] (left) and 6 . 4 M glyph[circledot] (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  \non 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 glyph[circledot] 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 [Z / H] > 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.  \nIn order to illustrate the range of metallicities available in the current MIST models, we show 10 Gyr isochrones at [Z / H] = -4, -3, -2, -1, 0, and 0 . 5 in Figure 11. For clarity,  \nFigure 12. Top: phase lifetimes as a function of initial mass for [Z / H] = -0 . 25, 0.0 and + 0 . 25. High mass stars are not included because they do not go through the same set of evolutionary phases featured here. 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 Γ = 20. Bottom: same as above but now showing the ratio of phase lifetimes to the MS lifetime instead.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 9. A series of 2 . 2 M glyph[circledot] (left) and 6 . 4 M glyph[circledot] (right) evolutionary tracks for three metallicities. The tracks become hotter and more luminous with decreasing metallicity due to a lower Rosseland mean opacity. In the bottom panels, we zoom in on the TPAGB (left) and CHeB (right) to further highlight the effects of metallicity on these evolutionary phases.  \n<!-- image -->  \n<!-- image -->  \n/circledot  \nFigure 10. The number of thermal pulses (TPs) executed by each model as a function of mass at a number of metallicities. There are two notable features: first, the maximum occurs at around 2 M glyph[circledot] independent of metallicity, and two, 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 [Z / H] > 0.  \nwe 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., [Z / H] glyph[lessorsimilar] -2.  \nIn the top panel of Figure 12 we show phase lifetimes as a function of initial mass for [Z / H] = -0 . 25, 0.0 and + 0 . 25 in solid, dotted-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 Γ = 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  \nFigure 11. Isochrones at 10 Gyr over a wide range of metallicities. For display purposes, we omit the phases beyond 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., [Z / H] glyph[lessorsimilar] -2.  \n<!-- image -->  \n11  \n<!-- image -->  \n/circledot  \nFigure 13. Left: similar to Figure 12 but now showing the cumulative age as a function of mass for [Z / H] = 0 . 0. The gray box marks the zoomed-in region shown in the right panel. Right: zooming in on a mass range that shows a critical transition from degenerate helium ignition to quiescent helium ignition. This produce a non-monotonic mass-age relationship for these phase. Note the linear scale in y . This figure is adapted from Figure 1 of Girardi et al. (2013).  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \nare 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 ∼ 2 M glyph[circledot] (see the discussion below).  \nThe 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 [Z / H] = 0 . 0. In the right panel, we zoom in on a particularly interesting mass range around 2 M glyph[circledot] , where there 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 RGBmass 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.", '6.2. Non-standard Models': '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.  \nIn the top left panel, we show two sets of isochrones which  \nare identical except for their surface boundary conditions. The three sets of ages correspond to log(Age) [yr] = 9 . 0 , 9 . 6, and 10.2. The solid lines are isochrones with our default implementation of boundary conditions; τ = 100 and photosphere tables from the ATLAS12/SYNTHE atmosphere models plus the simple Eddington gray T ( τ ) approximation for the hottest stars. The dashed lines are models with the Eddington gray T ( τ ) 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 eff on the RGB is smaller but important, amounting ∼ 50-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 [Z / H] = -1 and on the lower MS at [Z / H] = + 0 . 3.  \nIn the bottom left panel, we highlight the importance of core convective overshoot. The three sets of ages shown are log(Age) [yr] = 7 . 0 , 7 . 5, and 8.0. The solid lines are isochrones with the default core overshoot parameter f ov , core = 0 . 016 in the exponential diffusive overshoot formalism. The dashed and dotted lines represent decreased and increased  \nIn the top right panel, we show three sets of isochrones at log(Age) [yr] = 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 η R = 0 . 1 , η B = 0 . 2, and η 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 ˙ 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).  \n4  \nFigure 14. Solar metallicity isochrones at a number of ages showing the effects of varying input physics. The fiducial model is shown in black solid lines in all four panels. Note different x and y axes in each panel. Top left: variations in the adopted boundary conditions for log(Age) [yr] = 9 . 0 , 9 . 6, and 10.2. Top right: variations in the efficiency of mass loss for log(Age) [yr] = 6 . 7 , 7 . 4, and 9 . 0. The set of three numbers corresponds to η R, η B, and η Dutch. Bottom left: variations in the efficiency of overshoot in the core parameterized by f ov , core for log(Age) [yr] = 7 . 0 , 7 . 5, and 8.0. Bottom right: variations in the efficiency of convection parameterized by α MLT for log(Age) [yr] = 7 . 0 , 8 . 0 , 8 . 7, and 10.0.  \n<!-- image -->  \n( f ov , 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.  \nmissing from the majority of current stellar evolution models (but see Feiden & Chaboyer 2013). For example, α 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).', '7. COMPARISONS WITH EXISTING DATABASES': "In this section we compare the MIST models to several popular stellar evolution databases in the literature. Due to differences in the choice of input physics and their implementations in the codes, an apples-to-apples comparison is challenging. We stress that this is neither a comprehensive review of all published models in the literature nor a thorough and detailed comparison between different databases. Instead, we aim to provide the reader with a general impression of how the new MIST models compare to several widely used models. We refer the reader to the MESA instrument papers (Paxton et al. 2011, 2013, 2015) for closer comparisons at the level of the codes themselves and their evolutionary track outputs. Our goal here is to compare at the level of databases, which reflects the net effect of many different choices for input physics.  \nIn the bottom right panel, we show isochrones with different values of the mixing length parameter α MLT. The four ages shown are log(Age) [yr] = 7 . 0 , 8 . 0 , 8 . 7, and 10.0. The solid lines are isochrones with the solar-calibrated value α MLT = 1 . 82. The dashed and dotted lines correspond to isochrones with α MLT = 1 . 6 and α MLT = 2 . 0, respectively. For display purposes, we only show the evolution up through the RGBTip. The physical interpretation of a small value of α MLT is a fluid parcel traveling a short radial distance (in units of H 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, α MLT is used to mimic the effects of physical ingredients, e.g., the inhibition of convection by a magnetic field, that are  \nFigure 15. Acomparison between MIST (black; this work), PARSEC v1.2S (red; Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014), Y 2 (orange; Demarque et al. 2004), DSEP (green; Dotter et al. 2008), BaSTI 'noncanonical' (turquoise; Pietrinferni et al. 2004), and Lyon (navy; Baraffe et al. 1998, 2003, 2015) isochrones at Z = Z glyph[circledot] as defined by each model.  \n<!-- image -->", '7.1. Evolutionary Tracks and Isochrones': "We first compare isochrones at Z glyph[circledot] 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 2 (orange; Demarque et al. 2004), DSEP (green; Dotter et al. 2008), BaSTI 'non-canonical' with η = 0 . 4 (turquoise; Pietrinferni et al. 2004), and Lyon (navy; Baraffe et al. 1998, 2003, 2015) for log(Age) [yr] = 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. 19 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 glyph[circledot] , 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 -τ boundary conditions that have been empirically calibrated to match the observed mass-radius relations for cool dwarfs (Chen et al. 2014).  \nIn Figure 16, we now compare isochrones at fixed Z . Note  \n19 There is a version of Y 2 models with rotation for M < 1 . 25 M glyph[circledot] . See Spada et al. (2013).  \nFigure 16. The same as Fig 15, except now at Z = 0 . 0001 and Z = 0 . 03 for logAge = 10 . 0.  \n<!-- image -->  \nFigure 17. The evolutionary tracks for a 0 . 3 M glyph[circledot] star from MIST, PARSEC, and Lyon models at Z = Z glyph[circledot] .  \n<!-- image -->  \nthat 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 log(Age) [yr] = 10 . 0 isochrones at Z = 0 . 0001 and Z = 0 . 03 for PARSEC, Y 2 , 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 η R = 0 . 4, 0.2, and 0.1, respectively, while Y 2 does not include mass loss).  \nFigure 17 presents a comparison between the MIST, PARSEC, and Lyon models for a 0 . 3 M glyph[circledot] evolutionary track. As noted above, the PARSEC models now adopt a modified T -τ 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.  \nFigure 18. The evolution of integrated colors of a simple stellar population for MIST (black), PARSEC/COLIBRI (red), and BaSTI (blue) isochrones at solar metallicity. The colors were computed with the Flexible Stellar Population Synthesis code (FSPS, v2.6; Conroy et al. 2009; Conroy & Gunn 2010) assuming a Kroupa (2001) IMF and AGB circumstellar dust turned off. 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.  \n<!-- image -->", '7.2. Simple Stellar Population Colors': '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 20 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 log(Age) [yr] < 7 . 5 because they only go up to 10 M glyph[circledot] 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.  \nOverall, 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 log(Age) [yr] glyph[greaterorsimilar] 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 (log(Age) [yr] ∼ 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  \nFigure 19. A comparison of MIST and Geneva (Ekström et al. 2012) evolutionary tracks with rotation in black and pink, respectively.  \n<!-- image -->  \nto 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.', '7.3. The Effects of Rotation': '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 crit = 0 . 4 and adopt a similarly low-metallicity solar abundance scaleZ glyph[circledot] = 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.  \nDifferences in the efficiency of rotational mixing between the MIST and Geneva models is 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 crit = 0 . 4 and f µ = 0 . 05, so-  \nSome MIST models, such as the 120 M glyph[circledot] star in Figure 19, lose their hydrogen-rich envelope very promptly as they reach the so-called ΩΓ -limit (Maeder & Meynet 2000). As evident from Equations 21 and 22, massive stars with Γ = L / L Edd → 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 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.  \nFigure 20. The ratio of MS lifetimes for rotating and non-rotating models as a function of initial mass. Left: comparison at Z glyph[circledot] . The solid black, green, blue, and dashed pink lines correspond to the default model with v / v crit = 0 . 4 and f µ = 0 . 05, model with v / v crit = 0 . 6, model with f µ = 0 . 01, and the Geneva model, respectively. The default MIST model shows a modest ∼ 5% enhancement in the MS lifetime due to rotational mixing, whereas the Geneva model experiences a ∼ 20% increase due to more efficient rotational mixing. Right: comparison among MIST models at [Z / H] = -1 . 0, 0.0, and + 0 . 5 with default parameters v / v crit = 0 . 4 and f µ = 0 . 05. The efficiency of rotational mixing is larger in more metal-poor stars because line-driven mass loss-thus angular momentum loss efficiency-is lower.  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \nlar metallicity model with v / v crit = 0 . 6 and f µ = 0 . 05, solar metallicity model with v / v crit = 0 . 4 and f µ = 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. 21 The default MIST model experiences only a modest enhancement in the MS lifetime. In contrast, the Geneva model experiences an overall ∼ 25% increase in the MSlifetime for stars more massive than 2 M glyph[circledot] (Ekström et al. 2012; Georgy et al. 2013).  \nIn 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  \nAlthough 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 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 fc = 1 / 30 and f µ = 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 ( glyph[lessorsimilar] 1 . 5 M glyph[circledot] ), 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 nonrotating and rotating Geneva models.  \nmetallicity-dependent line-driven winds, rotational mixing becomes more important at low metallicities due to the lowered efficiency of angular momentum loss, as expected.', '8.1. Luminosity-Mass-Radius-Temperature Relations': '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).  \nIn Figure 21, we plot log( R ), log( L ), and log( T eff) as a function of stellar mass for the DEBCat 22 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 ∼ 1 M glyph[circledot] . We applied a log( g ) cut-log( g ) > 4 . 1 cm s -2 and 3 . 4 cm s -2 for M i > 1 . 2 M glyph[circledot] and < 1 . 2 M glyph[circledot] , 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 dotted-dashed lines demarcate the initial masses for which t MS ∼ t Hubble for [Z / 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 i glyph[lessorsimilar] 0 . 7 M glyph[circledot] . 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  \n1  \n0  \nFigure 21. log( R ), log( L ), and log( T eff) as a function of stellar mass measured for MS stars in detached eclipsing binaries (DEB). The Southworth (2015) sample in blue comes from DEBCat, an online catalog of DEBs with wellmeasured parameters gathered from the literature. The red points correspond to 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 ∼ 1 M glyph[circledot] . 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 dotted-dashed lines demarcate the initial masses for which t MS ∼ t Hubble for [Z / H] = -1 . 0, 0.0, and + 0 . 5, respectively. The insets highlight the well-known discrepancy for the low-mass stars ( < 0 . 7 M glyph[circledot] ).  \n<!-- image -->  \n/circledot  \ntemperatures by ∼ 0 . 05 dex, for a total deficit of ∼ 0 . 2 dex for the predicted luminosity.  \nThere 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). Largescale 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).', '8.2. Initial-Final Mass Relation': "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 IFMS (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. 2008b; Kalirai et al. 2008; Williams et al. 2009; Zhao et al. 2012).  \nWecompare the predicted IFMR to a sample of eight young open clusters, three older open clusters with ages > 1 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 turn-off and the  \nFigure 24 shows optical and near-infrared (NIR) CMDs and log(Age) [yr] = 9 . 58 (3.80 Gyr) isochrones with [Z / H] = 0 . 0, AV = 0 . 2, and µ = 9 . 7, where µ is the distance modulus. The Sandquist (2004) BV sample (blue points) was carefully selected using proper motion, radial velocity, variability, and  \n<!-- image -->  \n/circledot  \nFigure 22. The initial-final mass relation constructed using binned cluster data and Sirius B from Ferrario et al. (2005) (see references therein) and Kalirai et al. (2008). The predicted relations for [Z / H] = -0 . 5, 0.0, and + 0 . 5 are shown in blue, gray, and red. These metallicities roughly bracket the metallicity range of the systems in the sample.  \nWD 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.  \nFurthermore, the steep slope predicted around M i ∼ 3 M glyph[circledot] to 4 M glyph[circledot] (slope ∼ 0 . 16) is consistent with the empirical estimate from Cummings et al. (2015), who found a slope of M f = (0 . 163 ± 0 . 022) M i + (0 . 238 ± 0 . 071) M glyph[circledot] 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 glyph[circledot] . 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.  \nIn 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 i glyph[lessorsimilar] 2 M glyph[circledot] ( M f ∼ 0 . 6 M glyph[circledot] ) is marginally consistent with the peak of the galactic disk WD mass function near ∼ 0 . 6 M glyph[circledot] (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. 2008a).  \nIt is worth noting that calculations of M i and M f which rely on, e.g., WD spectral analysis and isochrone fitting, and the quality of the data vary between cluster to cluster. As noted in  \nKalirai 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 is 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 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. 2008a). Catalán et al. (2008a) found that reasonable variations in the envelope thickness can lead to differences as large as 1 M glyph[circledot] for progenitors with M i glyph[greaterorsimilar] 5 M glyph[circledot] but a more modest ∼ 0 . 1 M glyph[circledot] 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 i and M f.  \nWe 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 f -M 1stTP / M 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-3 M glyph[circledot] . 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.", '8.3. Optical and NIR Color Magnitude Diagrams of Clusters 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 V ≡ A 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.  \nM67 (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).  \nFigure 24. CMDs for M67 from 2MASS and BVI photometry (Sandquist 2004; Sarajedini et al. 2009). The Sandquist (2004) BVI sample (blue points) was carefully selected to exclude likely binaries, and the Sarajedini et al. (2009) sample (mauve points) reflects a membership probability cut of > 20%. MIST isochrones with [Z / H] = 0 . 0, log(Age) [yr] = 9 . 58 (3.80 Gyr), AV = 0 . 2, and µ = 9 . 7 are shown in black.  \n<!-- image -->  \n/circledot  \n/circledot  \nFigure 23. Left: a comparison highlighting the difference between the final mass (solid) from Figure 22 and the core mass at first thermal pulse (dashed) as a function of initial mass for three metallicities. The latter is devoid of large uncertainties in mass loss, third dredge-up, and hot bottom burning that strongly influence the TPAGB evolution, and consequently, the final remnant mass. Right: fractional growth in core mass during the TPAGB phase for the same three metallicities.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nCMD-location information to yield cluster members that are most likely to be single stars. The Sarajedini et al. (2009) 2MASS sample (mauve points) reflects a membership probability cut of > 20%.  \nPraesepe (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 glyph[greaterorsimilar] 1000; Kraus & Hillenbrand 2007), large proper motion, and proximity makes Prae-  \na favorable target for stellar population studies.  \nPleiades (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 ∼ 1400) and proximity make it a popular choice for test-  \nIn 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-motionselected cluster members that are likely to be single stars and includes stars with masses as low as ∼ 0 . 15 M glyph[circledot] . 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 i glyph[lessorsimilar] 0 . 8 M glyph[circledot] . We applied an additional cut to remove objects that were flagged as variable stars and/or had < 50% membership probability. We overplot log(Age) [yr] = 8 . 8 (630 Myr) isochrones with [Z / H] = + 0 . 15, AV = 0 . 08, and µ = 6 . 26. Overall, the MIST isochrones provide good fits in all three colors.  \nThe MS, MSTO morphology, as well as 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 ∼ 17 which corresponds to M i ∼ 0 . 7 M glyph[circledot] . This is a wellknown 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.  \nFigure 27. CMDs for NGC 6791 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. MIST isochrones with [Z / H] = + 0 . 47, log(Age) [yr] = 9 . 93 (8.5 Gyr), AV = 0 . 32 and µ = 13 . 1 are shown in black.  \n<!-- image -->  \n-  \n-  \n-  \nFigure 25. CMDs for Praesepe from 2MASS (Wang et al. 2014) and UKIDSS Galactic Clusters Survey photometry (Boudreault et al. 2012). The Wang et al. (2014) sample (mauve points) consists of proper-motion-selected cluster members that are likely to be single stars and the Boudreault et al. (2012) sample (blue points) was identified with an astrometric and five-band photometric cut. MIST isochrones with [Z / H] = + 0 . 15, log(Age) [yr] = 8 . 80 (630 Myr), AV = 0 . 08 and µ = 6 . 26 are shown in black.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \nFigure 26. CMDs for Pleiades from 2MASS and BVIC photometry Stauffer et al. (2007); Kamai et al. (2014). The mauve points represent the Stauffer et al. (2007) compilation of 2MASS and BVIC photometry from the literature, the green points correspond to a sample of newly identified 2MASS candidates by Stauffer et al. (2007), and the blue points are the Kamai et al. (2014) sample of proper motion members from the Stauffer catalog with updated BVIC photometry. MIST isochrones with [Z / H] = 0 . 0, log(Age) [yr] = 8 . 0 (100 Myr), AV = 0 . 1, and µ = 5 . 62 are shown in black. We omit the TPAGB and post-AGB phases for display purposes.  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \n<!-- image -->  \n-  \ning 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).  \nNGC 6791 is one of the most well-known and well-studied open clusters in the Milky Way. Its unusually old age ( ∼ 8 Gyr) and high metallicity ([Fe / H] ∼ 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.  \nIn Figure 26, we show optical and NIR CMDs constructed from the Stauffer et al. (2007) compilation of 2MASS and BVIC 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 BVIC photometry (blue points). We overplot log(Age) [yr] = 8 . 0 (100 Myr) isochrones with [Z / H] = 0 . 0, AV = 0 . 1, and µ = 5 . 62. Overall, the isochrones are wellmatched to the observed CMDs in all filters. However, as seen in M67, the isochrones depart blueward from the MS ridge line in V -IC and B -V at fainter than V ∼ 14, corresponding to M i ∼ 0 . 6 M glyph[circledot] .  \nIn 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 log(Age) [yr] = 9 . 93 (8.5 Gyr) isochrones with [Z / H] = + 0 . 47, AV = 0 . 32 and µ = 13 . 1.  \nRuprecht 106 is a relatively low-mass ( M ∼ 10 4 . 8 M glyph[circledot] ) globular cluster with a metallicity of [Fe / H] ∼ -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 anti-correlation, 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  \nThe 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 ∼ 21. We return to this point at the end of this section.  \nFigure 28 shows the optical CMD from the HST ACS observations in the F 606 W and F 814 broadband filters (Dotter et al. 2011). We plot a log(Age) [yr] = 10 . 08 (12.0 Gyr) isochrone with µ = 16 . 7, AV = 0 . 55, and [Z / 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.  \n<!-- image -->  \n-  \nFigure 28. CMD from the HST ACS observations in the F 606 W and F 814 broadband filters Dotter et al. (2011). We overplot a [Z / H] = -1 . 50, log(Age) [yr] = 10 . 08 (12.0 Gyr) isochrone with µ = 16 . 7 and AV = 0 . 55.  \n<!-- image -->  \n-  \nFigure 29. LkCa 3, a quadruple system of PMS stars in the Taurus-Auriga star-forming region (Torres et al. 2013) is compared to 1, 1.4, and 3 Myr isochrones at solar metallicity.  \nmore tests against a large sample of globular clusters with our α -enhanced models in Paper II.  \nWe 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 glyph[lessorsimilar] 0 . 6-0 . 7 M glyph[circledot] . These discrepancies are due to missing and/or  \ninaccurate 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 TaurusAuriga 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 HR 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.  \nIn 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 vs. MV plane. The observations as reported by Torres et al. (2013) were already corrected for interstellar extinction assuming A 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.', '8.4. The Age Discrepancy in Upper Scorpius': "Upper Scorpius (Upper Sco) is one of three subgroups (Upper Centaurus Lupus and Lower Centaurus Crux) in ScorpiusCentaurus (Sco-Cen), the nearest OB association from the Sun with d ∼ 145 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).  \nPecaut 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 datasets 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  \nFigure 30. HR diagram of Upper Sco stars from Preibisch & Zinnecker (1999), Preibisch et al. (2002), and Pecaut et al. (2012), in blue, green, and red symbols, respectively. Solar metallicity isochrones with ages of 2, 5, 10, and 20 Myr are shown in black. There is a well-known discrepancy in the ages inferred from stars with spectral types earlier than G-type versus those later than G-type (with the exception of Antares which has evolved off the MS).  \n<!-- image -->  \nall consistently older than 5 Myr albeit with a large scatter, and the authors concluded that Upper Sco has a median age of 11 ± 1 ± 2 (statistical, systematic) Myr.  \nRecent work by Kraus et al. (2015) on UScoCTIO 5, a lowmass 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 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. 23 A simple exercise where they horizontally shifted the evolutionary tracks in the log L -log T eff plane to bring the predicted T 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  \nFigure 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 glyph[lessorsimilar] 5 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).  \n23 The Lyon, Dartmouth, Pisa, and Siess models underpredict the radius, which is a well-known problem (see Section 8.1). On the other hand, the PARSEC model overpredicts the radius, which Kraus et al. (2015) attributes to a likely overcorrection in their new boundary conditions.  \n∼ 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.", '8.5. Asymptotic Giant Branch Stars': "Low- and intermediate-mass stars (1 M glyph[circledot] glyph[lessorsimilar] M i glyph[lessorsimilar] 8 M glyph[circledot] , 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 heliumrich 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 5 -10 6 years depending on mass and metallicity (see Figure 12), and the majority of that time is spent in the quiescent 'interpulse' state.  \nGiven 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  \nThe 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 ∼ 10 -3 M glyph[circledot] year -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 i glyph[greaterorsimilar] 4 -8 M glyph[circledot] ) superAGB 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).  \nFigure 31. Cumulative luminosity functions of AGB and C stars in the SMC. The observed LFs in black lines are constructed from the SAGE-SMC survey photometric catalogs (private communication, M. Boyer; see also Gordon et al. 2011; Boyer et al. 2011 for more details) and the gray bands reflect Poisson uncertainties. The predicted LFs are computed by convolving the isochrones with star formation and metallicity histories from Harris & Zaritsky (2004) to create composite CMDs in different filters. The stars are then selected using two methods: 'CMD cut' uses the same CMD-based criteria that were applied to the SAGE observed sample and 'phase cut' selects all stars that are phase-tagged as 'TPAGB' stars in the isochrones. The differences between the two methods emphasize the need for a careful analysis when comparing populations in the CMD.  \n<!-- image -->  \ncomputing 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.  \n5  \nFigure 32. The same as Figure 31, now for the LMC. The data come from the SAGE-LMC survey (private communication, M. Boyer; see also Meixner et al. 2006 for more details) and the star formation and metallicity histories are from Harris & Zaritsky (2009).  \n<!-- image -->  \n10", '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 (C / O > 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 bol ∼ -5 . 0 is sensitive to the efficiency of mass loss as well as TDU, and the bright end reflects the decreasing numbers both due to 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).  \nThe 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  \nFigure 33. The same as Figure 31, but now comparing MIST predictions to the PARSEC/COLIBRI (Bressan et al. 2012; Marigo et al. 2013; Rosenfield et al. 2014) and BaSTI (Pietrinferni et al. 2004) predictions in J and [3.6]. The predicted composite CMDs used to generate the LFs displayed here were selected with a CMD cut.  \n<!-- image -->  \n5  \nFigure 34. The same as Figure 33, now for the LMC.  \n<!-- image -->  \nMeixner et al. 2006; Gordon et al. 2011; Boyer et al. 2011 for more details). The very wide baseline in wavelengthoptical UBVI from the Magellanic Clouds Photometric Survey (MCPS; Zaritsky et al. 2002) all the way out to 160 µ 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  \nconstruct cumulative LFs.  \nIn 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.  \nTo create the predicted LFs, we utilize the Flexible Stellar Population Synthesis code (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.  \nFirst, 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.  \nFinally, 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  \nSecond, 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 (C / O < 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.  \nstars in the C star sample.  \nWe 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 ˙ M required for the computation of AGB circumstellar dust effects came directly from the isochrone files for MIST and PARSEC/COLIBRI, while for BaSTI ˙ 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.  \nTo 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 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, as well as 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.', '9.1. Width of the MS': "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 MSTOcalibrated 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 HR diagram.' It differs from an ordinary HR diagram in that it still has log( T eff) on the x -axis but a new quantity L ≡ T 4 eff / g on the y -axis. The main advantage of a spectroscopic HR diagram is that all of the relevant quantities can be obtained from spectroscopic analysis without having to worry about ambiguities in distance or extinction.  \nIn 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 log L / L glyph[circledot] glyph[greaterorsimilar] 4. This continuous distribu-  \nFigure 35. Spectroscopic HR diagram, with log T eff on the x -axis and the quantity L ≡ T 4 eff / g on the y -axis. Black dashed lines correspond to empirical fits to the probability density distribution constructed from a sample of more than 600 stars (Castro et al. 2014). The red and blue tracks show the MS in the solar metallicity MIST evolutionary tracks with and without rotation, respectively.  \n<!-- image -->  \ntion 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 L and suggested that core overshoot efficiency may be mass-dependent.  \nWe overplot a series of solar-metallicity MIST evolutionary tracks with masses ranging from 10 to 80 M glyph[circledot] 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.", '9.2. Locations of Red Supergiants on the HR Diagram': "When a high-mass star runs out of hydrogen in its core, it may migrate toward the Hayashi track and become a 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.  \nLevesque 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 HR 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.  \nFigure 36. A comparison between the MIST evolutionary tracks with rotation and observed RSGs (Levesque et al. 2005, 2006; Massey et al. 2009). The masses of the tracks displayed are 10, 16, 20, 26, and 30 M glyph[circledot] . The observed L bol is calculated from K -band photometry. Top left: M31 ([Z / H] = + 0 . 3 in black; [Z / H] = 0 . 0 in gray). Top right : Milky Way ([Z / H] = 0 . 0). Bottom left: LMC ([Z / H] = -0 . 5). Bottom right: SMC ([Z / H] = -0 . 75).  \n<!-- image -->  \n(2006) also found that RSGs in the SMC span a wide range of T eff for a given M bol, possibly due to the increased importance of rotational mixing at lower metallicities (see Section 7.3 for more details).  \nIn 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 [Z / 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 bol shown here is derived from MK rather than MV since the former is less sensitive to extinction. The typical measurement uncertainty in T eff ranges from ∼ 100 K for the warmest stars to ∼ 20 K for the coolest stars, and the uncertainty in 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 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.", '9.3. Relative Lifetimes': '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 obs , A / N obs , B:  \n∫ t A( M ) φ ( M ) dM ∫ t B( M ) φ ( M ) dM = N obs , A N obs , B , (24)  \n<!-- image -->  \n<!-- image -->  \nFigure 37. Left: the predicted WC to WN ratio as a function of metallicity. The red circles are observations from Neugent & Massey (2011) and Neugent et al. (2012), and the blue diamond point is the observed ratio from Georgy et al. (2012), estimated in the 3 kpc radius volume around the solar neighborhood. Middle: the same as the left panel except now showing the WR to O ratio. The red triangles are observed number ratios from Table 6 in Maeder & Meynet (1994) and the blue diamond point is the observed ratio from Georgy et al. (2012), estimated in the 2.5 kpc radius volume around the solar neighborhood. Right: the same as the left panel except now showing the BSG to RSG ratio. The gray points are observed ratios from young star clusters in the Milky Way and the Magellanic Clouds (Eggenberger et al. 2002). The blue hexagons and triangles are observed ratios computed with and without K-type stars in the RSG category (Massey 2002). The red squares are observed ratios computed from spectroscopically confirmed RSG stars (Massey & Olsen 2003).  \n<!-- image -->  \nwhere 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 i glyph[greaterorsimilar] 10 M glyph[circledot] and MS lifetimes glyph[lessorsimilar] 20 Myr (but see also e.g., Dohm-Palmer & Skillman 2002 where the authors examined the ratio of blue to red supergiants as a function of age in Sextans A).  \nHere we present three such tests: the ratio of WR subtypes WC to WN, the ratio of WR to O-type stars, and the ratio of blue supergiants (BSGs) to RSGs. We convert metallicities reported in the literature-log(O / 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 WR 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.  \nThere 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 WR 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  \nratios 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.  \nFirst 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.  \nWe identify WR 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).  \nIn the left panel of Figure 37, we plot the predicted WC to WNratio 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 to 3. We plan to improve this further in future work.', '9.3.2. WR/O': 'The predicted ratio of WR to O-type stars is computed from the models using the method outlined in the previous section and compared to observations. Here, the WR population is the sum of all WR subtypes. O-type MS stars are identified according to the Georgy et al. (2012) classification scheme.  \nIn the middle panel of Figure 37, we show the predicted WR to O ratio as a function of metallicity in a solid black  \nline. The red triangles are observed number ratios from Table 6 in Maeder & Meynet (1994) 24 and the blue diamond point is the observed ratio computed by Georgy et al. (2012), estimated in the 2.5 kpc radius volume around the solar neighborhood. Again, the model prediction is in qualitative agreement with the observations; WR stars become more abundant in higher metallicity environments. This is to be expected because efficient mass loss at high metallicities readily removes the hydrogen-rich outer layers and promotes the formation of WRstars.', '9.3.3. BSG/RSG': 'The number ratio of blue to red supergiants (BSG/RSG) has long been known to decrease with increasing galactocentric radius in the Milky Way, the Magellanic Clouds, and M33 (e.g., Walker 1964; Hartwick 1970; Cowley et al. 1979; Humphreys 1979; Meylan & Maeder 1982). This was explained by invoking radial metallicity gradients in the disks (van den Bergh 1968). The BSG/RSG ratio is an excellent diagnostic tool because whether a star becomes a RSG or a BSG hinges very sensitively on, for example, the details of semiconvection and convective overshoot (e.g, Langer 1991).  \nIn the right panel of Figure 37, we plot the observed BSG/RSG ratios in the Magellanic Clouds from Massey (2002) and Massey & Olsen (2003). The two Massey (2002) symbols correspond to ratios computed including (blue hexagon) and excluding (blue triangle) K-type stars in the RSG category from their UBVR -photometry-selected sample. Their sample was limited to M bol < -7 . 5 (roughly log( L / L glyph[circledot] ) > 4 . 9) in order to minimize contamination by the AGB stars and improve completeness. There was an estimated ∼ 10% contamination by foreground red dwarfs in their photometric sample, but this was dramatically improved in Massey & Olsen (2003) with spectroscopic radial velocity measurements. The updated ratios from spectroscopically confirmed RSG stars (Massey & Olsen 2003) are shown as red squares.  \nMost stellar evolution models have struggled to reproduce this radial/metallicity trend (see Langer & Maeder 1995, for an in-depth discussion of this topic). However, Maeder & Meynet (2001) found that the inclusion of rotation in their models produced more RSGs at low metallicities, which dramatically lowered the predicted BSG/RSG ratio in the SMC and brought the model prediction into better agreement with the observations. Eldridge et al. (2008) computed their model predictions with and without the effects of binarity and concluded that their single-star model underpredicted the BSG/RSG ratio. A mixed population model with a two-thirds binary fraction was required to reproduce the observations.  \nWe also show the observed BSG/RSG ratios computed from spectroscopically identified candidates in a sample of 45 young (6 . 8 < log(Age) [yr] < 7 . 5) clusters in the Milky Way and the Magellanic Clouds (gray; Eggenberger et al. 2002). The authors computed the observed BSG/RSG ratios at different metallicities by binning the stars according to their galactocentric distances and assuming a radial metallicity gradient to assign metallicities to each radius bin. Their results showed increasing BSG/RSG ratio with increasing metallicity.  \nWe overplot our prediction as a solid black line. RSGs and BSGs were identified using the selection criteria from Eldridge et al. (2008), which are consistent with observational cuts made by Massey (2002) and Massey & Olsen (2003).  \nFigure 38. Surface nitrogen enrichment midway through the MS ( X c ∼ 0 . 5; red) and at TAMS (blue) as a function of initial mass. The red and blue shaded boxes correspond to the average nitrogen excess observed for galactic MS B-type stars with M i < 20 M glyph[circledot] and the maximum observed excess (Gies &Lambert 1992; Kilian 1992; Morel et al. 2008; Hunter et al. 2009). Without rotation, the predicted nitrogen enrichment during the MS at these masses is zero. This figure is adapted from Figure 11 of Ekström et al. (2012).  \n<!-- image -->  \n/circledot  \nThe model predictions are bracketed by the two Massey (2002) points and marginally consistent with Massey & Olsen (2003). We do not reproduce the positive trend reported by Eggenberger et al. (2002), but Eldridge et al. (2008) raised the concern that clusters in the Eggenberger et al. (2002) sample generally have an age spread similar to the age of the cluster itself. Additional data over a wider range of environments would be valuable in assessing the quality of the massive star models.', '9.4. The Effects of Rotation on Surface Abundances': '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. 25 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 fc ; Pinsonneault et al. 1989; Heger et al. 2000; Brott et al. 2011a).  \nFollowing 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 c ∼ 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 glyph[circledot] 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  \n<!-- image -->  \nFigure 39. Left: evolutionary tracks in the log( g ) -log( T eff) plane for 7, 16, and 40 M glyph[circledot] stars at solar metallicity. Solid and dotted lines are models with and without rotation and open and filled symbols correspond to O- and B-type galactic dwarfs and supergiants, respectively. This figure is adapted from Figure 12 of Ekström et al. (2012). Right: surface nitrogen abundance evolution for the same three models.  \n<!-- image -->  \na 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 ∼ 2 M glyph[circledot] . This is due to the inclusion of magnetic braking effects in the Geneva models; stars with masses below ∼ 1 . 7 M glyph[circledot] 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 i < 1 . 2 M glyph[circledot] and gradually increase the rotation rate from v ZAMS / v crit = 0 . 0 to 0 . 4 over the mass range M i = 1 . 2 to 1 . 8 M glyph[circledot] . 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 i < 10 M glyph[circledot] would provide a valuable constraint on the models.  \n(2012) to aid comparison with their figure. As in Ekström et al. (2012), we do not compare surface rotation velocities because only v 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 glyph[circledot] and that the observed 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.  \nWe also note that for some samples consisting exclusively of slow rotators (e.g., the Morel et al. (2008) sample of B stars whose v sin i values range from 10 to 60 km 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 eff, 0 . 1-0 . 2 dex for log( g ), and 0 . 2 dex for the N/H ratio. This is meant to be a coarsegrained 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.', '10. CAVEATS AND FUTURE WORK': "In this paper we presented an overview of the MIST models, including a comprehensive discussion of the input physics  \nIn the left panel of Figure 39, we show 7, 16, and 40 M glyph[circledot] solar metallicity evolutionary tracks with (solid) and without (dashed) rotation in the log( g )-log( T 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.  \nand comparisons with existing databases and observational constraints. We conclude with a discussion of some of the caveats and plans for future work.  \nAnother major caveat is that binary interaction is not taken into account in these models. Multiplicity is extremely common among O- and B-type stars; binary mass exchange is believed to occur for glyph[greaterorsimilar] 70% of O-type stars and about a third of those stars will ultimately form a binary merger product (e.g., Chini et al. 2012; Sana et al. 2012, 2013; de Mink et al. 2014). These numbers have serious implications for the evolution of massive stars and their explosive final fates, including the expected frequency of different types of core-collapse supernovae (see the review by Smartt 2009). Since binarity dramatically expands the size of the parameter space to be explored (e.g., mass ratio, eccentricity, and separation), it is currently computationally infeasible to construct a multidimensional grid of binary models from full stellar evolution calculations. In order to make this a tractable problem, the standard approach has been to either couple detailed stellar evolution codes such as MESA to binary population synthesis codes (e.g., Eldridge et al. 2008) or to make use of fitting formulae to approximate models of single stars (e.g., Hurley et al. 2002; Izzard et al. 2006). Although binarity is not included in the current models, some of its evolutionary consequences must be at least partially mimicked through the effects of rotation and ordinary single-star mass loss, given that the models are in broad agreement with a number of observations. Moreover, it is possible to couple MIST to available  \nPerhaps the most significant shortcoming common to all stellar evolutionary tracks and isochrones of this generation is that they are computed within a 1D framework despite the inherently 3D nature of stellar astrophysical phenomena, e.g., mass-loss, convective mixing, rotation, and magnetic fields. There has been recent progress in 2D and 3D simulations of stellar interiors and atmospheres, but they are limited to small spatial scales and physical time durations of order only ∼ 10 hours (e.g., Woodward et al. 2015). Although full 3D simulations of stellar evolution are much too far beyond grasp today, we are nevertheless taking small steps forward. For example, it is becoming increasingly common to map 3D simulation results to a 1D formulation in order to incorporate the valuable insight we are gaining from these sophisticated simulations into standard 1D stellar evolution models (see e.g., Brown et al. 2013; Trampedach et al. 2014; Arnett et al. 2015; Magic et al. 2015). This is an active area of development in MESA.  \nbinary synthesis codes to investigate and model binary effects in detail in the future.  \nWethank the referee for his/her constructive comments. We thank Martha Boyer, Jason Kalirai, and Guillermo Torres for generously sharing their data with us and Francis Timmes for his assistance with nuclear reaction networks in MESA. We thank Leo Girardi and Dan Weisz for fruitful conversations and comments on the manuscript. We thank Conny Aerts, Lars Bildsten, Jonathan Bird, Phillip Cargile, Selma de Mink, Falk Herwig, Jessica Lu, Marc Pinsonneault, Philip Rosenfield, and Victor Silva-Aguirre for sharing their expertise. We also thank the entire MESA community for making this work possible, with special thanks to the MESA Council. Finally, we thank the participants of the KITP workshop 'Galactic Archaeology and Precision Stellar Astrophysics' for stimulating discussions. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915. JC acknowledges support from the National Science Foundation Graduate Research Fellowship Program. AD received support from the Australian Research Council under grant FL110100012. CC acknowledges support from NASA grant NNX13AI46G, NSF grant AST-1313280, and the Packard Foundation. This work was supported in part by NASA grant NNX15AK14G. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.", 'FEATURES IN THE EVOLUTIONARY TRACKS AND ISOCHRONES': "Here we identify and discuss some of the interesting and unusual features in the evolutionary tracks and isochrones. Figure A1 shows five such examples. For each row, the left panel shows the evolutionary track in its entirety, and the gray box marks the zoomed-in region shown in the middle panel. The right panel shows a series of relevant physical quantities plotted as a function of time.  \nThe first row displays an example of 'born-again' evolution during the post-AGB phase (e.g., Schoenberner 1979; Iben 1982; Iben et al. 1983). Following the end of the TPAGB phase, the star enters the post-AGB phase and rapidly evolves toward hotter temperatures. During this short-lived phase ( ∼ 10 3 -10 4 years), the luminosity remains roughly constant, powered primarily by hydrogen shell-burning with some contribution from gravitational contraction. In the most uneventful scenario (e.g., the 1 M glyph[circledot] model in the left panel, second row), the surface hydrogen envelope mass eventually falls below the critical value and hydrogen burning shuts off as a result. Having lost its nuclear energy source, the star then begins to cool and fade away as a WD. In contrast, more interesting things can happen if the star leaves the TPAGB in the middle of its last TP cycle, in which case it may subsequently undergo a late TP either as a post-AGB star (Late Thermal Pulse; LTP) or as a WD (Very Late Thermal Pulse; VLTP) (Kawaler et al. 1996). In either case, the star loops back around toward cooler temperatures as it puffs up from the sudden injection of energy into the envelope. The middle panel zooms in on the VLPT, where the beginning and the end of this episode are marked by the red circle and diamond, respectively. As the right panel shows, the sudden ignition of hydrogen is what triggers  \nIn Paper II, we will present models with non-solar-scaled abundance patterns for the same large range of masses, ages, metallicities, and evolutionary phases, and make extensive comparisons with the observed properties of globular clusters. Most of our intended future directions are also areas of known limitations. Nevertheless, other future projects we hope to pursue include: investigating the effects of varying α MLT across the HR diagram and as a function of metallicity (e.g., Trampedach et al. 2014; Magic et al. 2015), exploring the effects of magnetic fields in low-mass and PMS stars, improving the WD cooling models for joint fitting of entire CMDs of star clusters, implementing low temperature C/Ovariable molecular opacities for modeling the envelopes in AGB stars, performing more precise calibration of the AGB phase by jointly modeling the SFHs of the Magellanic Clouds and their AGB LFs, and developing non-standard and innovative isochrone construction methods.  \n)  \nL  \n/  \nL  \n(  \nlog  \nFigure A1. An illustration of interesting and unusual features in the evolutionary tracks. For each row, the left panel shows the evolutionary track in its entirety, and the gray box marks the zoomed-in region shown in the middle panel. The right panel shows the relevant physical quantities as a function of time. First row: 'born-again' evolution during the post-AGB phase. Second row: helium flashes following helium ignition at the tip of the RGB. Third row: TPAGB phase and post-AGB bump. Fourth row: a shift in T eff due to the initialization of rotation near ZAMS. Fifth row: 3 He-driven instability near the transition from fully convective to radiative core during the MS.  \n<!-- image -->  \n×  \nthe VLTP in the star just as it enters the WD cooling phase. The stellar bolometric luminosity and hydrogen burning luminosity are shown in blue and orange.  \nThe third row shows the TPAGB phase followed by a bump in the post-AGB track of a 2 . 6 M glyph[circledot] star. The TPAGB phase is famous both for its namesake thermal pulses triggered by alternating hydrogen and helium burning in shells (see Section 8.5 for a more detailed explanation) and high mass loss rates (see Herwig 2005 for an excellent review). The TPAGB phase is terminated once the star sheds almost all of its outer hydrogen-rich envelope to reveal the hot CO core, which launches the star leftward in the HR diagram. In the right panel, the alternating hydrogen and helium burning luminosities for the last ∼ two TPs are plotted in blue dashed and red dotted-dashed lines. The envelope mass shown in orange decreases rapidly until it reaches 0, marking the end of the TPAGB phase. Since the mass of the envelope is dramatically reduced but the stellar luminosity remains more or less the same, the Eddington ratio L / L Edd goes up and the stellar surface becomes unstable. The star glitches over a very short timescale ( < 1 yr) as it attempts to achieve hydrostatic equilibrium, which shows up as a sharp bump at around log( T eff) ∼ 4 . 2 K in the middle panel. We note that this feature only appears sometimes (generally for > 2 . 5 M glyph[circledot] at solar metallicity). Since it is unclear whether or not this behavior is real or a numerical artifact, and it has zero bearing on the evolution of the star due to its extremely short-lived nature, we post-process this feature out of the final evolutionary tracks in order to facilitate the construction of smooth isochrones. Specifically, we trim out any points with | d (log L ) / dt | < 0 . 1 during the post-AGB phase only. The resulting track is shown in black in the left and middle panels.  \nThe second row illustrates a series of helium flashes in the degenerate core of a low-mass star as it settles into quiescent helium burning. As the star climbs up the RGB, the inert helium core becomes denser and more degenerate as it undergoes gravitational contraction, continuously growing in mass as the hydrogen-burning shell above rains down freshly fused helium. For sufficiently large central densities, neutrino cooling becomes significant and the peak in temperature actually shifts away from the center. As a result, helium ignition occurs off-center at the RGBTip (Thomas 1967). 26 The initial helium flash is followed by a series of weaker helium flashes that move radially inward until helium burning reaches the center. This lifts the degeneracy of the core and the star commences quiescent helium burning. The oscillatory features in the middle panel correspond to these successive helium flashes. The right panel shows the helium burning luminosity, stellar bolometric luminosity, and temperature in orange, blue, and red, respectively. The sharp decline in stellar luminosity in concomitance with a sharp peak in helium burning luminosity marks the location of the RGBTip. The degeneracy parameter η = glyph[epsilon1] F / kBT in the center (shown in green, rescaled for display purposes) plummets to ∼ 0 once the helium flash reaches the innermost core.  \nThe fourth row illustrates the effect of turning on rotation at ZAMS. Once the star reaches ZAMS (defined to be L nuc / L ≥ 0 . 9 in MESA) solid-body rotation is established over 10 time steps. In the absence of efficient rotational mixing, rotation simply makes stars cooler (see Section 6.1). The middle panel shows that the rotating (black) and non-rotating (gray) tracks overlap up until ZAMS marked by an orange circle. Once rotation is established, the two curves diverge as the rotating model settles to a lower T eff. This jump in temperature is extremely short-lived and is purely a feature of our implementation of rotation; a real star is born from a birth cloud with nonzero initial angular momentum. A more realistic model of the PMS phase that includes the effects of rotation will be explored in the future. In the right panel, we plot v / v crit, central hydrogen abundance, and log L as a function of time. The small increase in velocity immediately following velocity initialization is due to the star experiencing additional contraction before it begins to steadily burn hydrogen. During the MS, v / v crit decreases as the stellar radius and luminosity gradually increase.  \nThe fifth and final row demonstrates the presence of a 3 He-driven instability in a 0 . 34 M glyph[circledot] star. This stellar mass represents the transition from stars that are fully convective (lower mass) on the MS to those with radiative cores and convective envelopes (higher mass). This instability was first reported by van Saders & Pinsonneault (2012) who found that these low-mass stars develop small convective cores-normally characteristic of CNO burning in more massive stars ( glyph[greaterorsimilar] 1 . 1 M glyph[circledot] )-due to nonequilibrium 3 He abundances at these low central temperatures. The net production of 3 He in the center leads to the development of a small convective core, which steadily grows in extent and eventually makes contact with the bottom of the deep convective envelope. As a result of vigorous convective mixing, the local 3 He enhancement in the center is erased as the central 3 He abundance dilutes back to the bulk abundance and the cycle resets. This event occurs as the star leaves the Hayashi track, as can be seen in the left and middle panels. As the right panel shows, a sudden drop in central 3 He abundance (orange) caused by the convective core coming into contact with the envelope coincides with a shrinking of the convective core (blue) and a drop in the stellar radius (red) and luminosity (not shown for display purposes).", 'THE EFFECTS OF VARYING TIME STEP AND MESH CONTROLS IN MESA EVOLUTION CALCULATIONS': 'Here we examine the temporal and spatial convergence in our MESA stellar evolution calculations following the methodology in Paxton et al. (2011). In MESA, the sizes of time steps and grid cells are controlled by varcontrol\\_target and mesh\\_delta\\_coeff , respectively, but there is additional flexibility to adjust the resolution specifically for certain zones or evolutionary phases. varcontrol\\_target is the target value for relative changes in the stellar structure variables, e.g., ρ ( m ), between two consecutive time steps and mesh\\_delta\\_coeff is the analogous parameter for differences between adjacent grid cells.  \nAdopting the notations introduced in Paxton et al. (2011), we vary temporal and spatial resolution in lockstep according to C , a numerical factor by which we multiply the default values for varcontrol\\_target and mesh\\_delta\\_coeff simultaneously. A small C demands small time steps and fine grid cells which increase the numerical accuracy but at the expense of longer computation times. It is thus advantageous to search for the largest C we can afford to use without sacrificing the quality of the models.  \nHere we consider five values of C : 1/4, 1/2, 1, 2, and 4. 27 Note that C = 1 corresponds to the models computed in this paper. For each value of C , we compute evolutionary tracks for four masses-1, 5, 20, and 100 M glyph[circledot] -that are representative of the range of stellar types in the MIST models. Each M i and C combination is computed until the end of hydrogen and helium burning in the core. Convergence of a given model is quantified using the ξ parameter, which is the difference between a variable at a given resolution and at the highest resolution in the study ( C = 1 / 4 in this case). For 100 M glyph[circledot] , however, the highest resolution is C = 1 / 2 because the C = 1 / 4 models did not finish running in a reasonable amount of time. Table B1 summarizes the results.  \nThe behavior is less clear at higher masses (20 and 100 M glyph[circledot] ), but this is not unexpected since massive star evolution is much more complex and sensitive to the choice of input physics. Although the answer can change as much as ∼ 12% between our default resolution and C = 1 / 2, this is not a huge effect given the uncertainties in the evolution calculations of massive stars. We regard agreement at the ∼ 10% level as satisfactory because at this level of precision nearly every other detail of the input physics matters, from MLT implementation to mass loss and rotation.  \nFor the low-mass stars, whose early evolution is relatively simple, the result at our default resolution C = 1 differs at most by a few percent compared to the model calculations at four times the resolution. We were also able to successfully run C = 1 / 8 and 1 / 16 models for 1 and 5 M glyph[circledot] models to TAMS. Even compared to model calculations at 16 times the resolution, our default C = 1 model is converged to within a few percent. As an additional test, we ran C = 2 and C = 4 models through to the WDCS to check if the final mass of the WD depends sensitively on the resolution; for example, improperly resolved convective boundaries may influence the growth of the core during the TPAGB phase. We find that decreasing the resolution changes the final mass of the WDby less than a percent, which is highly encouraging.  \n)  \n-0.057  \n-0.051  \n-0.02  \n-0.019  \n0.0  \nTable B1 Convergence test results.  \nC  \n: end of H burn  \n4  \n2  \n1  \n1  \n/  \n2  \n1  \n/  \n4  \nC  \n: end of He burn  \n4  \n2  \n1  \n1  \n/  \n2  \n1  \n/  \n4  \nN  \ncells  \ntimesteps  \nN  \nM  \n/  \nM  \nξ  \n(  \n×  \nlog(  \nL  \n/  \nL  \nξ  \n(  \n×  \nlog(  \nT  \neff) [K]  \n4.159  \n4.16  \n4.16  \n4.16  \n4.164  \nlog(  \nT  \neff [K])  \n3.614  \n3.624  \n3.63  \n3.63  \n3.631  \nξ  \n(  \n×  \nlog(  \nt  \n)  \n8.025  \n8.021  \n8.01  \n8.009  \n8.008  \nlog(  \nt  \n)  \n8.077  \n8.076  \n8.071  \n8.072  \n8.071  \nξ  \n(  \n×  \nlog(  \nT  \nc) [K]  \n7.663  \n7.658  \n7.657  \n7.657  \n7.657  \nlog(  \nT  \nc) [K]  \n8.379  \n8.371  \n8.368  \n8.371  \n8.369  \nξ  \n(  \n×  \n10  \nρ  \nlog(  \nξ  \nC  \n: end of H burn  \n4  \n2  \n1  \n1  \n/  \n2  \n1  \n/  \n4  \nC  \n: end of He burn  \n4  \n2  \n1  \n1  \n/  \n2  \n1  \n/  \n4  \nN  \ncells  \ntimesteps  \nN  \nM  \n/  \nM  \nξ  \n(  \n×  \nlog(  \nL  \n/  \nL  \nξ  \n(  \n×  \nlog(  \nT  \neff) [K]  \n4.468  \n4.442  \n4.437  \n4.434  \n4.435  \nlog(  \nT  \neff [K])  \n3.568  \n3.535  \n3.532  \n3.537  \n3.537  \nξ  \n(  \n×  \nlog(  \nt  \n)  \n6.998  \n6.961  \n6.953  \n6.949  \n6.947  \nlog(  \nt  \n)  \n7.034  \n7.002  \n6.995  \n6.993  \n6.991  \nξ  \n(  \n×  \nlog(  \nT  \nc) [K]  \n7.819  \n7.815  \n7.815  \n7.815  \n7.814  \nlog(  \nT  \nc) [K]  \n8.511  \n8.505  \n8.504  \n8.503  \n8.503  \nξ  \n(  \n×  \n10  \nρ  \nlog(  \nξ  \nC  \n: end of H burn  \n4  \n2  \n1  \n1  \n/  \n2  \n1  \n/  \n4  \nC  \n: end of He burn  \n4  \n2  \n1  \n1  \n/  \n2  \n1  \n/  \n4  \nN  \ncells  \ntimesteps  \nN  \nM  \n/  \nM  \nξ  \n(  \n×  \nlog(  \nL  \n/  \nL  \nξ  \n(  \n×  \nlog(  \nT  \neff) [K]  \n4.603  \n4.882  \n4.877  \n4.874  \n-  \nlog(  \nT  \neff [K])  \n5.228  \n5.2  \n5.199  \n5.198  \n-  \nξ  \n(  \n×  \nlog(  \nt  \n)  \n6.476  \n6.546  \n6.553  \n6.56  \n-  \nlog(  \nt  \n)  \n6.52  \n6.595  \n6.601  \n6.608  \n-  \nξ  \n(  \n×  \nlog(  \nT  \nc) [K]  \n7.894  \n7.87  \n7.868  \n7.867  \n-  \nlog(  \nT  \nc) [K]  \n8.56  \n8.529  \n8.53  \n8.528  \n-  \nξ  \n(  \n×  \n10  \nρ  \nlog(  \nξ  \n×  \n(  \n10  \n10  \n10  \nglyph[circledot]  \n10  \n-  \n3  \n)  \n1775.107  \n82.495  \n34.981  \n0.0  \n-  \nglyph[circledot]  \n)  \n6.274  \n5.805  \n5.774  \n5.749  \n-  \nlog(  \nL  \n/  \nL  \n)  \n2352.624  \n139.574  \n58.767  \n0.0  \n-  \nξ  \n×  \n(  \n10  \n1468  \n2727  \n5448  \n10818  \n-  \nN  \n275  \n323  \n323  \n328  \n-  \nN  \n65.916  \n25.712  \n24.583  \n23.752  \n-  \nM  \n/  \nM  \nξ  \n×  \n(  \ncells  \ntimesteps  \nglyph[circledot]  \n10  \n-  \n3  \n)  \n1952.055  \n40.126  \n18.353  \n0.0  \n-  \nglyph[circledot]  \n)  \n6.162  \n5.504  \n5.489  \n5.475  \n-  \n)  \n3867.197  \n67.71  \n31.926  \n0.0  \n-  \n-  \n3  \n-  \n3  \n10  \n-  \n3  \n10  \n-  \n3  \n10  \n×  \nρ  \n10  \n10  \n10  \nglyph[circledot]  \n10  \n-  \n3  \nglyph[circledot]  \n)  \n3.168  \n3.154  \n3.139  \n3.137  \n3.136  \nlog(  \nL  \n/  \nL  \n)  \n74.552  \n40.912  \n5.28  \n0.766  \n0.0  \nξ  \n×  \n(  \n10  \n345  \n547  \n1015  \n2003  \n3992  \nN  \n383  \n584  \n1019  \n1919  \n3699  \nN  \n4.997  \n4.997  \n4.997  \n4.997  \n4.997  \nM  \n/  \nM  \nξ  \n×  \n(  \ncells  \ntimesteps  \nglyph[circledot]  \n10  \n-  \n3  \n)  \n-0.355  \n-0.124  \n-0.033  \n-0.07  \n0.0  \nglyph[circledot]  \n)  \n3.294  \n3.23  \n3.173  \n3.175  \n3.17  \n)  \n330.852  \n147.317  \n7.776  \n12.006  \n0.0  \n-  \n3  \n-  \n3  \n10  \n-  \n3  \n10  \n-  \n3  \n10  \n×  \nρ  \nglyph[circledot]  \n]  \n2.001  \n2.001  \n2.006  \n2.007  \n2.007  \nlog(  \n5  \nM  \n-  \n3  \nglyph[circledot]  \n]  \n1.375  \n1.393  \n1.396  \n1.397  \n1.397  \nlog(  \n-  \n3  \n-  \n3  \n-  \n3  \n-  \n3  \n)  \n13.611  \n2.911  \n0.973  \n0.874  \n0.0  \nc) [g cm  \n10  \n)  \n-13.275  \n-14.014  \n-0.887  \n0.314  \n0.0  \n20  \nM  \n)  \n3.456  \n-0.532  \n-1.475  \n-1.287  \n0.0  \n-  \n3  \nglyph[circledot]  \n)  \n5.083  \n5.03  \n5.02  \n5.017  \n5.015  \nlog(  \nL  \n/  \nL  \n)  \n169.241  \n35.182  \n11.427  \n5.929  \n0.0  \nξ  \n×  \n(  \n10  \n1467  \n2683  \n5024  \n9912  \n17655  \nN  \n1803  \n2493  \n4460  \n8438  \n16409  \nN  \n19.279  \n19.202  \n19.184  \n19.188  \n19.212  \nM  \n/  \nM  \nξ  \n×  \n(  \ncells  \ntimesteps  \nglyph[circledot]  \n10  \n-  \n3  \n)  \n-181.746  \n-139.532  \n-117.17  \n-13.479  \n0.0  \nglyph[circledot]  \n)  \n5.174  \n5.131  \n5.121  \n5.101  \n5.099  \n)  \n190.493  \n77.659  \n53.169  \n4.741  \n0.0  \n-  \n3  \n-  \n3  \n10  \n-  \n3  \n10  \n-  \n3  \n10  \n×  \nρ  \n×  \n)  \n17.555  \n5.504  \n2.622  \n0.795  \n0.0  \nc) [g cm  \n10  \n)  \n-58.531  \n-13.623  \n-6.043  \n-1.536  \n0.0  \n-  \n3  \nglyph[circledot]  \n10  \n-  \n3  \n10  \n-  \n3  \n10  \n-  \n3  \n10  \n-  \n3  \n)  \n10.992  \n2.764  \n1.188  \n0.834  \n0.0  \nc) [g cm  \n10  \n)  \n-48.836  \n-10.201  \n-2.678  \n-0.329  \n0.0  \n-  \n3  \n100  \nM  \n-  \n3  \n-  \n3  \nglyph[circledot]  \n]  \n1.03  \n1.208  \n1.216  \n1.223  \n-  \nlog(  \nξ  \n×  \n(  \n-  \n3  \nξ  \n(  \nξ  \n×  \n(  \n)  \n21.884  \n3.362  \n-2.952  \n4.452  \n0.0  \nc) [g cm  \n10  \n)  \n2.19  \n6.022  \n7.495  \n9.271  \n0.0  \n-  \n3  \n)  \n78.495  \n3.91  \n4.801  \n0.0  \n-  \nc) [g cm  \n10  \n)  \n-426.172  \n-19.84  \n-0.399  \n0.0  \n-  \n-  \n3  \n-  \n3  \n-  \n3  \n-  \n3  \n-  \n3  \n)  \n65.67  \n8.314  \n3.609  \n0.0  \n-  \nc) [g cm  \n10  \n)  \n-359.052  \n-35.601  \n-15.776  \n0.0  \n-  \n-  \n3  \nξ  \n(  \nξ  \n×  \n(  \n)  \n-463.387  \n18.996  \n8.502  \n0.0  \n-  \n)  \n-174.392  \n-29.716  \n-15.581  \n0.0  \n-  \nξ  \n×  \n(  \n)  \n71.311  \n4.274  \n3.703  \n0.0  \n-  \n)  \n-183.629  \n-28.855  \n-14.912  \n0.0  \n-  \n1569  \n3152  \n6486  \n12465  \n-  \n696  \n572  \n559  \n584  \n-  \n33.813  \n11.914  \n11.664  \n11.454  \n-  \n×  \n(  \nξ  \n(  \nξ  \n×  \n(  \n)  \n78.496  \n17.202  \n4.979  \n-1.602  \n0.0  \n)  \n123.993  \n31.461  \n12.596  \n5.123  \n0.0  \nξ  \n×  \n(  \n)  \n74.031  \n-4.799  \n-10.998  \n-0.389  \n0.0  \n)  \n103.348  \n25.481  \n9.44  \n4.182  \n0.0  \n2325  \n4598  \n9133  \n18303  \n36666  \n2951  \n3635  \n6218  \n11591  \n22501  \n12.156  \n12.783  \n13.115  \n14.655  \n14.856  \n×  \n(  \nξ  \n(  \nξ  \n×  \n(  \n)  \n-12.433  \n-10.03  \n-7.975  \n-8.023  \n0.0  \n)  \n40.171  \n29.525  \n5.393  \n2.926  \n0.0  \nξ  \n×  \n(  \n)  \n-37.609  \n-15.963  \n-1.504  \n-1.559  \n0.0  \n)  \n13.82  \n12.067  \n-0.536  \n2.725  \n0.0  \n603  \n1248  \n2148  \n4582  \n7283  \n644  \n1117  \n2123  \n4262  \n8454  \n4.984  \n4.985  \n4.985  \n4.985  \n4.986  \n-  \n3  \n-  \n3  \n-  \n3  \n]  \n4.243  \n4.245  \n4.246  \n4.247  \n4.243  \n]  \n3.525  \n3.545  \n3.548  \n3.55  \n3.551  \n]  \n3.151  \n3.383  \n3.392  \n3.392  \n-', 'REFERENCES': "- Aerts, C., Molenberghs, G., Kenward, M. G., & Neiner, C. 2014, ApJ, 781, 88  \nAlastuey, A., & Jancovici, B. 1978, ApJ, 226, 1034 Alecian, G., Michaud, G., & Tully, J. 1993, ApJ, 411, 882 An, D., Pinsonneault, M. H., Masseron, T., et al. 2009, ApJ, 700, 523 Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A, 656, 3 Antia, H. M., & Basu, S. 2005, ApJL, 620, L129 Aringer, B., Girardi, L., Nowotny, W., Marigo, P., & Lederer, M. T. 2009,  \n- A&A, 503, 913  \nArnett, W. D., Meakin, C., Viallet, M., et al. 2015, ApJ, 809, 30  \n- Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 336, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, ed. T. G. Barnes, III & F. N. Bash, 25\n- Asplund, M., Grevesse, N., Sauval, A. J., Allende Prieto, C., & Kiselman, D. 2004, A&A, 417, 751\n- Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 Baglin, A., Auvergne, M., Barge, P., et al. 2006, in ESA Special Publication, Vol. 1306, ESA Special Publication, ed. M. Fridlund, A. Baglin, J. Lochard, & L. Conroy, 33\n- Bahcall, J. N., Basu, S., Pinsonneault, M., & Serenelli, A. M. 2005a, ApJ, 618, 1049\n- Bahcall, J. N., Basu, S., & Serenelli, A. M. 2005b, ApJ, 631, 1281\n- Bahcall, J. N., Serenelli, A. M., & Basu, S. 2006, ApJS, 165, 400\n- Bailey, J. E., Nagayama, T., Loisel, G. P., et al. 2015, Nature, 517, 56Ð59\n- Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1995, ApJL, 446, L35\n- -. 1997, A&A, 327, 1054\n- -. 1998, A&A, 337, 403\n- Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701\n- Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJL, 702, L27\n- Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42 Basu, S., & Antia, H. M. 2004, ApJL, 606, L85\n- Bedin, L. R., Salaris, M., Piotto, G., et al. 2005, ApJL, 624, L45\n- Bessell, M., Bloxham, G., Schmidt, B., et al. 2011, PASP, 123, 789\n- Bessell, M., & Murphy, S. 2012, PASP, 124, 140\n- Bessell, M. S. 2001, PASP, 113, 66\n- -. 2011, PASP, 123, 1442\n- Bird, J. C., & Pinsonneault, M. H. 2011, ApJ, 733, 81  \nBjorkman, J. E., & Cassinelli, J. P. 1993, ApJ, 409, 429  \nBlöcker, T. 1995, A&A, 297, 727  \nBöhm, K.-H. 1963, ApJ, 138, 297  \n- Böhm-Vitense, E. 1958, ZAp, 46, 108\n- Boudreault, S., Lodieu, N., Deacon, N. R., & Hambly, N. C. 2012, MNRAS, 426, 3419\n- Bouret, J.-C., Lanz, T., Martins, F., et al. 2013, A&A, 555, A1 Bouvier, J. 2008, A&A, 489, L53\n- Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112 Boyer, M. L., Srinivasan, S., van Loon, J. T., et al. 2011, AJ, 142, 103 Braithwaite, J. 2006, A&A, 449, 451\n- Braithwaite, J., & Spruit, H. C. 2004, Nature, 431, 819\n- Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1\n- Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127\n- Bressan, A. G., Chiosi, C., & Bertelli, G. 1981, A&A, 102, 25\n- Brogaard, K., VandenBerg, D. A., Bruntt, H., et al. 2012, A&A, 543, A106\n- Brott, I., de Mink, S. E., Cantiello, M., et al. 2011a, A&A, 530, A115\n- Brott, I., Evans, C. J., Hunter, I., et al. 2011b, A&A, 530, A116\n- Brown, J. A., Wallerstein, G., & Zucker, D. 1997, AJ, 114, 180\n- Brown, J. M., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 34\n- Brown, T. M., Christensen-Dalsgaard, J., Dziembowski, W. A., et al. 1989, ApJ, 343, 526\n- Brown, T. M., Gilliland, R. L., Noyes, R. W., & Ramsey, L. W. 1991, ApJ, 368, 599\n- Burbidge, E. M., Burbidge, G. R., Fowler, W. A., & Hoyle, F. 1957, Reviews of Modern Physics, 29, 547\n- Caffau, E., Ludwig, H.-G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255\n- Caloi, V., & D'Antona, F. 2011, MNRAS, 417, 228\n- Cantiello, M., & Langer, N. 2010, A&A, 521, A9\n- Cantiello, M., Mankovich, C., Bildsten, L., Christensen-Dalsgaard, J., & Paxton, B. 2014, ApJ, 788, 93\n- Canto Martins, B. L., Lèbre, A., Palacios, A., et al. 2011, A&A, 527, A94 Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245\n- Carney, B. W., Lee, J.-W., & Dodson, B. 2005, AJ, 129, 656  \nCarrera, R., & Pancino, E. 2011, A&A, 535, A30 Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2009, A&A, 505, 117 Casagrande, L., Flynn, C., & Bessell, M. 2008, MNRAS, 389, 585 Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094  \n- Cassisi, S., & Salaris, M. 2013, Old Stellar Populations: How to Study the Fossil Record of Galaxy Formation  \nCastelli, F., Gratton, R. G., & Kurucz, R. L. 1997, A&A, 318, 841 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 Castro, N., Fossati, L., Langer, N., et al. 2014, A&A, 570, L13 Catalán, S., Isern, J., García-Berro, E., & Ribas, I. 2008a, MNRAS, 387,  \n- 1693\n- Catalán, S., Isern, J., García-Berro, E., et al. 2008b, A&A, 477, 213 Chaboyer, B., Fenton, W. H., Nelan, J. E., Patnaude, D. J., & Simon, F. E. 2001, ApJ, 562, 521\n- Chaboyer, B., & Zahn, J.-P. 1992, A&A, 253, 173\n- Chabrier, G., Baraffe, I., & Plez, B. 1996, ApJL, 459, L91\n- Charbonneau, P., & MacGregor, K. B. 1993, ApJ, 417, 762\n- Charbonnel, C., & Lagarde, N. 2010, A&A, 522, A10\n- Charbonnel, C., & Talon, S. 2005, Science, 309, 2189\n- Charbonnel, C., & Zahn, J.-P. 2007, A&A, 467, L15\n- Chen, C. H., Mamajek, E. E., Bitner, M. A., et al. 2011, ApJ, 738, 122\n- Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525 Chieffi, A., & Limongi, M. 2013, ApJ, 764, 21\n- Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925\n- Christensen-Dalsgaard, J., di Mauro, M. P., Houdek, G., & Pijpers, F. 2009, A&A, 494, 205\n- Claver, C. F., Liebert, J., Bergeron, P., & Koester, D. 2001, ApJ, 563, 987 Cohen, M., Wheaton, W. A., & Megeath, S. T. 2003, AJ, 126, 1090\n- Conroy, C. 2013, ARA&A, 51, 393\n- Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833\n- Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486\n- Cowley, A. P., Dawson, P., & Hartwick, F. D. A. 1979, PASP, 91, 628  \nCox, A. N., & Stewart, J. N. 1965, ApJS, 11, 22  \n-. 1970, ApJS, 19, 243  \n- Cox, J. P. 1968, Principles of stellar structure - Vol.1: Physical principles; Vol.2: Applications to stars\n- Crowther, P. A., Lennon, D. J., & Walborn, N. R. 2006, A&A, 446, 279 Cummings, J. D., Kalirai, J. S., Tremblay, P.-E., & Ramirez-Ruiz, E. 2015, ApJ, 807, 90\n- Cunha, K., Hubeny, I., & Lanz, T. 2006, ApJL, 647, L143\n- Cunha, K., & Lambert, D. L. 1994, ApJ, 426, 170\n- Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240\n- Daflon, S., Cunha, K., Becker, S. R., & Smith, V. V. 2001, ApJ, 552, 309  \nDalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18  \n- D'Antona, F., & Caloi, V. 2008, MNRAS, 390, 693\n- D'Antona, F., & Mazzitelli, I. 1997, Mem. Soc. Astron. Italiana, 68, 807\n- de Geus, E. J., de Zeeuw, P. T., & Lub, J. 1989, A&A, 216, 44\n- de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259\n- de Mink, S. E., Cantiello, M., Langer, N., & Pols, O. R. 2010, in American Institute of Physics Conference Series, Vol. 1314, American Institute of Physics Conference Series, ed. V. Kalogera & M. van der Sluys, 291-296  \nde Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243  \n- de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764, 166\n- de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7\n- De Silva, G. M., Freeman, K. C., Bland-Hawthorn, J., et al. 2015, MNRAS, 449, 2604\n- de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354\n- Demarque, P., Woo, J.-H., Kim, Y.-C., & Yi, S. K. 2004, ApJS, 155, 667 Denissenkov, P. A. 2010, ApJ, 723, 563\n- Denissenkov, P. A., & Pinsonneault, M. 2007, ApJ, 655, 1157\n- Dobbie, P. D., Pinfield, D. J., Napiwotzki, R., et al. 2004, MNRAS, 355, L39 Doherty, C. L., Gil-Pons, P., Siess, L., Lattanzio, J. C., & Lau, H. H. B.\n- 2015, MNRAS, 446, 2599  \nDohm-Palmer, R. C., & Skillman, E. D. 2002, AJ, 123, 1433  \n- Dolphin, A. E. 2002, MNRAS, 332, 91\n- Dotter, A. 2016, ApJS, 222, 8\n- Dotter, A., Chaboyer, B., Jevremovi'c, D., et al. 2008, ApJS, 178, 89\n- Dotter, A., Sarajedini, A., & Anderson, J. 2011, ApJ, 738, 74\n- Drake, J. J., & Testa, P. 2005, Nature, 436, 525\n- Drout, M. R., Massey, P., & Meynet, G. 2012, ApJ, 750, 97\n- Eggenberger, P., Maeder, A., & Meynet, G. 2005, A&A, 440, L9\n- Eggenberger, P., Meynet, G., & Maeder, A. 2002, A&A, 386, 576\n- Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43  \nEggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4  \n- Eggleton, P. P. 1971, MNRAS, 151, 351\n- Eggleton, P. P., Dearborn, D. S. P., & Lattanzio, J. C. 2006, Science, 314, 1580  \nEkström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146  \nEldridge, J. J., Izzard, R. G., & Tout, C. A. 2008, MNRAS, 384, 1109  \nEldridge, J. J., & Tout, C. A. 2004, MNRAS, 348, 201  \nElias, J. H., Frogel, J. A., & Humphreys, R. M. 1985, ApJS, 57, 91  \nEndal, A. S., & Sofia, S. 1976, ApJ, 210, 184  \n-. 1978, ApJ, 220, 279  \n- Faulkner, D. J., & Cannon, R. D. 1973, ApJ, 180, 435\n- Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10\n- Feiden, G. A., & Chaboyer, B. 2012, ApJ, 757, 42\n- -. 2013, ApJ, 779, 183\n- Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 Ferrario, L., Wickramasinghe, D., Liebert, J., & Williams, K. A. 2005, MNRAS, 361, 1131\n- Fontaine, G., Brassard, P., Green, E. M., et al. 2008, A&A, 486, L39 Francois, P., Danziger, J., Buonanno, R., & Perrin, M. N. 1997, A&A, 327, 121  \nFreedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504  \nFreytag, B., Ludwig, H.-G., & Steffen, M. 1996, A&A, 313, 497  \n- Fricke, K. J., & Kippenhahn, R. 1972, ARA&A, 10, 45\n- Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701\n- Frogel, J. A., Mould, J., & Blanco, V. M. 1990, ApJ, 352, 96\n- Frommhold, L., Abel, M., Wang, F., et al. 2010, Molecular Physics, 108, 2265\n- Fynbo, H. O. U., Diget, C. A., Bergmann, U. C., et al. 2005, Nature, 433, 136  \nGallet, F., & Bouvier, J. 2013, A&A, 556, A36 Gáspár, A., Rieke, G. H., Su, K. Y. L., et al. 2009, ApJ, 697, 1578 Gathier, R., Lamers, H. J. G. L. M., & Snow, T. P. 1981, ApJ, 247, 173 Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29 Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 Gies, D. R., & Lambert, D. L. 1992, ApJ, 387, 673 Gilliland, R. L., Brown, T. M., Christensen-Dalsgaard, J., et al. 2010, PASP,  \n- 122, 131\n- Girardi, L. 1999, MNRAS, 308, 818\n- Girardi, L., Bertelli, G., Bressan, A., et al. 2002, A&A, 391, 195\n- Girardi, L., & Marigo, P. 2007, A&A, 462, 237\n- Girardi, L., Marigo, P., Bressan, A., & Rosenfield, P. 2013, ApJ, 777, 142\n- Girardi, L., Dalcanton, J., Williams, B., et al. 2008, PASP, 120, 583\n- Girardi, L., Williams, B. F., Gilbert, K. M., et al. 2010, ApJ, 724, 1030\n- Glebbeek, E., Gaburov, E., de Mink, S. E., Pols, O. R., & Portegies Zwart, S. F. 2009, A&A, 497, 255\n- Gordon, K. D., Meixner, M., Meade, M. R., et al. 2011, AJ, 142, 102 Graboske, H. C., Dewitt, H. E., Grossman, A. S., & Cooper, M. S. 1973, ApJ, 181, 457  \nGräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56  \nGratton, R., Sneden, C., & Carretta, E. 2004, ARA&A, 42, 385  \nGratton, R. G., Carretta, E., & Bragaglia, A. 2012, A&A Rev., 20, 50  \n- Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161  \nGroenewegen, M. A. T., & de Jong, T. 1993, A&A, 267, 410  \n- Grossman, S. A., & Taam, R. E. 1996, MNRAS, 283, 1165  \nHarris, J., & Zaritsky, D. 2004, AJ, 127, 1531  \n-. 2009, AJ, 138, 1243  \n- Hartwick, F. D. A. 1970, Astrophys. Lett., 7, 151\n- Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377\n- Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368\n- Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532\n- Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350\n- Henriques, B., Maraston, C., Monaco, P., et al. 2011, MNRAS, 415, 3571\n- Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841  \nHenyey, L. G., Forbes, J. E., & Gould, N. L. 1964, ApJ, 139, 306  \nHerbig, G. H. 1962, ApJ, 135, 736  \n- Herczeg, G. J., & Hillenbrand, L. A. 2015, ApJ, 808, 23\n- Herwig, F. 2000, A&A, 360, 952\n- -. 2004, ApJS, 155, 651\n- -. 2005, ARA&A, 43, 435\n- Herwig, F., & Austin, S. M. 2004, ApJL, 613, L73\n- Hewett, P. C., Warren, S. J., Leggett, S. K., & Hodgkin, S. T. 2006, MNRAS, 367, 454  \nHirschi, R., Meynet, G., & Maeder, A. 2004, A&A, 425, 649 Holtzman, J. A., Burrows, C. J., Casertano, S., et al. 1995, PASP, 107, 1065  \nHosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140 Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398 Hu, H., Tout, C. A., Glebbeek, E., & Dupret, M.-A. 2011, MNRAS, 418, 195 Huang, W., Gies, D. R., & McSwain, M. V. 2010, ApJ, 722, 605 Humphreys, R. M. 1979, ApJS, 39, 389  \n- Hunter, I., Dufton, P. L., Smartt, S. J., et al. 2007, A&A, 466, 277\n- Hunter, I., Brott, I., Lennon, D. J., et al. 2008, ApJL, 676, L29\n- Hunter, I., Brott, I., Langer, N., et al. 2009, A&A, 496, 841  \nHurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897  \nIben, Jr., I. 1982, ApJ, 260, 821  \n-. 1983, ApJL, 275, L65  \n- Iben, Jr., I., Kaler, J. B., Truran, J. W., & Renzini, A. 1983, ApJ, 264, 605 Iglesias, C. A., & Rogers, F. J. 1991, ApJ, 371, 408\n- -. 1993, ApJ, 412, 752\n- -. 1996, ApJ, 464, 943\n- Irwin, J. M., Quinn, S. N., Berta, Z. K., et al. 2011, ApJ, 742, 123\n- Itoh, N., Totsuji, H., Ichimaru, S., & Dewitt, H. E. 1979, ApJ, 234, 1079\n- Ivezic, Z., Tyson, J. A., Abel, B., et al. 2008, ArXiv e-prints, arXiv:0805.2366\n- Izzard, R. G., Dray, L. M., Karakas, A. I., Lugaro, M., & Tout, C. A. 2006, A&A, 460, 565\n- Jackson, R. J., & Jeffries, R. D. 2014, MNRAS, 441, 2111\n- Jeffries, R. D., & Oliveira, J. M. 2005, MNRAS, 358, 13\n- Jiang, Y.-F., Cantiello, M., Bildsten, L., Quataert, E., & Blaes, O. 2015, ApJ, 813, 74\n- Juarez, A. J., Cargile, P. A., James, D. J., & Stassun, K. G. 2014, ApJ, 795, 143\n- Kaiser, N., Burgett, W., Chambers, K., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7733, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 0\n- Kalirai, J. S., Bergeron, P., Hansen, B. M. S., et al. 2007, ApJ, 671, 748\n- Kalirai, J. S., Hansen, B. M. S., Kelson, D. D., et al. 2008, ApJ, 676, 594\n- Kalirai, J. S., Marigo, P., & Tremblay, P.-E. 2014, ApJ, 782, 17\n- Kalirai, J. S., Saul Davis, D., Richer, H. B., et al. 2009, ApJ, 705, 408\n- Kaluzny, J., Krzeminski, W., & Mazur, B. 1995, AJ, 110, 2206\n- Kamai, B. L., Vrba, F. J., Stauffer, J. R., & Stassun, K. G. 2014, AJ, 148, 30 Kawaler, S. D., Novikov, I., Srinivasan, G., Meynet, G., & Schaerer, D. 1996, Stellar Remnants\n- Kepler, S. O., Pelisoli, I., Koester, D., et al. 2015, MNRAS, 446, 4078 Kilian, J. 1992, A&A, 262, 171\n- Kim, Y.-C., Demarque, P., Yi, S. K., & Alexander, D. R. 2002, ApJS, 143, 499\n- Kippenhahn, R., Ruschenplatt, G., & Thomas, H.-C. 1980, A&A, 91, 175 Kippenhahn, R., & Thomas, H.-C. 1970, in IAU Colloq. 4: Stellar Rotation, ed. A. Slettebak, 20\n- Kjeldsen, H., & Bedding, T. R. 1995, A&A, 293, 87\n- Kleinman, S. J., Kepler, S. O., Koester, D., et al. 2013, ApJS, 204, 5\n- Koester, D. 2010, Mem. Soc. Astron. Italiana, 81, 921  \nKöhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71  \nKraft, R. P. 1967, ApJ, 150, 551  \n-. 1994, PASP, 106, 553  \n- Kraus, A. L., Cody, A. M., Covey, K. R., et al. 2015, ApJ, 807, 3\n- Kraus, A. L., & Hillenbrand, L. A. 2007, AJ, 134, 2340\n- Kraus, A. L., Tucker, R. A., Thompson, M. I., Craine, E. R., & Hillenbrand, L. A. 2011, ApJ, 728, 48\n- Kroupa, P. 2001, MNRAS, 322, 231\n- Kurucz, R. L. 1970, SAO Special Report, 309\n- -. 1993, SYNTHE spectrum synthesis programs and line data\n- Lafrenière, D., Jayawardhana, R., & van Kerkwijk, M. H. 2008, ApJL, 689, L153  \nLanger, N. 1991, A&A, 252, 669  \n- -. 1998, A&A, 329, 551\n- -. 2012, ARA&A, 50, 107\n- Langer, N., El Eid, M. F., & Baraffe, I. 1989, A&A, 224, L17\n- Langer, N., El Eid, M. F., & Fricke, K. J. 1985, A&A, 145, 179\n- Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207\n- Langer, N., & Kudritzki, R. P. 2014, A&A, 564, A52\n- Langer, N., & Maeder, A. 1995, A&A, 295, 685  \nLattanzio, J. C. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 378, Why Galaxies Care About AGB Stars: Their Importance as Actors and Probes, ed. F. Kerschbaum, C. Charbonnel, & R. F. Wing, 3 Lau, H. H. B., Stancliffe, R. J., & Tout, C. A. 2009, MNRAS, 396, 1046 Law, N. M., Kulkarni, S. R., Dekany, R. G., et al. 2009, PASP, 121, 1395 Levesque, E. M., Massey, P., Olsen, K. A. G., et al. 2005, ApJ, 628, 973 -. 2006, ApJ, 645, 1102  \nLiebert, J., Bergeron, P., & Holberg, J. B. 2005, ApJS, 156, 47  \n- Lodders, K., Palme, H., & Gail, H.-P. 2009, Landolt Börnstein, 44 Lodieu, N. 2013, MNRAS, 431, 3222\n- Lodieu, N., Hambly, N. C., Jameson, R. F., & Hodgkin, S. T. 2008, MNRAS, 383, 1385\n- Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879\n- MacDonald, J., & Mullan, D. J. 2012, MNRAS, 421, 3084\n- -. 2014, ApJ, 787, 70\n- Maeder, A. 1975, A&A, 40, 303\n- Maeder, A., & Conti, P. S. 1994, ARA&A, 32, 227\n- Maeder, A., & Meynet, G. 1994, A&A, 287, 803\n- -. 2000, ARA&A, 38, 143\n- -. 2001, A&A, 373, 555\n- -. 2012, Reviews of Modern Physics, 84, 25\n- Magic, Z., Serenelli, A., Weiss, A., & Chaboyer, B. 2010, ApJ, 718, 1378 Magic, Z., Weiss, A., & Asplund, M. 2015, A&A, 573, A89\n- Maraston, C. 2005, MNRAS, 362, 799\n- Marigo, P. 2002, A&A, 387, 507\n- Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488\n- Marigo, P., & Girardi, L. 2001, A&A, 377, 132\n- -. 2007, A&A, 469, 239\n- Marigo, P., Girardi, L., & Bressan, A. 1999, A&A, 344, 123\n- Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883\n- Marigo, P., Girardi, L., & Chiosi, C. 2003, A&A, 403, 225\n- Massey, P. 2002, ApJS, 141, 81\n- Massey, P., & Olsen, K. A. G. 2003, AJ, 126, 2867\n- Massey, P., Zangari, A. M., Morrell, N. I., et al. 2009, ApJ, 692, 618\n- Mauron, N., & Josselin, E. 2011, A&A, 526, A156\n- McQuinn, K. B. W., Skillman, E. D., Dalcanton, J. J., et al. 2011, ApJ, 740, 48\n- Meixner, M., Gordon, K. D., Indebetouw, R., et al. 2006, AJ, 132, 2268 Melbourne, J., Williams, B. F., Dalcanton, J. J., et al. 2012, ApJ, 748, 47 Melis, C., Reid, M. J., Mioduszewski, A. J., Stauffer, J. R., & Bower, G. C.\n- 2014, Science, 345, 1029\n- Mestel, L. 1968, MNRAS, 138, 359\n- Mestel, L., & Weiss, N. O. 1987, MNRAS, 226, 123\n- Meylan, G., & Maeder, A. 1982, A&A, 108, 148\n- Meynet, G., Eggenberger, P., & Maeder, A. 2011, A&A, 525, L11\n- Meynet, G., & Maeder, A. 1997, A&A, 321, 465\n- -. 2000, A&A, 361, 101\n- -. 2005, A&A, 429, 581\n- Michaud, G., Fontaine, G., & Beaudet, G. 1984, ApJ, 282, 206\n- Michaud, G., Richard, O., Richer, J., & VandenBerg, D. A. 2004, ApJ, 606, 452\n- Miglio, A., Brogaard, K., Stello, D., et al. 2012, MNRAS, 419, 2077 Miller Bertolami, M. M. 2016, A&A, 588, A25\n- Morales, J. C., Gallardo, J., Ribas, I., et al. 2010, ApJ, 718, 502\n- Morales, J. C., Ribas, I., & Jordi, C. 2008, A&A, 478, 507\n- Morel, P., & Baglin, A. 1999, A&A, 345, 156\n- Morel, P., & Thévenin, F. 2002, A&A, 390, 611\n- Morel, T., Hubrig, S., & Briquet, M. 2008, A&A, 481, 453\n- Neugent, K. F., & Massey, P. 2011, ApJ, 733, 123\n- Neugent, K. F., Massey, P., & Georgy, C. 2012, ApJ, 759, 11\n- Nieuwenhuijzen, H., & de Jager, C. 1988, A&A, 203, 355 -. 1990, A&A, 231, 134\n- Noël, N. E. D., Greggio, L., Renzini, A., Carollo, C. M., & Maraston, C. 2013, ApJ, 772, 58\n- Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227\n- Origlia, L., Valenti, E., Rich, R. M., & Ferraro, F. R. 2006, ApJ, 646, 499 Paczy'nski, B. 1970, Acta Astron., 20, 47\n- Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603\n- Paquette, C., Pelletier, C., Fontaine, G., & Michaud, G. 1986, ApJS, 61, 177\n- Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3\n- Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4\n- Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15\n- Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. 2012, ApJ, 746, 154\n- Petrovic, J., Langer, N., Yoon, S.-C., & Heger, A. 2005, A&A, 435, 247  \nPietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168  \nPinsonneault, M. 1997, ARA&A, 35, 557  \n- Pinsonneault, M. H., Kawaler, S. D., Sofia, S., & Demarque, P. 1989, ApJ, 338, 424\n- Piotto, G., Milone, A. P., Anderson, J., et al. 2012, ApJ, 760, 39\n- Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015, ArXiv e-prints, arXiv:1502.01589\n- Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964\n- Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82\n- Potter, A. T., Tout, C. A., & Eldridge, J. J. 2012, MNRAS, 419, 748 Prada Moroni, P. G., & Straniero, O. 2002, ApJ, 581, 585\n- Preibisch, T., Brown, A. G. A., Bridges, T., Guenther, E., & Zinnecker, H. 2002, AJ, 124, 404\n- Preibisch, T., & Zinnecker, H. 1999, AJ, 117, 2381\n- Przybilla, N., Firnstein, M., Nieva, M. F., Meynet, G., & Maeder, A. 2010, A&A, 517, A38\n- Rauscher, T., Heger, A., Hoffman, R. D., & Woosley, S. E. 2002, ApJ, 576, 323\n- Reimers, D. 1975, Memoires of the Societe Royale des Sciences de Liege, 8, 369\n- Renzini, A., & Fusi Pecci, F. 1988, ARA&A, 26, 199\n- Rhodes, Jr., E. J., Kosovichev, A. G., Schou, J., Scherrer, P. H., & Reiter, J. 1997, Sol. Phys., 175, 287\n- Richard, O., Vauclair, S., Charbonnel, C., & Dziembowski, W. A. 1996, A&A, 312, 1000\n- Rivero González, J. G., Puls, J., Najarro, F., & Brott, I. 2012, A&A, 537, A79\n- Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064\n- Romero, A. D., Campos, F., & Kepler, S. O. 2015, MNRAS, 450, 3708\n- Rosen, A. L., Krumholz, M. R., & Ramirez-Ruiz, E. 2012, ApJ, 748, 97  \nRosenfield, P., Marigo, P., Girardi, L., et al. 2014, ApJ, 790, 22  \nRoxburgh, I. W. 1978, A&A, 65, 281  \n- Rüdiger, G., Gellert, M., Spada, F., & Tereshin, I. 2015, A&A, 573, A80 Salaris, M., & Cassisi, S. 2015, A&A, 577, A60\n- Salaris, M., Groenewegen, M. A. T., & Weiss, A. 2000, A&A, 355, 299\n- Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 Sanders, N. E., Caldwell, N., McDowell, J., & Harding, P. 2012, ApJ, 758,\n- 133\n- Sandquist, E. L. 2004, MNRAS, 347, 101\n- Sarajedini, A., Dotter, A., & Kirkpatrick, A. 2009, ApJ, 698, 1872\n- Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713\n- Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269 Schoenberner, D. 1979, A&A, 79, 108\n- Schwarzschild, M. 1958, Structure and evolution of the stars.\n- Schwarzschild, M., & Härm, R. 1962, ApJ, 136, 158\n- -. 1965, ApJ, 142, 855\n- Searle, S. C., Prinja, R. K., Massa, D., & Ryans, R. 2008, A&A, 481, 777 Seaton, M. J. 2005, MNRAS, 362, L1\n- Serenelli, A. M., Basu, S., Ferguson, J. W., & Asplund, M. 2009, ApJL, 705, L123\n- Sestito, P., & Randich, S. 2005, A&A, 442, 615\n- Shaviv, G., & Salpeter, E. E. 1973, ApJ, 184, 191\n- Siess, L. 2009, A&A, 497, 463\n- Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593\n- Silva Aguirre, V., Ballot, J., Serenelli, A. M., & Weiss, A. 2011, A&A, 529, A63\n- Slesnick, C. L., Hillenbrand, L. A., & Carpenter, J. M. 2008, ApJ, 688, 377 Smartt, S. J. 2009, ARA&A, 47, 63\n- Smith, N. 2014, ARA&A, 52, 487\n- Smith, N., & Owocki, S. P. 2006, ApJL, 645, L45\n- Soderblom, D. R., Laskar, T., Valenti, J. A., Stauffer, J. R., & Rebull, L. M. 2009, AJ, 138, 1292\n- Soderblom, D. R., Nelan, E., Benedict, G. F., et al. 2005, AJ, 129, 1616 Söderlind, G., & Wang, L. 2006, Journal of Computational and Applied Mathematics, 185, 244 , special Issue: International Workshop on the Technological Aspects of Mathematics\n- Somers, G., & Pinsonneault, M. H. 2014, ApJ, 790, 72\n- Song, I., Zuckerman, B., & Bessell, M. S. 2012, AJ, 144, 8\n- Southworth, J. 2015, in Astronomical Society of the Pacific Conference Series, Vol. 496, Living Together: Planets, Host Stars and Binaries, ed. S. M. Rucinski, G. Torres, & M. Zejda, 164\n- Spada, F., Demarque, P., Kim, Y.-C., & Sills, A. 2013, ApJ, 776, 87 Spruit, H. C. 2002, A&A, 381, 923\n- -. 2013, A&A, 552, A76\n- Spruit, H. C., & Weiss, A. 1986, A&A, 166, 167\n- Stancliffe, R. J. 2010, MNRAS, 403, 505\n- Stancliffe, R. J., Dearborn, D. S. P., Lattanzio, J. C., Heap, S. A., &\n- Campbell, S. W. 2011, ApJ, 742, 121\n- Stancliffe, R. J., Glebbeek, E., Izzard, R. G., & Pols, O. R. 2007, A&A, 464, L57\n- Stauffer, J. R., Jones, B. F., Backman, D., et al. 2003, AJ, 126, 833\n- Stauffer, J. R., Hartmann, L. W., Fazio, G. G., et al. 2007, ApJS, 172, 663\n- Stetson, P. B., Bruntt, H., & Grundahl, F. 2003, PASP, 115, 413\n- Stothers, R., & Chin, C.-W. 1975, ApJ, 198, 407\n- Stothers, R., & Simon, N. R. 1969, ApJ, 157, 673\n- Stothers, R. B., & Chin, C.-W. 1991, ApJL, 381, L67\n- Strittmatter, P. A. 1969, ARA&A, 7, 665\n- Suijs, M. P. L., Langer, N., Poelarends, A.-J., et al. 2008, A&A, 481, L87 Sukhbold, T., & Woosley, S. E. 2014, ApJ, 783, 10\n- Takeda, Y., & Takada-Hidai, M. 2000, PASJ, 52, 113\n- Talon, S., & Charbonnel, C. 2005, A&A, 440, 981\n- Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287\n- Tang, J., Bressan, A., Slemer, A., et al. 2016, MNRAS, 455, 3393\n- Tassoul, J.-L. 1978, Theory of rotating stars\n- Taylor, B. J. 2006, AJ, 132, 2453\n- -. 2007, AJ, 133, 370\n- Thomas, H.-C. 1967, ZAp, 67, 420\n- Thoul, A. A., Bahcall, J. N., & Loeb, A. 1994, ApJ, 421, 828\n- Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501\n- Tognelli, E., Prada Moroni, P. G., & Degl'Innocenti, S. 2011, A&A, 533, A109\n- Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99\n- Torres, G. 2013, Astronomische Nachrichten, 334, 4\n- Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67\n- Torres, G., Ruíz-Rodríguez, D., Badenas, M., et al. 2013, ApJ, 773, 40\n- Townsend, R. H. D., & Teitler, S. A. 2013, MNRAS, 435, 3406\n- Trampedach, R. 2007, in American Institute of Physics Conference Series, Vol. 948, Unsolved Problems in Stellar Physics: A Conference in Honor of Douglas Gough, ed. R. J. Stancliffe, G. Houdek, R. G. Martin, & C. A. Tout, 141-148\n- Trampedach, R., Stein, R. F., Christensen-Dalsgaard, J., Nordlund, Å., & Asplund, M. 2014, MNRAS, 445, 4366\n- Traxler, A., Garaud, P., & Stellmach, S. 2011, ApJL, 728, L29\n- Trundle, C., Dufton, P. L., Lennon, D. J., Smartt, S. J., & Urbaneja, M. A. 2002, A&A, 395, 519\n- Turcotte, S., Richer, J., & Michaud, G. 1998, ApJ, 504, 559\n- ud-Doula, A., & Owocki, S. P. 2002, ApJ, 576, 413\n- Ulrich, R. K. 1972, ApJ, 172, 165\n- -. 1986, ApJL, 306, L37\n- van den Bergh, S. 1968, JRASC, 62, 145\n- van der Hucht, K. A. 2001, NewAR, 45, 135\n- van Saders, J. L., & Pinsonneault, M. H. 2012, ApJ, 751, 98\n- VandenBerg, D. A., Bergbusch, P. A., & Dowler, P. D. 2006, ApJS, 162, 375\n- Vardya, M. S. 1985, ApJ, 299, 255\n- Varenne, O., & Monier, R. 1999, A&A, 351, 247\n- Vassiliadis, E., & Wood, P. R. 1993, ApJ, 413, 641\n- -. 1994, ApJS, 92, 125\n- Vauclair, S. 1983, in Saas-Fee Advanced Course 13: Astrophysical Processes in Upper Main Sequence Stars, ed. A. N. Cox, S. Vauclair, & J. P. Zahn, 167  \nVenn, K. A., McCarthy, J. K., Lennon, D. J., et al. 2000, ApJ, 541, 610 Villamariz, M. R., & Herrero, A. 2005, A&A, 442, 263  \n- Villamariz, M. R., Herrero, A., Becker, S. R., & Butler, K. 2002, A&A, 388, 940\n- Villanova, S., Geisler, D., Carraro, G., Moni Bidin, C., & Muñoz, C. 2013, ApJ, 778, 186\n- Villaume, A., Conroy, C., & Johnson, B. D. 2015, ApJ, 806, 82 Vink, J. S., & de Koter, A. 2005, A&A, 442, 587\n- Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 -. 2001, A&A, 369, 574\n- Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132 Wachlin, F. C., Miller Bertolami, M. M., & Althaus, L. G. 2011, A&A, 533, A139\n- Wachlin, F. C., Vauclair, S., & Althaus, L. G. 2014, A&A, 570, A58 Wagenhuber, J., & Groenewegen, M. A. T. 1998, A&A, 340, 183 Walker, M. F. 1964, AJ, 69, 744\n- Wang, P. F., Chen, W. P., Lin, C. C., et al. 2014, ApJ, 784, 57\n- Weber, E. J., & Davis, Jr., L. 1967, ApJ, 148, 217\n- Weidemann, V. 1977, A&A, 59, 411\n- -. 2000, A&A, 363, 647\n- Weiss, A., & Ferguson, J. W. 2009, A&A, 508, 1343\n- Weiss, A., & Schlattl, H. 2008, Ap&SS, 316, 99\n- Wheeler, J. C., Kagan, D., & Chatzopoulos, E. 2015, ApJ, 799, 85\n- Williams, K. A., & Bolte, M. 2007, AJ, 133, 1490\n- Williams, K. A., Bolte, M., & Koester, D. 2009, ApJ, 693, 355\n- Willson, L. A. 2000, ARA&A, 38, 573\n- Wolff, S., & Simon, T. 1997, PASP, 109, 759\n- Wood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157\n- Woodward, P. R., Herwig, F., & Lin, P.-H. 2015, ApJ, 798, 49\n- Woosley, S. E., & Heger, A. 2006, ApJ, 637, 914\n- Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868\n- Yi, S., Demarque, P., Kim, Y.-C., et al. 2001, ApJS, 136, 417\n- Yi, S. K., Kim, Y.-C., & Demarque, P. 2003, ApJS, 144, 259\n- Yoon, S.-C., & Cantiello, M. 2010, ApJL, 717, L62\n- Yoon, S.-C., & Langer, N. 2005, A&A, 443, 643\n- Yoon, S.-C., Langer, N., & Norman, C. 2006, A&A, 460, 199\n- Young, P. A., & Arnett, D. 2005, ApJ, 618, 908\n- Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, MNRAS, 413, 1828 Zahn, J.-P. 1983, in Saas-Fee Advanced Course 13: Astrophysical Processes in Upper Main Sequence Stars, ed. A. N. Cox, S. Vauclair, & J. P. Zahn, 253\n- Zahn, J.-P. 1992, A&A, 265, 115\n- Zahn, J.-P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145\n- Zaritsky, D., Harris, J., Thompson, I. B., Grebel, E. K., & Massey, P. 2002, AJ, 123, 855\n- Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, AJ, 146, 81 Zhang, F., Li, L., Han, Z., Zhuang, Y., & Kang, X. 2013, MNRAS, 428, 3390 Zhao, J. K., Oswalt, T. D., Willson, L. A., Wang, Q., & Zhao, G. 2012, ApJ, 746, 144  \nZurita, A., & Bresolin, F. 2012, MNRAS, 427, 1463"}