We substantially update the capabilities of the open source software package Modules for Experiments in Stellar Astrophysics MESA and its onedimensional stellar evolution module MESA star. Improvements in MESA stars ability to model the evolution of giant planets now extends its applicability down to masses as low as onetenth that of Jupiter. The dramatic improvement in asteroseismology enabled by the spacebased Kepler and CoRoT missions motivates our full coupling of the ADIPLS adiabatic pulsation code with MESA star. This also motivates a numerical recasting of the Ledoux criterion that is more easily implemented when many nuclei are present at nonnegligible abundances. This impacts the way in which MESA star calculates semiconvective and thermohaline mixing. We exhibit the evolution of 38 M SUBSUB stars through the end of core He burning the onset of He thermal pulses and arrival on the white dwarf cooling sequence. We implement diffusion of angular momentum and chemical abundances that enable calculations of rotatingstar models which we compare thoroughly with earlier work. We introduce a new treatment of radiationdominated envelopes that allows the uninterrupted evolution of massive stars to core collapse. This enables the generation of new sets of supernovae long gammaray burst and pairinstability progenitor models. We substantially modify the way in which MESA star solves the fully coupled stellar structure and composition equations and we show how this has improved the scaling of MESAs calculational speed on multicore processors. Updates to the modules for equation of state opacity nuclear reaction rates and atmospheric boundary conditions are also provided. We describe the MESA Software Development Kit that packages all the required components needed to form a unified maintained and wellvalidated build environment for MESA. We also highlight a few tools developed by the community for rapid visualization of MESA star results.
['asteroseismology', 'methods: numerical', 'planets and satellites: physical evolution', 'stars: evolution', 'stars: massive', 'stars: rotation', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Modules for Experiments in Stellar Astrophysics MESA Planets Oscillations Rotation and Massive Stars
{'MODULES FOR EXPERIMENTS IN STELLAR ASTROPHYSICS (MESA): PLANETS, OSCILLATIONS, ROTATION, AND MASSIVE STARS': 'B ill P axton 1 , M atteo C antiello 1 , P hil A rras 2 , L ars B ildsten 1,3 , E dward F. B rown 4 , A aron D otter 5 , C hristopher M ankovich 3 , M. H. M ontgomery 6 , D ennis S tello 7 , F. X. T immes 8 , and R ichard T ownsend 9  \nto appear in The Astrophysical Journal Supplement Series', 'ABSTRACT': "We substantially update the capabilities of the open source software package Modules for Experiments in Stellar Astrophysics ( MESA ), and its one-dimensional stellar evolution module, MESA star . Improvements in MESA star 's ability to model the evolution of giant planets now extends its applicability down to masses as low as one-tenth that of Jupiter. The dramatic improvement in asteroseismology enabled by the space-based Kepler and CoRoT missions motivates our full coupling of the ADIPLS adiabatic pulsation code with MESA star . This also motivates a numerical recasting of the Ledoux criterion that is more easily implemented when many nuclei are present at non-negligible abundances. This impacts the way in which MESA star calculates semi-convective and thermohaline mixing. We exhibit the evolution of 3-8 M GLYPH<12> stars through the end of core He burning, the onset of He thermal pulses, and arrival on the white dwarf cooling sequence. We implement di GLYPH<11> usion of angular momentum and chemical abundances that enable calculations of rotating-star models, which we compare thoroughly with earlier work. We introduce a new treatment of radiation-dominated envelopes that allows the uninterrupted evolution of massive stars to core collapse. This enables the generation of new sets of supernovae, long gamma-ray burst, and pair-instability progenitor models. We substantially modify the way in which MESA star solves the fully coupled stellar structure and composition equations, and we show how this has improved the scaling of MESA 's calculational speed on multi-core processors. Updates to the modules for equation of state, opacity, nuclear reaction rates, and atmospheric boundary conditions are also provided. We describe the MESA Software Development Kit ( SDK ) that packages all the required components needed to form a unified, maintained, and well-validated build environment for MESA . We also highlight a few tools developed by the community for rapid visualization of MESA star results.  \nSubject headings: asteroseismology - methods: numerical - planets and satellites: physical evolution stars: evolution - stars: massive - stars: rotation", '1. INTRODUCTION': "As the most commonly observed objects, stars remain at the forefront of astrophysical research. Advances in optical detector technology, computer processing power, and data storage capability have enabled new sky surveys (e.g., the Sloan Digital Sky Survey; York et al. 2000); triggered many new optical transient surveys, such as the Palomar Transient Factory (Law et al. 2009) and Pan-STARRS1 (Kaiser et al. 2010); and allowed for space missions (e.g., Kepler ; Koch et al. 2010) that continuously monitor more than 100,000 stars. The stellar discoveries from these surveys include revelations about rare stars, unusual explosive outcomes, and remarkably complex binaries. The immediate future holds tremendous promise, as both the space-based survey Gaia (de Bruijne 2012; Liu et al. 2012) and the ground based Large Synoptic Survey Telescope (LSST; Ivezic et al. 2008) come to fruition.  \nThese developments have created a new demand for a reliable and publicly available research and education tool in computational stellar astrophysics. We introduced the open source community tool MESA (Paxton et al. 2010, hereafter Paper I) to meet these new demands. This first 'instrument' paper described the design, implementation, and realm of validity of MESA modules for numerics, microphysics, and macrophysics, and introduced the stellar evolution module, MESA star . We presented a multitude of tests and code comparisons that served as our initial verification and demon-  \nstrated MESA star 's initial capabilities. Since Paper I, MESA has attracted over 500 registered users, witnessed over 5,000 downloads from http://mesa.sourceforge.net/ , started an annual Summer School program, and provided a portal ( http://mesastar.org ) for the community to openly share knowledge (e.g., the specific settings for a published MESA star run), codes, and publications.  \nThis paper describes the major new MESA capabilities for modeling giant planets, asteroseismology, and the treatment of rotation and evolution of massive stars. We also describe numerous advances since Paper I. These include the incorporation of composition gradients in the determination of convective mixing and additional verification for evolution of intermediate mass stars and the white dwarfs they create.  \nOur improvements to MESA star for gas giant planets were motivated by the dramatic growth in this field. Over 800 exoplanets have been confirmed, and their study has prompted enormous progress in our understanding of the formation and migration of giant planets, and of the importance of factors such as stellar mass (Laughlin et al. 2004; Alibert et al. 2011; Boss 2011), composition (Fischer & Valenti 2005; Young et al. 2012), and binarity (Patience et al. 2002; Mugrauer & Neuhauser 2009; Roell et al. 2012). Puzzles remain, though, both in our solar system and in the studies of the plethora of these newly discovered exoplanets, including the characteristics of the planet-hosting stars and the interiors, atmospheres, surface gravities, temperatures, and compositions of the planets (e.g., Udry & Santos 2007; Seager & Deming 2010). Many of these variations can now be numerically explored, as can the incorporation of an inert core in an otherwise regular gas giant and the impact of irradiation.  \nThe ability to infer stellar properties (e.g., mass, radius, internal state, and rotation) from measurements of the radial and non-radial oscillation modes has been dramatically improved by two space-based optical telescopes (Convection Rotation and Planetary Transits, CoRoT ; Baglin et al. 2009 and Kepler ; Borucki et al. 2009). The high cadences and precision (often better than ten parts per million) reveal and accurately measure multitudes of oscillation frequencies for over 10,000 stars, substantially raising the need for accurate and e GLYPH<14> cient computations of stellar mode frequencies and the resulting eigenfunctions. The intrinsic flexibility of MESA star allows for the exploration of model-space required to precisely infer stellar properties from the observed frequencies.  \nAn important new addition to MESA is the incorporation of stellar rotation and magnetic fields in radiative regions. As stars are not solid bodies, they undergo radial di GLYPH<11> erential rotation (Thompson et al. 2003; Balbus et al. 2012) and also rotate at di GLYPH<11> erent angular velocities at di GLYPH<11> erent latitudes (Ruediger et al. 1998; Bonanno et al. 2007; Kuker et al. 2011). These rotational shears have a significant impact on the evolution of the stellar magnetic field. Despite the resulting 3D nature of magnetism and rotation, the stellar evolution community has come a long way in understanding stars with 1D simulations (Meynet & Maeder 1997; Langer et al. 1999; Maeder & Meynet 2000b; Heger & Langer 2000; Hirschi et al. 2004; Cantiello & Langer 2010), thus motivating our need to fully incorporate rotation within MESA . The new flexibility in angular momentum transport mechanisms allows for numerical exploration of alternate rotational outcomes should the observations (e.g., asteroseismology) require it.  \nThe paper is outlined as follows. Section 2 describes the new capability of MESA to evolve models of giant planets, while x 3 discusses the new asteroseismology capabilities.  \nThe MESA implementation of composition gradients in stellar interiors and their impact on convective mixing is described in x 4. The status of the evolution of intermediate mass stars and the MESA star construction and evolution of white dwarfs is described in x 5. The new capabilities for evolving rotating stars is described in x 6. The onset of near Eddington luminosities and radiation pressure dominance in the envelopes of evolving massive stars has been a challenge for many stellar evolution codes ever since the realization of the iron opacity bump at log T GLYPH<25> 5 : 3 (Iglesias et al. 1992). We discuss in x 7 the resulting improvements for evolving massive stars. This allows for the uninterrupted evolution of rotating massive stars to the onset of core collapse. We conclude in x 8 by highlighting where additional improvements to MESA are likely to occur in the near future. Appendix A describes the many improvements to the physics modules since Paper I; Appendix B presents 'nuts and bolts' information on the primary components of evolution calculations; and Appendix C presents the MESA Software Development Kit ( SDK ). All of our symbols are defined in Table 1. We denote components of MESA , such as modules and routines, in Courier font, e.g., evolve star .  \nTable 1 Variable Index.Table 1 -Continued", '2. GIANT PLANETS AND LOW-MASS STARS': "Evolutionary models of giant planets and low-mass stars di GLYPH<11> er from their higher-mass stellar counterparts in both the microphysics needed to describe the interior and the role of stellar irradiation in the outer boundary condition. For masses M . 84 MJ, hydrogen burning is insu GLYPH<14> cient to prevent cooling and contraction. Deuterium burning can briefly slow the cooling for M & 13 MJ, where MJ = 9 : 54 GLYPH<2> 10 GLYPH<0> 4 M GLYPH<12> is Jupiter's mass, but has a negligible influence on the cooling for smaller masses. Hence nuclear burning can be ignored in the planetary mass regime.  \nFor hydrogen-helium rich objects with M GLYPH<29> MJ, an ideal gas equation of state (EOS), with arbitrary degeneracy, is a good approximation while for M . MJ particle interactions play an important role. Specifically, pressure ionization of hydrogen at GLYPH<26> ' 1 g cm GLYPH<0> 3 and T ' 10 4 K causes a sudden change from a H2-dominated phase to an ionized phase. MESA employs the Saumon et al. (1995) equation of state (SCVH EOS), smoothly interpolated from the low to high pressure phase, for this complicated region of parameter space where thermal, Fermi, and electrostatic energies may all be comparable. The SCVH EOS includes pressure ionization of hydro-  \ngen, but not helium. The temperature range covered by the tables is 2 : 10 < log T (K) < 7 : 06, and the pressure ranges from log P (dyne cm GLYPH<0> 2 ) = 4 to a maximum value 19 dependent on the temperature. Smooth interpolation to other EOS occurs near the SCVH boundaries (for more details see Paper I). At the low temperatures in planetary atmospheres, abundant species such as CNO atoms will be in molecular form, and may condense into clouds. MESA does not follow the transition from atomic to molecular form for these species in the EOS-they are currently included by increasing the helium abundance from Y to Y + Z when calling the SCVH EOS. MESA does, however, include the e GLYPH<11> ect of molecules in the Rosseland opacities. Currently, the Ferguson et al. (2005) and Freedman et al. (2008) tables, which include the opacity from molecules, but ignore condensates, are available.  \nLastly, for planets in close-in orbits about their parent star, the external irradiation flux may be orders of magnitude larger than the cooling flux from the planet's interior. This may dramatically increase the surface temperature and a GLYPH<11> ect the outer boundary condition. MESA now implements several options for this surface heating, including the flexibility to include user-supplied prescriptions.  \nIn the following subsections, we discuss a new MESA module that creates initial models in the planetary mass range M ' 0 : 1-10 MJ, and present a suite of evolutionary calculations. We discuss how surface irradiation may be included, as well as an inert core at the center of the planet. We also show what MESA star yields for the mass-radius relation for sub-solar mass stars in x 2.5.", '2.1. Construction of Starting Models': "For stellar mass objects, the pre ms model routine constructs pre-main-sequence (PMS) models assuming L ( r ) / m , where L ( r ) is the luminosity at radius r , by iterating on the starting conditions at the center to find a model with a given M and central temperature T c. This PMS routine works well for M & 0 : 03 M GLYPH<12> , but lower masses may not converge when the guess for central density GLYPH<26> c and luminosity are not close enough to the (unknown) true values. As a result, it is di GLYPH<14> -cult and time consuming to create models with M < 0 : 03 M GLYPH<12> using the same routine for giant planets as for stars.  \nA new routine called create initial model builds a model of given M and radius R using an adiabatic temperature profile. Given the central pressure P c and specific entropy S , the equation of hydrostatic balance is integrated outward, and the temperature at each step determined from the equation of state using T = T ( P ; S ). The values of P c and S are iterated to attain the desired M and R . The luminosity profile is then derived treating S as constant in space for the fully convective planet (e.g., Ushomirsky et al. 1998), so  \nZ m 0 d m 0 T ( m 0 ) d S d t ' d S d t Z m 0 d m 0 T ( m 0 ) = GLYPH<0> L ( m ) : (1)  \nThe luminosity at the surface, L ( M ), is estimated using the radius R and temperature T e GLYPH<11> at the GLYPH<28> = GLYPH<20> P = g = 2 = 3 point as L ( M ) = 4 GLYPH<25> R 2 GLYPH<27> SB T 4 e GLYPH<11> . Given L ( M ), the luminosity at interior points is found by  \nL ( m ) = L ( M ) 0 B B B B B B B @ R m 0 d m 0 T ( m 0 ) R M 0 d m 0 T ( m 0 ) 1 C C C C C C C A : (2)  \nThis procedure works well for M down to GLYPH<24> 0 : 1 MJ and over a range of initial radii, allowing the user to choose either a  \nGLYPH<24> 1 RJ radius appropriate for a cold planet, to radii GLYPH<24> 2-3 RJ appropriate for young or inflated planets (e.g., Marley et al. 2007). Here RJ = 7 : 192 GLYPH<2> 10 9 cm is the equatorial radius of Jupiter.", '2.2. Evolutionary Calculations': "Figure 1. The solid black lines show T c versus GLYPH<26> c during the evolution. Each line is labeled on the left by the mass in units of MJ. The dotted red lines show constant values of log(age[yr]), labeled at the base of each line. The blue dashed lines show fixed values of S = ( N A k B), labeled at the top of each line. The large black dots show the position of maximum T c along the evolutionary track.  \n<!-- image -->  \nFigures 1 and 2 show evolutionary calculations for models with masses M = 0 : 2-20 MJ. All models were evolved for 20 Gyr. The initial models from create initial model had a large radius R = 5 RJ. The other parameters used are Y = 0 : 27, Z = 0 : 02 and GLYPH<11> MLT = 2. The opacity and EOS tables used are eos file prefix = mesa , kappa file prefix = gs98 and kappa lowT prefix = lowT Freedman11 . The atmosphere model is which atm option = simple photosphere .  \nFigure 1 is a low mass extension of Figure 16 from Paper I, showing evolution in the GLYPH<26> cT c plane. Each track (solid black curve) is labeled on the left by the planet's mass, and evolution goes from left to right. Initially the planet is nondegenerate and contraction increases both GLYPH<26> c / R GLYPH<0> 3 and T c / R GLYPH<0> 1 / GLYPH<26> 1 = 3 c . A maximum T c is reached when k B T c GLYPH<24> E F ; c, where E F ; c is the electron Fermi energy at the center, beyond which GLYPH<26> c approaches a constant as T c decreases further. Ignoring Coulomb interactions in the EOS, S is a function of the electron degeneracy parameter GLYPH<22> e = k B T , where GLYPH<22> e is the electron chemical potential and all models should have maximum k B T c GLYPH<24> E F ; c at the same S . The line labeled S = ( N A k B) = 10 : 3 indeed coincides with maximum T c down to M ' 1 MJ, but at smaller masses where non-ideal e GLYPH<11> ects are more important, maximum T c occurs when S = ( N A k B) < 10 : 3. Also shown in Figure 1 are lines of constant age, shown as dotted red lines, and labeled on the bottom of the plot.  \nThe same evolutionary calculations are used in Figure 2 to  \nFigure 2. Radius versus mass iso-contours from a suite of evolutionary calculations. The solid red lines show R = RJ versus M = MJ at fixed values of log(age[yr]), labeled on the left of each curve. The dashed blue curves are for fixed entropy, with each curve labeled by S = ( N A k B) on the right. The dotted black curves are for fixed luminosity, with each curve labeled by log( L [erg s GLYPH<0> 1 ]) above M = 1 MJ. The green curve at the bottom is the T = 0 M -R relation from Zapolsky & Salpeter (1969) for a solar mixture of H and He.  \n<!-- image -->  \nshow radius versus mass at fixed values of age, entropy or luminosity. At late times, or low entropy and luminosity, the radius approaches the zero-temperature value (green curve; Zapolsky & Salpeter 1969) for which thermal support is insignificant. The maximum radius occurs where gravitational and Coulomb energies, per ion, are comparable. The solid red lines, labeled by age on the left, show that contraction down to R ' 1 : 5 RJ is rapid, taking less than 10 Myr for M . 10 MJ. This initial rapid cooling phase occurs because the initial luminosity is orders of magnitude higher than the luminosity around one Gyr. This can been seen in the black dotted contours of constant log( L ), where L is larger by a factor of 100 for R = 1 : 3 RJ and 10 4 for R = 1 : 7 RJ, as compared to R = 1 : 1 RJ. The blue dashed lines show contours of constant entropy, labeled on the right by S = ( N A k B).", '2.3. Implementation of Inert Cores': "In the core accretion model of planet formation (e.g., Pollack et al. 1996; Hubickyj et al. 2005), a rock / ice core is first assembled. Once this core grows to GLYPH<24> 10 M GLYPH<8> , where M GLYPH<8> denotes an Earth mass, it can initiate rapid accretion of nebular gas, which could then dominate the mass of the planet. For studies of planetary radii, a central core composed of high mean molecular weight material can decrease the radius of the planet by a significant amount ( ' 0 : 1-0 : 2 RJ). The MESA star inert core feature allows one to add a core of specified mass M c and radius R c, or more conveniently, density GLYPH<26> c. A luminosity L c may also be specified, although the high mean molecular weight of the core, as compared to the overlying H / He envelope, implies that even large cores will tend to have small heat content (Fortney et al. 2006). This inert core is not presently evolved in any way, and changes in P c during evolution are neglected as R changes. While cores of mass  \n. 10-20 M GLYPH<8> are commonly used for modeling solar system giants (e.g. Guillot 2005), the large masses and small radii of some exoplanets may imply far larger core masses (e.g. HD 149026; Sato et al. 2005). In addition, Neptune-like planets with smaller ratios of envelope to core masses may be modeled with MESA (Owen & Wu 2013).", '2.4. Irradiation': "Figure 3. Radius versus age for the planet HD 209458b. The solid red lines are for MESA , using the grey irradiated atmosphere model. The dotted black lines show the CEPAM code results. The dashed blue lines show the MESA calculation using the F ? -GLYPH<6> ? surface heat source. The data point with error bars is the observed value of the radius for HD 209458b quoted in Guillot (2010). The two sets of curves are deep heating (upper three curves) and shallow heating (lower three curves).  \n<!-- image -->  \nSurface heating by stellar irradiation changes the boundary condition for the planet's cooling and contraction. This modifies the planetary radius versus age for exoplanets at orbital separation . 0 : 1 AU. MESA provides several ways to implement surface heating with varying degrees of fidelity to the true solution. These presently include:  \na) An energy generation rate GLYPH<15> = F ?= 4 GLYPH<6> ? applied in the outer mass column GLYPH<6> GLYPH<20> GLYPH<6> ? . Here F ? is the day-side flux from the star, and GLYPH<6> ( r ) = R R r d r 0 GLYPH<26> ( r 0 ) is the mass column. In steady-state, this generates an outward flux F ?= 4, which is meant to simulate the angle-averaged flux over the planetary surface. This model implicitly assumes that day-night heat transport is e GLYPH<14> cient, and at the depths of interest the temperature is uniform over the surface. The parameters F ? and GLYPH<6> ? are specified through the user-specified variables irradiation flux and column depth for irradiation , making this the simplest method to use. This heating mechanism represents absorption of stellar optical radiation well below the photosphere of the planet's thermal radiation and gives rise to greenhouse heating of the atmosphere where GLYPH<15> , 0.  \nb) MESA 's grey irradiated atmosphere model (see also x A.5) implements the angle-averaged temperature profile of Guillot (2010). This approximate solution to the transfer equation assumes two frequency bands: optical radiation from the star (with user-specified opacity GLYPH<20> v) and thermal radiation originating in the planet (with user-specified opacity GLYPH<20> th). The temperature profile is derived using the Eddington approximation, assuming an external flux from the star as well as a flux from the planetary interior. While the Guillot (2010) model implemented in MESA uses a single temperature as a function of depth, it is derived allowing for local temperature variations over the surface which are then averaged over angle. This temperature profile is shown to be valid in the presence of horizontal heat transport by fluid motions. This is the only MESA atmosphere model that uses pressure instead of optical depth to determine the surface boundary condition. As this pressure may be relatively deep in the atmosphere, a correction to the radius may be required to give either the vertical thermal photosphere, or the optical photosphere in transit along a chord. Lastly, the relax irradiation routine improves initial convergence by providing a starting model closer to the irradiated one.  \nc) Finally, MESA allows user-specified heating functions (e.g., F ? -GLYPH<6> ? surface heating) or atmosphere models (e.g., grey irradiated ). User-supplied routines may be easily implemented by using the other energy module.  \nFigure 3 shows radius versus age for the planet HD 209458b (Guillot 2010). The two groupings of lines are for di GLYPH<11> erent heating depths, and within each grouping of lines, there are three calculations: MESA using grey irradiated surface boundary condition (solid red line), MESA using the F ? -GLYPH<6> ? surface heating profile (dashed blue line), and CEPAM (Guillot & Morel 1995) using the same grey irradiated boundary condition (dotted black line; kindly provided by Tristan Guillot). The lower curves , corresponding to shallow heating, use fiducial values ( GLYPH<20> th ; GLYPH<20> v) = (10 GLYPH<0> 2 ; 6 GLYPH<2> 10 GLYPH<0> 3 ) cm 2 g GLYPH<0> 1 and give a model radius significantly smaller than the observed radius. The upper curves , corresponding to deep heating, use ( GLYPH<20> th ; GLYPH<20> v) = (10 GLYPH<0> 2 ; 6 GLYPH<2> 10 GLYPH<0> 4 ) cm 2 g GLYPH<0> 1 , yielding significantly hotter temperatures deep in the surface radiative zone, which slow the cooling enough to agree with the observed radius. The choice GLYPH<6> ? = 2 =GLYPH<20> v gives agreement between the grey irradiated and F ? -GLYPH<6> ? methods, where the factor of 2 accounts for the fact that the grey irradiated boundary condition has some heating below GLYPH<6> = 1 =GLYPH<20> v. The radii are at the GLYPH<28> th = 2 = 3 photosphere for a vertical path into the atmosphere.  \nThe agreement between all three methods is excellent, at the 1-2% level after 100 Myr. The remaining discrepancy between the MESA and CEPAM grey irradiated results are likely due to di GLYPH<11> erent opacity tables, with the MESA result using an update of Freedman et al. (2008) (Freedman 2011, priv. comm.) while the CEPAM run uses the Allard et al. (2001) COND table. The di GLYPH<11> erences at ages . 100 Myr are due to di GLYPH<11> erent starting conditions. The CEPAM calculation started with initial radius 2 RJ, whereas the MESA calculations started with 5 RJ. The MESA grey irradiated and F ? -GLYPH<6> ? calculations di GLYPH<11> er at . 100 Myr, likely because the former has a fixed thermal opacity while the latter allows the opacity to change.", '2.5. Low-Mass Main Sequence Stars': "Most of MESA star 's capability to evolve low-mass ( M < 2 M GLYPH<12> ) stars was demonstrated in Section 7.1 of Paper I. MESA has seen use in the asteroseismology of helium core flashing stars (Bildsten et al. 2012) and the discovery of a new instability from the onset of 3 He burning (van Saders & Pinsonneault 2012). We expect the future use of MESA star for asteroseismic investigations of these stars to be substantial (see x 3).  \nFigure 4. Stellar isochrones at solar composition spanning 0.1 to 1 M GLYPH<12> from MESA star (solid lines) and Dotter et al. (2008, dashed lines) in the mass-radius plane. The data points plotted are the same as shown by Bass et al. (2012).  \n<!-- image -->  \nglyph[circledot]  \nThe derivation of accurate planetary radii based on transits requires accurate radii of the host stars; this motivates MESA star investigations of low-mass stars (Lloyd 2011). Figure 4 shows 1 and 5 Gyr isochrones at solar composition ( Y = 0 : 27 ; Z = 0 : 019) from MESA star (solid lines) and Dotter et al. (2008, dashed lines) in the mass-radius diagram. Data points shown in Figure 4 are taken from Torres et al. (2010), Carter et al. (2011), Irwin et al. (2011), and Bass et al. (2012). This figure is a reproduction of the upper panel of Figure 11 from Bass et al. (2012). Figure 4 indicates that MESA star is capable of producing mass-radius relations for main sequence stars that are consistent with other widely-used models as well as observational data. The MESA star models were computed using, as much as possible, the same physical assumptions as the models used by Dotter et al. (2008). The main di GLYPH<11> erence is the equation of state, for which Dotter et al. (2008) used FreeEOS 10 and MESA star uses a combination of the OPAL (Rogers & Nayfonov 2002) and SCVH EOS for thermodynamic parameters relevant to this diagram.", '3. ASTEROSEISMOLOGY': 'With its highly configurable output options, and its ability to calculate asteroseismic variables, MESA star can readily produce models suitable for use with a range of oscillation codes. In addition to its own text output files, MESA can  \nproduce outputs in formats widely used by stellar oscillation codes, such as fgong and osc (Monteiro 2009).  \nIn Figure 5 we show the evolution of a 1 M GLYPH<12> model in the Hertzsprung-Russell Diagram (HRD) and in T cGLYPH<26> c space. These were evolved following the test case found in 1M pre ms to wd , which was modified to include di GLYPH<11> usion. This runs without user intervention from pre-main sequence to white dwarf. To demonstrate the changing stellar structure as the model evolves from the main sequence to post helium-core burning on the Asymptotic Giant Branch (AGB), we show in Figure 6 some of the fundamental quantities extracted from the corresponding profile.data files for the models marked in Figure 5. These include the Lamb and Brunt-Vaisala frequencies defined respectively by  \nS 2 \' = \' ( \' + 1) c 2 s r 2 ; (3)  \nN 2 = g r " 1 GLYPH<0> 1 d ln P d ln r GLYPH<0> d ln GLYPH<26> d ln r # ; (4)  \nwhere c s is the adiabatic sound speed and \' is the spherical harmonic degree.', '3.1. The Solar Sound Speed Profile': "The seismic properties of the Sun provide a test of stellar evolution models, and an opportunity to calibrate GLYPH<11> MLT for any particular set of input physics and other assumptions. The MESA star test case solar calibration produces a calibrated Standard Solar Model. Figure 7 shows the di GLYPH<11> erence between the helioseismically-inferred solar sound speed profile and this model. We also show 'Model S' from Christensen-Dalsgaard et al. (1996). Both models employ comparable input physics and assume solar abundances from Grevesse & Noels (1993) and Grevesse & Sauval (1998). One clear improvement since Paper I is a smoother sound speed profile at small r = R , which is primarily due to improvements in the diffusion module. This is particularly important for asteroseismology, where sharp features in the sound speed profile can influence the stellar oscillation frequencies. The results are based on the solar calibration test case compiled with the GNU Fortran compiler version 4.7.2 on Mac OS X 10.7.5; Appendix B.11 provides information about how the solar calibration results may depend on di GLYPH<11> erent operating systems and compilers.", '3.2. New Asteroseismic Capabilities in MESA': "The 'astero' extension to MESA star implements an integrated approach that passes results automatically between MESA star and the new MESA module based on the adiabatic code ADIPLS (Christensen-Dalsgaard 2008a, June 2011 release). The MESA module ADIPLS also supports independent use for post-processing, including the calculation of pulsation frequencies.  \nThis astero extension enables calculation of selected pulsation frequencies by MESA star during the evolution of the model. This allows fitting to the observations that can include spectroscopic constraints (e.g., [Fe / H] and T e GLYPH<11> ), asteroseismic constraints, such as the large frequency separation ( GLYPH<1> GLYPH<23> ) and the frequency of maximum power ( GLYPH<23> max), and even individual frequencies. A variety of approaches for finding a best-fitting model are available, including grid searches and automatic GLYPH<31> 2 minimization by the Hooke-Jeeves algorithm (Hooke & Jeeves 1961) or by the 'Bound Optimization BY Quadratic  \nApproximation' (BOBYQA; Powell 2009) technique. These searches are user controlled through a number of parameter bounds and step sizes. Users also have full control over the relative weight assigned to the seismic and spectroscopic parts of the GLYPH<31> 2 statistic.  \nFor the automated GLYPH<31> 2 minimization, astero will evolve a pre-main sequence model from a user defined starting point, and find the best match along that single evolutionary track. The code then recalculates the track, again initiated at the premain sequence, with di GLYPH<11> erent initial parameters such as mass, composition, mixing length parameter and overshoot, and repeats until the lowest GLYPH<31> 2 has been found.  \nCalculating specific mode frequencies is computationally intensive. Hence, a number of options exist to improve the e GLYPH<14> ciency of the minimization when individual frequencies are included. Bounds can be established on stellar parameters (e.g., T e GLYPH<11> , central H mass fraction, GLYPH<1> GLYPH<23> ), so that ADIPLS is invoked only when the model falls within these bounds. This enables certain evolutionary stages to be skipped when other observational diagnostics rule them out-if a star is known to be a red giant, for instance, there is no sense in invoking ADIPLS when models are on the main sequence. The large frequency separation, GLYPH<1> GLYPH<23> , of the model is calculated as the inverse of the sound travel time through the star, GLYPH<1> GLYPH<23> = [2 R d r = c s] GLYPH<0> 1 (Tassoul 1980; Gough 1986). There is also the option to derive GLYPH<1> GLYPH<23> using simple solar scaling: GLYPH<1> GLYPH<23> / ( M = R 3 ) 0 : 5 (Kjeldsen & Bedding 1995). To obtain GLYPH<23> max, MESA scales the solar value with the acoustic cut-o GLYPH<11> frequency: GLYPH<23> max / g = p T e GLYPH<11> (Brown et al. 1991; Kjeldsen & Bedding 1995).  \nMoreover, hierarchical approaches to the frequency fitting can be selected, saving large amounts of computational time. In one case the radial modes are first calculated, and only when they match reasonably well are the non-radial mode frequencies derived and included in the GLYPH<31> 2 . This is particularly beneficial for red giants where the calculation of the non-radial frequencies is extremely time consuming. Another example is when the time steps in the stellar evolution calculations are too large to find an accurate minimum of GLYPH<31> 2 . Hence, as a further option to increase e GLYPH<14> ciency while attaining accuracy, the time steps can be set to automatically reduce when the model comes close to the 'target box' of the observational constraints. As for other modules used in MESA star , astero o GLYPH<11> ers a range of graphical outputs including an 'echelle diagram where the fitting process can be followed in real time.  \nThere is also an option for including corrections to the model frequencies on-the-fly to compensate for the inadequate modelling of the near surface layers of the star. The e GLYPH<11> ect, known as the 'surface term,' is seen as a frequency dependent o GLYPH<11> set between the modelled and observed acoustic frequencies of the Sun (e.g. Christensen-Dalsgaard & Thompson 1997). The o GLYPH<11> set increases towards higher frequencies and is well described by a power law (Kjeldsen et al. 2008). MESA star follows the approach described by Kjeldsen et al. (2008) for correcting the surface term.  \nTo illustrate the performance of astero, we show here a fit to the star HD49385. The input frequencies and the spectroscopic constraints are from Deheuvels et al. (2010). We first ran a wide-range grid search over M , GLYPH<11> MLT, [Fe / H], and Y , including only [Fe / H], T e GLYPH<11> , and GLYPH<1> GLYPH<23> as observational constraints. The results of this initial search guided our starting parameters and ranges for the next automatic GLYPH<31> 2 minimization. We first compare our grid results with those of the RADIUS grid  \nFigure 5. Hertzsprung-Russell diagram and T cGLYPH<26> c evolution of a 1 M GLYPH<12> model evolving from pre-main sequence to the white dwarf cooling sequence. The number labels denote selected models, for which we show internal profiles in Fig. 6.  \n<!-- image -->  \nsearch routine (Stello et al. 2009), which is based on a grid of ASTEC models (Christensen-Dalsgaard 2008b) and find agreement within uncertainties.  \nWe then include the individual oscillation frequencies and use the Hooke-Jeeves algorithm for the GLYPH<31> 2 minimization. Model frequencies were corrected for the surface term, and the part of the GLYPH<31> 2 coming from the frequencies was given 2 / 3 of the weight in the final GLYPH<31> 2 , similar to that used by Metcalfe et al. (2012). To ensure we adequately sample the parameter space, we initiate the search at several initial values within a broad range. By starting the search from multiple initial values, we aim to reduce the chance of ending up in a local minimum, which could potentially provide unphysical results, such as the spuriously low helium abundances reported by Mathur et al. (2012). Current developments in astero further seeks to overcome such problems and improve the robustness of the results by including frequency ratios (Roxburgh & Vorontsov 2003; Silva Aguirre et al. 2013) in the GLYPH<31> 2 minimization.  \nEach 'Hooke' search generates several stellar evolution tracks, each with a best GLYPH<31> 2 value. We then combine the data from about 1400 tracks to estimate the 1GLYPH<27> uncertainties in the varied parameters following the approach by Deheuvels et al. (2010). The lowest (reduced) GLYPH<31> 2 value we obtained was 2.4 with a few tens of models in the 2.4-4.0 range, which all fit the frequencies similarly well. Among these models there are two families of results, one of which has slightly lower [Fe / H] and Y , and a slightly increased value for the spectroscopic part of the GLYPH<31> 2 .  \nThe comparison of the observed and modeled frequencies for the realization with the lowest GLYPH<31> 2 is shown in the 'echelle diagram format in Figure 8. A plot of the internal structure including the Brunt-Vaisala and Lamb frequencies is shown in Figure 9, and the parameters of the model are listed in Table 2. We set f ov = 0 : 015 and use the GN98 solar abundances. Our results can be best compared to those listed as 'low GLYPH<11> ov' and 'GN93' in Table 4 of Deheuvels & Michel (2011) and agree within the uncertainties.", '3.3. The E GLYPH<11> ect of Composition Gradients on the Brunt-Vaisala Frequency': "Including the e GLYPH<11> ect of composition gradients in the calculation of the Brunt-Vaisala frequency is important for two rea-a [  \nTable 2 Properties of best fitting model to HD49385  \n= H] s is the log of the ratio of the surface ( Z / X ) relative to the solar value of 0.02293.  \nsons. First, it is necessary for implementing the Ledoux criterion for convection, which is used to determine the chemical mixing and convective heat transport in a region (see x 4.1). Second, a smooth and accurate method for calculating N 2 is crucial for studies of g-mode pulsation in stars. In a highly degenerate environment, the pressure is nearly independent of temperature, and P / GLYPH<26> GLYPH<0> 1 , so from eq. (4) we see that N 2 depends on the di GLYPH<11> erence of two large and nearly equal quantities. This can lead to a loss of precision and a noisy N 2 . To eliminate this problem, N 2 is re-written into a form that depends on the di GLYPH<11> erence of the adiabatic and true temperature gradients and on the composition gradient:  \nN 2 = g 2 GLYPH<26> P GLYPH<31> T GLYPH<31>GLYPH<26> ( r ad GLYPH<0> r T + B ) : (5)  \nThe term B explicitly takes into account the e GLYPH<11> ect of composition gradients and is commonly called the Ledoux term (e.g., Unno et al. 1989; Brassard et al. 1991). For the general case of an N -component plasma with mass fractions f Xi g , the  \nFigure 6. Internal structure of the five points (indicated by the numbers in each panel) marked in Figure 5. The left panel for each point shows N and S ' for harmonic degrees ' = 1 and 2. The dashed line indicates the frequency of maximum power GLYPH<23> max of the stochastically excited solar-like modes. The vertical dotted lines mark the radius of the model. Right panels show temperature, hydrogen and helium mass fractions, mass, luminosity, and the nuclear energy generation rate. Grey areas mark convective regions according to the Schwarzschild criterion.  \n<!-- image -->  \n=  \n<!-- image -->  \n=  \nGLYPH<12>  \n1  \nFigure 7. Comparison of the di GLYPH<11> erence between the helioseismicallyinferred sound speed profile (Bahcall et al. 1998) of a MESA star model and Model S (Christensen-Dalsgaard et al. 1996).  \n<!-- image -->  \nFigure 8. ' Echelle diagram of the oscillation frequencies of the subgiant HD49385. Observed frequencies are shown with filled symbols as blue squares ( ' = 2), black circles ( ' = 0), and red triangles ( ' = 1), and the matched model frequencies are shown with open symbols. Black horizontal lines indicates 1GLYPH<27> error bars.  \n<!-- image -->  \nFigure 9. Same format as in Figure 6, but for the best-fitting model of HD49385 (see also Table 2).  \n<!-- image -->  \n=  \nstandard formula for B is (e.g., Unno et al. 1989)  \nB = GLYPH<0> 1 GLYPH<31> T N GLYPH<0> 1 X i = 1 @ ln P @ ln Xi ! GLYPH<26>; T ; f Xj , i g d ln Xi d ln P : (6)  \nSince P N i = 1 Xi = 1, one of the mass fractions can be eliminated, so that the sum in eq. (6) runs from 1 to N GLYPH<0> 1. We note that the partial derivatives in eq. (6) hold all the f Xj g constant except for Xi and XN , where XN is varied so as to maintain P N i = 1 Xi = 1.  \nAlthough eq. (6) is correct as written, we have developed a new, formally-equivalent prescription that is both numerically robust and simpler to implement. We define a new Ledoux term by taking a directional derivative along the radial composition gradient in the stellar model,  \nB GLYPH<17> GLYPH<0> 1 GLYPH<31> T lim GLYPH<14> ln P ! 0 ln P ( GLYPH<26>; T ; ~ X + (d ~ X = d ln P ) GLYPH<14> ln P ) GLYPH<0> ln P ( GLYPH<26>; T ; ~ X ) GLYPH<14> ln P : (7)  \nThe implementation of the above derivative typically involves the use of quantities on neighboring mesh points. Using the subscript k to denote the value of a given quantity on the k th mesh point, we therefore have  \nB = GLYPH<0> 1 GLYPH<31> T ln P ( GLYPH<26> k ; Tk ; ~ Xk + 1) GLYPH<0> ln P ( GLYPH<26> k ; Tk ; ~ Xk ) ln Pk + 1 GLYPH<0> ln Pk : (8)  \nThis is the form of the Ledoux term that is implemented in MESA and we term it the 'New Ledoux' formulation. Since MESA ensures that P N i = 1 Xi = 1 at each mesh point, this condition does not have to be separately enforced. This formulation requires just one numerical di GLYPH<11> erence along ~ X that is consistent with the stellar model and equation of state. Because MESA 's EOS does not directly supply the partial derivatives required for the formulation in eq. 6, an implementation of that method would su GLYPH<11> er in both accuracy and e GLYPH<14> ciency from having to do a large number of numerical di GLYPH<11> erences. Brassard et al. (1991) dealt with a similar problem by using a restricted form of eq. 6 that included only the helium composition gradient. They showed that for cases where their restricted form applied, it gave significantly better numerical results than an implementation of eq. 4 based on finite differences. Figure 10 shows that our New Ledoux prescription (grey heavy curve) retains their good results compared to eq. 4 (thin black curve) while extending the applicability to cases that cannot be dealt with using only helium gradients.", '4. MIXING MECHANISMS INVOLVING COMPOSITION GRADIENTS': 'We described the implementation of mixing-length theory (MLT) in Paper I, including the allowance for overshoot beyond the boundaries of the convective zones as determined by the standard Schwarzschild condition, r rad > r ad. Overshooting is implemented via an exponential decay of the convective di GLYPH<11> usion coe GLYPH<14> cient beyond the boundary of full convection, following Herwig (2000a):  \nD ov = D conv ; 0 exp GLYPH<0> 2 GLYPH<1> r f ov GLYPH<21> P ! ; (9)  \nwhere D conv ; 0 is the di GLYPH<11> usion coe GLYPH<14> cient at the convective border, GLYPH<1> r is the distance from the start of overshoot, and GLYPH<21> P is the local pressure scale height. The user-adjusted dimensionless parameter f ov then determines the extent of the overshooting  \nFigure 10. Acomparison of the new Ledoux prescription for N 2 (eqs. 5 and 8) versus the direct numerical calculation (eq. 4). This calculation is for a 0 : 535 M GLYPH<12> white dwarf model at T e GLYPH<11> = 12 ; 300 K.  \n<!-- image -->  \nregion. MESA also allows for the adoption of a step-function overshooting model, where the mixing region extends a distance f ov GLYPH<21> P beyond the convective boundary with a constant specified di GLYPH<11> usion coe GLYPH<14> cient.  \nIn Paper I we did not implement the influence of composition gradients on mixing and the resulting di GLYPH<11> usion coefficients when instabilities are operative. The description of how MESA star calculates the Ledoux criterion is in x 3.3. In this section, we describe the implementation of mixing due to composition gradients in stellar interiors.  \nWe refer to x B.7 for a discussion of the free parameters involved in the implementation of these mixing mechanisms.', '4.1. Semiconvection': 'Semiconvection refers to mixing in regions unstable to Schwarzschild but stable to Ledoux, that is  \nr ad < r T < r L ; (10)  \nwhere r L is the sum of the adiabatic gradient and the Brunt composition gradient term (see eqs. [5] and [8]),  \nr L GLYPH<17> r ad + B : (11)  \nOnce r L is calculated, regions satisfying equation (10) undergo mixing via a time-dependent di GLYPH<11> usive process with a di GLYPH<11> usion coe GLYPH<14> cient calculated by the mlt module following Langer et al. (1983),  \nDsc = GLYPH<11> sc K 6 CP GLYPH<26> ! r T GLYPH<0> r ad r L GLYPH<0> r T ; (12)  \nwhere K = 4 acT 3 = (3 GLYPH<20>GLYPH<26> ) is the radiative conductivity, CP is the specific heat at constant pressure, and GLYPH<11> sc a dimensionless e GLYPH<14> ciency parameter. See x B.7 for a discussion of the range of values for GLYPH<11> sc.  \nWe stress that semiconvection and overshooting have distinct implementations in MESA . Both are time-dependent diffusive processes. As an example, in Figure 11 we display profiles of thermodynamic gradients and their resulting di GLYPH<11> usion coe GLYPH<14> cients during core helium burning in a semiconvective model with GLYPH<11> sc = 0 : 01 and in an exponentially overshooting model with f ov = 10 GLYPH<0> 5 .  \nFigure 11. Sample profiles of semiconvective (left) and exponentially overshooting (right) 3M GLYPH<12> models undergoing core helium burning. Top panels show the radiative, adiabatic, temperature, and Ledoux gradients that determine mixing boundaries and di GLYPH<11> usion coe GLYPH<14> cients. Bottom panels show the resulting di GLYPH<11> usion coe GLYPH<14> cients for energy and chemical transport. In either case, a thin dotted line spanning a single intermediate cell joins the convective and semiconvective / overshoot curves. This is intended merely as a guide for the eye, as di GLYPH<11> usion coe GLYPH<14> cients are defined only at the two boundaries of a cell. In particular, di GLYPH<11> usion for this intermediate cell is governed by convection at its interior boundary and semiconvection / overshoot at the exterior. The semiconvective model shown here was run with GLYPH<11> sc = 0 : 01; the exponentially overshooting model with f ov = 10 GLYPH<0> 5 . The profiles are taken at the points marked in Figure 15.  \n<!-- image -->', '4.2. Thermohaline Mixing': 'Thermohaline mixing arises in the presence of an inversion of the mean molecular weight in regions that are formally stable against convection according to the Ledoux criterion,  \nr T GLYPH<0> r ad GLYPH<20> B GLYPH<20> 0 ; (13)  \nIn MESA thermohaline mixing is treated in a di GLYPH<11> usion approximation, with a di GLYPH<11> usion coe GLYPH<14> cient motivated by the linear stability analysis of Ulrich (1972) and Kippenhahn et al. (1980)  \nD th = GLYPH<11> th 3 K 2 GLYPH<26> CP B ( r T GLYPH<0> r ad) : (14)  \nThe quantity GLYPH<11> th is a dimensionless e GLYPH<14> ciency parameter. In the linear analysis it depends on the aspect ratio of the blobs / fingers arising from the instability. In the case of salt fingers such a value is calibrated using laboratory experiments in water (e.g. Krishnamurti 2003), where the fingers have an aspect ratio of GLYPH<25> 5. In the stellar case the value of this parameter is vexatious (e.g. Charbonnel & Zahn 2007; Denissenkov &Pinsonneault 2008; Cantiello & Langer 2010; Wachlin et al. 2011), with recent 2D and 3D hydrodynamical calculations pointing toward a much reduced value of GLYPH<11> th relative to the salt fingers case (Denissenkov 2010; Traxler et al. 2011; Brown et al. 2013). Figure 12 shows a calculation including the effects of thermohaline mixing during the RGB phase of a 1M GLYPH<12> star after the luminosity bump (e.g. Charbonnel & Zahn 2007; Cantiello & Langer 2010). For this calculation a value GLYPH<11> th = 2  \nhas been adopted, but see x B.7 for a discussion of the range of options.', '4.3. Impact of Mixing on Convective Core Hydrogen and Helium Burning': 'The duration of the hydrogen and helium core burning depends on the extent of the convective core, so we focus here on exhibiting the MESA capabilities during these phases. As we noted above, there are many physical e GLYPH<11> ects that change the size of the convective core, such as semiconvection, overshooting, and rotation-induced mixing. For example, the Schwarzschild criterion implies larger cores than the Ledoux criterion, but when using Ledoux alone, the region above the convective boundary is overstable and so semiconvection occurs (see x 4.1).  \nWe evolved a non-rotating 1 : 5 M GLYPH<12> star with ( Y ; Z ) = (0 : 23 ; 0 : 02) through central hydrogen burning using Ledoux, Ledoux plus semiconvection, Schwarzschild, and Schwarzschild plus overshoot. As is evident in Figure 13, this set of physical processes leads to a large range of convective core masses and thereby main sequence lifetimes. For the parameters explored we found that overshooting increases the lifetime by a factor . 1.2 for Schwarzschild and . 2.5 for Ledoux. Figure 14 shows an HR diagram for each of the 1 : 5 M GLYPH<12> models undergoing core hydrogen burning, showing the impact of convective core extent on main-sequence turno GLYPH<11> morphology.  \nWe also evolved a non-rotating 3 M GLYPH<12> star with ( Y ; Z ) = (0 : 25 ; 0 : 02) through central helium burning. Overshooting extends the burning lifetime by a factor . 1.6 for Schwarzschild and . 2.8 for Ledoux. Although this lengthening of the core burning phase is always true of convective overshoot, we find that the extension of the overshoot and convective regions is sensitive to the temporal resolution adopted. With su GLYPH<14> ciently large values of f ov the upper boundary develops oscillatory behavior which can also a GLYPH<11> ect the lifetime. This behavior also occurs with the step-function implementation of overshoot. This instability is not seen in overshoot during hydrogen burning and has yet to be studied in detail.', '5. EVOLUTION BEYOND THE MAIN SEQUENCE AND WHITE DWARFS': 'Extending the verification of Paper I, we now compare to other available codes for intermediate-mass stars, 3-8 M GLYPH<12> . Wedescribe the techniques used by MESA star to evolve stars through the AGB phase to the white dwarf cooling sequence. We also demonstrate how MESA star incorporates compressional heating from accretion.', '5.1. Code Comparisons during Helium Core Burning': 'We start by comparing the results of MESA star to those from the Dartmouth Stellar Evolution Program (DSEP; Dotter et al. 2008) for stars with M = 3-8 M GLYPH<12> . In both cases, the models were evolved from the pre-main sequence to the depletion of helium in their cores. For completeness, the MESA star models were further evolved to the occurrence of the first helium thermal pulse.  \nAll models have an initial composition Y = 0 : 272, Z = 0 : 02, and no mass loss or rotation was included. The boundaries of mixing zones are determined by the Schwarzschild criterion with GLYPH<11> MLT = 2. In order to compare the codes, we do not allow overshooting or semiconvection. We adopt the Kunz et al. (2002) rate for 12 C( GLYPH<11>; GLYPH<13> ) 16 Oand the Imbriani et al.  \n(2004) rate for 14 N( p ; GLYPH<13> ) 15 O; for all other rates we use the NACRE compilation (Angulo et al. 1999). We use the OPAL Type 2 opacity tables (Iglesias & Rogers 1993) to account for the carbon- and oxygen-enhanced opacities during helium burning.  \nThe resulting tracks in the HR diagram of Figure 16 and the evolution in the T cGLYPH<26> c plane of Figure 17 show excellent agreement between the codes. Figures 18 and 19 show the hydrogen-burning luminosity, the helium-burning luminosity, and the extent of the convective core during convective helium core burning for a 4 M GLYPH<12> model (Fig. 18) and a 6 M GLYPH<12> model (Fig. 19). Table 3 gives a summary of the core hydrogen burning lifetime, the core helium burning lifetime, the final extent of the convective core during central helium burning, and the final carbon mass fraction X C in the core for each model. For the MESA models, we also show the maximum extent of the convective core during central hydrogen burning, the mass of the helium core before helium ignition, and the mass of the C / O core at the time of the first helium thermal pulse.  \nWe close with an additional comparison of the helium core burning phase of a M = 3 M GLYPH<12> , Z = 0 : 02 model computed by MESA to that of Straniero et al. (2003). Both models were evolved using the Kunz et al. (2002) rate for 12 C( GLYPH<11>; GLYPH<13> ) 16 O. The results for MESA star are a helium core burning lifetime of 83.6 Myr and final C / O mass fractions of X C = 0 : 43, X O = 0 : 55; Straniero et al. (2003) find a lifetime of 88 Myr and X C = 0 : 42, X O = 0 : 56.', '5.2. Making and Cooling White Dwarfs': "In the previous section, we discussed the evolution of 3-8 M GLYPH<12> stars up to the occurrence of the first He thermal pulse. In Paper I we showed detailed comparisons of the evolution of a 2 M GLYPH<12> star to the EVOL code (Herwig 2004), exhibiting the ability of MESA star to calculate multiple helium shell pulses. We now illustrate the final evolution of intermediate-mass stars, and how to construct white dwarfs (WDs) by using winds.  \nWe evolve 3 M GLYPH<12> , 5 M GLYPH<12> , and 7 M GLYPH<12> stars from the ZAMS using the test suite case make co wd . This makes use of RGB mass loss following Reimers (1975) with an e GLYPH<14> ciency parameter GLYPH<17> = 0 : 5 and AGB mass loss following Bloecker (1995) using GLYPH<17> = 0 : 1 until the occurrence of the first helium shell flash. At that time, an increased Bloecker GLYPH<17> = 5 is adopted to allow only a small number of thermal pulses before the wind mass loss eliminates the envelope. Such intervention allows MESA star to make a high-mass WD. To avoid shortening of timesteps due to radiation-dominated envelopes, these cases also use the MLT ++ capability described in x 7.2.  \nFigure 20 shows the resulting tracks on the HR diagram. The 3 M GLYPH<12> star underwent eight thermal pulses after the enhancement of Bloecker winds, while the 5 M GLYPH<12> and 7 M GLYPH<12> stars lost their envelopes so quickly that thermal pulses were immediately halted. The 5 M GLYPH<12> star ended up as an M = 0 : 844 M GLYPH<12> C / O WD with a helium shell of thickness M He = 0 : 009 M GLYPH<12> and a hydrogen envelope of M H = 2 : 3 GLYPH<2> 10 GLYPH<0> 5 M GLYPH<12> . Note that the C / O WD mass is only slightly larger than the C / O mass at the first thermal pulse (0 : 827 M GLYPH<12> ) reported in Table 3.  \nAfter removal of the envelope, the evolution of the white dwarf is continued through its cooling phase past solidification. We include gravitational settling and chemical di GLYPH<11> usion of the outermost layers. Figure 21 shows T -GLYPH<26> profiles taken at various e GLYPH<11> ective temperatures during the cooling of the M = 0 : 844 M GLYPH<12> C / O WD made from the 5 M GLYPH<12> star. The growing depth of the convection zone is shown by the dashed  \nFigure 12. Thermohaline mixing during the RGB phase of a Z = 0 : 02, 1 M GLYPH<12> model, initially rotating with an equatorial velocity of 10 km s GLYPH<0> 1 and adopting GLYPH<11> th = 2. In the left panel a Kippenhahn diagram shows, in mass coordinate and as function of model number, the locations of the retreating convective envelope (blue), of the H-burning shell (red) and of the thermohaline mixing region (magenta). The right panel shows di GLYPH<11> usion coe GLYPH<14> cient profiles extracted at model number 1849, which is the last model shown in the Kippenhahn plot. The H-burning shell and the convective envelope are shaded in red and blue, respectively. Thermohaline mixing (magenta line) transports chemicals between the burning shell and the convective envelope. Also shown are the di GLYPH<11> usion coe GLYPH<14> cients resulting from Eddington-Sweet circulation (ES), magnetic torques by dynamo generated fields (ST), Dynamical Shear (DSI), Secular Shear (SSI) and Goldreich-Schubert-Fricke (GSF) instability (see x 6 for details).  \n<!-- image -->  \n<!-- image -->  \n/circledot  \nTable 3  \nProperties of the 3-8 M GLYPH<12> evolution (masses in solar units). Selected quantities are also shown from DSEP for comparison.  \n- (1) Central H burning lifetime (Myr)\n- (2) Maximum extent of the convective core during core H burning\n- (3) Mass of the He core before central He ignition\n- (4) Central He burning lifetime (Myr)\n- (5) Stable final extent of the Schwarzschild convective core during core He burning\n- (6) Central mass fraction of 12 C at the end of core He burning\n- (7) Mass of the C / O core at the time of the first thermal pulse  \nline, and the open circles designate the H / He transition, while the filled circles denote the He / CO transition. Figure 22 illustrates the resulting L -T c relation as these models cool.  \nThe test suite case wd diffusion uses the implementation of di GLYPH<11> usion described in Paper I to evolve a WD of mass 0 : 535 M GLYPH<12> until the M H = 5 : 9 GLYPH<2> 10 GLYPH<0> 5 M GLYPH<12> hydrogen layer and the M He = 1 : 0 GLYPH<2> 10 GLYPH<0> 2 M GLYPH<12> helium layer approach di GLYPH<11> usive equilibrium. At this point, the WD has an e GLYPH<11> ective temperature of T e GLYPH<11> GLYPH<25> 5 ; 000 K. We show the resulting abundance profiles in Figure 23, and, for comparison, the abundance profiles derived from the analytic form for di GLYPH<11> usive equilibrium (eq. (22) of Althaus et al. 2003). This formula is obtained by integrating equation (A.5) of Arcoragi & Fontaine (1980) and assuming an ideal gas equation of state and complete ionization of both species.  \nThe specific treatment of convection can also impact WD evolution. In Paper I, MESA used the Cox & Giuli (1968) pre-  \nscription for convection as its default convective MLT, with the optional extension of Henyey et al. (1965). Since Paper I, we have added support for the formulations of Bohm-Vitense (1958), Bohm & Cassinelli (1971), and Mihalas (1978). In particular, the Bohm & Cassinelli prescription, often referred to as 'ML2,' is frequently employed in WD studies (e.g., Bergeron et al. 1995). In Figure 24 we show a comparison of the Brunt-Vaisala frequency calculated with MESA to that using the Warsaw envelope code (Paczy'nski 1969, 1970; Pamyatnykh 1999), assuming the ML2 prescription. This is the same WD as in Figure 23, but now at a lower T e GLYPH<11> = 11 ; 354 K. To more accurately integrate these opaque but thin layers, we reduce GLYPH<28> at the boundary of the model by a factor of 1000 from its photospheric value of 2 = 3. This calculation is a sensitive test of the envelope integrations because N 2 is a derivative of the envelope structure. The two codes give indistinguishable results for this case and all other cases that we have calculated.  \nFigure 13. History of convective core extent during the main sequence for a non-rotating 1 : 5 M GLYPH<12> star with various mixing options. The plot shows the boundary of convection not including the extent of semiconvection or overshooting.  \n<!-- image -->  \nFigure 14. The HR diagram for the non-rotating 1 : 5 M GLYPH<12> star with various mixing options. Tracks are displayed from ZAMS until depletion of core hydrogen to X = 10 GLYPH<0> 5 .  \n<!-- image -->  \nMESA now includes atmospheric tables based on the nongrey model atmospheres for hydrogen-atmosphere WDs (Rohrmann 2001; Rohrmann et al. 2012), spanning the following range of parameters: 2 ; 000 K GLYPH<20> T e GLYPH<11> GLYPH<20> 40 ; 000 K and 5 : 5 GLYPH<20> log g GLYPH<20> 9 : 5. Such an approach is necessary at T e GLYPH<11> . 6000 K, where WDs develop deeper convection zones. When the convection zone comes in contact with the degenerate, nearly isothermal core, energy is able to flow out of the core much more e GLYPH<14> ciently. The use of non-grey atmosphere models results in shallower convection zones, so this convective coupling of the core and envelope is delayed. For reliable cooling ages, we therefore recommend using non-grey atmospheres when T e GLYPH<11> . 6000 K. Figure 25 demonstrates the impact of non-grey atmospheres with the 0 : 535 M GLYPH<12> WD, which was cooled with and without the non-grey atmosphere.  \nMESA currently treats crystallization by employing the Potekhin & Chabrier (2010) EOS (PC EOS). The PC EOS is callable for arbitrary mixtures of chemical species and for  \nFigure 15. History of convective core extent during the core helium burning phase for a non-rotating 3 M GLYPH<12> star with various mixing options, as in Figure 13. Time is measured relative to the onset of the convective core burning. E GLYPH<14> cient semiconvection ( GLYPH<11> sc = 0 : 01) and ine GLYPH<14> cient overshooting ( f ov = 10 GLYPH<0> 5 ) coincide with the pure Schwarzschild model. The filled (open) circle indicates the time for which we display a profile detailing semiconvection (overshooting) in Figure 11.  \n<!-- image -->  \nFigure 16. Hertzsprung-Russell diagram for evolution of 3-8 M GLYPH<12> stars from the pre-main sequence through core helium depletion. Models are from MESA (thick grey lines) and DSEP (dashed black lines). Each curve is labeled with its corresponding initial mass in solar units.  \n<!-- image -->  \ndensities with log GLYPH<26> GLYPH<21> 2 : 8; it is applicable in the domains of non-degenerate and degenerate, non-relativistic and relativistic electrons, weakly and strongly coupled Coulomb liquids, and classical and quantum Coulomb crystals. The phase transition is first-order, so the PC EOS exhibits a latent heat between the solid and liquid phases, i.e., the entropy and internal energy both experience finite jumps. This energy is included in MESA star models of cooling white dwarfs through the  \nFigure 17. Same as Fig. 16, but in the T cGLYPH<26> c plane. The MESA models (thick grey lines) are evolved until the occurrence of the first thermal pulse.  \n<!-- image -->  \nFigure 18. History of hydrogen burning luminosity (top), helium-burning luminosity (center), and convective core extent (bottom) during the core helium burning phase for the 4 M GLYPH<12> models. Time is measured relative to the onset of the convective core.  \n<!-- image -->  \ngravitational source term in the energy equation,  \nGLYPH<15> grav GLYPH<17> GLYPH<0> T d S d t : (15)  \nThis form for GLYPH<15> grav replaces the default one (see eq. (16) below) in cells where GLYPH<0> GLYPH<21> 160 ( GLYPH<0> is the Coulomb coupling parameter), and is smoothly interpolated with the default form in cells where 130 GLYPH<20> GLYPH<0> < 160. The PC EOS uses the criterion GLYPH<0> = 175 to determine crystallization, but it is straightforward to include explicit crystallization curves for C / O and  \nFigure 19. History of hydrogen burning luminosity (top), helium-burning luminosity (center), and convective core extent (bottom) during the core helium burning phase for the 6 M GLYPH<12> models. Time is measured relative to the onset of the convective core.  \n<!-- image -->  \nFigure 20. Evolution of 3, 5 and 7 M GLYPH<12> models from zero-age main sequence to cooling white dwarfs. A Bloecker mass loss strips the stars of their envelopes on the thermally pulsing AGB to make the three C / O white dwarfs. The single 0 : 32 M GLYPH<12> He white dwarf was made with mass loss after the hydrogen main sequence for the 3 M GLYPH<12> model was completed.  \n<!-- image -->  \nother mixtures (e.g., Schneider et al. 2012; Medin & Cumming 2010). For example, using the parameters of the model in Figure 23, the age di GLYPH<11> erence at late times ( T e GLYPH<11> < 3 ; 500 K) between a model with and without the latent heat of crystallization is GLYPH<25> 0 : 8 Gyr; a slightly larger value would be obtained using the phase diagram of Schneider et al. (2012). MESA does not currently treat phase separation of di GLYPH<11> erent chemical species upon crystallization.  \nFigure 21. Profiles in log T -log GLYPH<26> space of the cooling 0 : 844 M GLYPH<12> C / O white dwarf evolved from a 5 M GLYPH<12> progenitor. Each model is labeled to the right by T e GLYPH<11> . The outermost point of the model is at GLYPH<28> = 25. Dotted curves denote convective regions. Going toward the interior, open circles designate the transition into the helium-rich shell, and filled circles designate the transition into the C / O core.  \n<!-- image -->  \nFigure 22. Surface luminosity as a function of central temperature for the cooling 0 : 32, 0 : 574, 0 : 844, and 0 : 928 M GLYPH<12> WDs evolved from 3, 5, and 7 M GLYPH<12> progenitors.  \n<!-- image -->  \nLow mass WDs ( M . 0 : 4 M GLYPH<12> ) with helium cores and hydrogen envelopes may be produced in binary systems when the envelope is stripped by the companion as the primary evolves up the giant branch (Iben 1991, and references therein). He-core WDs of mass M ' 0 : 4-0 : 5 M GLYPH<12> may also be produced through strong RGB winds (D'Cruz et al. 1996), although we do not discuss this possibility further here.  \nHere we discuss the prescription for stripping the envelope used in the test case make he wd . The first step is to evolve a star, M = 3 : 0 M GLYPH<12> in this example, from the PMS until a He core of the correct size has been made. The remnant total mass is determined by the mass interior to where the H  \nFigure 23. A comparison of time-dependent di GLYPH<11> usion calculations for a M = 0 : 535 M GLYPH<12> WD with MH = 5 : 9 GLYPH<2> 10 GLYPH<0> 5 M GLYPH<12> and M He = 1 : 0 GLYPH<2> 10 GLYPH<0> 2 M GLYPH<12> with MESA star (solid lines) to those assuming di GLYPH<11> usive equilibrium and an ideal gas equation of state (dashed lines).  \n<!-- image -->  \nFigure 24. A comparison of the Brunt-Vaisala frequency calculated with MESA (solid grey line) to that using the Warsaw envelope code pig35.f (dashed line) for the same WD in Figure 23, but at a cooler T e GLYPH<11> = 11 ; 354 K.  \n<!-- image -->  \nabundance has dropped below a preset value, for example, X H = 0 : 1, moving in from the surface. Next, the routine relax mass is used to remove mass from the model until it has the desired remnant mass. After the initial remnant has been constructed, di GLYPH<11> usion can then be turned on to allow an outer H layer to form. After this stage, normal evolution of the WD occurs, as shown in Figures 20 and 22.", '5.3. Compressional Heating and Accretion': 'Accretion onto stars occurs in many contexts and requires special treatment for the outermost layers added in each timestep. In particular, a special evaluation of the GLYPH<15> grav = GLYPH<0> T d S = d t term is required for fluid parcels that were not present in the previous timestep. Prior to addressing that subtlety, we restate (as discussed in x 6.2 of Paper I) that MESA star calculates GLYPH<15> grav of eq. (15) in terms of the local  \nFigure 25. The age di GLYPH<11> erence (non-grey minus grey) in Gyr as a function of T e GLYPH<11> .  \n<!-- image -->  \nthermodynamic variables ( T and GLYPH<26> ) used by MESA ,  \nGLYPH<15> grav = GLYPH<0> CPT " GLYPH<16> 1 GLYPH<0> r ad GLYPH<31> T GLYPH<17> d ln T d t GLYPH<0> r ad GLYPH<31>GLYPH<26> d ln GLYPH<26> d t # : (16)  \nMESA star takes the quantities in this equation as provided by eos , and computes the Lagrangian time derivatives to find GLYPH<15> grav. MESA star can alternatively work under the assumption that P = P gas + P rad, in which case MESA star treats P gas rather than GLYPH<26> as its basic variable (see x B.2 for a discussion). In that case,  \nGLYPH<15> grav = GLYPH<0> CPT " GLYPH<18> 1 GLYPH<0> 4 r ad P rad P GLYPH<19> d ln T d t GLYPH<0> r ad P gas P d ln P gas d t # : (17)  \nEither formulation can be used deep within the star, as long as the location is safely removed from any phase transition. Paper I described the validation of these formulations.  \nWe now turn to the complication which arises when GLYPH<15> grav needs to be evaluated in material that was not present in the previous timestep. Defining the envelope mass coordinate GLYPH<1> M GLYPH<17> M GLYPH<0> m , we need to resolve the entropy for GLYPH<1> M < GLYPH<14> M = ˙ M GLYPH<14> t , as the explicit Lagrangian time derivatives of eqs. (16) and (17) cannot be numerically evaluated. Since there can be important physics that needs to be resolved for these mass shells for GLYPH<1> M GLYPH<28> GLYPH<14> M , an approximation must be derived that allows for accurate modeling of the star\'s outermost layers without having to result to a dramatic shortening of GLYPH<14> t .  \nThe luminosity L acc = GM ˙ M = R from the accretion shock (or boundary layer) goes outwards and does not determine the entropy of the material as it becomes part of the hydrostatically adjusting star. Rather, the entropy of the material at GLYPH<1> M GLYPH<28> GLYPH<14> M is determined by the the transport of L (Nomoto & Sugimoto 1977; Nomoto 1982; Townsley & Bildsten 2004). Consider such an outermost layer, where there are two relevant timescales, the thermal time, t th = CPT GLYPH<1> M = L , and the local accretion time, t acc = GLYPH<1> M = ˙ M . In nearly all relevant cases, the ratio t th = t acc = CPT ˙ M = L GLYPH<28> 1; this implies that the fluid element adjusts its temperature to that needed to transport the stellar luminosity from deep within. This simplifies GLYPH<15> grav in  \nthat part of the star (following Townsley & Bildsten 2004) to  \nGLYPH<15> grav = CPTGm ˙ M 4 GLYPH<25> r 4 P ( r ad GLYPH<0> r T ) ; (18)  \nenabling accurate modeling within MESA star of nearly all fluid elements that become part of the star during each timestep, many of which have envelope mass coordinates GLYPH<1> M GLYPH<28> ˙ M GLYPH<14> t .  \nWe give an explicit example of this thin-shell radiative calculation of GLYPH<15> grav in a C / O white dwarf accreting hydrogenrich material and undergoing classical nova (CN) cycles. We present two models accreting at rates of ˙ M = 10 GLYPH<0> 11 M GLYPH<12> yr GLYPH<0> 1 and 10 GLYPH<0> 10 M GLYPH<12> yr GLYPH<0> 1 . Both cases were evolved from a 0 : 6 M GLYPH<12> starting model with T c = 10 7 K which had undergone a few flashes while accreting at ˙ M = 10 GLYPH<0> 11 M GLYPH<12> yr GLYPH<0> 1 . The accreted material has solar-like composition X = 0 : 70, Y = 0 : 26, and Z = 0 : 04 where the metal mass fractions are taken from Lodders (2003).  \nProfiles of the envelope during the mass accumulation phase between CN outbursts for the two accretion rates are displayed in Figures 26 and 27. Each line represents a di GLYPH<11> erent time in the accumulation cycle up to the unstable ignition, when the hydrogen mass reaches M H = M ign. All material at pressures smaller than that shown by the open circle is new to the model in that timestep (e.g., it has GLYPH<1> M < GLYPH<14> M ) and employs the modified GLYPH<15> grav of eq. (18). This highlights the significance of this approximation as it allows MESA star to calculate material properties at GLYPH<1> M GLYPH<24> 10 GLYPH<0> 8 GLYPH<14> M . The solid points show where GLYPH<15> grav switches to the explicit form employing the Lagrangian time derivatives, such as eq. (16).  \nThe middle panel shows GLYPH<15> grav P / GLYPH<15> grav GLYPH<1> M , which reflects the contribution of GLYPH<15> grav to the outward luminosity. The discontinuity of GLYPH<15> grav at the solid point reflects the error associated with the abrupt transition in the calculational approach. The substantially larger luminosity of the early ( M H = M ign = 0 : 22) stages is due to the ongoing transfer of heat from the previous outburst. The near-discontinuous drop in GLYPH<15> grav occurs at the base of the hydrogen-rich envelope, and reflects the jump in composition from the accreted material to the nearly pure 4 He layer. The expected amplitude of the jump in GLYPH<15> grav depends on both the composition jump and the local degree of electron degeneracy (see Appendix B of Townsley & Bildsten 2004 for a discussion).', '6. ROTATION': "A star's rotational energy is usually a small fraction of the gravitational energy: for the Sun it is GLYPH<24> 10 GLYPH<0> 5 and for a 25 M GLYPH<12> star rotating with a typical equatorial velocity 3 eq = 200 km s GLYPH<0> 1 on the main sequence it is GLYPH<24> 0 : 04. Therefore the e GLYPH<11> ects on the stellar hydrostatic equilibrium are marginal, with the exception of stars close to critical rotation (see x 6.4). Even in the case of a small perturbation to hydrostatic equilibrium, rotation induces a modification to the star's thermal equilibrium (von Zeipel 1924). Together with the emergence of rotationally-induced dynamical and secular instabilities, this can significantly a GLYPH<11> ect the evolution of stars (Maeder & Meynet 2000b). Due to the destabilizing e GLYPH<11> ect of increasing radiation pressure, rotation is particularly important in massive stars (see, e.g., Heger et al. 2000; Meynet & Maeder 2000). Moreover, the final fate of a massive star depends chiefly on the relative importance of rotation during its evolution (e.g., Heger et al. 2000; Hirschi et al. 2004; Heger et al. 2005; Yoon et al. 2006; Woosley & Heger 2006; Ekstrom et al.  \nFigure 26. Envelope profiles as a function of pressure of the accreting white dwarf for three instants during the mass accumulation phase; ˙ M = 10 GLYPH<0> 11 M GLYPH<12> yr GLYPH<0> 1 model. The top panel shows temperature, the central panel shows the gravitational energy release rate, and the bottom shows the luminosity. Material to the right of the open circle is newly accreted. The code treats material to the right of the filled circle using the thin-shell radiative calculation of GLYPH<15> grav.  \n<!-- image -->  \nFigure 27. Same as Fig. 26, but for model accreting 10 GLYPH<0> 10 M GLYPH<12> yr GLYPH<0> 1 .  \n<!-- image -->  \n2012; Georgy et al. 2012; Langer 2012).  \nHere we describe the implementation of rotation in MESA star . We briefly discuss the modification to the stellar structure equations and the inclusion of rotationally- and magnetically-induced mixing. Magnetic fields generated by di GLYPH<11> erential rotation in radiative regions have been implemented following the work of Spruit (2002) and in the same  \nfashion as in Petrovic et al. (2005) and Heger et al. (2005). Rotationally enhanced mass loss is also discussed.  \nWe compare rotating massive-star models calculated with MESA star to previous calculations performed with KEPLER (Heger et al. 2005). We also directly compare runs from MESA star and STERN (Petrovic et al. 2005; Yoon & Langer 2005; Brott et al. 2011). The purpose of these tests is to verify our implementation of rotation, which is derived from STERN. We do not compare to codes that have a di GLYPH<11> erent implementation of rotation (e.g., Hirschi et al. 2004; Ekstrom et al. 2012; Georgy et al. 2012; Potter et al. 2012b,a). Although beyond the scope of this paper, such comparisons are critical when coupled to observations of the e GLYPH<11> ects of rotation in stars (e.g., Hunter et al. 2007; Evans et al. 2011) including asteroseismology (Beck et al. 2012; Mosser et al. 2012).", '6.1. Implementation of Shellular Rotation': 'Stellar structure deviates from spherical symmetry in the presence of rotation. While the structure is inherently threedimensional, it su GLYPH<14> ces to solve the stellar structure equations in one dimension if the angular velocity, ! , is constant over isobars (the so-called shellular approximation; see, e.g., Meynet & Maeder 1997). This is expected in the presence of strong anisotropic turbulence acting along isobars. In radiative regions such turbulence is a consequence of di GLYPH<11> erential rotation (Zahn 1992) and e GLYPH<14> ciently erases gradients along isobars and enforces shellular rotation (Meynet & Maeder 1997). Turbulence in the vertical direction (i.e., perpendicular to the isobars) is much weaker due to the stabilizing e GLYPH<11> ect of stratification. In MESA star we adopt the shellular approximation (Meynet & Maeder 1997) and calculate the modification to the stellar equations due to centrifugal acceleration (Kippenhahn & Thomas 1970; Endal & Sofia 1976).', '6.1.1. Stellar Structure': 'Anisobar with volume V P and surface area S P deviates from spherical symmetry in the presence of rotation. However one can retain a 1D approximation by re-defining the radius coordinate as the radius of a sphere containing the same volume V P = 4 GLYPH<25> r P 3 = 3, allowing an equation of continuity in the usual form  \n@ m P @ r P ! t = 4 GLYPH<25> r P 2 GLYPH<26> ; (19)  \nwith GLYPH<26> being the density and m P the mass enclosed by S P. The energy equation also retains its usual, non-rotating form  \n@ L P @ m P ! t = GLYPH<15> nuc GLYPH<0> GLYPH<15>GLYPH<23> + GLYPH<15> grav ; (20)  \nwhere L P is the rate of energy flow through the equipotential surface S P. Then the next step is to define mean values for the quantities varying on isobars,  \nh GLYPH<1> i GLYPH<17> 1 S P I S P GLYPH<1> d GLYPH<27> ; (21)  \nwhere d GLYPH<27> is an isobaric surface area element. The equation of momentum balance can be written as  \n@ P @ m P ! t = GLYPH<0> Gm P 4 GLYPH<25> r P 4 f P GLYPH<0> 1 4 GLYPH<25> r P 2 @ 2 r P @ t 2 ! m P ; (22)  \nwhere P is the pressure, G is the gravitational constant and t the time. The last term in the equation is the inertia term.  \nRotation enters the momentum equation through the quantity  \nf  \nP  \nf P GLYPH<17> 4 GLYPH<25> r P 4 Gm P S P h g GLYPH<0> 1 i GLYPH<0> 1 ; (23)  \nwhere g GLYPH<17> GLYPH<12> GLYPH<12> GLYPH<12> ~ g GLYPH<12> GLYPH<12> GLYPH<12> , with g the e GLYPH<11> ective gravitational acceleration ( ~ g is normal to S P). Then the radiative temperature gradient becomes  \n@ ln T @ ln P ! t = 3 GLYPH<20> 16 GLYPH<25> acG P T 4 L P m P f T f P 2 6 6 6 6 4 1 + r P 2 Gm P f P @ 2 r P @ t 2 ! m P 3 7 7 7 7 5 GLYPH<0> 1 ; (24)  \nwith a the radiation constant, GLYPH<20> is the opacity, T the temperature and L P, the energy flux through S P. The last factor on the right hand side accounts for inertia, and  \nf T GLYPH<17> 4 GLYPH<25> r P 2 S P ! 2 GLYPH<16> h g ih g GLYPH<0> 1 i GLYPH<17> GLYPH<0> 1 : (25)  \nIn rotating models the values of f T and f P di GLYPH<11> er from 1 mostly in the outer stellar layers. Limits to the minimum values of f T and f P are set in the code (default values are 0.95 and 0.75, respectively). This prevents numerical instabilities in models approaching critical rotation ( GLYPH<10> = GLYPH<10> crit = 1, see 6.4). In such cases the outer layers greatly deviates from spherical symmetry and the results from 1D calculations should be considered particularly uncertain.', '6.1.2. Mixing and angular momentum transport': 'Transport of angular momentum and chemicals due to rotationally-induced instabilities is implemented in a di GLYPH<11> usion approximation (e.g., Endal & Sofia 1978; Pinsonneault et al. 1989; Heger et al. 2000). This choice has also been adopted by other stellar evolution codes (e.g., KEPLER, Heger et al. 2000; STERN, Yoon & Langer 2005). We stress that this is not the only possibility, and other groups have implemented a di GLYPH<11> usion-advection approach (e.g., GENEVA, Eggenberger et al. 2008; RoSE, Potter et al. 2012b). The RoSE code can switch between the two di GLYPH<11> erent implementations. The two approaches are equivalent for the transport of chemicals. Potentially large di GLYPH<11> erences can arise, however, for the transport of angular momentum. A detailed description of the advection-di GLYPH<11> usion equation for angular momentum is given in Zahn (1992) and Maeder & Zahn (1998).  \nIn MESA star the turbulent viscosity GLYPH<23> is determined as the sum of the di GLYPH<11> usion coe GLYPH<14> cients for convection, semiconvection, and rotationally-induced instabilities. In convective regions, the very large di GLYPH<11> usion coe GLYPH<14> cient implies that the rotation law is not far from solid body. This is a very common assumption in stellar evolution codes (e.g. Pinsonneault et al. 1989; Heger et al. 2000; Eggenberger et al. 2008); note however that helioseismology has clearly shown this is not the case for the solar convection zone (e.g. Brown et al. 1989; Thompson et al. 1996; Schou et al. 1998). MESA star calculates di GLYPH<11> usion coe GLYPH<14> cients for five di GLYPH<11> erent rotationally-induced mixing processes: dynamical shear instability, Solberg-Høiland instability, secular shear instability, Eddington-Sweet circulation, and the Goldreich-SchubertFricke instability. See Heger et al. (2000) for a detailed description of the physics of the di GLYPH<11> erent instabilities and the calculation of the respective di GLYPH<11> usion coe GLYPH<14> cients. These enter the angular momentum and abundance di GLYPH<11> usion equations that are solved at each timestep (see x B.6).', '6.2. Internal Magnetic Fields': 'It has been suggested that di GLYPH<11> erential rotation in the radiative layers of a star can amplify a seed magnetic field. Such a dynamo process has been proposed by Spruit (2002, Spruit-Tayler dynamo); a theoretical debate on this is still ongoing (Braithwaite 2006; Zahn et al. 2007; Denissenkov & Pinsonneault 2007). From the observational point of view, pure hydrodynamic models fail to predict the solar core rotation (e.g., Pinsonneault et al. 1989), with the exception of models that include transport of angular momentum by gravity waves (Charbonnel & Talon 2005). Models that include the Spruit-Tayler dynamo can reproduce the flat rotation profile of the Sun. Note however that these have di GLYPH<14> culty explaining the core-envelope decoupling observed in low-mass, young cluster stars (Denissenkov et al. 2010). On the other hand, observations of the final spins of both WDs and neutron stars (Heger et al. 2005; Suijs et al. 2008) suggest that angular momentum transport with an e GLYPH<14> ciency similar to the torques provided by the Spruit-Tayler dynamo operates. Models that only include angular momentum transport through rotational instabilities do not produce the core-envelope ratio of angular velocity observed through the splitting of mixed modes in red giant stars (Eggenberger et al. 2012).  \nMESA star accounts for transport by magnetic fields of angular momentum and chemicals due to the Spruit-Tayler dynamo. We refer to Spruit (2002) for a description of the physics of the dynamo loop and to Maeder & Meynet (2003), Maeder & Meynet (2004) and Heger et al. (2005) for a discussion of its inclusion in stellar evolution codes. We implement the Spruit-Tayler dynamo in MESA star following KEPLER (Heger et al. 2005) and STERN (Petrovic et al. 2005).', '6.3. Surface Magnetic Fields': 'Rotating stars that have a significant outer convective zone can produce surface magnetic fields through a dynamo (See e.g., Brandenburg & Subramanian 2005, for a review on astrophysical dynamos). This is the case for low-mass main sequence stars below about 1.5 M GLYPH<12> , and observationally the break in the rotation properties around this mass is attributed to the presence of magnetized stellar winds (e.g., Schatzman 1962; Kawaler 1988). Note that dynamo action in a subsurface convective layer is in principle possible also in early-type stars (Cantiello et al. 2009; Cantiello & Braithwaite 2011). Surface magnetic fields can also be of fossil origin, as is usually discussed in the context of Ap stars (Braithwaite & Spruit 2004). Whatever the origin of surface magnetic fields, these are expected to couple to the wind mass-loss and, if strong enough, produce magnetic braking (e.g., Weber & Davis 1967; ud-Doula & Owocki 2002; Meynet et al. 2011). Such magnetic braking has been directly observed in the case of the main sequence massive star GLYPH<27> -Ori E (Townsend et al. 2010). Here we do not include the physics of magnetic braking, as we only consider the evolution of stars without surface magnetic fields.', '6.4. Rotationally-Enhanced Mass Loss': 'We include the rotational modification to the wind mass loss rate (Friend & Abbott 1986; Langer 1998; Heger & Langer 1998; Maeder & Meynet 2000a). Similar to other codes (e.g., Heger et al. 2000; Brott et al. 2011; Potter et al. 2012a), in MESA star the stellar mass loss is enhanced as the  \nrotation rate increases according to the prescription  \n˙ M ( GLYPH<10> ) = ˙ M (0) 1 1 GLYPH<0> GLYPH<10> = GLYPH<10> crit ! GLYPH<24> ; (26)  \nwhere GLYPH<10> is the value of the surface angular velocity and GLYPH<10> crit is the critical angular velocity at the surface. This last quantity is defined as GLYPH<10> 2 crit = (1 GLYPH<0> L = L Edd) GM = R 3 , where L Edd = 4 GLYPH<25> cGM =GLYPH<20> is calculated as a mass-weighted average in a userspecified optical depth range (default value GLYPH<28> 2 [1 GLYPH<0> 100]). In MESA star the default value for the exponent GLYPH<24> is 0.43 (Langer 1998). Other implementations of rotationally enhanced mass loss can be found in Maeder & Meynet (2000a) and Georgy et al. (2011).  \nFor stars approaching GLYPH<10> = GLYPH<10> crit = 1, the mass loss calculated using equation (26) diverges. Notice that luminous stars can approach this limit without having to rotate very rapidly as GLYPH<10> crit ! 0 when L = L Edd ! 1. Following Yoon et al. (2010) we limit the mass loss timescale to the thermal timescale of the star GLYPH<28> KH  \n˙ M = min " ˙ M ( GLYPH<10> ) ; f M GLYPH<28> KH # (27)  \nwhere f is an e GLYPH<14> ciency factor of order unity (default value is f = 0 : 3).', '6.5. Initial models': 'In all the rotating models presented in this paper, rotation is initialized by imposing a solid body rotation law on the zero-age main sequence (ZAMS, L = L nuc). In these massive stars this is motivated by the presence of rotationally-induced angular momentum transport during the pre-main sequence evolution. This alone is able to enforce a state of close-torigid rotation by the time the star reaches the ZAMS (Heger et al. 2000). Overall initial solid body rotation is a common choice in stellar evolution codes, but other rotational laws are certainly possible.  \n6.6. Test Cases: 15 M GLYPH<12> and 25 M GLYPH<12>  \nAs a first test we initialize a 15M GLYPH<12> model with Z = 0 : 02 and initial equatorial rotational velocity 3 eq = 200 km s GLYPH<0> 1 and run two calculations:  \n- GLYPH<15> 15MAG includes the e GLYPH<11> ects of rotation and SpruitTayler magnetic fields on both the transport of chemicals and angular momentum.\n- GLYPH<15> 15ROT includes only the e GLYPH<11> ect of rotation on both the transport of chemicals and angular momentum;  \nThe initial conditions have been calibrated to match as closely as possible the KEPLER 15 M GLYPH<12> models (Heger et al. 2005). Moreover, we directly compare the MESA star models with calculations from STERN (see e.g., Yoon & Langer 2005; Yoon et al. 2006). In particular we adopt a value of fc = 1 = 30 for the ratio of the turbulent viscosity to the di GLYPH<11> usion coe GLYPH<14> -cient and a value f GLYPH<22> = 0 : 1 for the sensitivity to GLYPH<22> -gradients (see Heger et al. 2000, for a discussion of these calibration parameters). The Ledoux criterion is used for the treatment of convective boundaries, together with semiconvection ( GLYPH<11> sc = 1). We use GLYPH<11> MLT = 1 : 6, mass loss as in Yoon et al. (2006) with rotational enhancement as described in x 6.4.  \nIn Fig. 28 we show the evolutionary track and the evolution of surface equatorial rotational velocity for the 15MAG model. Results of a similar calculation using STERN are  \nshown as a dashed curve. The two results are in excellent agreement. Small di GLYPH<11> erences in luminosity and lifetimes are not unexpected, as we have only matched the physics of rotation between the two calculations and not other ingredients. Values for the di GLYPH<11> usion coe GLYPH<14> cients for rotationally induced mixing and magnetic torques during the main sequence of 15MAG are shown in Fig. 29. The comparison reveals a very good agreement. Both stars are kept in solid-body rotation during the main sequence by the e GLYPH<14> cient transport of angular momentum provided by the Eddington-Sweet circulation and Spruit-Tayler magnetic fields.  \nThe amplitude and location of the azimuthal ( B GLYPH<30> ) and radial ( Br ) components of the magnetic fields during di GLYPH<11> erent phases of the evolution of 15MAG are shown in Fig. 30. As expected, these fields are generated only in radiative regions of the star and B GLYPH<30> > Br (Spruit 2002). As the star evolves away from the main sequence its structure departs from solidbody rotation with the core rotating faster than the envelope. During this stage the role of magnetic fields is very important in transporting angular momentum from the core to the envelope. The e GLYPH<11> ect can be seen in Fig. 31, which shows the evolution of the internal specific angular momentum in models 15ROT and 15MAG. The presence of magnetic torques results in a dramatic spin-down of the core of 15MAG with respect to 15ROT (see also Table 4). These results are in very good agreement with the ones obtained by STERN and KEPLER.  \nAs a second test, we now evolve a 25 M GLYPH<12> model (25MAG) with the same physics as in 15MAG. Figure 32 directly compare results with calculations performed with STERN. In Fig. 33 we show a detailed comparison of the evolution of the internal specific angular momentum profile. We find a very good quantitative agreement between MESA star and STERN down to He depletion in the core. The timescale for nuclear burning decreases substantially after He-burning and becomes shorter than the angular momentum transport timescale after C depletion. Thus only minor changes in the final angular momentum content of the stellar core are expected after this stage. Figure 34 shows the full evolution of the specific angular momentum profile of the MESA star calculation from ZAMS to Si exhaustion.', '6.7. Rapidly Rotating Massive Stars': 'MESA star can calculate the evolution to core collapse of rapidly rotating massive stars. Rotational instabilities can be e GLYPH<14> cient enough to erase the compositional gradients built by nuclear burning. In such cases the model never develops a compositional stratification and remains almost completely mixed throughout its evolution (Maeder 1987b; Yoon & Langer 2005; Woosley & Heger 2006). This process leads to a bifurcation in the HR-diagram, with stars above a certain mass and rotation rate becoming more luminous and hotter. The threshold required for this bifurcation depends mostly on the initial mass of the star (Yoon et al. 2006). Metallicity also plays an important role, as angular momentum is lost through line-driven stellar winds, with mass-loss rates depending on the metallicity at the stellar surface (Vink et al. 2001). For the calculations in this section, we adopt the same mass-loss prescription of Yoon et al. (2006).  \nFigure 35 shows the evolution of two 16 M GLYPH<12> models at metallicity Z = 0 : 0002 with rotation initialized at the ZAMS. One model is rotating very rapidly, with 3 eq = 450 km s GLYPH<0> 1 (corresponding to GLYPH<10> = GLYPH<10> crit = 0 : 55 and J = 3 : 23 GLYPH<2> 10 52 ergs s),  \n<!-- image -->  \nFigure 28. Comparison of evolutionary tracks (left) and equatorial rotational velocities (right) for a 15 M GLYPH<12> model with Z = 0 : 02 rotating initially with 3 eq = 200 km s GLYPH<0> 1 (15MAG). The solid black line shows MESA star results, and the dashed gold line shows the STERN calculations.  \n<!-- image -->  \n2  \n<!-- image -->  \nFigure 29. Same as Fig. 28. As function of mass coordinate we plot the values of the di GLYPH<11> usion coe GLYPH<14> cient for convection (MLT), Eddington-Sweet circulation (ES), magnetic torques by dynamo generated fields (ST), Dynamical Shear (DSI), Secular Shear (SSI) and Goldreich-Schubert-Fricke (GSF) instability. Following STERN, we turn o GLYPH<11> the Solberg-Høiland instability (SH) for this comparison. This does not a GLYPH<11> ect the results, as the di GLYPH<11> usion coe GLYPH<14> cient for SH is usually smaller than the ones for ES and ST. The values of the specific angular momentum j and the angular velocity ! are also plotted. Left panel shows the results using MESA star , while the right panel shows analogous STERN calculations.  \n<!-- image -->  \nwhile the other rotates at 3 eq = 280 km s GLYPH<0> 1 (corresponding to GLYPH<10> = GLYPH<10> crit = 0 : 39 and J = 2 : 52 GLYPH<2> 10 52 ergs s). The model with GLYPH<10> = GLYPH<10> crit = 0 : 55 avoids the core-envelope structure and becomes a compact Wolf-Rayet star. The absence of a RSG phase eliminates the large magnetic torques from an extended envelope. The evolution of the internal profile of specific angular momentum in the two models clarify this point: the model with GLYPH<10> = GLYPH<10> crit = 0 : 39 becomes a RSG, and the core spins down rapidly. When it reaches core-collapse its structure is extended, as implied by the large free-fall timescale shown in the left panel of Fig. 35. As a consequence, there is not enough angular momentum in its core to build an accretion disk around a newly formed compact object. This model is expected to produce a Type IIP supernova. On the contrary, the model with GLYPH<10> = GLYPH<10> crit = 0 : 55 is compact (the free-fall timescale is on the order of seconds, right panel of Fig. 35) and has enough angular momentum to produce an  \naccretion disk around the central compact object. Therefore this model is a candidate progenitor for a long gamma-ray burst (Woosley 1993). This last calculation can be directly compared to the KEPLER model 16TI in Woosley & Heger (2006).  \nWe further test MESA capabilities by evolving two rotating 40 M GLYPH<12> models at Z = 10 GLYPH<0> 5 . One model is initialized at the ZAMS with 3 eq = 260 km s GLYPH<0> 1 , while the other has 3 eq = 630 km s GLYPH<0> 1 . The results of these calculations can be compared with the models shown in Yoon & Langer (2005). Figure 36 shows that for the more rapidly rotating model, rotational mixing (mainly due to the Eddington-Sweet circulation) is large enough that the star evolves blueward in the HRdiagram. This evolution results in a compact configuration and enough angular momentum to fulfill the requirements of the collapsar scenario for long gamma-ray bursts, as shown in Fig. 37 (right panel). On the other hand, the slower rotat-  \nFigure 30. Magnetic field structure and angular momentum distribution for model 15MAG at di GLYPH<11> erent evolutionary stages (see Table 4). The curves show profiles for specific angular momentum ( j ), angular velocity ( ! ), azimuthal and radial components of magnetic field ( B GLYPH<30> and Br ). The shaded regions represent convective parts of the star. Compare with Fig. 1 in Heger et al. (2005).  \n<!-- image -->  \nFigure 31. Specific angular momentum distribution at di GLYPH<11> erent evolutionary stage for 15MAG and 15ROT. See Table 4 for the definitions of these times. Compare with Fig. 2 in Heger et al. (2005).  \n<!-- image -->  \nEvolution of Angular Momentum at Fiducial Mass Coordinates for a Z = 0 : 02, 15 M GLYPH<12> star initially rotating with 3 eq = 200 km s with (15MAG) and without (15ROT) the inclusion of magnetic fields.  \nTable 4 GLYPH<0> 1  \nNOTE: a Results from Table 1 of Heger et al. (2005); b See e.g., Petrovic et al. (2005); Yoon & Langer (2005); Yoon et al. (2006); c 40 % central hydrogen mass fraction; d 1 % hydrogen left in the core; e 1 % helium burnt; f 50% central helium mass fraction; g 1 % helium left in the core; h central temperature of 5 GLYPH<2> 10 8 K; i central temperature of 1.2 GLYPH<2> 10 9 K; j central oxygen mass fraction drops below 5 %; k central Si mass fraction drops below 10 GLYPH<0> 4 ;  \n<!-- image -->  \nFigure 32. Same as Fig. 28, except for the 25MAG model.  \n<!-- image -->  \nFigure 34. Evolution to Si-depletion of the internal specific angular momentum for the 25MAG model.  \n<!-- image -->  \nenclosed mass/M  \nFigure 33. Evolution of internal specific angular momentum for the 25MAG model. Solid lines show MESA star result, while dashed lines refer to STERN. Left panel shows the evolution from zero age main sequence to He ignition. Right panel shows the evolution during core He-burning (from 50% of He in the core to He depletion). Notice the di GLYPH<11> erent axis range in the two plots.  \n<!-- image -->  \n<!-- image -->  \n2  \ning model becomes a RSG and loses most of its core angular momentum, as shown in Fig. 37 (left panel).', '7. MASSIVE STELLAR EVOLUTION': "Modeling massive stars is numerically di GLYPH<14> cult. One problem is they develop loosely bound, radiation pressure dominated envelopes that can cause density and gas pressure inversions. Indeed, very massive stars are observed to su GLYPH<11> er sporadic 'eruptions' of extreme mass loss (i.e., the Luminous Blue Variables), and the tendency to form inversions has been speculatively mentioned as playing a role in such episodes (see Humphreys & Davidson 1994, and references therein). This environment poses a physical and numerical challenge that all stellar evolution codes must address to evolve mas-  \nsive stars past the main sequence. In this section we discuss MESA star 's capability to evolve rotating massive stars from their zero age main sequence to core-collapse.", '7.1. Evolution of Massive Stars with MESA': 'Previous computations with MESA star found these envelopes to be numerically (and probably physically) unstable. This is a known issue in the literature (e.g., Maeder 1987a), which reveals the limitations of the 1D treatment of late phases of evolution of massive stars. The evolution of stars with radiation-dominated envelopes can require prohibitively short timesteps in MESA star if the standard mixing length theory is adopted. This problem usually appears during the evolution of high mass and / or high metallicity stars after hydrogen-core burning and prevents evolution to core collapse. We discuss in 7.2 our treatment of superadiabatic convection in these envelopes, which allows uninterrupted evolution, from ZAMS to core collapse.  \nSince it is relevant to later discussions we start with a plot of the OPAL opacity data (Iglesias & Rogers 1996) and 60 M GLYPH<12> ZAMS models in Figure 38. The plot is inspired by Figure 1 of Cantiello et al. (2009). The left-hand panel of Figure 38 shows the OPAL data for five di GLYPH<11> erent Z values at constant X = 0 : 7 and log( GLYPH<26>= T 6 3 ) = GLYPH<0> 5, where T 6 is the temperature in units of 10 6 K. The right-hand panel shows the opacity profiles of five 60 M GLYPH<12> ZAMS models for the same five Z values. The model profiles exhibit the same general behavior in the opacity-temperature profile as the raw opacity data. Of particular importance are the iron opacity bumps that occur at log T GLYPH<25> 5 : 3 and 6.3. These bumps cause both the local radiation pressure to dominate and the luminosity to approach the Eddington luminosity L Edd.  \nWhere both the pressure is dominated by radiation and L rad approaches L Edd, specific conditions can be reached that cause convection and inversions in density and gas pressure. To define the conditions under which these occur, we follow the discussion of Joss et al. (1973), going from high to low L rad. We assume that d T = d r < 0, d P = d r < 0, and that the inertial terms in the momentum equation are small. First, we establish a condition for the occurrence of an inversion in the gas pressure P gas. Recasting the equation for the temperature gradient  \n/s)  \n2  \nlog j (cm  \n<!-- image -->  \nFigure 35. Specific angular momentum distribution for the two 16 M GLYPH<12> models. In the top panels, the solid curves show the distribution of specific angular momentum at di GLYPH<11> erent evolutionary stages. The other curves in the top panel show the specific angular momentum of the last stable orbit around a Schwarzschild black hole, a maximally rotating Kerr black hole ( a = 1), and a black hole with a Kerr parameter corresponding to the angular momentum content of the stellar progenitor at that mass coordinate. The bottom panels show the free-fall time at the relative mass coordinate at the end of Si-burning. Notice the di GLYPH<11> erent ranges of the y-axis. These models can be compared to the calculations of Woosley & Heger (2006), in particular their models 16SG and 16TI respectively.  \n<!-- image -->  \nFigure 36. Evolution in the HR diagram for two rotating 40 M GLYPH<12> models at Z = 10 GLYPH<0> 5 . The slower rotating model evolves toward the red part of the HRD; the other model evolves toward the blue part of the HRD. The internal evolution of the angular momentum is shown in Fig. 37. This can be compared to Fig. 2 of Yoon & Langer (2005).  \n<!-- image -->  \ngives  \nL rad = GLYPH<0> 4 GLYPH<25> r 2 c GLYPH<26>GLYPH<20> d P rad d r ; (28)  \nand using the equation of hydrostatic equilibrium, one obtains  \nd P rad d P = L rad L Edd : (29)  \nWriting d P gas = d r = d( P GLYPH<0> P rad) = d r and using equation (29) and the fact that both P rad and P monotonically decrease with r ,  \none obtains  \nd P gas d r = d P rad d r ! " L Edd L rad GLYPH<0> 1 # : (30)  \nSince d P rad = d r < 0, equation (30) implies that for L rad > L Edd, the gas pressure gradient will increase outward, d P gas = d r > 0, as shown by Joss et al. (1973).  \nThe next step is to establish the condition for a density inversion to occur. Writing the gas equation of state as P gas = P gas( GLYPH<26>; P rad) gives  \nd P gas d r = @ P gas @GLYPH<26> ! P rad d GLYPH<26> d r + @ P gas @ P rad ! GLYPH<26> d P rad d r : (31)  \nSolving eq. (31) for d GLYPH<26>= d r and using eq. (30) eliminates d P gas = d r . Gas equations of state have ( @ P gas =@GLYPH<26> ) T > 0, so that for d GLYPH<26>= d r > 0 (a density inversion), one must have  \nd P rad d r ! 2 6 6 6 6 4 L Edd L rad GLYPH<0> 1 GLYPH<0> @ P gas @ P rad ! GLYPH<26> 3 7 7 7 7 5 > 0 :  \nRecognizing that d P rad = d r < 0, we find that a density inversion occurs when  \nL rad L Edd > L inv L Edd GLYPH<17> 2 6 6 6 6 4 1 + @ P gas @ P rad ! GLYPH<26> 3 7 7 7 7 5 GLYPH<0> 1 : (32)  \nThis equation is identical to eq. (8) of Joss et al. (1973). Since under conditions of interest ( @ P gas =@ P rad) GLYPH<26> > 0, we have L inv < L Edd. For L inv < L rad < L Edd, a density inversion will occur even though d P gas = d r < 0.  \nNext, we shall consider the luminosity L onset at which convection occurs. In a convective region, the entropy is either constant or declining with radius. Hence, convection will occur once  \nd ln P rad d ln P > @ ln P rad @ ln P ! s ; (33)  \nusing equation (29) and solving for the luminosity, we find that convection starts once  \nL rad L Edd > L onset L Edd GLYPH<17> 1 GLYPH<0> P gas P ! @ ln P rad @ ln P ! s : (34)', 'Paxton et al.': '<!-- image -->  \nFigure 37. Specific angular momentum distribution for two 40 M GLYPH<12> models at Z = 10 GLYPH<0> 5 with 3 eq = 260 km s GLYPH<0> 1 (top panel) and 630 km s GLYPH<0> 1 (bottom panel). Lines are showing the distribution of specific angular momentum at di GLYPH<11> erent stages of the evolution, together with the specific angular momentum of the last stable orbit around a Schwarzschild black hole, a maximally rotating Kerr black hole ( a = 1) and a black hole with a Kerr parameter corresponding to the angular momentum content of the stellar progenitor at that mass coordinate. Note the di GLYPH<11> erent ranges of the y-axis. The evolutionary tracks for these models are shown in Fig. 36. These calculations should be compared to Fig. 5 of Yoon & Langer (2005).  \n<!-- image -->  \nFigure 38. Left: A plot of the OPAL opacity data for five Z values at X = 0 : 7 and log( GLYPH<26>= T 6 3 ) = GLYPH<0> 5. These curves show the increase in the iron opacity bumps at log T GLYPH<25> 5 : 3 and 6.3 as Z increases from 10 GLYPH<0> 4 to 0.02. Right: The opacity-temperature profiles of 60 M GLYPH<12> ZAMS models for the same Z values.  \n<!-- image -->  \nEquation (34) corresponds to eq. (9) of Joss et al. (1973). As argued in that paper, entropy decreases as density increases; therefore a density inversion implies a superadiabatic gradient, and as a result, L onset < L inv. This can be shown explicitly for a chemically homogenous mixture of an ideal gas and radiation. For such a mixture, equation (32) becomes  \nL rad L Edd > L inv L Edd = " 1 GLYPH<0> P gas = P 1 GLYPH<0> 3 P gas = 4 P # ; (35)  \nand equation (34) becomes  \nL rad L Edd > L onset L Edd = 8(1 GLYPH<0> P gas = P )(4 GLYPH<0> 3 P gas = P ) 32 GLYPH<0> 24 P gas = P + 3( P gas = P ) 2 ; (36)  \nallowing one to show that L onset < L inv. At high luminosities where the gas becomes radiation-dominated, however, the di GLYPH<11> erence between L onset and L inv becomes small.  \nExpanding equations (35) and (36) for P gas = P GLYPH<28> 1 gives L inv GLYPH<0> L onset GLYPH<25> (3 = 4) GLYPH<2> ( P gas = P ) GLYPH<2> L Edd. For such high-luminosity, radiation-dominated stars, a small ine GLYPH<14> ciency in convection is su GLYPH<14> cient to drive a density inversion.  \nWe now demonstrate that such ine GLYPH<14> cient convection can arise in the convective, radiation-dominated, envelopes of massive stars. In order of magnitude the convective and radiative fluxes are, respectively, F conv GLYPH<24> GLYPH<26> c 3 s ( r T GLYPH<0> r ad) 3 = 2 and F rad GLYPH<24> cP rad =GLYPH<28> . To carry the flux, we need F conv GLYPH<24> F rad; equating and substituting GLYPH<26> c 2 s GLYPH<24> P GLYPH<24> P gas, we arrive at an expression that sets the level of superadiabaticity,  \n( r T GLYPH<0> r ad) 3 = 2 GLYPH<24> c c s P rad P gas GLYPH<28> GLYPH<0> 1 : (37)  \nUnder typical conditions in massive star envelopes, c = c s GLYPH<24> 10 4 at the iron opacity bump, but at this location, GLYPH<28> is not  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \nlarge enough to prevent the superadiabaticity from triggering a density inversion. ● ● ● ● ● ● ● ● ● ● ●  \nThe lines in Fig. 39 show these luminosity conditions as a function of P gas = P , and reveal that as the stellar conditions become radiation dominated, there is only a small gap between a convective model that is adiabatically stratified and a model with a density inversion. This corresponds to the region between the curves L rad = L onset (dot-dashed line) and L rad = L inv (dashed line). The gas pressure does not invert until L > L Edd, which in Fig. 39 is the region above the solid horizontal line. We show profiles from a 30 M GLYPH<12> (left panel) and a 70 M GLYPH<12> model (right panel). These are from the first crossing of the Hertzsprung gap when T e GLYPH<11> = 5000 K. Each dot corresponds to a zone in the calculation; as the profile moves outward from center to surface the traces go from bottom to top in the plot. The blue dots indicate zones where the star is radiative; red indicates convection; a black border denotes a density inversion, d GLYPH<26>= d r > 0; and yellow indicates a gas pressure inversion, d P gas = d r > 0. There is excellent agreement between the detailed MESA evolutionary calculations and the analytical conditions (eq. [35] and [36]). The 70 M GLYPH<12> profile goes into the low P gas = P , high L rad = L Edd regime. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \nFigure 40. The panels display, from top to bottom, the density, gas pressure, total pressure, and entropy as functions of radius for the 70 M GLYPH<12> model shown in Fig. 39. The range of radii is restricted to the region where density and gas pressure inversions develop. Each zone is marked by a dot; a small red dot indicates convection with no predicted gas or gas pressure inversion ( L rad < L inv); a large red dot with black border indicates a predicted density inversion but no gas pressure inversion ( L inv < L rad < L Edd); and a yellow dot with black border indicates a convective region with a predicted gas pressure inversion ( L rad > L Edd). The total pressure (panel c) is well-behaved at all radii. Note also the decrease in entropy (panel d): the region is superadiabatic.  \n<!-- image -->  \n●●●●  \n●●  \n●●●  \n●●●●●  \n●●●●●●●●●●●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●●●●●●●●●●●●●●●●●●●●●●●●●●●●  \n●  \n●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●  \n● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ●●● ●●●● ●●●● ●●● ●●● ●●●●●● ●●●● ●●● ●●●●● ●●●●● ●●●●●● ●●●●●● ●●●●●●● ●●●●●●●● ●●●●●●●●●●● ●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● Figure 39. The critical luminosities L rad = L onset (eq. [36], dot-dashed line), L rad = L inv (eq. [35], dashed line), and L rad = L Edd (solid line) as a function of P gas = P for an ideal gas-radiation mixture. Compare this with Fig. 1 of Joss et al. (1973). For L rad < L onset, the gas is convectively stable; for L onset < L rad < L inv, the gas is convective; for L inv < L rad < L Edd, the density is inverted, d GLYPH<26>= d r > 0; and for L Edd < L rad, the gas pressure is inverted, d P gas = d r > 0. Overlaid on the plots are the profiles from a 30 M GLYPH<12> (left panel) and a 70 M GLYPH<12> (right panel) model with Z = 0 : 02: blue dots indicate zones that are radiative; red dots indicate r rad > r ad; dots with a black border have a density inversion; and the yellow dots with black borders indicate a gas pressure inversion. As the profile moves out from the stellar center it traces out the points on the plot from bottom to top. Only a part of the model profiles are visible in the plot. The calculations correspond to the first crossing of the Hertzsprung gap when T e GLYPH<11> = 5000 K.  \nFigure 40 displays the physical conditions in the 70 M GLYPH<12> model where the density and gas pressure inversions develop. The panels display, from top to bottom, density, gas pressure, total pressure, and entropy, all as functions of radius. The total radius is R = 1330 R GLYPH<12> . Regions with r > r ad and L rad < L inv < L Edd are marked with a small red dot. Regions where L inv < L rad < L Edd (cf. eq. [35]) are marked with a large red dot with a black border. Regions where L rad > L Edd are marked with a large yellow dot with black border. Although the pressure (panel c) is well-behaved in this superadiabatic (panel d) region, a density inversion does develop  \n●●●●●●●●●●●●●●●●●●●●●●●●●●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \nwhere L Edd > L rad > L inv (panel a) and a gas pressure inversion develops (panel b) where L rad > L Edd, as predicted. In this region the superadiabaticity r T GLYPH<0> r ad > 10 GLYPH<0> 2 and is greater than unity for r = R GLYPH<12> & 1300. This is much larger than a typical value ( GLYPH<24> 10 GLYPH<0> 6 ) where convection is e GLYPH<14> cient and results in the entropy decreasing with r as shown in panel c. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n<!-- image -->  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●  \n●', '7.2. Treatment of Superadiabatic Convection in Radiation-Dominated Regions': "In MESA star the superadiabatic gradient arising in radiation-dominated envelopes can force the adoption of prohibitively short timesteps. Energy is mostly transported by radiation, and the convective velocities resulting from MLT approach the sound speed. The stability of such radiationdominated envelopes has been discussed in the past, and is still a matter of debate (see, e.g., Langer 1997; BisnovatyiKogan & Dorodnitsyn 1999; Maeder 2009; Su'arez-Madrigal et al. 2013). In this regime, the treatment of convective energy transport by MLT is admittedly out of its domain of applicability. Hydrodynamical instabilities and the transport of energy from waves excited by near-sonic turbulent convection are important for energy transport, and three-dimensional hydrodynamical calculations are required to capture fully the complex physics occurring in this regime.  \nHere we develop a treatment of convection, known as MLT ++ , that reduces the superadiabaticity in some  \nradiation-dominated convective regions. This treatment allows MESA star to calculate models of massive stars up to core collapse. For every model, MESA star computes the values of  \nGLYPH<21> max GLYPH<17> max L rad L Edd ! and GLYPH<12> min GLYPH<17> min P gas P ! : (38)  \nWhen GLYPH<12> min is small and GLYPH<21> max is large, and MLT yields a GLYPH<14> r > GLYPH<14> r ; thresh, we artificially decrease the superadiabaticity, GLYPH<14> r GLYPH<17> r T GLYPH<0> r ad, implied by MLT. The default of the userspecified parameter GLYPH<14> r ; thresh is su GLYPH<14> ciently large, GLYPH<24> 10 GLYPH<0> 3 , so that convection is still ine GLYPH<14> cient.  \nMESA star sets r T to reduce the GLYPH<14> r GLYPH<0> GLYPH<14> r ; thresh by a factor GLYPH<11> r f r , where f r is specified by the user, and GLYPH<11> r is updated at each timestep to a linear combination of its previous value and a value f GLYPH<11> r ( GLYPH<21> max ; GLYPH<12> min). For large values of GLYPH<21> max and small values of GLYPH<12> min, f GLYPH<11> r ! 1; in typical usage, the transition happens where GLYPH<21> max GLYPH<25> 0 : 5 and GLYPH<12> min GLYPH<25> 0 : 3. For small values of GLYPH<21> max and large values of GLYPH<12> min, f GLYPH<11> r ! 0. Thus f r sets the maximum reduction of GLYPH<14> r GLYPH<0> GLYPH<14> r ; thresh. Figure 41 shows how MESA star turns on the reduction in GLYPH<14> r as a star evolves. Tracks in the HR diagram are shown for four stellar models: 15, 25, 30, and 70 M GLYPH<12> . The color of each line indicates the value of GLYPH<11> r at each point.  \nFigure 41. HRdiagram of 15, 25, 30, and 70 M GLYPH<12> models. The color indicates the value of GLYPH<11> r at that point in the star's evolution. For the 25 M GLYPH<12> and 30 M GLYPH<12> stars, there is a sharp spike in GLYPH<11> r as the star crosses the Hertzsprung gap followed by a sharp drop at the base of the red giant branch. The 70 M GLYPH<12> model has GLYPH<11> r > 0 : 9 for its entire evolution.  \n<!-- image -->  \nSuch a decrease of the temperature gradient reduces L rad and implies additional physical transport. Potential agents for the excess transport include waves excited by turbulent convection (see, e.g., Maeder 1987a) and radiative di GLYPH<11> usion enhanced by porous clumping of the envelope (e.g., Owocki et al. 2004). As these radiation-dominated envelopes might be physically unstable, with a resulting strong enhancement of mass loss, we caution that the results of any 1D stellar evolution calculation for the late evolutionary phases of massive stars should be considered highly uncertain.  \nWe now show a comparison of MESA star calculations of rotating massive stars done with and without MLT ++ . We used the 25 M GLYPH<12> model described in x 6.6, which at Z = 0 : 02 is around the upper mass limit that can converge using a reasonably short timestep without having to rely on the MLT ++ . The most prominent di GLYPH<11> erence between the calculations is the evolutionary track in the HR-diagram (Fig. 42). This is not surprising, as MLT and MLT ++ result in di GLYPH<11> erent e GLYPH<14> ciencies of energy transport in radiation-dominated stellar envelopes. The sharp drop in L for the MLT ++ case is the result of a brief period of enhanced mass loss due to super-critical rotation. The structure and the angular momentum content of the collapsing core are weakly dependent, however, on the choice of MLT vs. MLT ++ (Fig. 43).", '7.3. Core-Collapse Progenitor Models': 'We evolve a grid of massive stars initially rotating with GLYPH<10> = GLYPH<10> crit = 0 : 2. The models have been initialized using solid body rotation. Models with initial M = M GLYPH<12> = 30, 40, 50, 60, 70, 80, 90 and 100 have initial Z = 0 : 02, while models with initial M = M GLYPH<12> = 120, 150, 250, 500 and 1000 have been initialized with Z = 0 : 001. To calculate convective boundaries we adopt the Ledoux criterion including the impact of semiconvection (with GLYPH<11> sc = 0 : 02, see x 4.1). The transport of angular momentum and chemicals by rotational instabilities and magnetic torques is included and calibrated following Heger et al. (2000, 2005) and Yoon & Langer (2005). Wind massloss is been implemented following the recipe of Glebbeek et al. (2009). For T e GLYPH<11> > 10 4 K and H-surface fraction > 0 : 4, the mass-loss prescription of Vink et al. (2001) is used. In the same temperature range, but when the H-surface fraction decreases below 0.4, Nugis & Lamers (2000) determine the mass-loss rate. At low temperatures ( T e GLYPH<11> < 10 4 K) the massloss rate of de Jager et al. (1988) is used.  \nFigure 44 shows the central conditions of these massive rotating models. For each model the calculation stops when any part of the collapsing core reaches an infall velocity of 1000 km s GLYPH<0> 1 . Some of the initial and final properties are summarized in Table 5. These calculations are performed to reveal the new capabilities of MESA star . The values of the parameters for these calculations have not been calibrated against existing calculations or observations.', '7.4. Radial Instability of Red Supergiants': 'Massive red supergiants (RSG) are unstable to radial pulsations driven by the GLYPH<20> -mechanism in the hydrogen ionization zone. Both linear and non-linear calculations show the occurrence of oscillations with the period and growth rate of the dominant fundamental mode increasing with L = M (Li & Gong 1994; Heger et al. 1997; Yoon & Cantiello 2010). The periods are of the order of years. As discussed by Yoon & Cantiello (2010) the occurrence of RSG pulsations can impact stellar mass-loss rates and modify the evolution of massive stars above a certain mass. We study the occurrence of RSG pulsations with MESA star and compare results with existing non-linear calculations.  \nIn Fig. 45 we show the capability of MESA star to exhibit radial oscillations in luminous RSGs. We use the same 25 M GLYPH<12> rotating model discussed in x 6.6, and we restart the calculation when the He mass fraction in the core is Yc = 0 : 7. For non-rotating RSG with Z = 0 : 02, Yoon & Cantiello (2010) found pulsation periods in the range 1-8 yr. To resolve the RSG pulsations we force the timestep to < 0 : 01 yr,  \n<!-- image -->  \nFigure 42. Comparison of evolutionary tracks (left) and equatorial rotational velocity (right) for a 25 M GLYPH<12> model with Z = 0 : 02 and 3 eq = 200 km s GLYPH<0> 1 . The solid black lines show MESA star results with MLT (black) and MLT ++ (orange), while the dashed blue line refers to STERN calculations. The star symbol shows the location where we started the calculation for the RSG pulsations discussed in x 7.4.  \n<!-- image -->  \nTable 5 Initial (ZAMS) and final (core-collapse) properties.  \n- a Initial rotation rate, see definition in x . 6.4.\n- b Initial equatorial rotational velocity.\n- c Total initial angular momentum.\n- d Stellar lifetime.\n- e Main sequence and core He-burning lifetimes. These are defined as the interval between onset of core burning and depletion of central hydrogen (or helium) to 1% by mass.  \n- g Mass of the Iron core (if present).\n- h Final total angular momentum.  \ni Final total angular momentum of the iron-core.  \nmuch shorter than the usual timestep during He-burning ( GLYPH<14> t & 10 2 yr, see Appendix B.3). This explains why RSG pulsations are usually not found during the evolution of massive stars. Before the code stops due to the emergence of supersonic radial velocities in the envelope, we find a pulsational period GLYPH<25> 4 yr, in good agreement with the results of Yoon & Cantiello (2010).', '8. SUMMARY AND CONCLUSIONS': "We have explained and, where possible, verified the improvements and major new capabilities implemented in MESA since the publication of Paper I. These advancements include evolutionary modeling for giant planets ( x 2), tools for aster-  \noseismology ( x 3), implementation of composition gradients in stellar interiors and their impact on convective mixing ( x 4) , the evolution of intermediate mass stars and white dwarfs ( x 5) the treatment of rotation during stellar evolution ( x 6), addressing the onset of radiation pressure dominance in the envelopes of evolving massive stars due to the iron opacity bump, and evolving massive stars to the onset of core collapse ( x 7). The enhancements include the physics modules (Appendix A), the algorithms (Appendix B), and the addition of a MESA Software Development Kit (Appendix C). MESA star input files and related materials for all the figures are avaliable at http://mesastar.org .  \nThese hitherto unpublished advancements have already en-  \nFigure 43. Evolution of internal specific angular momentum for the two 25 M GLYPH<12> models of Fig. 42. The dashed lines show models calculated with MLT ++ . Due to di GLYPH<11> erent excursions in the HRD (see Fig. 28) calculations with MLT and MLT ++ end with di GLYPH<11> erent final masses. There are no substantial changes, however, in the specific angular momentum content of the stellar cores.  \n<!-- image -->  \n10.0  \nFigure 45. Surface properties of a pulsating RSG. This is the same 25M GLYPH<12> model discussed in Sec. 6.6, evolved from t 0 = 6 : 851 Myr (corresponding to Yc = 0 : 7, star symbol in Fig. 32) with timesteps GLYPH<14> t GLYPH<20> 0 : 01 yr. The black line shows the evolution of the stellar radius, while the orange line shows the value of the surface radial velocity (in units of the local sound speed). The inset shows the corresponding evolution in the HR-diagram.  \n<!-- image -->  \nc  \nFigure 44. Evolution of T c and GLYPH<26> c in the massive rotating models. The locations of core helium, carbon, neon, oxygen, and silicon burning are labeled. A dashed curve marks the electron-positron pair-instability region where GLYPH<0> 1 < 4 = 3. All models are rotating initially at 20% of critical rotation. The calculations include the e GLYPH<11> ects of rotation and Spruit-Tayler magnetic fields as discussed in x 6. Models with initial mass GLYPH<20> 100 M GLYPH<12> have initial metallicity Z = 0 : 02, while models with mass GLYPH<21> 120 M GLYPH<12> have initial metallicity Z = 0 : 001. The end of the line for each mass corresponds to the time of core-collapse, defined as when any part of the collapsing-core exceeds an in-fall velocity of 1000 km s GLYPH<0> 1 . The tracks for the 60 M GLYPH<12> and 70 M GLYPH<12> overlap in this plot.  \n/R  \nR  \n<!-- image -->  \nMach #  \nabled a number of studies in planets (e.g., Passy et al. 2012; Huang & Cumming 2012; Carlberg et al. 2012), classical novae (Denissenkov et al. 2013), asteroseismology (e.g., Yang et al. 2012; Burkart et al. 2012b; Moravveji et al. 2012), rotationally induced mixing (e.g., Denissenkov 2010; Chatzopoulos et al. 2012; Chatzopoulos & Wheeler 2012) and enabled the discovery of new features in the evolution of low-mass stars (Denissenkov 2012). In addition, these enhanced capabilities have allowed for applications of MESA star that were not initially envisioned, such as explorations of stars under modified gravity (Chang & Hui 2011; Davis et al. 2012), and stellar oscillations induced by tidal disturbances in double white dwarf binaries (Fuller & Lai 2012a,b; Burkart et al. 2012a).  \nAs an open source 'instrument' for stellar astrophysics, it is di GLYPH<14> cult to predict all the ways in which future development of MESA will occur. We do know, however, that future versions of MESA will include advances in physics modules, features driven by the MESA user community, and architectural refinements. For example, the plethora of asteroseismological data is driving new initiatives to incorporate nonadiabatic pulsation codes, where possible, into MESA . The prevalence of interacting binary star systems, especially for massive stars, has increased the pressure for MESA development e GLYPH<11> orts that would yield the capability to simultaneously evolve two interacting stellar models. Physics module developments will likely include general relativistic corrections to the stellar structure equations (e.g., di GLYPH<11> erence between gravitational and baryonic mass), the mass di GLYPH<11> usion coe GLYPH<14> cients in electron degenerate environments, phase separation in cooling white dwarfs, and nuclear statistical equilibrium solvers. We also expect the transition from multicore systems (with order 10 cores) to many-core architectures (with order 100 cores) to drive new directions in MESA 's algorithmic and architectural  \ndevelopment.  \nIt is a pleasure to thank Falk Herwig for significant contributions to the MESA project and Jørgen Christensen-Dalsgaard for kindly providing the ADIPLS code for inclusion in MESA and assisting with its integration. We likewise thank Jared Brooks for documenting the MESA test suite, Tristan Guillot for providing the CEPAM evolutionary tracks, Alexander Potekhin for an update to his EOS code, Didier Saumon and Jim MacDonald for providing EOS tables, Ren'e Rohrman for providing atmosphere tables, Richard Freedman for providing opacity tables, Evert Glebbeek, Alex Heger, and Norbert Langer for providing code for implementing rotation, and Haili Hu for providing code for implementing di GLYPH<11> usion.  \nWe also thank David Arnett, Tim Bedding, Kent Budge, Phil Chang, Pieter DeGroote, Pavel Denisenkov, Jonathan Fortney, Chris Fryer, Gustavo Hime, Raphael Hirschi, Sam Jones, Steve Kawaler, Phillip Macias, Pablo Marchant, Travis Metcalfe, Kevin Moore, Ehsan Moravveji, Jean-Claude Passy, Hideyuki Saio, Josiah Schwab, Aldo Serenelli, Josh Shiode, Steinn Sigurdsson, Anne Thoul, Roni Waldman, Achim Weiss, Stan Woosley, Sung-Chul Yoon, and Patrick Young for providing valuable discussions and correspondence. Some of the simulations for this work were made possible by the Triton Resource, a high performance research computing system operated by San Diego Supercomputer Center at UC San Diego.  \nWe thank the participants of the 2012 MESA Summer School for their willingness to experiment with the new capabilities: Je GLYPH<11> Andrews, Umberto Battino, Keaton Bell, Harshal Bhadkamkar, Kristen Boydstun, Emmanouil Chatzopou-  \nlos, Eugene Chen, Jieun Choi, Alex Deibel, Luc Dessart, Ian Dobbs-Dixon, Tassos Fragos, Samuel Harrold, Daniel Huber, Joe Hughto, Max Katz, Agnes Kim, Io Kleiser, Shri Kulkarni, Gongjie Li, Christopher Lindner, Jing Luan, Mia Lundkvist, Morgan MacLeod, Jo˜ao Marques, Grant Newsham, Rachel Olson, Richard O'Shaughnessy, Kuo-Chuan Pan, Ilka Petermann, Theodore Sande, Ken Shen, Natalia Shabaltas, Dave Spiegel, Jie Su, Tuguldur Sukhbold, David Tsang, Bill Wolf, Angie Wolfgang, Tsing Wai Wong, and Alexey Zinger.  \nThis project was broadly supported by the NSF under grants PHY 11-25915 and AST 11-09174. M.C. acknowledges partial support from the 'Alberto Barlettani' Prize 2012. P.A. acknowledges support by NSF AST-0908873 and NASA NNX09AF98G. L.B. acknowledges support from the Wayne Rosing, Simon and Diana Raab Chair in Theoretical Astrophysics at KITP. E.F.B acknowledges support by the Joint Institute for Nuclear Astrophysics under NSF PHY grant 08-22648. A.L.D received support from the Australian Research Council under grant FL110100012. M.H.M acknowledges support by the NSF grant AST-0909107, the NASA grant NNX12AC96G, the Norman Hackerman Advanced Research Program under grant 003658-0252-2009, and the Delaware Asteroseismic Research Center. D.S. acknowledges support by the Australian Research Council. F.X.T acknowledges support from the NSF under grants AST 08-06720, AST 6736821, AST 09-07919, AST 10-07977, PHY 08-22648, and from NASA under grants 08-NAI5-0018 and NNX11AD31G. R.T. acknowledges support by NSF grants AST-0908688 and AST-0904607, and NASA grant NNX12AC72G.", 'UPDATES TO INPUT PHYSICS MODULES': 'There have been many updates and improvements to the physics modules since Paper I. In this appendix, we describe the changes that have been made to the microphysics modules chem ( x A.1), eos ( x A.2), kap ( x A.3), and net ( x A.4). We conclude by listing updates to the atmosphere boundary conditions ( x A.5).', 'Atomic and Nuclear Data': 'The chem module now has the latest version (v2.0) of the JINA reaclib nuclide data (Cyburt et al. 2010). This contains updated mass evaluations, and now includes 7853 nuclides up to 337 Cn. For precision work, the chem module now distinguishes between the atomic mass number Ai -the number of nucleons in a given isotope-and the atomic mass Wi . The abundance of a species i is defined as  \nYi GLYPH<17> ni n B ; (A1)  \nwhere n B is the baryonic number density. The baryon fraction Xi is then  \nXi = YiAi = niAi n B ; (A2)  \nNote that P i Xi = n B = n B = 1 and is invariant under nuclear reactions. We then define the baryon density (in mass units) as  \nGLYPH<26> = n B m u ; (A3)  \nwhere m u = 1 : 660538782 GLYPH<2> 10 GLYPH<0> 24 g is the atomic mass unit (CODATA 2006 value; Mohr et al. 2008). Note that the numerical value m u, along with other physics constants, are defined in the const module. The atomic mass of isotope i is defined in MESA as  \nWi = Ai + GLYPH<1> i m u c 2 ; (A4)  \nwhere GLYPH<1> i = c 2 is the mass excess of isotope i . This treatment neglects the electronic binding energy, and GLYPH<1> is therefore independent of the ionization state of a given species. The electron rest masses are, however, included in this definition, since the Wi are atomic masses.  \nThe MESA microphysics moduleskap , eos , neu , and net -use GLYPH<26> , T , and f Xi g as inputs. MESA star multiplies GLYPH<26> by a mass correction factor ¯ W = ¯ A = P i WiYi = P i AiYi to distinguishes between Ai and Wi before starting the calculation for a timestep. A  \ncall to the routine composition info in the chem module returns the following averaged quantities: the mean atomic mass number, ¯ A GLYPH<17> P i YiAi = P i Yi , mean atomic charge number, ¯ Z GLYPH<17> P i ZiYi = P i Yi , mean square atomic charge number, P i Z 2 i Yi = P i Yi , the electron abundance, Ye = ¯ Z = ¯ A , and the mass correction term, ¯ W = ¯ A . In addition, the routine returns the derivatives of ¯ A , ¯ Z , and ¯ W = ¯ A with respect to the baryon fractions Xi :  \n@ ¯ A @ Xi GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<26>; Xj , i = ¯ A Ai GLYPH<16> Ai GLYPH<0> ¯ A GLYPH<17> 1 P i Xi ; (A5)  \n@ ¯ Z @ Xi GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<26>; Xj , i = ¯ A Ai GLYPH<16> Zi GLYPH<0> ¯ Z GLYPH<17> 1 P i Xi ; (A6)  \n@ ( ¯ W = ¯ A ) @ Xi GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<26>; Xj , i = Wi Ai GLYPH<0> ¯ W ¯ A ! 1 P i Xi : (A7)  \nNote that the routine does not make any assumption in these derivatives that P i Xi GLYPH<17> P i AiYi = 1; in this formulation, P i Xi is not explicitly set to unity.  \nAt the beginning of each Newton iteration, the abundances are checked. A mass fraction is considered good if its value exceeds min xa hard limit . If all mass fractions meet this standard, then the mass fractions are clipped to range from 0 to 1, and the mass fractions are summed. If the sum di GLYPH<11> ers from unity by less than a value sum xa tolerance , then the mass fractions are renormalized to sum to unity; otherwise, the code reports an error. Currently composition derivatives are ignored in the eos and kap routines. Equations (A5)-(A7) allow, however, future additions to these routines to compute these derivatives analytically.', 'Equation of State': "The only significant change to the eos module since Paper I is the addition of tables for Z > 0 : 04, where Z is the mass fraction of all elements heavier than He. The eos module as described in Paper I supplied equation of state (EOS) tables for Z = 0 : 0, 0.02, and 0.04 at temperatures and densities for which neutral and partially-ionized species are present (see Paper I, Figure 1). For Z > 0 : 04 MESA switched to the HELM EOS (Timmes & Swesty 2000), which assumes full ionization. In order to rectify the inconsistent treatment of the partially-ionized region at high Z , new EOS tables have been computed (J. MacDonald, priv. comm.) using the MacDonald EOS code (MacDonald & Mullan 2012) for Z = 0 : 2 (scaled-solar), and two Z = 1 : 0 compositions: one with 49.5% C, 49.5% O, and 1% scaled-solar by mass; and one with 50% C and 50% O by mass. Here 'scaled-solar' refers to the Grevesse & Noels (1993) solar heavy element distribution adopted in the OPAL EOS tables (Rogers & Nayfonov 2002).", 'Opacities': 'The kap module now divides the opacity tables into a high-temperature domain, log( T = K) & 4, and a low-temperature domain, log( T = K) . 4; the exact range of log T over which the tables are blended can be adjusted at runtime. This treatment di GLYPH<11> ers from the opacity tables described in Paper I, which combined high- and low-temperature opacities into a single set of tables. The motivation for separating the tables is to facilitate using di GLYPH<11> erent sources of lowT opacity data. The kap module now supports lowT opacities from either Ferguson et al. (2005) or Freedman et al. (2008) with updates to the molecular hydrogen pressureinduced opacity (Frommhold et al. 2010) and the ammonia opacity (Yurchenko et al. 2011). Either set may be selected at run time. The electron conduction opacity tables, based on Cassisi et al. (2007), have been expanded (Potekhin 2011, priv. comm.) to cover higher temperatures (up to 10 10 K, originally 10 9 K) and densities (up to 10 11 : 5 g cm GLYPH<0> 3 , originally 10 9 : 75 g cm GLYPH<0> 3 ).', 'Nuclear Reactions': "A reaction network is defined by a set of isotopes and a set of reactions; these sets are specified in a reaction network definition file. MESA comes with many predefined reaction networks in data/net data/nets and can also incorporate user-defined networks. To use a custom network, a user creates a reaction network definition file containing the command  \nadd isos and reactions(isos list) , which will automatically add all reactions linking the isotopes in isos list . The sequence of isotopes in isos list may be specified by the name of the isotope: for example, add isos and reactions(he4) adds 4 He. Alternatively, one can specify the name of element followed by the desired minimum and maximum nucleon number. For example, the command add isos and reactions(o 16 18) adds 16 O, 17 O, and 18 O. Note that because many of the predefined networks may use e GLYPH<11> ective rates-that is, using one reaction rate to represent a reaction sequence or group of reaction sequences-it is not recommended that the user extend one of the pre-existing networks with this command.  \nMESA creates and stores reaction rate tables for each reaction whose entries are derived from evaluating standard analytic fitting formulas (see x A.4), but these reaction rates may be replaced with user-specified values. To change a rate for a given reaction,  \n- 1. create a file with two columns: the temperature in units of 10 8 K and the rate N A h GLYPH<27> v i in units of cm 3 g GLYPH<0> 1 s GLYPH<0> 1 ;\n- 2. list the file name in a local file rate list.txt along with its 'handle' for the reaction rate in question (see discussion below); and\n- 3. set the parameter rate tables dir in the namelist star job to the name of the directory in which rate list.txt is located; by default this is data/rates/rate tables .  \nThe handle for a reaction is derived from the input and output channel isotopes according to a few rules. Capture reactions, such as x ( p ; GLYPH<13> ) y , have handles of the form r x pg y and exchange reactions, such as x ( GLYPH<11>; p ) y , have handles of the form r x ap y . Other arbitrary reactions may be added by listing them in a form r inputs to outputs where inputs and outputs are isotopes separated by ' '. If the same isotope appears two or more times, the isotope name may be repeated. For example, the tripleGLYPH<11> reaction is specified as r he4 he4 he4 to c12 . Isotopes are ordered by increasing Z and N , e.g., r h3 be7 to neut h1 he4 he4 . To see a list of reactions used, the parameter show net reactions info in namelist star job should be set to .true. .", 'Atmosphere Boundary Conditions': 'The atm module provides the surface boundary condition for the interior model. A collection of four new options that extend the set described in Paper I are now available.', '1. solar Hopf grey : Implements the solar-calibrated Hopf T ( GLYPH<28> ) relation, where': "T 4 ( GLYPH<28> ) = 3 4 T 4 e GLYPH<11> GLYPH<2> GLYPH<28> + q ( GLYPH<28> ) GLYPH<3> ; (A8)  \nand  \nq ( GLYPH<28> ) = q 1 + q 2 exp( GLYPH<0> q 3 GLYPH<28> ) + q 4 exp( GLYPH<0> q 5 GLYPH<28> ) : (A9)  \nThe qi are fit to the solar atmosphere with resulting values q 1 = 1 : 0361, q 2 = GLYPH<0> 0 : 3134, q 3 = 2 : 448, q 4 = GLYPH<0> 0 : 2959, and q 5 = 30 : 0 (J. Christensen-Dalsgaard, 2011, priv. comm.).  \n- 2. grey and kap : Expands on the simple assumption that P ' GLYPH<28> g =GLYPH<20> by iterating to find a consistent solution among P , T , and GLYPH<20> ( GLYPH<26>; T ).\n- 3. grey irradiated : Implements the Guillot (2010) T ( GLYPH<28> ) relation that includes both external irradiation by the star and cooling flux from the interior; see Guillot (2010, eq. 49) along with the discussion and results in x 2.4. In addition to the external and internal fluxes, this boundary condition requires two constant opacity values: GLYPH<20> v for the external radiation, and GLYPH<20> th for the thermal radiation generated within the atmosphere. This boundary condition is unique in that it is applied at a specified pressure level, as opposed to optical depth. This pressure must be chosen su GLYPH<14> ciently high to capture any heating of the atmosphere by the irradiation.  \n- 4. WD tau 25 tables : Provides as outer boundary conditions the values of P gas and T at log( GLYPH<28> ) = 1 : 4 as extracted from pure hydrogen model atmospheres of WDs (Rohrmann et al. 2012; Rohrmann 2001). The tables span a range of e GLYPH<11> ective temperatures and surface gravities: 2 ; 000 K GLYPH<20> T e GLYPH<11> GLYPH<20> 40 ; 000 K and 5 : 5 GLYPH<20> log g GLYPH<20> 9 : 5. See x 5.2 for an example of the use of these tables.", 'NUTS AND BOLTS': 'We now briefly describe the primary components of evolution calculations. MESA star first reads the input files and initializes the physics modules to create a nuclear reaction network and access the EOS and opacity data. The specified starting model is then loaded into memory and the evolution loop is entered.', 'Evolve a Step': "The top level routine for evolving a star for a single timestep is do evolve step . If this is the first attempt to do a step starting from the current state, the model is remeshed (see x B.4), and information for MLT ++ is prepared by the routine set gradT excess alpha (see x 7.2). Su GLYPH<14> cient information is saved so that if necessary it will be possible to make other attempts (i.e., after a redo, a retry, or a backup). In addition to the current state, we keep the previous state (called 'old'), and the one that came before 'old' (called 'older'). During the step, the current state is modified, and the old one holds the state at the start of the step. If we do a redo or a retry, we copy old to current to restore the starting state. If we do a backup, we copy older to old before copying old to current, making us start at the state prior to the current one. Note that the duration of the timestep is determined before the call on do evolve step by the process described in x B.3.  \nAfter remeshing and the other initial preparations, do evolve step begins the operations that are done on every attempt. It first calls the routine do winds which sets ˙ M based on the current radius, luminosity, mass, metallicity, and other properties as needed. During the evalution of do winds there is a call on the user-defined other wind routine giving users an easy way to define di GLYPH<11> erent schemes for setting ˙ M .  \nInformation for evaluating the Lagrangian time derivatives is stored by a call to save for d dt . The ensuing call to do adjust mass adds or removes mass without changing the number of grid points (see x B.5). Information for evaluating the Lagrangian time derivatives is updated at this point. Variables for the model are evaluated to reflect the changes made by remeshing and changing mass. This includes evalution of the Brunt-Vaisala frequency (see x 3.3), and the di GLYPH<11> usion coe GLYPH<14> cients for the mixing of composition (see x 4.1 and 4.2). The user-defined routines other brunt , other mlt , and other mixing are called as part of this. If rotation is enabled, there is a call to set rotation mixing info (see x 6) which in turn calls other am mixing . If element di GLYPH<11> usion from gravitational settling and chemical di GLYPH<11> usion is active, the routine do element diffusion adjusts the composition and includes a call on other diffusion The ensuing call to do struct burn mix solves for the new structure and composition of the star through repeated Newton iterations (see x B.2). Non-convergence causes do evolve step to return with a result indicating a failure. Convergence is followed by a call to the routine do solve omega mix which adjusts the total angular momentum by solving a di GLYPH<11> usion equation (see x B.6); it calls other torque . There is an option to repeat the operations described in this paragraph in case rotationally enhanced mass loss (see x 6.4) has not been su GLYPH<14> cient to eliminate super-critical surface velocities. In such a situation, the mass loss is adjusted iteratively until slightly sub-critical velocities result. In e GLYPH<11> ect, this is an implicit solution for the appropriate ˙ M when super-critical rotation occurs.  \nNext, if specified by the user, smooth convective bdy is called to smooth abundances behind retreating convection boundaries. Finally, a call to do report gathers information and metrics about the timestep for the user. This information will then be available to the user's extras check model routine.", 'Solving the Coupled Structure, Burn, and Mix Equations': "A call to do struct burn mix invokes a Newton method-an N-dimensional root find-to solve a system of N nonlinear di GLYPH<11> erential-algebraic equations for the new structure and composition of the stellar model. Here N is the number of zones in the current model times the number of basic variables per zone and can exceed 100,000.  \nThe equations to be solved are written as the relation F ( basic vars ) = 0, where F is the vector-valued function of the residuals. If the basic vars were a perfect solution to the equations, we would have F = 0; in practice, the solution is never perfect. The solution strategy is to iteratively adjust the values of the basic vars to reduce F towards zero. An approximate solution is accepted depending on both the magnitude of F and the relative size of the adjustments to basic vars . Adjustments are chosen using the Jacobian matrix of partial derivatives of all the F equations with respect to all the basic vars .  \nFigure 47 shows the three blocks making up the row of the block tridiagonal Jacobian matrix for the tenth from the center cell of a non-rotating 2 : 5 M GLYPH<12> ZAMSmodel, with black dots showing non-zero entries. The partial derivatives of the equations for cell k form the rows of the blocks. In this case, we have 4 equations for the structure of the model ( P , T , L , and r ) and 8 equations for the chemical abundances ( 1 H through 24 Mg). Each block of the tridiagonal matrix is demarcated by dashed black vertical lines. The block matrix on the left shows the dependencies of the equations for cell k on the variables of cell k GLYPH<0> 1, the one in the middle shows the dependencies of the equations for cell k on the variables of cell k , and the one on the right shows the dependencies of the equations for cell k on the variables of cell k + 1. The dashed lines partition each each block into four sub-blocks to highlight the structure and abundance portions of each block.  \nThe structure of the lower-right subblocks in the left and right blocks shows that the chemical abundance of a particular species in cell k depends on the chemical abundances of that species in cells k GLYPH<0> 1 and k + 1; this is because of mixing between neighboring cells. In this specific case of a non-rotating 2 : 5 M GLYPH<12> ZAMS model the mixing of chemical elements between cells is only due to the treatment of convection. The lower-right subblock of the center block also shows the interdependencies of abundances due to nuclear reactions in the cell. The bottom-left subblocks are zero in the left and right blocks but show dependencies on the P  \nand T variables of the center block. This is because the nuclear reactions that change the abundances depend of P and T of that cell but do not depend on P and T in the neighboring cells. The columns for L and r are zero in the center lower-left subblock because the equations for the abundances do not directly depend on those variables. The upper-right subblocks are zero in the left and right blocks but show that the equation for L depends on the abundance variables in the center block. This is because the L equation includes results from nuclear burning, and that depends directly on the composition of cell k but not on the composition in neighboring cells. The other rows in the center upper-right subblock are zero because the equations for P , T , and r do not directly depend on composition.  \nFinally, consider the upper-left subblocks that show the dependencies of structure equations on structure variables. The upperleft subblock in the center shows that each structure equation in k depends on 3 or 4 of the structure variables in k , The P and T equations for cell k also depend on both of the variables P and T in k GLYPH<0> 1. while the L and r equations for k depend on the corresponding variables in k GLYPH<0> 1. This pattern reflects the form of the finite di GLYPH<11> erences in the implementation of the structure equations: P and T di GLYPH<11> erences use the outer neighbor ( k GLYPH<0> 1) while L and r di GLYPH<11> erences use the inner neighbor ( k + 1). The L and r equations for innermost cell k = n use L center and R center; the P and T equations for the outermost cell k = 1 use the surface boundary conditions.  \nFigure 47. One row of the block tridiagonal Jacobian matrix for a 2 : 5 M GLYPH<12> ZAMS model, with black dots showing the locations of non-zero entries.  \n<!-- image -->  \nThe structure variables for each zone always include the zone average of the natural logarithm of the temperature, ln T , the luminosity at the outer edge of the zone, L , the natural logarithm of the radius at the outer edge of the zone, ln r , and a second thermodynamic variable-either the zone average of the natural logarithm of the mass density, ln GLYPH<26> , or the the zone average of the natural logarithm of the gas pressure, ln P gas. Ideally it would not matter whether ln GLYPH<26> or ln P gas was used as the second thermodynamic variable-for a given temperature and composition the equation of state permits going back and forth between the two. Microphysics packages tend to use mass density as a primary input (i.e., they use a Helmholtz free energy basis) leading to the common choice of ln GLYPH<26> . However, the structure equations are solved only to within a finite but non-zero residual (see above). Approximately correct values for the density and temperature can then lead to anomalous pressure profiles, with tiny violations of hydrostatic balance. These local violations tend to appear near large jumps in density, such as at a sharp H / He boundary. Using P gas as the second thermodynamic variable (e GLYPH<11> ectively using a Gibbs free energy basis) removes these anomalous pressure profiles. For example, in stellar models without overshooting or semiconvection, the H / He boundary is extremely sharp. Using the gas pressure as the second thermodynamic variable results in single zone step function transitions in the abundances and in the density, while the temperature and pressure are smooth across the transition. Applications that demand smooth pressure profiles, such as pulsation analysis (see x 3), should generally specify the gas pressure as the second thermodynamic variable.  \nMESA star treats convective mixing as a time-dependent, di GLYPH<11> usive process with a di GLYPH<11> usion coe GLYPH<14> cient, D , determined by the MLT module. See Paper I for the implementaton details of standard mixing length treatment. In addition to this standard MLT treatment, the MLT module includes the option to use the modified MLT ++ prescription described in x 7.1 and x 7.2. After the convective mixing calculations calculations have been performed, MESA star calculates the overshoot mixing di GLYPH<11> usion coe GLYPH<14> cient as described in Paper I. During the solution of the coupled structure, burning, and mixing equations the equation for mass fraction Xi ; k of species i in cell k is determined by  \nXi ; k ( t + GLYPH<14> t ) GLYPH<0> Xi ; k ( t ) = dX burn + dX mix = dXi ; k dt GLYPH<14> t + ( Fi ; k + 1 GLYPH<0> Fi ; k ) GLYPH<14> t dmk ; (B1)  \nwhere dXi ; k = dt is the rate of change from nuclear reactions, Fi ; k is the mass of species i flowing across face k  \nFi ; k = GLYPH<0> Xi ; k GLYPH<0> Xi ; k GLYPH<0> 1 GLYPH<1> GLYPH<27> k dmk ; (B2)  \nwhere GLYPH<27> k is the Lagrangian di GLYPH<11> usion coe GLYPH<14> cient from the combined e GLYPH<11> ects of convection and overshoot mixing and dmk = 0 : 5( dmk GLYPH<0> 1 + dmk ). For numerical stability, GLYPH<27> k is calculated at the beginning of the timestep and held constant during the implicit solver iterations. This assumption accommodates the non-local overshooting algorithm and significantly improves the numerical convergence. The structure of the lower-right subblocks in the left and right blocks in Figure 47 shows that the dependency of the chemcial abundances in cell k depends on the chemical abundances of that species in cells k GLYPH<0> 1 and k + 1 as a result of convective mixing.  \nIf the optional hydrodynamic mode is activated, then the radial velocity at the outer edge of the cell, v , is added to the structure variables. Figure 47 shows the order of the model variables in the Jacobian: each cell includes the structure variables followed by the mass fraction Xi of each isotope. Mass and the local angular velocity ! are not treated as structure variables because they are held constant during the Newton iterations. The mass is set before the iterations, while ! is determined after convergence. This is computed taking into account loss / gain of angular momentum during the time step, the new stellar structure and internal transport of angular momentum calculated by a di GLYPH<11> usion equation (see x B.6).  \nThe program flow to solve the coupled structure, burning, and mixing equations is to first create the matrix of partial derivatives using the current candidate solution, solve the block tridiagonal system of linear equations for the corrections to the basic variables, apply the possibly damped corrections (see next paragraph) to update the candidate solution, and calculate the residual F . Then, if the residual is small enough, we declare victory, otherwise we repeat the general flow with the updated candidate solution.  \nEach iteration of the Newton solver uses a linear approximation to create a vector of corrections to the model. These corrections do not include the physical requirement that the abundance mass fractions need to remain positive. To reduce the possible occurrence of negative abundances MESA star now uses a damped Newton scheme. This checks for proposed corrections that would produce negative abundances and multiplies the entire correction vector by a factor less than one, so that only part of the the full correction is applied. In many cases, this is su GLYPH<14> cient to significantly improve the convergence properties of a model. In other cases, the damped correction scheme may force so many small corrections that the Newton solver cannot converge within the user-specified maximum number of iterations, forcing the previous model to be attempted again with a smaller timestep (termed 'a backup'). On balance, this is usually a small price to pay for an improved conservation of species and more accurate solutions.  \nThe modules in star provide routines to evaluate the residual equations and create the Jacobian matrix. Given a candidate solution (i.e., the set of basic variables for each cell), the microphysics for each cell (EOS, thermal neutrino loss, opacity, nuclear reaction rates) are evaluated in parallel (see x B.9). The Jacobian matrix is then further populated with elements from rotation, artificial viscosity, and mixing length theory for the temperature gradient, and these are also evaluated in parallel. Each of the routines that evaluate these components returns output values and partial derivatives of the output values with respect to the input values. Analytic partial derivaties are used whenever feasible, otherwise numerical partial derivatives are supplied.", 'Timestep Controls': 'Control of the timestep is a critical part of stellar evolution and requires careful trade-o GLYPH<11> s. The timesteps must be small enough to allow convergence in comparatively few iterations but large enough to allow su GLYPH<14> ciently e GLYPH<14> cient evolutions. Changes to the timestep must respond rapidly to varying structure or composition conditions, but they need to be controlled to avoid large jumps that can reduce the convergence rate or the accuracy of the results. The routine pick next timestep performs timestep control as a two-stage process. The first stage proposes a new timestep using the H211B low-pass filter (Soderlind & Wang 2006), a scheme based on digital control theory. The second stage implements a wide range of tests that can reduce the proposed timestep if certain selected properties of the model are changing too much in a single timestep.  \nFor the first stage, routine hydro timestep sets the variable for the next timestep, dt next , according to the relative magnitude of changes to the basic vars . The variable reflecting the size of these changes is called varcontrol and is calculated by the routine eval varcontrol . For improved stability and response, the low-pass controller uses previous and current values of varcontrol to make the next timestep match the varcontrol target , wt , which is 10 GLYPH<0> 4 by default. To make this explicit, let GLYPH<14> ti GLYPH<0> 1, GLYPH<14> ti , and GLYPH<14> ti + 1 be the previous, current, and next timestep, respectively, while wc ; i GLYPH<0> 1 and wc ; i are the previous and current values of varcontrol . The maximum timestep for model i + 1 is then determined by  \nGLYPH<14> ti + 1 = GLYPH<14> ti f " f ( wt = wc ; i ) f ( wt = wc ; i GLYPH<0> 1) f ( GLYPH<14> ti =GLYPH<14> ti GLYPH<0> 1) # 1 = 4 ; (B3)  \nwhere f ( x ) = 1 + 2 tan GLYPH<0> 1 [0 : 5( x GLYPH<0> 1)]. This control scheme allows rapid changes in the timestep without undesirable fluctuations. The timestep proposed by this low-pass filtering scheme can be reduced in the second stage according to a variety of special tests that have hard and soft limits. If a change exceeds its specified hard limit, the current trial solution for the new step is rejected, and the code is forced to do a retry or a backup. If a change exceeds its specified soft limit, the next timestep is reduced proportionally. The current classes of special cases that can reduce the next timestep are limits based on:  \n- 1. Number of Newton iterations required to converge.\n- 2. Maximum absolute change in the mass fraction of hydrogen or helium in any cell.  \n- 3. Maximum relative change in any mass fraction at any cell.\n- 4. Magnitude in the relative change in the structure variables in each cell.\n- 5. Nuclear energy generated in each cell for several categories of nuclear reactions.\n- 6. Changes in the luminosity resulting from nuclear burning.\n- 7. Changes at the photosphere in ln L and T e GLYPH<11> .\n- 8. Changes in ln GLYPH<26> center, ln T center, X (H)center, X (He)center.\n- 9. Magnitude of the change in log( M = M GLYPH<12> ) due to winds or accretion.\n- 10. Mass accreted so that compressional heating is correct (see x 5.3).\n- 11. Changes in the logarithm of the total angular momentum.\n- 12. Distance moved in the HR diagram.\n- 13. Maximum allowed timestep under any circumstance.\n- 14. Any user specifed timestep limit, accomplished by setting max years for timestep , in the optional routine extras check model .  \nFor convergence studies with respect to the timestep it is vital to change the control parameters that are actually setting the timestep. Often, this is just varcontrol target , but in many situations the timestep will be set by one of the special timestep control parameters.', 'Mesh Controls': "Control of the spatial mesh is a key ingredient of a stellar evolution instrument, and requires careful trade-o GLYPH<11> s. The mesh must respond to gradients in the structure, chemical composition, and energy generation, in order to give an accurate result, but it should not be overly dense since that will unnecessarily increase the cost of the calculation.  \nSince MESA star allows for simulations with a fixed inner core mass, M c, the total mass M is M c + M m where M m is the modeled mass. For cell k , MESA star stores the relative cell mass d qk = d mk = M m where d mk is the mass contained in cell k . The relative mass interior to the outer cell face is qk = 1 GLYPH<0> P i = k GLYPH<0> 1 i = 1 d qi , and the total mass interior to the outer cell face is mk = qk GLYPH<3> M m + M c. In all cases, m 1 = M and q 1 = 1. We explicitly keep d qk in addition to q and define q in terms of d qk to avoid the need for evaluating qk GLYPH<0> qk + 1 since that can involve the subtraction of almost equal numbers leading to an undesirable loss of precision (Lesa GLYPH<11> re et al. 2006). For example, in the outer envelope of a star where the qk approach 1, the d qk can be 10 GLYPH<0> 12 or smaller. By storing d qk we have 16 digit precision, whereas in this case, qk GLYPH<0> qk + 1 would only give us 4 digits at best for the relative cell mass.  \nMESA star checks the structure and composition profiles of the model at the beginning of each timestep and, if necessary, adjusts the mesh. A single cell can be split into two or more cells and two or more adjacent cells can be merged. In practice, only a small fraction of the cells are changed during a remesh. This minimizes numerical di GLYPH<11> usion, aids convergence, and keeps the cost of remeshing relatively small. Remeshing is divided into a planning stage and an adjustment stage.  \nThe planning stage determines which cells to split or merge based on the magnitude of allowed cell-to-cell changes in a variety of mesh functions. Built-in mesh functions include gradients of the mass, radius, pressure, temperature, adiabatic gradient, angular velocity and mass fractions above some threshold. Users can add others by defining their own other mesh functions routine.  \nOther controls are provided to increase the sensitivity in regions selected by the user. Examples include increasing the spatial resolution in regions with changes in user-specified abundances with respect to pressure, changes in the energy generation rate with respect to pressure for di GLYPH<11> erent types of burning (e.g., the pp chains, CNO cycles, tripleGLYPH<11> , and others), for regions near burning or non-burning convective boundaries, and others.  \nAfter the mesh functions are evaluated, the relative magnitude of the changes between adjacent cells are determined. The magnitude of change is multiplied by mesh delta coeff to obtain a weighted mesh function. Cells where the weighted changes are 'too large' are marked for splitting, and cells where the changes are 'too small' are marked for merging. For example, if the weighted changes in all mesh functions from cells k to k + n are less than 1, the series of cells from k to k + n are marked for merging. If any weighted mesh function changes from cell k to k + 1 by an amount greater than 1, the larger of cell k and cell k + 1 is marked for splitting. Finally, if adjacent cells have too large of a relative size di GLYPH<11> erence (as defined by mesh max allowed ratio which defaults to 2.5), the larger cell is marked for splitting and the check for excessive ratios is repeated. This can lead to a cascade of splitting in order to ensure that cells sizes do not have excessive jumps.  \nThe adjustment stage executes the remesh plan by performing the merge and split operations to calculate new values for basic variables. Special care is taken to use physical knowledge whenever possible when setting new values. For example, conservation of mass is accounted for when determining new densities, and species conservation is used when setting new mass fractions. Energy conservation is used when setting the temperature (see Paper I), and conservation of angular momentum plays a role in determining the angular velocity. Cells to be split are constructed by first performing a monotonicity-preserving cubic  \ninterpolation (Ste GLYPH<11> en 1990) in mass to obtain the luminosities and ln r values at the new cell boundaries. The new densities are then calculated from the new cell masses and volumes. Next, new composition mass fraction vectors are calculated. For cells being merged, the mass averaged abundances are used. For cells being split, neighboring cells are used to form a linear approximation of mass fraction for each species as a function of mass coordinate within the cell. The slopes are adjusted so that the mass fractions sum to one everywhere, and the functions are integrated over the new cell mass to determine the abundances.", 'Mass Adjustment': 'Mass adjustment for mass loss or accretion is performed at each timestep when do evolve step calls the routine do adjust mass (see x B.1). MESA star o GLYPH<11> ers several ways to set the rate of mass change ˙ M . A constant mass accretion rate (positive ˙ M ) or mass loss rate (negative ˙ M ) can be specified in the input files, a wind can produce a mass loss, the user can set ˙ M in an other wind routine or in an other check model routine. When do adjust mass is called, the timestep dt and the rate of mass change ˙ M are known, and thus the change in mass, GLYPH<14> M = ˙ M GLYPH<14> t .  \nWhen there is a change in mass, instead of adding or removing cells, the total mass is changed by modifying the modeled mass M m, and cell mass sizes are changed by revisions to d qk which in turn changes cell mass locations qk (see x B.4). The mass structure is divided into an inner region where the mk and d mk are unchanged but the qk and d qk change, an outer region where the qk and d qk are unchanged but the mk and d mk change, and an intermediate blending region where all of these change. The selection of the region boundaries is discussed in detail in Paper I. The implementation of GLYPH<15> grav in the newly accreted matter is described in x 5.3.  \nOnce the three regions have been defined, the d qk are updated. In the inner region they are rescaled by M = ( M + GLYPH<14> M ). Thus, d mk , mk , and Xk have the same value before and after a change in mass. This eliminates the possibility of unwanted numerical di GLYPH<11> usion causing unphysical mixing in the center region. In the outer region, cells retain the same value of d qk to improve convergence in the high entropy parts of the star (Sugimoto et al. 1981). In the intermediate region, the d qk are uniformly scaled to make P d qk = 1.  \nThe chemical mass fractions of cells in the intermediate and outer regions are then updated by summing the abundances between the new cell mass boundaries. This step is not necessary for the inner region since those cells have not changed mass location. In the case of mass accretion, the composition of the outermost cells whose enclosed mass totals GLYPH<14> M is set to match the specified accretion abundances. The single cell that is part old material and part newly accreted material is given an appropriately mixed composition.  \nFinally, to create a somewhat better starting model for the Newton iterations (see x B.2), the ln T and ln GLYPH<26> and ln P gas values are revised by monotonic cubic interpolating to the cell center by mass from the values prior to mass adjustment. The ln r and material speed v are also set by monotonic cubic interpolation to the value at the new outer mass boundary. The angular velocity is set by integrating the angular momentum between the new cell mass boundaries and using the new ln r values, conserving the total angular momentum to the floating point limit of the arithmetic.', 'Evolving the Angular Velocity': 'Initialization of rotation in MESA star begins from a non-rotating model. The angular velocity ! is added to the set of model variables and initialized to a constant value throughout the model (i.e., solid body rotation). The initial value of ! can be specified as a surface rotational velocity (in km / s) or as a fraction of the surface critical rotation rate (see x 6). During the subsequent evolution, ! is changed at each timestep by remeshing, mass adjustment, radius adjustment (as part of the structure evolution), optional extra angular momentum removal in the outer layers, and the transport of angular momentum optionally with user-defined source terms for external torques.  \nThe angular velocity ! is defined at cell boundaries. Thus omega(k) is at the outer boundary of cell k , which is the same location as the radius, r(k) , the specific moment of inertia, i rot(k) , and the specific angular momentum, j rot(k) . The mass associated with omega(k) spans the range from the center by mass of cell k outward to the center by mass of cell k-1 and is referred to as dm bar(k) to distinguish it from the cell mass dm(k) .  \nThe remeshing operation splits and merges cells but does not change the physical stellar structure (see x B.4). For regions where there has been a change in the mesh, the values of ! are adjusted to give the same angular momentum as before. More specifically, the angular momentum from the original model is summed over the mass range encompassed by the new dm bar(k) , and omega(k) is adjusted to give the same total for the new model.  \nDuring the mass adjustment operation, when mass is added or removed from the model, cells in the outer layers are moved to new mass locations (see x B.5). As part of this process, the angular velocity values are updated to conserve angular momentum using the same scheme as for remeshing: sum the angular momentum in the original model and set omega(k) in the new model to conserve it. Newly added material from accretion is given the current surface angular velocity. In the case of mass loss, this operation removes the amount of angular momentum contained in the lost mass at the start of the timestep; it does not deal with possible transport of angular momentum into the lost mass during the timestep. That is dealt with by an optional, user-specified removal prior to the angular momentum transport.  \nMESA star performs the transport of angular momentum as a separate operation from the evolution of structure and composition. This is done in order to obtain high accuracy in the angular momentum transport by using substeps and quad-precision linear algebra. It does not introduce additional operator splitting errors since ! is not used in the structure and abundance equations. So we solve for the new structure and composition after any mass change and before the transport of angular momentum. Calculation of the new stellar structure changes the radii but does not change the mass partitioning of the model (see x B.2). Given the new radius r(k) , we calculate the new i rot(k) . Then using the unchanged j rot(k) , omega(k) is set to j rot(k)/i rot(k) to conserve specific angular momentum. Since dm bar(k) has not changed, this also conserves total angular momentum.  \nNext, MESA star applies an optional, user-specified amount of angular momentum loss in the outer surface layers. This is to account for possible transport of angular momentum during the timestep from these outer layers into the mass removed by the mass adjustment operation.  \nThe final operation is the transport of angular momentum within the star, which is treated with a di GLYPH<11> usion approximation (Endal &Sofia 1978; Pinsonneault et al. 1989; Heger et al. 2000)  \n@! @ t ! m = 1 i @ @ m ! t " (4 GLYPH<25> r 2 GLYPH<26> ) 2 i GLYPH<23> @! @ m ! t # GLYPH<0> 2 ! r @ r @ t ! m 1 2 d ln i d ln r ! ; (B4)  \nwhere i is the specific moment of inertia of a shell at mass coordinate m , and GLYPH<23> is the turbulent viscosity determined as the sum of the di GLYPH<11> usion coe GLYPH<14> cients for convection, double di GLYPH<11> usion, overshooting and rotationally-induced instabilities (see x 6). The di GLYPH<11> usive transport is carefully implemented to accurately conserve angular momentum. The angular momentum associated with location k is dm bar(k)*i rot(k)*omega(k) . The change in angular momentum for k is determined by the flux in angular momentum from k GLYPH<0> 1 to k and from k + 1 to k . The flux from k GLYPH<0> 1 to k is set by GLYPH<23> ( k GLYPH<0> 1) and the di GLYPH<11> erence between omega(k) and omega(k-1) . The flux from k + 1 to k is found similarly using GLYPH<23> ( k ) and the di GLYPH<11> erence between omega(k) and omega(k+1) . Source terms for location k are applied by user-supplied values for extra jdot(k) or extra omegadot(k) . The finite di GLYPH<11> erence equation for the e GLYPH<11> ects of the transport and source terms is solved over the stellar timestep with an implicit time integration that uses multiple smaller timesteps. The sizes of these substeps are determined by the timescale set by the di GLYPH<11> usion coe GLYPH<14> cients and the di GLYPH<11> erences in ! . It is not unusual to use 10 or more substeps to evolve omega(k) over the stellar timestep. Each implicit substep is solved using a quad-precision tridiagonal matrix routine. The conservation of total angular momentum is monitored and the stellar timestep is rejected if there is any deviation from conservation by more than a user-specified factor. In practice, we find the total angular momentum is conserved over the stellar timestep to within a few digits of the floating point limit of the arithmetic.', 'Free Parameters': 'Stellar evolution calculations involve the choice of a number of free parameters. The values of these parameters are not determined by first principles, and in the literature one can find a range of possibilities. In some cases the parameters can be constrained by matching a restricted set of observations; in other cases they represent common choices. Users need to be aware that their results will depend on these values, and that in some cases the sensitivity can be large. Below we illustrate this by discussing some of the main parameters involved in the mixing of stellar interiors.', 'Convection': 'In the literature the value of the mixing length parameter GLYPH<11> MLT (see e.g. Paper I for a definition) is usually found to vary within the range 1 : 0 . GLYPH<11> MLT . 2 : 0. E GLYPH<11> orts are ongoing to eliminate this free parameter (e.g., Arnett et al. 2010).', 'Overshooting': 'In the literature the adopted value for the convective core overshooting parameter is in the range 0 : 1 . f ov . 0 : 6, in units of the pressure scale height GLYPH<21> P , when the overshoot zone is considered to be fully mixed (Maeder & Meynet 1987; Dupret et al. 2004; Straka et al. 2005; Claret 2007; Briquet et al. 2007). When overshoot mixing is treated as an exponential decay process the free parameter should be smaller, f ov GLYPH<24> 0 : 016, (see the discussion by Herwig 2000b). MESA has the ability to treat overshoot mixing zones as either fully mixed or in the exponential decay formalism.  \nIt has been suggested that the overshooting parameter is a function of both mass and metal abundance, in that it should transition smoothly from zero to a maximum value over a small range of stellar mass where a convective core is present on the main sequence (Woo & Demarque 2001; VandenBerg et al. 2006), but see also Claret (2007). A dependency on the evolutionary stage seems also likely (Herwig et al. 1997; Meakin & Arnett 2007; Tian et al. 2009). See Fig. 13, 14 and 15 in x 4 for an example of the sensitivity of the calculations to changes in the f ov parameter.', 'Semiconvection': 'Semiconvection, as implemented in MESA , requires a choice of the free parameter GLYPH<11> sc (see x 4.1). In the literature this spans the range 0 : 001 . GLYPH<11> sc . 1 : 0 (Langer 1991; Yoon et al. 2006). See Fig. 13, 14 and 15 in x 4 for an example of the sensitivity of the calculations to changes in the GLYPH<11> sc parameter. Research is ongoing to eliminate this free parameter (Wood et al. 2013; Spruit 2013).', 'Thermohaline Mixing': 'The implemented formulation for thermohaline mixing requires the adoption of the free parameter GLYPH<11> th (see x 4.2). In the literature this parameter can be usually found within the range 1 . GLYPH<11> th . 667 (Kippenhahn et al. 1980; Charbonnel & Zahn 2007; Cantiello & Langer 2010; Stancli GLYPH<11> e 2010; Wachlin et al. 2011). Research is ongoing to eliminate this free parameter (Traxler et al. 2011; Brown et al. 2013).', 'Multicore Performance': "MESA implements shared memory multiprocessing via OpenMP 11 . Paper I explored the runtime scaling of MESA star which at that time used a banded matrix linear algebra solver that did not benefit from multiple cores. A large part of the performance improvement in MESA star since Paper I comes from converting to a parallel block tridiagonal linear algebra solver derived from BCYCLIC (Hirshman et al. 2010). This improved solver is particularly important since linear algebra is typically the largest part of the runtime in MESA star . In addition, the new algorithm has the desirable property of producing numerically identical results independent of the number of cores, an attribute that is not generally true of parallel matrix solvers.  \nOur test case is a 1 : 5 M GLYPH<12> model with Z = 0 : 02 that is evolved from the ZAMS until the central H mass fraction falls to 0.35. This model includes 25 isotopes and 4 structure variables per cell with a variable number of zones typically exceeding 1700. The test takes GLYPH<24> 55 time steps to cover GLYPH<24> 1 : 4 Gyr and uses the default amount of I / O.  \nFigure 48 shows the scaling behavior of some key components of MESA star under GFORTRAN 4.7.2. on a 12 core 2010 Apple MacPro. The dotted line shows the ideal scaling relation where doubling the number of cores cuts the run time in half. The linear algebra, labeled 'mtx,' dominates the total run time as the number of cores increases. For example, in the case of 12 cores it accounts for about half the total and is 2.5 times larger than 'net', the evaluation of the nuclear reaction network. The net evaluations closely approach the ideal scaling behavior because they can be done in parallel, each cell independent of the others, with one core working on one cell at a time. The equation of state component, labeled 'eos,' also closely approaches the ideal scaling law while consuming less than a third of the run time for the net. The component labeled 'eqns,' which includes the evaluation of the structure equations and the creation of the block tridiagonal matrix, also is close to the ideal scaling law and costs about the same as the eos. The 'other' component is everything else. It is dominated by processes that currently are not e GLYPH<14> cient to parallelize because of the relatively large overhead for OpenMP operations. Consequently it remains at roughly a constant run time independent of the number of cores. When a significantly larger number of cores per processor becomes available, the larger operations in this category will have to be reworked or they will dominate the total run time.  \nThe run time also depends on the hardware, the quality of the compiled code, and the e GLYPH<14> ciency of the OpenMP implementation. For example, we ran the test case under GFORTRAN 4.7.2 on a 40-core server. While we obtained a speedup of 5.6 in going from 1 to 12 cores for the 12 core machine (see Figure 48), we find a speedup of only 4.8 on the 40-core server in going from 1 to 12 cores. Moreover, the speedup per core dropped steeply beyond 8-12 cores on the 40-core server, confirming the expectation that much work will be required to make full use of machines with many cores.", 'Visualization': "MESA star provides alphanumeric output at user-specified regular intervals. In addition, the routines in module star/public/pgstar.f provide an option for concurrent graphical output with the PGPLOT 12 library to create on-screen plots that can be saved for post-processing into animations of an evolutionary sequence. A variety of options are provided and are all configurable through the PGstar inlist. For example, a PGstar X11 window can simultaneously hold an H-R diagram, a T cGLYPH<26> c diagram, and interior profiles of physical variables, such as nuclear energy generation and composition. The PGstar inlist is read at each timestep, so the display options can be changed without have to stop MESA star .  \nSince Paper I, a number of MESA star users have developed and released toolkits 13 to visualize the alphanumeric output with common graphical packages including: Mathematica scripts (contributed by Richard O'Shaughnessy) and the intuitive and e GLYPH<14> cient graphical user interface MESAFace (Giannotti et al. 2012); MatLab utilities (contributed by Dave Spiegel and  \nFigure 48. Scaling behavior of various components of MESA star using 1, 2, 3, 6, and 12 cores. The top curve shows the total run time, and the lower solid curves show the run times for the components of the total. The dotted line shows the ideal scaling relation.  \n<!-- image -->  \nGongjie Li); IDL functions (contributed by Rich Townsend); Python scripts (contributions from Falk Herwig and the NuGrid collaboration, David Kaplan, Alfred Gautschy, William Wolf); and Tioga scripts (contributed by Christopher Mankovich and Bill Paxton).", 'Operating System and Compiler Considerations': 'We next consider the implications of running MESA compiled with di GLYPH<11> erent compilers and on di GLYPH<11> erent operating systems. The operating systems examined are Linux (Gentoo 2.1; kernel 3.6.11) and Apple OS X (10.7.5), both 64-bit; the compilers are GNU gfortran 4.9.0 on Linux and OS X and Intel ifort 13.1.0 on Linux. We used OpenMP in all cases. For optimization we used -O2 with gfortran and -O1 with ifort. The models described here were computed with MESA revision 4942. The comparison case is the example solar model from the test suite . It evolves a calibrated solar model from the pre-main sequence to the solar age, 4.57 Gyr.  \nFigure 49. Relative di GLYPH<11> erence in the solar sound speed profile from three models run using the same inlist but on di GLYPH<11> erent operating systems and / or compilers.  \n<!-- image -->  \nglyph[circledot]  \nIn Figure 49 we compare the sound speed profiles of these solar models in the same way that standard solar models are compared with the solar sound speed profile in Section 3.1, Figure 7. The reference was chosen to be the Linux / gfortran model; the other two models were compared with it by interpolating their sound speed profiles onto the radial grid of the reference model and then taking the relative sound speed di GLYPH<11> erence with respect to that reference model. Figure 49 indicates that models computed on di GLYPH<11> erent platforms are consistent with one another at the level of the numerical tolerances with which the equations are solved. The di GLYPH<11> erences shown here are about one part per million or less, whereas the di GLYPH<11> erences between standard solar models and the helioseismic data are in the parts per thousand (Figure 7)-a factor of 1000 di GLYPH<11> erence. While this result is reassuring, it is worthwhile to close this discussion with the comment that the consistency found in a low-mass, non-rotating model evolved about half way through the main sequence phase will not be representative of other cases dealing with di GLYPH<11> erent physics, stellar masses, and evolutionary phases. Projects using MESA with heterogeneous architectures should perform their own consistency checks.  \nTable 6 Principal components of the MESA Software Development Kit', 'THE MESA SOFTWARE DEVELOPMENT KIT': "MESA is provided as source code, allowing users access to all of the implementation details. Installation necessarily involves building the code from source, which is a non-trivial task. A successful build requires cooperation between the operating system, compiler, libraries, and utilities.  \nTo address this issue we have created the MESA Software Development Kit ( SDK ), which packages everything necessary to establish a unified and maintained build environment. 14 The principal components of the SDK are summarized in Table 6; all of these are distributed under an open-source license (detailed in the table), permitting their redistribution without financial or copyright encumbrances. Perhaps the most important component is the GFORTRAN compiler, part of the GNU Compiler Collection. GFORTRAN implements almost all of the Fortran 2003 (F2003) standard, and benefits from a high level of community support.  \nThe SDK is available for Intel x86 and x86-64 CPU architectures running the Linux and Mac OS X operating systems (these platforms comprise most of the MESA user base). Installation of the kit is straightforward, requiring a tar archive to be unpacked (Linux) or an application folder to be copied (OS X), followed by the initialization of a few environment variables. By default, MESA is configured to compile 'out-of-the-box' with the SDK . MESA can also be compiled without the SDK , using any alternate compiler which supports the F2003 standard. In this respect, GFORTRAN should not be viewed as the MESA compiler (nor the full SDK as the MESA build environment). MESA will adhere to Fortran standards rather than rely on vendor-specific extensions.  \nUptake of the SDK has been very rapid: at the time of writing, we estimate over 90% of the MESA community (over 500 users) are using the SDK . This growth has been matched by a significant decline in the number of installation support requests, and a corresponding reduction in the time taken to resolve these requests. With these maintenance overheads curbed, the MESA developers are able to devote more of their time to refining and extending the code.", 'REFERENCES': "Alibert, Y., Mordasini, C., & Benz, W. 2011, A&A, 526, A63  \nAllard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357  \nAlthaus, L. G., Serenelli, A. M., C'orsico, A. H., & Montgomery, M. H. 2003, A&A, 404, 593  \nAngulo, C., et al. 1999, Nuclear Physics A, 656, 3  \nArcoragi, J.-P., & Fontaine, G. 1980, ApJ, 242, 1208  \nArnett, D., Meakin, C., & Young, P. A. 2010, ApJ, 710, 1619  \nBohm-Vitense, E. 1958, Zeitschrift fur Astrophysik, 46, 108  \nBaglin, A., Auvergne, M., Barge, P., Deleuil, M., Michel, E., & CoRoT Exoplanet Science Team. 2009, in IAU Symposium, Vol. 253, IAU Symposium, 71-81  \nBahcall, J. N., Basu, S., & Pinsonneault, M. H. 1998, Physics Letters B, 433, 1  \nBalbus, S. A., Latter, H., & Weiss, N. 2012, MNRAS, 420, 2457  \nBass, G., Orosz, J. A., Welsh, W. F., Windmiller, G., Ames Gregg, T., Fetherolf, T., Wade, R. A., & Quinn, S. N. 2012, ArXiv e-prints  \nBeck, P. G., et al. 2012, Nature, 481, 55  \nBergeron, P., Wesemael, F., Lamontagne, R., Fontaine, G., Sa GLYPH<11> er, R. A., & Allard, N. F. 1995, ApJ, 449, 258  \nBildsten, L., Paxton, B., Moore, K., & Macias, P. J. 2012, ApJ, 744, L6  \nBisnovatyi-Kogan, G. S., & Dorodnitsyn, A. V. 1999, A&A, 344, 647  \nBloecker, T. 1995, A&A, 297, 727  \nBohm, K. H., & Cassinelli, J. 1971, A&A, 12, 21  \nBonanno, A., Kuker, M., & Patern'o, L. 2007, A&A, 462, 1031  \nBorucki, W., et al. 2009, in IAU Symposium, Vol. 253, IAU Symposium, 289-299  \nBoss, A. P. 2011, ApJ, 731, 74  \nBraithwaite, J. 2006, A&A, 449, 451  \nBraithwaite, J., & Spruit, H. C. 2004, Nature, 431, 819  \nBrandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1  \nBrassard, P., Fontaine, G., Wesemael, F., Kawaler, S. D., & Tassoul, M. 1991, ApJ, 367, 601  \nBriquet, M., Morel, T., Thoul, A., Scuflaire, R., Miglio, A., Montalb'an, J., Dupret, M.-A., & Aerts, C. 2007, MNRAS, 381, 1482  \nBrott, I., et al. 2011, A&A, 530, A115  \nBrown, J. M., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 34  \nBrown, T. M., Christensen-Dalsgaard, J., Dziembowski, W. A., Goode, P., Gough, D. O., & Morrow, C. A. 1989, ApJ, 343, 526  \nBrown, T. M., Gilliland, R. L., Noyes, R. W., & Ramsey, L. W. 1991, ApJ, 368, 599  \nBurkart, J., Quataert, E., Arras, P., & Weinberg, N. N. 2012a, ArXiv e-prints  \n-. 2012b, MNRAS, 421, 983  \nCalder, A. C., et al. 2007, ApJ, 656, 313  \nCantiello, M., & Braithwaite, J. 2011, A&A, 534, A140  \nCantiello, M., & Langer, N. 2010, A&A, 521, A9  \nCantiello, M., et al. 2009, A&A, 499, 279  \nCarlberg, J. K., Cunha, K., Smith, V. V., & Majewski, S. R. 2012, ApJ, 757, 109  \nCarter, J. A., et al. 2011, Science, 331, 562  \nCassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094  \nChang, P., & Hui, L. 2011, ApJ, 732, 25  \nCharbonnel, C., & Talon, S. 2005, Science, 309, 2189  \nCharbonnel, C., & Zahn, J.-P. 2007, A&A, 467, L15  \nChatzopoulos, E., Robinson, E. L., & Wheeler, J. C. 2012, ApJ, 755, 95  \nChatzopoulos, E., & Wheeler, J. C. 2012, ApJ, 748, 42  \nChristensen-Dalsgaard, J. 2008a, Ap&SS, 316, 113  \n-. 2008b, Ap&SS, 316, 13  \nChristensen-Dalsgaard, J., & Thompson, M. J. 1997, MNRAS, 284, 527  \nChristensen-Dalsgaard, J., et al. 1996, Science, 272, 1286  \nClaret, A. 2007, A&A, 475, 1019  \nCox, J. P., & Giuli, R. T. 1968, Principles of stellar structure (New York: Gordon and Breach)  \nCyburt, R. H., et al. 2010, ApJS, 189, 240  \nDavis, A.-C., Lim, E. A., Sakstein, J., & Shaw, D. J. 2012, Phys. Rev. D, 85, 123006  \nde Bruijne, J. H. J. 2012, Ap&SS, 341, 31  \nde Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259  \nDeheuvels, S., & Michel, E. 2011, A&A, 535, A91  \nDeheuvels, S., et al. 2010, A&A, 515, 87  \nDenissenkov, P. A. 2010, ApJ, 723, 563  \n-. 2012, ApJ, 753, L3  \nDenissenkov, P. A., Herwig, F., Bildsten, L., & Paxton, B. 2013, ApJ, 762, 8  \nDenissenkov, P. A., & Pinsonneault, M. 2007, ApJ, 655, 1157  \n-. 2008, ApJ, 684, 626  \nDenissenkov, P. A., Pinsonneault, M., Terndrup, D. M., & Newsham, G. 2010, ApJ, 716, 1269  \nDotter, A., Chaboyer, B., Jevremovi'c, D., Kostov, V., Baron, E., & Ferguson, J. W. 2008, ApJS, 178, 89  \nDupret, M.-A., Thoul, A., Scuflaire, R., Daszy'nska-Daszkiewicz, J., Aerts, C., Bourge, P.-O., Waelkens, C., & Noels, A. 2004, A&A, 415, 251  \nEggenberger, P., Meynet, G., Maeder, A., Hirschi, R., Charbonnel, C., Talon, S., & Ekstrom, S. 2008, Ap&SS, 316, 43  \nEggenberger, P., Montalb'an, J., & Miglio, A. 2012, A&A, 544, L4  \nEkstrom, S., et al. 2012, A&A, 537, A146  \nEndal, A. S., & Sofia, S. 1976, ApJ, 210, 184  \n-. 1978, ApJ, 220, 279  \nEvans, C. J., et al. 2011, A&A, 530, A108  \nFerguson, J. W., Alexander, D. R., Allard, F., Barman, T., Bodnarik, J. G., Hauschildt, P. H., He GLYPH<11> ner-Wong, A., & Tamanai, A. 2005, ApJ, 623, 585 Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102  \nFortney, J. J., Saumon, D., Marley, M. S., Lodders, K., & Freedman, R. S. 2006, ApJ, 642, 495  \nFreedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504  \nFriend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701  \nFrommhold, L., Abel, M., Wang, F., Gustafsson, M., Li, X., & Hunt, K. 2010, Molecular Physics, 108, 2265  \nFuller, J., & Lai, D. 2012a, ArXiv e-prints  \n-. 2012b, ApJ, 756, L17  \nGeorgy, C., Ekstrom, S., Meynet, G., Massey, P., Levesque, E. M., Hirschi, R., Eggenberger, P., & Maeder, A. 2012, A&A, 542, A29  \nGeorgy, C., Meynet, G., & Maeder, A. 2011, A&A, 527, A52  \nGiannotti, M., Wise, M., & Mohammed, A. 2012, ArXiv e-prints  \nGlebbeek, E., Gaburov, E., de Mink, S. E., Pols, O. R., & Portegies Zwart, S. F. 2009, A&A, 497, 255  \nGough, D. O. 1986, in Hydrodynamic and Magnetodynamic Problems in the Sun and Stars, ed. Y. Osaki, 117  \nGrevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, ed. N. Prantzos, E. Vangioni-Flam, & M. Casse, 15-25  \nGrevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161  \nGuillot, T. 2005, Annual Review of Earth and Planetary Sciences, 33, 493  \n-. 2010, A&A, 520, A27  \nGuillot, T., & Morel, P. 1995, A&AS, 109, 109  \nHairer, E., & Wanner, G. 1996, Solving Ordinary Di GLYPH<11> erential Equations. II. Sti GLYPH<11> and Di GLYPH<11> erential-Algebraic Problems, 2nd edn. (Berlin: Springer)  \nHeger, A., Jeannin, L., Langer, N., & Bara GLYPH<11> e, I. 1997, A&A, 327, 224  \nHeger, A., & Langer, N. 1998, A&A, 334, 210  \n-. 2000, ApJ, 544, 1016  \nHeger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368  \nHeger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350  \nHenyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841  \nHerwig, F. 2000a, A&A, 360, 952  \n-. 2000b, A&A, 360, 952  \n-. 2004, ApJS, 155, 651  \nHerwig, F., Bloecker, T., Schoenberner, D., & El Eid, M. 1997, A&A, 324, L81  \nHirschi, R., Meynet, G., & Maeder, A. 2004, A&A, 425, 649  \nHirshman, S. P., Perumalla, K. S., Lynch, V. E., & Sanchez, R. 2010, J. Comput. Phys., 229, 6392  \nHooke, R., & Jeeves, T. A. 1961, J. ACM, 8, 212  \nHuang, X., & Cumming, A. 2012, ApJ, 757, 47  \nHubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415  \nHumphreys, R. M., & Davidson, K. 1994, Publ.Astron.Soc.Pac., 106, 1025  \nHunter, I., et al. 2007, A&A, 466, 277  \nIben, Jr., I. 1991, ApJS, 76, 55  \nIglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752  \n-. 1996, ApJ, 464, 943  \nIglesias, C. A., Rogers, F. J., & Wilson, B. G. 1992, ApJ, 397, 717  \nImbriani, G., et al. 2004, A&A, 420, 625  \nIrwin, J. M., et al. 2011, ApJ, 742, 123  \nIvezic, Z., et al. 2008, ArXiv e-prints  \nJoss, P. C., Salpeter, E. E., & Ostriker, J. P. 1973, ApJ, 181, 429  \nKaiser, N., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7733, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series  \nKawaler, S. D. 1988, ApJ, 333, 236  \nKippenhahn, R., Ruschenplatt, G., & Thomas, H.-C. 1980, A&A, 91, 175  \nKippenhahn, R., & Thomas, H.-C. 1970, in IAU Colloq. 4: Stellar Rotation, ed. A. Slettebak, 20  \nKjeldsen, H., & Bedding, T. R. 1995, A&A, 293, 87  \nKjeldsen, H., Bedding, T. R., & Christensen-Dalsgaard, J. 2008, ApJ, 683, L175  \nKoch, D. G., et al. 2010, ApJ, 713, L79  \nKrishnamurti, R. 2003, Journal of Fluid Mechanics, 483, 287  \nKuker, M., Rudiger, G., & Kitchatinov, L. L. 2011, A&A, 530, A48  \nKunz, R., Fey, M., Jaeger, M., Mayer, A., Hammer, J. W., Staudt, G., Harissopulos, S., & Paradellis, T. 2002, ApJ, 567, 643  \nLanger, N. 1991, A&A, 252, 669  \nLanger, N. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 120, Luminous Blue Variables: Massive Stars in Transition, ed. A. Nota & H. Lamers, 83  \n-. 1998, A&A, 329, 551  \n-. 2012, ARA&A, 50, 107  \nLanger, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207  \nLanger, N., Heger, A., Wellstein, S., & Herwig, F. 1999, A&A, 346, L37  \nLaughlin, G., Bodenheimer, P., & Adams, F. C. 2004, ApJ, 612, L73  \nLaw, N. M., et al. 2009, PASP, 121, 1395  \nLesa GLYPH<11> re, P., Han, Z., Tout, C. A., Podsiadlowski, P., & Martin, R. G. 2006, MNRAS, 368, 187  \nLi, Y., & Gong, Z. G. 1994, A&A, 289, 449  \nLiu, C., Bailer-Jones, C. A. L., Sordo, R., Vallenari, A., Borrachero, R., Luri, X., & Sartoretti, P. 2012, MNRAS, 426, 2463  \nLloyd, J. P. 2011, ApJ, 739, L49  \nLodders, K. 2003, ApJ, 591, 1220  \nMacDonald, J., & Mullan, D. J. 2012, MNRAS, 421, 3084  \nMaeder, A. 1987a, A&A, 173, 247  \n-. 1987b, A&A, 178, 159  \n-. 2009, Physics, Formation and Evolution of Rotating Stars (Berlin: Springer)  \nMaeder, A., & Meynet, G. 1987, A&A, 182, 243  \n-. 2000a, A&A, 361, 159  \n- -. 2000b, ARA&A, 38, 143  \n-. 2003, A&A, 411, 543  \n-. 2004, A&A, 422, 225  \nMaeder, A., & Zahn, J.-P. 1998, A&A, 334, 1000  \nMarley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541  \nMathur, S., et al. 2012, ApJ, 749, 152  \nMeakin, C. A., & Arnett, D. 2007, ApJ, 667, 448  \nMedin, Z., & Cumming, A. 2010, Phys. Rev. E, 81, 036107  \nMetcalfe, T. S., et al. 2012, ApJ, 748, L10  \nMeynet, G., Eggenberger, P., & Maeder, A. 2011, A&A, 525, L11  \nMeynet, G., & Maeder, A. 1997, A&A, 321, 465  \n-. 2000, A&A, 361, 101  \nMihalas, D. 1978, Stellar atmospheres / 2nd edition / (San Francisco, W. H. Freeman and Co.)  \nMohr, P. J., Taylor, B. N., & Newell, D. B. 2008, Reviews of Modern Physics, 80, 633  \nMonteiro, M. J. P. F. G., ed. 2009, Astrophysics and Space Science, Vol. 316, Evolution and Seismic Tools for Stellar Astrophysics (Springer)  \nMoravveji, E., Moya, A., & Guinan, E. F. 2012, ApJ, 749, 74  \nMosser, B., et al. 2012, A&A, 548, A10  \nMugrauer, M., & Neuhauser, R. 2009, A&A, 494, 373  \nNomoto, K. 1982, ApJ, 253, 798  \nNomoto, K., & Sugimoto, D. 1977, PASJ, 29, 765  \nNugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227  \nOwen, J. E., & Wu, Y. 2013, ArXiv e-prints  \nOwocki, S. P., Gayley, K. G., & Shaviv, N. J. 2004, ApJ, 616, 525  \nPaczy'nski, B. 1969, Acta Astron., 19, 1  \n-. 1970, Acta Astron., 20, 47  \nPamyatnykh, A. A. 1999, Acta Astron., 49, 119  \nPassy, J.-C., Mac Low, M.-M., & De Marco, O. 2012, ApJ, 759, L30  \nPatience, J., et al. 2002, ApJ, 581, 654  \nPaxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesa GLYPH<11> re, P., & Timmes, F. 2010, ApJS, 192, 3  \nPetrovic, J., Langer, N., Yoon, S.-C., & Heger, A. 2005, A&A, 435, 247  \nPinsonneault, M. H., Kawaler, S. D., Sofia, S., & Demarque, P. 1989, ApJ, 338, 424  \nPollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., & Greenzweig, Y. 1996, Icarus, 124, 62  \nPotekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82  \nPotter, A. T., Tout, C. A., & Brott, I. 2012a, MNRAS, 2883  \nPotter, A. T., Tout, C. A., & Eldridge, J. J. 2012b, MNRAS, 419, 748  \nPowell, M. 2009, The BOBYQA algorithm for bound optimization without derivatives, Tech. rep., Department of Applied Mathematics and Theoretical Physics, University of Cambridge  \nRauscher, T., & Thielemann, F. 2000, At. Data Nucl. Data Tables, 75, 1  \nReimers, D. 1975, in Problems in stellar atmospheres and envelopes., ed. B. Baschek, W. H. Kegel, & G. Traving (New York: Spring-Verlag), 229-256  \nRoell, T., Neuhauser, R., Seifahrt, A., & Mugrauer, M. 2012, A&A, 542, A92  \nRogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064  \nRohrmann, R. D. 2001, MNRAS, 323, 699  \nRohrmann, R. D., Althaus, L. G., Garc'ıa-Berro, E., C'orsico, A. H., & Miller Bertolami, M. M. 2012, ArXiv e-prints  \nRoxburgh, I. W., & Vorontsov, S. V. 2003, A&A, 411, 215  \nRuediger, G., von Rekowski, B., Donahue, R. A., & Baliunas, S. L. 1998, ApJ, 494, 691  \nSato, B., et al. 2005, ApJ, 633, 465  \nSaumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713  \nSchatzman, E. 1962, Annales d'Astrophysique, 25, 18  \nSchneider, A. S., Hughto, J., Horowitz, C. J., & Berry, D. K. 2012, Phys. Rev. E, 85, 066405  \nSchou, J., et al. 1998, ApJ, 505, 390  \nSeager, S., & Deming, D. 2010, ARA&A, 48, 631  \nSilva Aguirre, V., et al. 2013, ArXiv e-prints  \nSoderlind, G., & Wang, L. 2006, J. Comput. Appl. Math., 185, 225  \nSpruit, H. C. 2002, A&A, 381, 923  \n-. 2013, A&A, 552, A76  \nStancli GLYPH<11> e, R. J. 2010, MNRAS, 403, 505  \nSte GLYPH<11> en, M. 1990, A&A, 239, 443  \nStello, D., et al. 2009, ApJ, 700, 1589  \nStraka, C. W., Demarque, P., & Guenther, D. B. 2005, ApJ, 629, 1075  \nStraniero, O., Dom'ınguez, I., Imbriani, G., & Piersanti, L. 2003, ApJ, 583, 878  \nSu'arez-Madrigal, A., Krumholz, M., & Ramirez-Ruiz, E. 2013, ArXiv e-prints  \nSugimoto, D., Nomoto, K., & Eriguchi, Y. 1981, Progress of Theoretical Physics Supplement, 70, 115  \nSuijs, M. P. L., Langer, N., Poelarends, A.-J., Yoon, S.-C., Heger, A., & Herwig, F. 2008, A&A, 481, L87  \nTassoul, M. 1980, ApJS, 43, 469", 'Modules for Experiments in Stellar Astrophysics (MESA)': 'Thompson, M. J., Christensen-Dalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599  \nThompson, M. J., et al. 1996, Science, 272, 1300  \nTian, C.-L., Deng, L.-C., & Chan, K.-L. 2009, MNRAS, 398, 1011  \nTimmes, F. X. 1999, ApJS, 124, 241  \nTimmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501  \nTorres, G., Andersen, J., & Gim´enez, A. 2010, A&A Rev., 18, 67  \nTownsend, R. H. D., Oksala, M. E., Cohen, D. H., Owocki, S. P., & ud-Doula, A. 2010, ApJ, 714, L318  \nTownsley, D. M., & Bildsten, L. 2004, ApJ, 600, 390  \nTraxler, A., Garaud, P., & Stellmach, S. 2011, ApJ, 728, L29  \nud-Doula, A., & Owocki, S. P. 2002, ApJ, 576, 413  \nUdry, S., & Santos, N. C. 2007, ARA&A, 45, 397  \nUlrich, R. K. 1972, ApJ, 172, 165  \nUnno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: University of Tokyo Press)  \nUshomirsky, G., Matzner, C. D., Brown, E. F., Bildsten, L., Hilliard, V. G., & Schroeder, P. C. 1998, ApJ, 497, 253  \nvan Saders, J. L., & Pinsonneault, M. H. 2012, ApJ, 751, 98  \nVandenBerg, D. A., Bergbusch, P. A., & Dowler, P. D. 2006, ApJS, 162, 375  \nVink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574  \nvon Zeipel, H. 1924, MNRAS, 84, 665  \nWachlin, F. C., Miller Bertolami, M. M., & Althaus, L. G. 2011, A&A, 533, A139  \nWeber, E. J., & Davis, Jr., L. 1967, ApJ, 148, 217  \nWoo, J.-H., & Demarque, P. 2001, AJ, 122, 1602  \nWood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157  \nWoosley, S. E. 1993, ApJ, 405, 273  \nWoosley, S. E., & Heger, A. 2006, ApJ, 637, 914  \nYang, X. H., Fu, J. N., & Zha, Q. 2012, AJ, 144, 92  \nYoon, S.-C., & Cantiello, M. 2010, ApJ, 717, L62  \nYoon, S.-C., & Langer, N. 2005, A&A, 443, 643  \nYoon, S.-C., Langer, N., & Norman, C. 2006, A&A, 460, 199  \nYoon, S.-C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940  \nYork, D. G., et al. 2000, AJ, 120, 1579  \nYoung, P. A., Liebst, K., & Pagano, M. 2012, ApJ, 755, L31  \nYurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, MNRAS, 413, 1828  \nZahn, J.-P. 1992, A&A, 265, 115  \nZahn, J.-P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145  \nZapolsky, H. S., & Salpeter, E. E. 1969, ApJ, 158, 809'}