We substantially update the capabilities of the opensource software instrument Modules for Experiments in Stellar Astrophysics MESA. MESA can now simultaneously evolve an interacting pair of differentially rotating stars undergoing transfer and loss of mass and angular momentum greatly enhancing the prior ability to model binary evolution. New MESA capabilities in fully coupled calculation of nuclear networks with hundreds of isotopes now allow MESA to accurately simulate the advanced burning stages needed to construct supernova progenitor models. Implicit hydrodynamics with shocks can now be treated with MESA enabling modeling of the entire massive star lifecycle from premainsequence evolution to the onset of core collapse and nucleosynthesis from the resulting explosion. Coupling of the GYRE nonadiabatic pulsation instrument with MESA allows for new explorations of the instability strips for massive stars while also accelerating the astrophysical use of asteroseismology data. We improve the treatment of mass accretion giving more accurate and robust nearsurface profiles. A new MESA capability to calculate weak reaction rates onthefly from input nuclear data allows better simulation of accretion induced collapse of massive white dwarfs and the fate of some massive stars. We discuss the ongoing challenge of chemical diffusion in the strongly coupled plasma regime and exhibit improvements in MESA that now allow for the simulation of radiative levitation of heavy elements in hot stars. We close by noting that the MESA software infrastructure provides bitforbit consistency for all results across all the supported platforms a profound enabling capability for accelerating MESAs development.
['binaries: general', 'methods: numerical', 'nuclear reactions', 'nucleosynthesis', 'abundances', 'shock waves', 'stars: evolution', 'stars: oscillations', 'Astrophysics - Solar and Stellar Astrophysics']
Modules for Experiments in Stellar Astrophysics MESA Binaries Pulsations and Explosions
{'MODULES FOR EXPERIMENTS IN STELLAR ASTROPHYSICS (MESA): BINARIES, PULSATIONS, AND EXPLOSIONS': 'B ill P axton , 1 P ablo M archant , 2 J osiah S chwab , 3,4 E van B. B auer , 5 L ars B ildsten , 1,5 M atteo C antiello , 1 L uc D essart , 6 R. F armer , 7 H. H u , 8 N. L anger , 2 R.H.D. T ownsend , 9 D ean M. T ownsley , 10 and F.X. T immes 7  \nDraft version January 10, 2017', 'ABSTRACT': "We substantially update the capabilities of the open-source software instrument Modules for Experiments in Stellar Astrophysics ( MESA ). MESA can now simultaneously evolve an interacting pair of di GLYPH<11> erentially rotating stars undergoing transfer and loss of mass and angular momentum, greatly enhancing the prior ability to model binary evolution. New MESA capabilities in fully coupled calculation of nuclear networks with hundreds of isotopes now allow MESA to accurately simulate advanced burning stages needed to construct supernova progenitor models. Implicit hydrodynamics with shocks can now be treated with MESA , enabling modeling of the entire massive star lifecycle, from pre-main sequence evolution to the onset of core collapse and nucleosynthesis from the resulting explosion. Coupling of the GYRE non-adiabatic pulsation instrument with MESA allows for new explorations of the instability strips for massive stars while also accelerating the astrophysical use of asteroseismology data. We improve treatment of mass accretion, giving more accurate and robust near-surface profiles. A new MESA capability to calculate weak reaction rates 'on-the-fly' from input nuclear data allows better simulation of accretion induced collapse of massive white dwarfs and the fate of some massive stars. We discuss the ongoing challenge of chemical di GLYPH<11> usion in the strongly coupled plasma regime, and exhibit improvements in MESA that now allow for the simulation of radiative levitation of heavy elements in hot stars. We close by noting that the MESA software infrastructure provides bit-for-bit consistency for all results across all the supported platforms, a profound enabling capability for accelerating MESA 's development.  \nKeywords: stars: evolution - methods: numerical - binaries: general - stars: oscillations - nuclear reactions - shock waves - di GLYPH<11> usion", 'Contents': '5.1. More Accurate Summations Yield More', '11. Summary and Conclusions': '41', '1. INTRODUCTION': "The development of a relatively complete and quantitative picture of stellar evolution is one of the great drivers of astrophysics. On the observational side of this impetus, the Kepler (Borucki et al. 2010) and the CoRoT (Baglin et al. 2009) missions continuously monitored more than 100,000 stars in a 100 deg 2 window with a dynamic range of apparent stellar brightness of 10 6 . Highlights include the discoveries that nearly all GLYPH<13> Doradus and GLYPH<14> Scuti stars are hybrid pulsators, and the detection of solar-like oscillations in a large sample of red giants (Auvergne et al. 2009; De Ridder et al. 2009; Grigahc'ene et al. 2010; Bedding et al. 2010; ChristensenDalsgaard & Thompson 2011; Chaplin & Miglio 2013). The Dark Energy Survey is scanning 5000 deg 2 of the southern sky in 5 optical filters every few days to discover and study thousands of supernovae (e.g., Papadopoulos et al. 2015; Yuan et al. 2015). Building upon the legacy of the Palomar Transient Factory (Law et al. 2009), the intermediate Palomar Transient Factory conducts a fully-automated, wide-field sur-  \nvey that systematically explores the transient sky with a 90 second to 5 day cadence (Vreeswijk et al. 2014). The forthcoming Zwicky Transient Facility will enable a survey more than an order of magnitude faster at the same depth as its predecessors. In its unique orbit the Transiting Exoplanet Survey Satellite will have an unobstructed view to scrutinize the light curves of the brightest 100,000 stars with a 1 minute cadence (Ricker et al. 2015). The Gaia mission aims to provide unprecedented distance and radial velocity measurements with the accuracies needed to reveal the evolutionary state, composition, and kinematics of about one billion stars in our Galaxy (e.g., Creevey et al. 2015; Sacco et al. 2015). The Large Synoptic Survey Telescope will image the entire Southern Hemisphere deeply in multiple optical colors every week with its three billion pixel digital camera, thus opening a new window on transient objects such as interacting close binary systems. (e.g., Oluseyi et al. 2012).  \nIt has been a little more than 200 years since Herschel (1802) announced, after 25 years of observation, that certain pairs of stars displayed evidence of orbital motion around their common center of mass. Binary systems allow the masses of their component stars to be directly determined, which in turn allows stellar radii to be indirectly estimated. This allows the calibration of an empirical mass-luminosity relationship from which the masses of single stars can be estimated (Torres et al. 2010). Recent surveys such as Raghavan et al. (2010) suggest 30% to 50% of solar-like systems in the Galactic disk are composed of binaries, where the binary fraction is higher for more massive stars (Sana et al. 2012; Kobulnicky et al. 2014). As argued by de Mink et al. (2013), the most rapidly rotating massive stars are expected in binary systems as a consequence of accretion-induced spin-up. Di GLYPH<11> erential rotation has a major impact on the evolution of massive stars (Maeder & Meynet 2000; Heger et al. 2000, 2005) and for single stars the corresponding physics has been included into MESA as described in Paper II. On the other hand very few works that include the physics of di GLYPH<11> erential rotation in binaries have been published (Petrovic et al. 2005; Cantiello et al. 2007). Our improvements to MESA now allow for the calculation of di GLYPH<11> erentially-rotating binary stars.  \nInterpreting these new observations and predicting new stellar phenomena propels the theoretical side, in particular the evolution of the community software instrument Modules for Experiments in Stellar Astrophysics ( MESA ) for research and education. We introduced MESA in Paxton et al. (2011, hereafter Paper I) and expanded its range of capabilities in Paxton et al. (2013, hereafter Paper II). This paper describes the major new advances for MESA modeling of binary systems, shock hydrodynamics, explosions of massive stars and x-ray bursts with large, in situ reaction networks. Moreover it details the coupling of MESA with the non-adiabatic pulsation software instrument GYRE (Townsend & Teitler 2013). We also describe advances made to existing modules since Paper II, including improved treatments of mass accretion, weak reaction rates, and particle di GLYPH<11> usion.  \nThe rapid expansion of extra-solar planet research has led to a revival of interest in the detailed properties of stars probed through space-based brightness variability studies and radial velocity measurements. Stellar properties can be derived from measurements of the radial and non-radial oscillation modes of a star, but this requires the accurate and e GLYPH<14> cient computations of mode frequencies and their eigenfunctions enabled by the coupling of GYRE with MESA .  \nThere are many ways M & 8 M GLYPH<12> stars can end their lives  \n(e.g., Woosley et al. 2002; Smartt 2009; Meynet et al. 2010; Langer 2012; Nomoto et al. 2013; Smith 2014). Some become electron capture supernovae; others collapse with most of their extended envelope intact and yield a Type II supernova; others can lead to pair instability; and some have envelopes thin enough to allow a jet to break through and appear as a long gamma-ray burst (MacFadyen & Woosley 1999; Woosley & Bloom 2006; Gehrels et al. 2009). There is a pressing need throughout the stellar community to routinely explore this entire mass range with new supernova progenitor and explosion models. The observational facilities discussed above have found explosions that indicate large amounts of mass are lost within a few years of explosion (Smith 2014); some show evidence of optically thick winds present at the moment of explosion (Ofek et al. 2014), while others have yet to be securely identified with a specific core collapse scenario. These mysteries, coupled with the community's call for new yields from massive stars for galactic chemical evolution studies motivate the development of implicit shock hydrodynamics and explosions with large, in situ reactions networks in MESA .  \nThe paper is outlined as follows. Section 2 describes the new capability of MESA to evolve binary systems. Section 3 discusses the new non-adiabatic pulsation capabilities resulting from fully coupling to the GYRE software instrument. Section 4 describes the improvements to accommodate implicit hydrodynamics with shocks. New capabilities for advanced burning and X-ray bursts with large, in situ reaction networks are described in Section 5. In Section 6 we model the pre-supernova evolution of massive stars and combine the implicit hydrodynamics module and the new capabilities for advanced burning to probe the nucleosynthesis and yields of core-collapse supernovae. Section 7 discusses the improvements for a more robust and e GLYPH<14> cient treatment of mass accretion. Section 8 presents a new option for an on-the-fly calculation of weak reaction rates and their application to the Urca process and accretion-induced-collapse models. Section 9 presents improvements in the physics implementation of particle di GLYPH<11> usion by including radiative levitation and pushing di GLYPH<11> usion methods into the strongly coupled, electron degenerate regime. In Section 10 we discuss improvements to the MESA software infrastructure, highlighting bit-for-bit consistency across operating systems and compilers. We conclude in Section 11 by noting additional improvements to MESA are likely to occur in the near future. Important 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.  \nTable 1 -Continued", '2. BINARIES': 'MESAbinary is a MESA module that uses MESAstar to evolve binary systems. It can be used to evolve a full stellar model plus a point mass companion or to simultaneously evolve the structure of two stars. It optionally allows the modeling of systems including stellar rotation, assuming the axis of rotation of each star to be perpendicular to the orbital plane, accounting for the e GLYPH<11> ects of tidal interaction and spin-up through accretion. The implementation of MESAbinary benefits from early contributions by Madhusudhan et al. (2006) and Lin et al. (2011) who modeled mass transfer from a star to a point mass.  \nHere we provide an overview of the modelled physical processes for circular binary systems and describe the tests against which we validate MESAbinary .', '2.1. Initialization of a Circular Binary System': 'Abinary system is initialized by specifying the components and either the orbital period P orb or separation a . Each component can be a point mass or a stellar model. The initial model(s) are provided by a saved MESA model or a zero-age main-sequence (ZAMS) specification. For stellar models including rotation, the initial rotational velocities of each component can be explicitly defined, or set such that the star is synchronized to the orbit at the beginning of the evolution. The orbital angular momentum of the system is  \nJ orb = M 1 M 2 r Ga M 1 + M 2 ; (1)  \nwhere M 1 and M 2 are the stellar masses. Evolution of M 1, M 2, and J orb is used to update a using Equation (1). Masses can be modified both by Roche lobe overflow (RLOF) and winds. The total time derivatives of the component masses are given by  \n˙ M 1 = ˙ M 1 ; w + ˙ M RLOF ; ˙ M 2 = ˙ M 2 ; w GLYPH<0> f mt ˙ M RLOF ; (2)  \nwhere M 1 is the donor mass and M 2 the accretor mass. The stellar wind mass loss rates are ˙ M 1 ; w and ˙ M 2 ; w (see Paper I and Paper II) and ˙ M RLOF is the mass transfer rate from RLOF, all defined as negative. The factor f mt represents the e GLYPH<14> ciency of accretion and can be used to limit accretion to the Eddington rate ˙ M Edd.', '2.2. Evolution of Orbital Angular Momentum': 'To compute the rate of change of orbital angular momentum, we consider the contribution of gravitational waves, mass loss, magnetic braking, and spin orbit (LS) coupling  \n˙ J orb = ˙ J gr + ˙ J ml + ˙ J mb + ˙ J ls ; (3)  \nfrom which the change in orbital angular momentum in one step is calculated as GLYPH<1> J orb = ˙ J orb GLYPH<14> t , where GLYPH<14> t is the timestep. Unless models with stellar rotation are being used, the ˙ J ls term is equal to zero, and the contribution of the individual spins of each star is not directly considered. On the other hand, the ˙ J mb term implicitly assumes a strong tide that keeps the orbit synchronized. The simultaneous usage of ˙ J mb with stellar rotation is not consistent (see Section 2.2.4). We now describe how these terms are computed.', '2.2.1. Gravitational Wave Radiation': 'Very compact binaries can experience significant orbital decay due to the emission of gravitational waves. Observations  \nof the Hulse-Taylor binary pulsar over three decades (Weisberg & Taylor 2005) and of the double pulsar (Kramer et al. 2006) have tested the predicted e GLYPH<11> ect from general relativity to a high precision. The angular momentum loss from gravitational waves is  \n˙ J gr = GLYPH<0> 32 5 c 5 2 GLYPH<25> G P orb ! 7 = 3 ( M 1 M 2) 2 ( M 1 + M 2) 2 = 3 : (4)  \n2.2.2. Mass Loss  \nWe assume the mass lost in a stellar wind has the specific orbital angular momentum of its star. For the case of ine GLYPH<14> -cient mass transfer, angular momentum loss follows Soberman et al. (1997), where fixed fractions of the transferred mass are lost either as a fast isotropic wind from each star or a circumbinary toroid with a given radius:  \n˙ J ml = h ( ˙ M 1 ; w + GLYPH<11> mt ˙ M RLOF) M 2 1 + ( ˙ M 2 ; w + GLYPH<12> mt ˙ M RLOF) M 2 2 i GLYPH<2> a 2 ( M 1 + M 2) 2 2 GLYPH<25> P orb + GLYPH<13> mt GLYPH<14> mt ˙ M RLOF p G ( M 1 + M 2) a ; (5)  \nwhere GLYPH<11> mt, GLYPH<12> mt, and GLYPH<14> mt are respectively the fractions of mass transferred that is lost from the vicinity of the donor, accretor and circumbinary toroid, and GLYPH<13> 2 mt a is the radius of the toroid. Ignoring winds, the e GLYPH<14> ciency of mass transfer is then given by f mt = 1 GLYPH<0> GLYPH<11> mt GLYPH<0> GLYPH<12> mt GLYPH<0> GLYPH<14> mt. When accretion is limited to ˙ M Edd, e GLYPH<14> ciency of accretion is given by  \nf mt = min(1 GLYPH<0> GLYPH<11> mt GLYPH<0> GLYPH<12> mt GLYPH<0> GLYPH<14> mt ; j ˙ M Edd = ˙ M RLOF j ) ; (6)  \nand the additional mass being lost is added to the GLYPH<12> mt ˙ M RLOF term in Equation (5), i.e., it is assumed to leave the system carrying the specific orbital angular momentum of the accretor.', '2.2.3. Spin Orbit Coupling': 'Tidal interaction and mass transfer can significantly modify the spin angular momentum of the stars in a binary system, acting as both sources and sinks for orbital angular momentum. The impact spin-orbit interactions have on orbital evolution depends on the orbital separation and the mass ratio, with the e GLYPH<11> ect being greater for tighter orbits and uneven masses. The corresponding contribution to ˙ J orb is computed by demanding conservation of the total angular momentum, accounting for losses due to the other ˙ J orb mechanisms and loss of stellar angular momentum due to winds.  \nIn a fully conservative system, the change in orbital angular momentum in one timestep is GLYPH<14> J orb = GLYPH<0> GLYPH<14> S 1 GLYPH<0> GLYPH<14> S 2, where GLYPH<14> S 1 and GLYPH<14> S 2 are the changes in spin angular momenta. This needs to be corrected if mass loss is included, as winds take away angular momentum from the system. If S 1 ; lost and S 2 ; lost are the amounts of spin angular momentum removed in a step from each star due to mass loss (including winds and RLOF),  \n˙ J ls = GLYPH<0> 1 GLYPH<14> t GLYPH<14> S 1 GLYPH<0> S 1 ; lost ˙ M 1 ; w ˙ M 1 + GLYPH<14> S 2 GLYPH<0> S 2 ; lost ! ; (7)  \nwhere the additional factor for the donor accounts for mass lost from the system, ignoring mass loss due to mass transfer. In the absence of RLOF this equation becomes symmetric between both stars, as then ˙ M 1 ; w = ˙ M 1 = 1.  \nThe form of Equation (7) is independent of how tides and angular momentum accretion work, as it is merely a statement on angular momentum conservation. The details of how we model these processes and their impact on the spin of each component are described in Section 2.4.', '2.2.4. Magnetic Braking': 'The rotational velocities of low mass stars are strongly correlated with their ages (Skumanich 1972). This spin-down arises from the coupling of the stellar wind to a magnetic field. If the star is in a binary system and tidally coupled to the orbit, magnetic braking can provide a very e GLYPH<14> cient sink for orbital angular momentum (Mestel 1968; Verbunt & Zwaan 1981). We implement this e GLYPH<11> ect following Rappaport et al. (1983), who assumed the star being braked is tidally synchronized:  \n˙ J mb = GLYPH<0> 6 : 82 GLYPH<2> 10 34 M 1 M GLYPH<12> ! R 1 R GLYPH<12> ! GLYPH<13> mb 1 d P orb ! 3 [dyn cm] ; (8)  \nwhere in the simplest approximation GLYPH<13> mb = 4 (Verbunt & Zwaan 1981). A similar contribution from the accretor can be included. As tidal synchronization is assumed, this formulation is incompatible with the use of LS coupling.  \nIt is normally assumed that once a star becomes fully convective, the dynamo process that regenerates the field will stop working or at least behave in a qualitatively di GLYPH<11> erent manner. Similarly, magnetic fields in stars with radiative envelopes are of a significantly di GLYPH<11> erent nature than those seen in stars with convective envelopes, and there is no simple way to predict even the presence of magnetism (Donati & Landstreet 2009). By default MESAbinary only accounts for magnetic braking as long as the star being braked has a convective envelope and a radiative core, though the process might still operate outside of these conditions.', '2.3. Mass Transfer from RLOF': 'Close binary stars are defined as systems tight enough to interact through mass transfer, with the most important mechanism being RLOF. This process is commonly modeled in 1D by considering the spherical-equivalent Roche lobe radius R RL of each object (Eggleton 1983)  \nR RL ; j = 0 : 49 q 2 = 3 j 0 : 6 q 2 = 3 j + ln GLYPH<16> 1 + q 1 = 3 j GLYPH<17> a ; (9)  \nwhere j is the index identifying each star, q 1 = M 1 = M 2 and q 2 = M 2 = M 1. This fit is correct up to a few percent for the full range of mass ratios, 0 < qj < 1 . Mass transfer occurs then when the radius of a star approaches or exceeds R RL. Depending on several factors, once a star begins RLOF the ensuing mass transfer phase can proceed on a nuclear, thermal, or dynamical timescale.  \nThe stability of mass transfer is normally understood in terms of mass-radius relationships (e.g., Tout et al. 1997; Soberman et al. 1997),  \nGLYPH<16> eq = d ln R 1 d ln M 1 ! eq ; (10)  \nGLYPH<16> ad = d ln R 1 d ln M 1 ! ad ; (11)  \nGLYPH<16> RL = d ln R RL ; 1 d ln M 1 : (12)  \nHere, GLYPH<16> eq gives the radial response of the donor to mass loss when it happens slowly enough for the star to remain in thermal equilibrium. When mass loss proceeds on a timescale much shorter than the thermal timescale of the star, but still slow enough for the star to retain hydrostatic equilibrium then the radial response will be given by GLYPH<16> ad. The dependency of the Roche lobe radius on mass transfer is encoded in GLYPH<16> RL. In general GLYPH<16> = d ln R 1 = d ln M 1 is a function of ˙ M RLOF, so requiring GLYPH<16> = GLYPH<16> RL will determine the value of ˙ M RLOF. If an overflowing star satisfies GLYPH<16> eq > GLYPH<16> RL, then it can remain inside its Roche lobe by transferring mass while retaining thermal equilibrium. If on the contrary GLYPH<16> ad > GLYPH<16> RL > GLYPH<16> eq, mass transfer will proceed on a thermal timescale, while for the extreme case GLYPH<16> RL > GLYPH<16> ad the star will depart from hydrostatic equilibrium and the process will be dynamical. MESA cannot model common envelope or contact binaries.  \nMESAbinary provides both explicit and implicit methods to compute mass transfer rates. An explicit computation sets the value of ˙ M RLOF at the start of a step, while an implicit one begins with a guess for ˙ M RLOF and iterates until the required tolerance is reached. The composition of accreted material is set to that of the donor surface, and the specific entropy of accreted material is the same as the surface of the accretor. In models with rotation the specific angular momentum of accreted material is described in Section 2.4.', '2.3.1. Explicit Methods': 'MESAbinary implements two mass transfer schemes: the model of Ritter (1988) which we refer to as the Ritter scheme and Kolb & Ritter (1990) which we refer to as the Kolb scheme. We use the mass ratio q 2 consistent with the Ritter scheme.  \nRitter scheme: - Stars have extended atmospheres therefore RLOF can take place through the L1 point even when R 1 < R RL ; 1. Ritter (1988) estimated the mass transfer rate for this case as  \n˙ M RLOF = GLYPH<0> ˙ M 0 exp R 1 GLYPH<0> R RL ; 1 HP ; 1 =GLYPH<13> ( q 2) ! ; (13)  \nwhere HP ; 1 is the pressure scale height at the photosphere of the donor and  \n˙ M 0 = 2 GLYPH<25> exp(1 = 2) F 1( q 2) R 3 RL ; 1 GM 1 k B T e GLYPH<11> m p GLYPH<22> ph ! 3 = 2 GLYPH<26> ph ; (14)  \nwhere m p is the proton mass, T e GLYPH<11> is the e GLYPH<11> ective temperature of the donor, and GLYPH<22> ph and GLYPH<26> ph are the mean molecular weight and density at its photosphere. The two fitting functions are  \nF 1( q 2) = 1 : 23 + 0 : 5 log q 2 ; 0 : 5 . q 2 . 10 : (15)  \nand  \nGLYPH<13> ( q 2) = ( 0 : 954 + 0 : 025 log q 2 GLYPH<0> 0 : 038(log q 2) 2 ; 0 : 04 . q 2 GLYPH<20> 1 0 : 954 + 0 : 039 log q 2 + 0 : 114(log q 2) 2 ; 1 GLYPH<20> q 2 . 20 : (16)  \nOutside the ranges of validity F 1( q 2) and GLYPH<13> ( q 2) are evaluated using the value of q 2 at the edge of their respective ranges.  \nKolb scheme: - Kolb & Ritter (1990) extended the Ritter scheme in order to cover the case R 1 > R RL ; 1 according to  \n˙ M RLOF = GLYPH<0> ˙ M 0 GLYPH<0> 2 GLYPH<25> F 1( q 2) R 3 RL ; 1 GM 1 GLYPH<2> Z P RL P ph GLYPH<0> 1 = 2 1 2 GLYPH<0> 1 + 1 ! ( GLYPH<0> 1 + 1) = (2 GLYPH<0> 1 GLYPH<0> 2) k B T m p GLYPH<22> ! 1 = 2 dP (17)  \nwhere GLYPH<0> 1 is the first adiabatic exponent, and P ph and P RL are respectively the pressures at the photosphere and at the radius for which r 1 = R RL ; 1.', '2.3.2. Implicit Methods': 'Explicit schemes exhibit large jumps in ˙ M RLOF unless the timestep is severely restricted. Therefore, if one needs accurate values of ˙ M RLOF and stellar radius, this requires use of an implicit scheme. Implicit schemes also allow the calculation these quantities when there is no general closed form formula for ˙ M RLOF.  \nThese implicit methods use a bisection-based root solve to satisfy j f ( ˙ M RLOF) j < GLYPH<24> at the end of the step, where GLYPH<24> is a given tolerance. The implicit schemes are then defined by the choice of the function f ( ˙ M RLOF). For the Ritter and the Kolb scheme the function is chosen as  \nf ( ˙ M RLOF) = ˙ M end GLYPH<0> ˙ M RLOF ˙ M end ; (18)  \nwith ˙ M end being the mass transfer rate computed at the end of each iteration.  \nA di GLYPH<11> erent implicit method is also provided. In this case, whenever the donor star overflows its Roche lobe the implicit solver will adjust the mass transfer rate until R 1 = R RL ; 1 within some tolerance (see e.g., Whyte & Eggleton 1980; Rappaport et al. 1982, 1983). In this case  \nf ( ˙ M RLOF) = 2( R 1 GLYPH<0> R RL ; 1) R RL ; 1 + GLYPH<24> ; (19)  \nand if ˙ M RLOF is below a certain threshold and f ( ˙ M RLOF) < GLYPH<0> GLYPH<24> then the system is assumed to detach and ˙ M RLOF is set to zero.', '2.4. E GLYPH<11> ect of Tides and Accretion on Stellar Spin': "To model tidal interaction we adjusted the model of Hut (1981) to include the case of di GLYPH<11> erentially rotating stars. The time evolution of the angular frequency for each component is  \nd GLYPH<10> i ; j dt = GLYPH<10> orb GLYPH<0> GLYPH<10> i ; j GLYPH<28> sync ; j ; 1 GLYPH<28> sync ; j = 3 ( qjr g ; j ) 2 k T ! c ; j Rj a ! 6 ; (20)  \nwhere j = 1 ; 2 is the index of each star, GLYPH<10> i ; j is the angular frequency at the face of cell i towards the surface, r 2 g ; j = Ij = ( MjR 2 j ) is the radius of gyration (with Ij being the moment of inertia of each star), and the ratio of the apsidal motion constant to the viscous dissipation timescale, ( k = T )c ; j , is computed as in Hurley et al. (2002). Similarly to Detmers et al. (2008), we assume constant GLYPH<28> sync ; j and GLYPH<10> orb through a step and therefore GLYPH<1>GLYPH<10> i ; j = [1 GLYPH<0> exp( GLYPH<0> GLYPH<14> t =GLYPH<28> sync ; j )]( GLYPH<10> orb GLYPH<0> GLYPH<10> i ; j ). This extension of Hut's work to di GLYPH<11> erentially rotating stars is not formally derived but merely applies his result for solid body rotators independently to each shell. The formulation of Hut (1981) can be recovered from Equation (20), by forcing  \nsolid body rotation with a large di GLYPH<11> usion coe GLYPH<14> cient for angular momentum throughout the star. In reality tides would act mostly on the outer layers, and whether the core synchronizes or not depends on the coupling between the core and the envelope.  \nTo compute the specific angular momentum carried by accreted material, we consider the possibility of both ballistic and Keplerian disk mass transfer (e.g., Marsh et al. 2004; de Mink et al. 2013). To distinguish which occurs, we compare the minimum distance of approach of the accretion stream (Lubow & Shu 1975; Ulrich & Burger 1976) 11  \nR min = 0 : 0425 a GLYPH<16> q 2 + q 2 2 GLYPH<17> 1 = 4 ; 0 : 0667 GLYPH<20> q 2 GLYPH<20> 15 (21)  \nto the radius of the accreting star. When outside the range of validity, R min is computed using the value of q 2 at the respective edge. Accretion is assumed to be ballistic whenever R 2 > R min and the specific angular momentum is (1 : 7 GM 2 R min) 1 = 2 . When R 2 < R min the specific angular momentum is taken as that of a Keplerian orbit at the surface ( GM 2 R 2) 1 = 2 .", '2.5. Treatment of Thermohaline Mixing in Accreting Models': 'In stars with radiative envelopes accreted material with a high mean molecular weight is expected to mix inwards due to thermohaline mixing, a process that is very sensitive to the GLYPH<22> -gradient (see e.g., Kippenhahn et al. 1980; Cantiello & Langer 2010). Thermohaline mixing is included in MESA (see Paper I). However, as mass with homogeneous composition is added during the accretion process, a jump is produced at the boundary between new and old material. MESAstar computes mixing coe GLYPH<14> cients explicitly at the start of each step, so this results in thermohaline mixing only operating near this boundary, leading to unphysical compositional staircases. To avoid this issue, we artificially soften the composition gradient in the outer ( GLYPH<1> q )large fraction of the star by mass. We do this starting at the surface and homogeneously mixing inwards a region of size ( GLYPH<1> q )small. Then, moving towards the center, the process is repeated at each cell while linearly (with respect to mass) reducing the size of the small mixed region such that it is zero after going ( GLYPH<1> q )large inwards. All the binary models where the accretor is not a point mass are calculated using ( GLYPH<1> q )large = 0 : 05 and ( GLYPH<1> q )small = 0 : 03.', '2.6. Numerical Tests': 'Here we describe tests designed to validate the implementation of the physics described in Section 2.2. We check orbital evolution in the presence of gravitational waves and mass loss by comparing to analytical solutions. We also verify total angular momentum conservation in calculations that include the physics of tides and spinup by accretion. To test for the thermal response of stellar models undergoing mass transfer, we compare MESAbinary results to those from the STARS code (Eggleton 1971; Pols et al. 1995; Stancli GLYPH<11> e & Eldridge 2009).', '2.6.1. Gravitational Wave Radiation': 'If gravitational waves are the only source of angular momentum loss and the masses of each component remain constant, Equation (3) can be integrated to obtain the time evolution of orbital separation (Peters 1964). We model a system consisting of a 0 : 5 M GLYPH<12> star and a 0 : 8 M GLYPH<12> point mass with  \n11 Note that there is a small typo in the fit given by Ulrich & Burger (1976). The corrected fit given here fits the results of Lubow & Shu (1975) to the 4% accuracy claimed by Ulrich & Burger (1976).  \nFigure 1. Evolution in the HR diagram for a 2 : 5 M GLYPH<12> star transferring mass to a 1 : 4 M GLYPH<12> point mass, assuming a mass transfer e GLYPH<14> ciency of 1%. Symbols are shown at zero-age main sequence (ZAMS) and terminal-age main sequence (TAMS), together with parts of the track where RLOF is occurring. The inset shows evolution from ZAMS up to the beginning of the first phase of mass transfer.  \n<!-- image -->  \na = 2 R GLYPH<12> . We ignore all e GLYPH<11> ects on the evolution of orbital angular momentum except its loss due to gravitational waves. In 3 : 5 Gyr the orbital separation of this system reduces to a = 1 : 3 R GLYPH<12> , at which point the 0 : 5 M GLYPH<12> star begins mass transfer. We terminate the run at the onset of RLOF. The maximum error in a is 0 : 35% relative to the analytical result.', '2.6.2. Ine GLYPH<14> cient Mass Transfer': 'An analytical expression for the evolution of orbital separation can be derived if ine GLYPH<14> cient mass transfer is the only contribution to the angular momentum evolution (Tauris & van den Heuvel 2006). We model a 2 : 5 M GLYPH<12> main sequence star together with a 1 : 4 M GLYPH<12> point mass with an initial orbital separation of 10 R GLYPH<12> . We choose GLYPH<11> mt = 0 : 03, GLYPH<12> mt = 0 : 95, GLYPH<14> mt = 0 : 01 and GLYPH<13> 2 mt = 2, which give a low mass transfer e GLYPH<14> -ciency of f mt = 0 : 01. Such a system is representative of the evolution of an intermediate mass X-ray binary (IMXB). The model initiates mass transfer just after the end of the main sequence, interrupting the evolution of the star through the Hertzsprung gap and producing a low mass white dwarf (WD) ( M He = 0 : 289 M GLYPH<12> ) with a small amount of hydrogen on its surface. As the WD evolves to the cooling track, it experiences several hydrogen flashes, one of them strong enough to produce an additional phase of RLOF (see Figure 1).  \nFigure 2 shows that MESAbinary computes the orbital evolution to a precision of a few parts in 10 4 . We run this system using both the Ritter and the Kolb implicit schemes to display that under some circumstances the precise choice of mass transfer scheme does not play a big role in the evolution.', '2.6.3. Spin Orbit Coupling': 'We now test angular momentum conservation by ignoring all the mechanisms that remove angular momentum from the  \nFigure 3 shows that spin angular momentum in both components increases during the pre-interaction phase, which is due to both stars expanding on the main sequence while remaining tidally locked. During Roche lobe overflow, the secondary is rapidly spun-up, reaching nearly 80% of critical rotation before contact. The calculation of total angular momentum requires the summation of di GLYPH<11> erent contributions (orbital angular momentum and spin of both components). Therefore the maximum accuracy to which we can conserve angular momentum is limited by rounding errors. Figure 3 shows that conservation of angular momentum in the run is very close to machine precision.  \n<!-- image -->  \n/circledot  \nFigure 2. Evolution of mass transfer rate from a 2 : 5 M GLYPH<12> to a 1 : 4 M GLYPH<12> point mass, assuming a mass transfer e GLYPH<14> ciency of 1%. The upper panel shows the di GLYPH<11> erence between the computed orbital separation and the analytical solution while the bottom one displays the evolution of the mass transfer rate, using two di GLYPH<11> erent schemes.  \nbinary system. For this purpose we model an 8 M GLYPH<12> + 6 M GLYPH<12> binary with rotating components and an initial orbital period of 1 : 5 days. Due to the short orbital separation we assume the initial spin periods of the two stars are equal to the orbital period. The primary undergoes RLOF during the main sequence, initiating a phase of mass transfer on a thermal timescale. After transferring just 0 : 3 M GLYPH<12> the accretor also fills its Roche lobe, producing a contact system. At this point we terminate the evolution.', '2.6.4. Thermal Response to Mass Loss': 'The fate of binary systems depends largely on the precise value of ˙ M during an interaction phase, which depends on the thermal response of the donor star to mass loss. For WDs there is a limited range of accretion rates for stable hydrogen burning (Nomoto et al. 2007; Shen & Bildsten 2007). In main sequence binaries the evolution of the accretor radius depends on the mass transfer rate, and expansion during the interaction phase can lead to contact or even a merger (Wellstein et al. 2001).  \nWe calculated an 8 M GLYPH<12> + 6 : 5 M GLYPH<12> binary system with an initial orbital period of 1 : 5 days using both MESAbinary and  \nFigure 3. Angular momentum evolution in an 8 M GLYPH<12> + 6 M GLYPH<12> binary with an initial orbital period of 1 : 5 days. Left panels show the evolution before the onset of RLOF, while right panels display evolution from the beginning of RLOF until contact, when both components fill their Roche lobe. The fractional error in the total angular momentum is plotted in the bottom panel and is of order machine-precision.  \n<!-- image -->  \nSTARS . To minimize the modeling di GLYPH<11> erences and focus on the thermal response of both components, we use an extremely simplified model that ignores internal mixing (including convective mixing). Under these conditions, the more massive star quickly depletes its central hydrogen and begins shell hydrogen burning, reaching RLOF and undergoing a phase of mass transfer on the thermal timescale. The resulting mass transfer rates are shown in Figure 4. The agreement is very good, despite mass transfer rates being computed in slightly di GLYPH<11> erent ways. Masses at detachment show a small di GLYPH<11> erence, with the MESAbinary model terminating mass transfer when M 1 = 0 : 952 M GLYPH<12> while the STARS calculation when M 1 = 0 : 935 M GLYPH<12> . The figure also shows the change in radius of the accreting star, with two prominent peaks at R 2 = R GLYPH<12> = 4 : 84 ; 5 : 34 for MESAbinary and 4 : 82 ; 5 : 28 for STARS . The larger radius of the MESAbinary model is likely associated to the slightly higher mass transfer rates.', '2.7. Period Gap of Cataclysmic Variables': "Although cataclysmic variables (CVs) span a wide range of periods, observations show a lack of systems in the range 2 h < P orb < 3 h (see, for instance, Gansicke et al. 2009). Such a feature is commonly explained by having an angular momentum loss mechanism 'turn o GLYPH<11> ' or become ine GLYPH<14> cient at some point. The most popular model for such a mechanism is magnetic braking (Rappaport et al. 1983), as the magnetic field of the donor is assumed to change quickly when the star loses enough mass to become fully convective.  \nWe now compare to the results of Howell et al. (2001), who performed a population synthesis study to explore in detail the standard scenario involving magnetic braking. In Figure 5 we show the evolution of mass transfer rates and orbital periods for a set of CV models with di GLYPH<11> erent component masses and  \n<!-- image -->  \n/circledot  \nFigure 4. Mass transfer rate and accretor radius as computed by MESA and STARS for an 8 M GLYPH<12> + 6 : 5 M GLYPH<12> binary with an initial orbital period of 3 days. All internal mixing processes (including convective mixing) are turned o GLYPH<11> in the calculations.Figure 5. Evolution of CV models under the e GLYPH<11> ect of magnetic braking and gravitational wave radiation. For each track the label gives the donor mass, the WD mass, and the initial orbital period respectively. The grey band shows the observed period gap for CVs. These results reproduce figure 1 in Howell et al. (2001).  \n<!-- image -->  \norbital periods. We run all models using GLYPH<12> mt = 1 and GLYPH<13> mb = 3 and magnetic braking is turned o GLYPH<11> when the donor star becomes fully convective. As an example the system with a 0 : 9 M GLYPH<12> donor (left panel in Figure 5) experiences a first phase of mass transfer induced by magnetic braking between 10 7 : 1 and 10 8 : 3 yr, a non-interacting phase (the gap) between 10 8 : 3 and 10 8 : 8 yr, and a subsequent phase of mass transfer dominated by gravitational wave radiation, reaching a minimum orbital period of about 1 hour at 10 9 : 6 yr. As a comparison, for the same model Howell et al. (2001) obtain a first phase of mass transfer between 10 7 : 3 and 10 8 : 4 yr, the gap occurs between 10 8 : 4 and 10 8 : 8 yr and a period minimum is reached at 10 9 : 4 yr. Figure 5 shows that our CV models spend most time away from the observed period gap.", '2.8. Evolution of Massive Binaries': "/circledot  \nFigure 6. Evolution of a 16 M GLYPH<12> + 14 M GLYPH<12> system with a 3 day initial orbital period. MESAbinary models are compared to the results of Wellstein et al. (2001), which were calculated using the STERN code. The terms primary and secondary are used throughout the evolution to describe the initially more massive and the less massive components, respectively. For each component in the MESAbinary model, squares mark the ZAMS and the depletion of the indicated nuclear fuel in the core.  \n<!-- image -->  \nIn massive stars, binary interactions have dramatic e GLYPH<11> ects on the evolution of both components. Kippenhahn & Weigert (1967) introduced the term 'case A' to refer to a mass transfer phase occurring in systems tight enough such that RLOF starts during the main sequence. This results in a large amount of mass being transferred on a thermal timescale, followed by a phase of mass transfer that proceeds on the nuclear timescale until the end of core H-burning. An additional phase of thermal-timescale mass transfer then follows (the so-called 'case AB'), which strips the donor and produces an almostnaked helium star.  \nFigure 7 illustrates the prevalence of both thermohaline mixing and semiconvection in the accreting star. Newly accreted material is e GLYPH<14> ciently mixed inwards by thermohaline mixing. On the other hand the GLYPH<22> -gradient formed before interaction prevents the convective core from growing, with the e GLYPH<14> ciency of semiconvection controlling whether or not the star rejuvenates. Due to the choice of ine GLYPH<14> cient semiconvection, the core remains small, preventing the star from becoming a red supergiant. The star accretes a large amount of CNO-processed and helium-rich material. After being mixed through the envelope this material results in the surface being nitrogen rich and carbon depleted, with a slight enhancement in helium.  \nHere we show that MESAbinary can calculate the evolution of massive interacting binaries. We reproduce one of the models from Wellstein et al. (2001), a 16 M GLYPH<12> + 14 M GLYPH<12> system with an initial period of 3 days, using the same semiconvection e GLYPH<14> -ciency of GLYPH<11> sc = 0 : 01. As shown in Figures 6 and 7 this system experiences case A and AB mass transfer, and the accretor becomes a blue supergiant after core hydrogen depletion. The accretor depletes carbon before its donor.  \nFigure 7. Kippenhahn diagram for the evolution of a 16 M GLYPH<12> + 14 M GLYPH<12> system with a 3 day initial orbital period. Most of the pre-interaction phase is not shown in this figure. The upper plot shows the evolution of the donor, while the lower plot displays that of the accretor.  \n<!-- image -->", '2.9. Rotating Binaries and the E GLYPH<14> ciency of Mass Transfer': 'The e GLYPH<14> ciency of mass transfer plays a key role in close binary systems, but the processes by which material is lost from the system are not well-understood. In particular, whenever an accreting star approaches GLYPH<10> = GLYPH<10> crit = 1, it is uncertain whether accretion can continue, one option being the development of a strong wind that prevents accretion (e.g., Petrovic et al. 2005; Cantiello et al. 2007). Whenever GLYPH<10> = GLYPH<10> crit approaches one, we use an implicit method to iteratively reduce ˙ M 2 until this ratio falls below a threshold.  \nFigure 8 shows that MESAbinary models using both the Zahn (1977) and Hurley et al. (2002) timescales for tidal coupling. These models experience highly non-conservative phases of mass transfer, corresponding to the accreting star evolving very close to critical rotation. In particular during case AB mass transfer the accretor needs to switch from mass accretion to mass loss in order to remain sub-critical. As expected, the system with the tidal timescale from Zahn (1977) has a significantly higher mass transfer e GLYPH<14> ciency, and during the first phase of RLOF it only experiences a brief period in which the accretor reaches critical rotation. This is in broad agreement with the model by Langer et al. (2003).  \nTides counteract the e GLYPH<11> ect of spin-up from accretion. Whether or not an accreting object reaches critical rotation depends on the e GLYPH<14> ciency of tidal coupling. Here we model a 16 M GLYPH<12> + 15 M GLYPH<12> binary system including di GLYPH<11> erential stellar rotation, with an initial orbital period of 3 days and assuming initial orbital synchronization. Langer et al. (2003) argue that turbulent processes in the radiative envelope can significantly enhance tidal strength. They model the same system using the simple estimate for the synchronization timescale for a star with a convective envelope given by Zahn (1977), GLYPH<28> sync ; j = 1 yr GLYPH<2> q 2 j ( a = R ) 6 . For our implicit modeling of stellar winds we use a threshold of ( GLYPH<10> = GLYPH<10> crit) = 0 : 99.  \nFigure 8. E GLYPH<14> ciency of mass transfer in a 16 M GLYPH<12> + 15 M GLYPH<12> binary system including di GLYPH<11> erential rotation. The system is modeled with tides as described by Hurley et al. (2002) for radiative envelopes, and also with the simple tidal timescale given by Zahn (1977). The upper panel shows the e GLYPH<14> ciency of mass transfer, the middle panel the angular frequency of each star in terms of its critical value, while the lower panel shows the evolution of ˙ M for both components.  \n<!-- image -->  \n/circledot', '2.10. Description of a Binary Run': 'MESAbinary performs each evolution step by independently solving the structure of each component and the orbital parameters, using the same timestep GLYPH<14> t for each. This approach di GLYPH<11> ers from STARS , which simultaneously solves for the structure of both stars and the orbit in a single NewtonRaphson solver. Our choice to solve for each star separately gives a significant amount of flexibility and simplicity, as the examples in this paper demonstrate.  \nThe top-level algorithm for evolving a star is described in appendix B1 of Paper II. We modified this algorithm to support the new implementation of binary interactions, which is described in detail in the MESA documentation. Additional timestep limits are imposed in MESAbinary that consider relative changes between the radius and Roche lobe radius of both components, the total orbital angular momentum, the orbital separation, and the envelope mass in the donor.', '3. PULSATIONS': 'The study of stellar pulsations (also termed oscillations) o GLYPH<11> ers unique insights into the interiors of stars (Aerts et al. 2010). In some classes of star (e.g., solar-type, red giant), the stochastic excitation of hundreds of oscillation modes, typically by convective motions, allows remarkably detailed measurements to be made of the interior, including nuclear burning state (Bedding et al. 2011) and internal rotation (Beck  \net al. 2012). In other classes (e.g., classical Cepheid, GLYPH<12> Cephei, GLYPH<14> Scuti and GLYPH<13> Doradus pulsators), modes are instead excited by linear instabilities, most often linked to opacity variations in the envelope (the GLYPH<20> mechanism). In these latter objects, typically too few modes are excited for detailed asteroseismic analysis to be feasible; nevertheless, mapping out the regions of the theoretical HR diagram where the instabilities are expected to operate, and then comparing these instability strips against observational surveys, can often lead to new science.  \nGYRE calculates the normal-mode eigenfrequencies GLYPH<27> of a stellar model by solving the system of linearized equations and boundary conditions governing small periodic perturbations ( / exp[i GLYPH<27> t ]) to the equilibrium state. It is based on a novel Magnus Multiple Shooting (MMS) numerical scheme which is robust and accurate, and makes full use of all available processors on multicore computer architectures. The MMS scheme and the initial release of the code, which focuses on adiabatic pulsations, is described in Townsend & Teitler (2013); extensions to the code to support non-adiabatic pulsations are described in Goldstein & Townsend (2015).  \nPaper II introduced the astero extension to MESAstar , which permits on-the-fly refinement of stellar model parameters in order to fit a set of observed oscillation frequencies and spectroscopic constraints. Subsequent improvements to the astero capabilities include frequency correction recipes from Ball & Gizon (2014); implementation of the downhill simplex (Nelder & Mead 1965) and NEWYUO (Powell 2004) algorithms for GLYPH<31> 2 minimization; parameter optimization using only spectroscopic constraints (e.g., T e GLYPH<11> and surface gravity); and coupling to the GYRE oscillation code, as an alternative to the ADIPLS code (Christensen-Dalsgaard 2008) used in the original implementation.  \nMESAstar couples to GYRE via two mechanisms. Loose coupling is achieved simply by MESAstar writing models out to disk, and GYRE subsequently reading these models in; we use this process below to map out massive-star instability strips. Tight coupling removes the intermediate disk usage, by handling all communication between MESAstar and GYRE in-memory; this permits fully closed-loop calculations, where the changes in the oscillation eigenfrequencies of an evolving stellar model are used to guide the further evolution of the model. Tight coupling allows GYRE to function as an alternative to ADIPLS in the astero extension, and opens up the possibility of other kinds of novel calculations, such as the automated location of instability-strip boundaries.', '3.1. Massive-Star Instability Strips': 'As an illustration of a large-scale calculation using MESAstar and GYRE loosely coupled, Figure 9 plots the instability strips for massive stars on and near the upper main sequence, for oscillation modes with harmonic degrees \' = 0-3. These strips are based on a set of 182 evolutionary tracks, each extending from the ZAMS across to a red limit at log( T e GLYPH<11> = K) = 3 : 75, with 101 tracks spanning the initial mass range 2 : 5 M GLYPH<12> GLYPH<20> M GLYPH<20> 25 M GLYPH<12> in uniform logarithmic increments, and the remaining 81 tracks spanning the mass range 6 M GLYPH<12> GLYPH<20> M GLYPH<20> 10 M GLYPH<12> in uniform linear increments (the latter set is designed to adequately resolve the \'fingers\' discussed below). OPAL opacity tables are used with the proto-solar abundances from Asplund et al. (2009), and for simplicity we neglect any rotation or mass loss. Convection is modeled with a mixing-length parameter GLYPH<11> MLT = 1.5 and an exponential overshoot parameter f ov = 0 : 024, and the Schwarzschild  \nFigure 9. Instability strips for \' = 0-3 oscillation modes in the upper part of the Hertzsprung-Russell diagram. Separate strips are shown for the GLYPH<12> Cephei ( ! > 1) and slowly pulsating B-type (SPB; ! < 1) classes of pulsating stars. The ZAMS and red edge of the main sequence (REMS) are shown for reference, as are evolutionary tracks for models with selected masses (labeled in solar units along the ZAMS). The red edges of the post-MS SPB strips are drawn with a dotted line, indicating that the positioning of these edges is an artifact of our numerical procedure.  \n<!-- image -->  \nstability criterion is assumed.  \nWeselect points i = i 1 ; i 2 ; i 3 ; : : : along each of the 182 tracks (where i is the timestep index; see section 6.4 of Paper I), chosen so that i 1 corresponds to the ZAMS,  \ni 2 GLYPH<0> 1 X i = i 1 8 > > < > > : " log( T e GLYPH<11> ; i + 1 = T e GLYPH<11> ; i ) GLYPH<1> T # 2 + " log( Li + 1 = Li ) GLYPH<1> L # 2 9 > > = > > ; 1 = 2 GLYPH<25> 1 (22)  \nacross the ( i 1 ; i 2) pair, and similarly for subsequent pairs. Here, GLYPH<1> T and GLYPH<1> L are dimensionless weights which control the spacing of points in e GLYPH<11> ective temperature and luminos-  \nity; we adopt the values 0.004 and 0.011, respectively, for these weights. At the selected points, GYRE searches for unstable oscillation modes with the harmonic degrees considered. First, GYRE solves the adiabatic oscillation equations to find eigenfrequencies GLYPH<27> ad falling in the range extending from the asymptotic frequency of the gravity (g) mode with radial order n = 400, up to the asymptotic frequency of the pressure (p) mode with radial order n = 10. Each GLYPH<27> ad is then used as an initial guess in finding a corresponding eigenfrequency GLYPH<27> of the full non-adiabatic oscillation equations. The real and imaginary parts of GLYPH<27> give the linear frequency GLYPH<23> osc and the  \nFigure 10. The \' = 1 dimensionless frequency spectrum of an 8.5 M GLYPH<12> stellar model as it evolves from the ZAMS to the REMS. Blue (orange) dots indicate which modes are stable (unstable); selected modes are labeled along the left / bottom edge using their classification.  \n<!-- image -->  \nFigure 11. Instability strips for dipole ( \' = 1) oscillation modes in the upper part of the HR diagram, but calculated using OP rather than OPAL opacities (cf. Figure 9).  \n<!-- image -->  \ngrowth e -folding time GLYPH<28> osc of a mode:  \nGLYPH<23> osc = Re( GLYPH<27> ) 2 GLYPH<25> ; GLYPH<28> osc = GLYPH<0> 1 Im( GLYPH<27> ) : (23)  \nIf GLYPH<28> osc is negative, the mode is damped. Separate strips are shown in Figure 9 for regions exhibiting  \nunstable modes with Re( ! ) > 1 and Re( ! ) < 1, where  \n! = GLYPH<27> r R 3 GM (24)  \nis the dimensionless eigenfrequency; these correspond, respectively, to the GLYPH<12> Cephei and slowly pulsating B-type (SPB) classes of pulsating stars. In GLYPH<12> Cephei stars during the MS phase, p and g modes with periods of a few hours and radial orders n GLYPH<25> 1-3 are excited by a GLYPH<20> mechanism operating on the iron opacity bump situated in the outer envelope at log( T = K) GLYPH<25> 5 : 3 (Cox et al. 1992; Dziembowski & Pamiatnykh 1993). In SPB stars during the MS phase, g modes with periods of a few days and radial orders n GLYPH<25> 20-50 are excited by the same mechanism (Dziembowski et al. 1993). For masses M & 9 M GLYPH<12> the strips for both classes of stars extend into the post-MS domain. During this phase, unstable modes couple with g modes trapped near the boundary of the inert helium core. In the case of the SPB stars this leads to very high overall radial orders, n & 100, and ultimately limits our ability to follow the instability strips all the way to the red edge (our calculations are restricted to n . 400 for computational e GLYPH<14> ciency reasons). Hence, in Figure 9 we plot the red edges of the post-MS SPB strips with dotted lines, to highlight that these are not the true red edges.  \nThis alternation between instability and stability, seen as fingers in Figure 9, stems from the fact that the GLYPH<20> mechanism only excites modes whose eigenfrequencies fall in a narrow range [ GLYPH<27> lo ; GLYPH<27> hi]. At frequencies Re( GLYPH<27> ) > GLYPH<27> hi, the pulsation period becomes comparable to the local thermal timescale in the envelope region above the iron opacity peak, and this region behaves as a damping zone, stabilizing the modes. Conversely, at frequencies Re( GLYPH<27> ) < GLYPH<27> lo, modes couple with gravity waves trapped in the GLYPH<22> -gradient zone developing at the core boundary, and are likewise damped. The intermediate stable phase in Figure 10, between log( T e GLYPH<11> = K) = 4 : 301 and log( T e GLYPH<11> = K) = 4 : 288 occurs when there are no modes in the [ GLYPH<27> lo ; GLYPH<27> hi] range. As the star evolves, the unstable range narrows: GLYPH<27> hi decreases due to lower T e GLYPH<11> , while GLYPH<27> lo increases due to the growth of the GLYPH<22> -gradient zone.  \nAllowing for di GLYPH<11> erences in adopted abundances and other modeling parameters, the instability strips plotted in Figure 9 are in general agreement with those published in the literature (e.g., Pamyatnykh 1999; Zdravkov & Pamyatnykh 2008; Saio 2011). The notable di GLYPH<11> erence is the presence of fingers in the lower boundaries of our GLYPH<12> Cephei strips for \' GLYPH<21> 1. Their appearance here is due to the unprecedented resolution in HR-diagram space of our stability calculations. To elucidate their origin, Figure 10 plots part of the \' = 1 frequency spectrum of an 8 : 5 M GLYPH<12> stellar model as it evolves from the ZAMS to the red edge of the main sequence (REMS), showing which modes are stable and which are unstable. The p1 mode is unstable over the e GLYPH<11> ective temperature range 4 : 358 GLYPH<21> log( T e GLYPH<11> = K) GLYPH<21> 4 : 317, and the g1 mode over the cooler but overlapping range 4 : 341 GLYPH<21> log( T e GLYPH<11> = K) GLYPH<21> 4 : 301. The star then passes through a phase with no unstable modes, before the instability reappears in the range 4 : 288 GLYPH<21> log( T e GLYPH<11> = K) GLYPH<21> 4 : 278 for the g2 mode.  \nFigure 11 shows a version of the \' = 1 panel calculated using OP opacity tables rather than OPAL tables. There is an overall shift of the instability strips toward higher luminosities, an e GLYPH<11> ect already noted by Pamyatnykh (1999). The fingers persist with much the same structure, supporting the fact that they are physical e GLYPH<11> ects rather than numerical artifacts.  \nReturning now to Figure 9, the post-MS extension of the SPB strips has been attributed in the literature to features in the Brunt-Vaisala frequency which reflect gravity waves at the boundary of the helium core, preventing them from penetrating into the core and being dissipated by strong radiative damping. Saio et al. (2006) and Godart et al. (2009) argue that the necessary feature is an intermediate convection zone (ICZ) associated with the hydrogen-burning shell, but more recently Daszy\'nska-Daszkiewicz et al. (2013) have shown that even a local minimum in the Brunt-Vaisala frequency is su GLYPH<14> cient to reflect modes. In the present case, the empirical mass threshold M & 9 M GLYPH<12> required for formation of an ICZ coincides with the lower boundaries of the SPB strip extensions. In the lowest-mass models above this threshold, the ICZ vanishes shortly after its appearance, but it leaves behind a narrow region with a steep molecular weight gradient. This gradient causes a spike in the Brunt-Vaisala frequency, which serves in a similar manner to prevent gravity waves from entering into the core and being dissipated.  \nThe corresponding post-MS extension of the GLYPH<12> Cephei strips was first noted by Dziembowski & Pamiatnykh (1993), but has not received much attention in the literature. Figure 9 shows that this extension has a well defined lower boundary, much like the SPB stars although situated at slightly higher masses, M & 10 : 5 M GLYPH<12> . We have determined that the extension is also a consequence of ICZ formation; the shift to higher masses arises because it appears that multiple convection zones, rather than a single one, are necessary to reflect waves at the core boundary in the case of GLYPH<12> Cephei pulsators.', '3.2. Asteroseismic Optimization': "Figure 12. Echelle diagram for the subgiant star HD 49385. Observed frequencies are shown as filled circles ( ' = 0), triangles ( ' = 1) and squares ( ' = 2); black horizontal lines indicate the 1 GLYPH<27> error bars. Calculated frequencies of the best-fit model are overplotted as the corresponding open symbols.  \n<!-- image -->  \nTo illustrate the updated asteroseismic capabilities of MESA , Figure 12 plots the echelle diagram for the subgiant star HD49385, showing both the frequencies of ' = 0 GLYPH<0> 2 modes  \nmeasured by Deheuvels et al. (2010), and the corresponding frequencies of the best-fit model determined using the astero extension. The calculations follow the same procedure detailed in section 3.2 of Paper II; the only significant di GLYPH<11> erences are that the initial mass, helium abundance, metal abundance and mixing length parameter are refined using the downhill simplex algorithm rather than the Hooke-Jeeves algorithm; oscillation frequencies are calculated using GYRE rather than ADIPLS ; and the surface corrections to frequencies are evaluated using equation 4 of Ball & Gizon (2014) rather than with the Kjeldsen et al. (2008) scheme.  \nComparing Figure 12 against figure 8 of Paper II reveals only small di GLYPH<11> erences between the two. The GLYPH<31> 2 of the best-fit models reported by astero is 2.3 in the former case, compared to 2.4 in the latter (cf. table 2 of Paper II).", '3.3. Automated Strip Location': "Figure 13. The growth timescale GLYPH<28> osc (left axis) and oscillation period P osc (right axis) of the fundamental and first-overtone radial modes of a 8.5 M GLYPH<12> model, plotted as a function of T e GLYPH<11> as the star evolves away from the ZAMS. The vertical dashed lines, determined automatically, show the points where the modes switch from stable ( GLYPH<28> osc < 0) to unstable ( GLYPH<28> osc > 0), and vice versa.  \n<!-- image -->  \nThe instability strips presented above involved the examination of GLYPH<24> 11 million modes of GLYPH<24> 40,000 stellar models. To partially automate the process we can leverage tight coupling between GYRE and MESAstar . This is achieved by making small modifications to the extras check model callback routine in MESAstar (see Appendix B.1 of Paper II), so that  \nFigure 14. The calculated blue edge of the classical instability strip, for fundamental and first-overtone radial modes. The corresponding dashed lines show the predictions from set B of Smolec & Moskalik (2008, their figure 1).  \n<!-- image -->  \nGYRE is run after each time step to determine the set of eigenfrequencies f GLYPH<27> g of a user-specified group of modes. When Im( GLYPH<27> ) changes sign from one time step to the next for any of these modes, indicating that an instability-strip boundary has been crossed, a search is performed to find Im( GLYPH<27> ) GLYPH<25> 0.  \nAs a further demonstration of automated instability strip location, Figure 14 plots the blue edges of the classical instability strip in the HR diagram, for fundamental and first-overtone radial modes. The edges are calculated for 51 evolutionary tracks spanning the initial mass range 1 : 25 M GLYPH<12> GLYPH<20> M GLYPH<20> 12 : 5 M GLYPH<12> in uniform logarithmic increments. At luminosities log( L = L GLYPH<12> ) & 2 : 5 corresponding to classical Cepheid pulsators, these edges show good agreement with the set B results published by Smolec & Moskalik (2008, their figure 1). At lu-  \nFigure 13 presents an application of the tight coupling to the fundamental and first-overtone radial modes of the 8.5 M GLYPH<12> model considered in Section 3.1, showing how the growth timescales GLYPH<28> osc and oscillation periods P osc ˜ = ˜ 1 / GLYPH<23> osc of the modes change as the star evolves from the ZAMS into the post-MS. The second-overtone radial mode remains stable, GLYPH<28> osc < 0, over the range plotted. The vertical lines show where GLYPH<28> osc changes sign. The blue and red edges of the ( ' = 0) GLYPH<12> Cephei instability strip can be seen in the upper panel of Figure 13 at log( T e GLYPH<11> = K) = 4 : 36 and log( T e GLYPH<11> = K) = 4 : 30, respectively. The blue edge 12 of the classical instability strip can likewise be seen in both panels at log( T e GLYPH<11> = K) = 3 : 75. The corresponding red edge is not found because GYRE does not include a treatment of the pulsation-convection interaction -anecessary ingredient for modeling the classical red edge (see, e.g, section 3.7.3 of Aerts et al. 2010, and references therein).  \n12 For classical ( GLYPH<14> ) Cepheid pulsators, the observational blue edge of the classical instability strip is in fact established by stars evolving to higher T e GLYPH<11> on their first blue loop, rather than stars on their first crossing of the Hertzsprung gap. However, the purpose of the present section is to demonstrate the capability of tightly coupling, and in this context the distinction between the blue edges from multiple crossings is unimportant.  \nminosities log( L = L GLYPH<12> ) . 1 : 6 corresponding to GLYPH<14> Scuti stars, the edges are somewhat cooler than results published in the literature; however, this is because we consider only fundamental and first-overtone modes, whereas the blue edge is typically set by higher overtones which are displaced toward hotter T e GLYPH<11> (e.g., Dupret et al. 2004, their figure 1).  \nIdeally, the same automated approach could be used to locate the boundaries of the non-radial ( ' > 0) instability strips plotted in Figure 9. In practice it is very challenging to devise a robust algorithm that can unambiguously interpret the eigenfrequencies produced by GYRE . Sometimes, acoustic glitches in a model can trap modes in surface layers, where they are strongly excited; however, these modes are very sensitive to model parameters, and it is unclear whether they are physically meaningful or not.", '4. IMPLICIT HYDRODYNAMICS': 'Shocks happen in stars, such as after a massive star collapses, or cyclically in the outer envelopes of stars pulsating at su GLYPH<14> ciently large amplitude. Previous versions of MESAstar allowed large velocities such as those encountered in the last few seconds leading to a core collapse ( GLYPH<25> 1000 km s GLYPH<0> 1 ), but there was no provision for large jumps in velocities leading to shocks. In this section we describe the changes that have been made to support an implicit treatment of hydrodynamic shocks that includes careful attention to conservation of energy. We demonstrate that the revised equations are intrinsically conservative in the sense that deviations from exact energy balance can only arise from residual numerical errors in the approximate solutions rather than from the form of the equations themselves. Following the description of the changes, we show a series of envelope shocks as a test of the implementation. The form of the equations and the demonstration of intrinsic conservation closely follow Fraley (1968) and Grott et al. (2005). The treatment of artificial viscosity is based on Weaver et al. (1978).  \nFigure 15. Schematic of relevant cell and face variables relevant for hydrodynamics in MESAstar .  \n<!-- image -->', '4.1. Mass Continuity': 'The specific volume of cell k is  \n1 GLYPH<26> k = 4 GLYPH<25> 3 0 B B B B @ r 3 k GLYPH<0> r 3 k + 1 dmk 1 C C C C A ; (25)  \nwhere rk is the outer face radius, rk + 1 is the inner face radius, dmk is the cell mass, and GLYPH<26> k is the cell average density (see Figure 15 for the layout of cells in MESAstar ). We create an initial algebraic form of the continuity di GLYPH<11> erential equation by dividing the change in the specific volume in a step by the length of time GLYPH<14> t , using step start and end values for rk , rk + 1, and GLYPH<26> k and the value dmk which is constant during the step:  \n1 =GLYPH<26> k GLYPH<0> 1 =GLYPH<26> start ; k GLYPH<14> t = 4 GLYPH<25> 3 dmk hGLYPH<16> r 3 k GLYPH<0> r 3 start ; k GLYPH<17> GLYPH<0> GLYPH<16> r 3 k + 1 GLYPH<0> r 3 start ; k + 1 GLYPH<17>i GLYPH<14> t : (26)  \nNext, we rewrite the right hand side introducing variables for the time centered velocity ˆ v and the e GLYPH<11> ective area A to get the final form of the mass continuity equation as used in MESAstar :  \n1 =GLYPH<26> k GLYPH<0> 1 =GLYPH<26> start ; k GLYPH<14> t = 1 dmk ( Ak ˆ vk GLYPH<0> Ak + 1ˆ vk + 1) ; (27)  \nˆ vk = ( vk + v start ; k ) = 2 (28)  \nwhere  \nand rk is evaluated as  \nrk = r start ; k + ˆ vk GLYPH<14> t : (29)  \nAlgebraic simplification then shows that  \nAk = 4 GLYPH<25> 3 GLYPH<16> r 2 k + rkr start ; k + r 2 start ; k GLYPH<17> : (30)  \nTo be consistent with the mass continuity equation, we use these expressions for e GLYPH<11> ective area and time centered velocity in the following momentum and energy equations. It will be shown below that to get intrinsic energy conservation, we must time center the velocity and use special combinations of starting and ending radius in a couple of places, but all other terms in the equations can remain fully implicit to avoid degrading the numerical stability as would happen in a uniformly time-centered scheme.', '4.2. Artificial Viscosity': 'In MESAstar , the artificial dynamic viscosity coe GLYPH<14> cient GLYPH<17> visc (which has the dimensions g cm GLYPH<0> 1 s GLYPH<0> 1 )  \nGLYPH<17> visc ; k = GLYPH<17> visc ; linear ; k + GLYPH<17> visc ; quad ; k ; (31)  \nwhere the linear term is  \nGLYPH<17> visc ; linear ; k = 3 4 l 1 GLYPH<26> kr mid ; kc s ; k (32)  \nand the quadratic term is  \nGLYPH<17> visc ; quad ; k = 3 4 l 2 2 GLYPH<26> 2 k 4 GLYPH<25> r 2 mid ; k dmk max GLYPH<16> 0 ; r 2 k + 1 ˆ vk + 1 GLYPH<0> r 2 k ˆ vk GLYPH<17> ; (33)  \nwith r mid ; k = ( rk + 1 + rk ) = 2, c s ; k the sound speed in cell k , and l 1 ( l 2) is a dimensionless coe GLYPH<14> cient for the linear (quadratic) term. The linear term is rarely used; it provides for a general damping of pulsations. The quadratic term is only nonzero in regions of compression and is the primary control for the strength of artificial viscosity. Assuming the usual case of l 1 = 0, the shock front is spread over a distance GLYPH<24> l 2 rk . We follow Dorfi (1998) in opting for a shock spread proportional to the local radius r rather than the local cell width. This choice is dictated by the fact that step-by-step adjustments to the mesh resolution lead to dynamically changing, nonmonotonic variations in cell widths of up to a factor of 2 or  \nmore between neighboring cells. Making the shock spread directly dependent on the local cell widths would produce numerically intolerable dynamically changing, non-monotonic variations in cell-to-cell values for the shock spread. Use of a local running average cell width is also ruled out by the need to keep algebraic equations dependent only on nearest neighbors to allow a block tridiagonal matrix solution. The use of a small fraction of the local radius gives a smoothly varying shock spread that avoids the numerical problems associated with using the cell width.  \nWe define the quantity Q visc ; k (having dimensions of energy), in cell k as  \nQ visc ; k = GLYPH<17> visc ; k GLYPH<26> k 4 GLYPH<25> r 6 mid ; k dmk ˆ vk rk GLYPH<0> ˆ vk + 1 rk + 1 ! : (34)  \nThe momentum equation uses Q visc in an expression that defines an artificial acceleration analogous to the pressure gradient term, and the energy equation uses it to define an artificial viscous heating analogous to the mechanical work term.', '4.3. Specific Linear Momentum Equation': 'The local linear momentum conservation equation at face k between inner cell k and outer cell k GLYPH<0> 1 is  \nvk GLYPH<0> v start ; k GLYPH<14> t = GLYPH<0> Gmk rkr start ; k GLYPH<0> Ak Pk GLYPH<0> 1 GLYPH<0> Pk dmk ! + g visc ; k : (35)  \nwhere dmk = ( dmk + dmk GLYPH<0> 1) = 2 is the mass associated with face k , and the viscous acceleration term at face k is  \ng visc ; k = 4 GLYPH<25> rk Q visc ; k GLYPH<0> 1 GLYPH<0> Q visc ; k dmk ! : (36)  \nThe use of the product rkr start ; k in the denominator of the gravitation term is necessary for intrinsic energy conservation as will be shown below.', '4.4. Specific Energy Equation': 'The local energy conservation equation for cell k between outer face k and inner face k + 1 is  \nek GLYPH<0> e start ; k GLYPH<14> t = GLYPH<0> Lk GLYPH<0> Lk + 1 dmk GLYPH<0> Pk Ak ˆ vk GLYPH<0> Ak + 1 ˆ vk + 1 dmk ! (37)  \n+ GLYPH<15> visc ; k + GLYPH<15> nuc ; k GLYPH<0> GLYPH<15>GLYPH<23>; k + GLYPH<15> extra ; k ;  \nwhere ek is the specific thermal energy for cell k . The viscous heating rate for cell k is  \nGLYPH<15> visc ; k = 4 GLYPH<25> Q visc ; k dmk ˆ vk rk GLYPH<0> ˆ vk + 1 rk + 1 ! : (38)  \nEnergy loss from weak reaction neutrinos is already subtracted from the nuclear burning term, GLYPH<15> nuc ; k , so only the neutrino energy loss rate from thermal processes, GLYPH<15>GLYPH<23>; k , is explicitly accounted for in Equation (37). An example of GLYPH<15> extra ; k is artificial injection of energy to trigger a shock.  \nAn alternative form of the energy equation equates the model dL = dm to the expected value  \nLk GLYPH<0> Lk + 1 dmk = dL dm ! expected ; k ; (39)  \nwhere  \ndL dm ! expected ; k = GLYPH<15> grav ; k + GLYPH<15> visc ; k + GLYPH<15> nuc ; k GLYPH<0> GLYPH<15>GLYPH<23>; k + GLYPH<15> extra ; k ; (40)  \nand  \nGLYPH<15> grav ; k = GLYPH<0> ek GLYPH<0> e start ; k GLYPH<14> t GLYPH<0> Pk Ak ˆ vk GLYPH<0> Ak + 1 ˆ vk + 1 dmk ! : (41)  \nUsing Equation (27), the expression for GLYPH<15> grav can be rewritten  \nGLYPH<15> grav ; k = GLYPH<0> ek GLYPH<0> e start ; k GLYPH<14> t GLYPH<0> Pk (1 =GLYPH<26> k GLYPH<0> 1 =GLYPH<26> start ; k ) GLYPH<14> t (42)  \nthereby avoiding the use of velocities and thus be appropriate for hydrostatic cases.', '4.5. Intrinsic Energy Conservation': 'The summed kinetic, potential, internal energies are  \nKE = X k 1 2 dmk v 2 k ; (43)  \nPE = X k GLYPH<0> Gmkdmk rk ; (44)  \nIE = X k ek dmk ; (45)  \nand thus the total energy of the star is E = KE + PE + IE. We now explicitly demonstrate that the equations we solve are formulated in such a way that the rate of change of total energy exactly equals the combined energy sources and sinks.  \nMultiplying Equation (35) by ˆ vkdmk gives an equation with units of luminosity:  \n1 2 dmk 0 B B B B B @ v 2 k GLYPH<0> v 2 start ; k GLYPH<14> t 1 C C C C C A = GLYPH<0> Gmkdmk ˆ vk rkr start ; k ! GLYPH<0> Ak ˆ vk ( Pk GLYPH<0> 1 GLYPH<0> Pk ) + 4 GLYPH<25> ˆ vk rk ! ( Q visc ; k GLYPH<0> 1 GLYPH<0> Q visc ; k ) : (46)  \nUsing Equation (29) to eliminate ˆ vk in the first term on the right,  \nGLYPH<0> Gmkdmk ˆ vk rkr start ; k ! = Gmkdmk GLYPH<14> t 1 rk GLYPH<0> 1 r start ; k ! ; (47)  \nshows that this term is the negative of the rate of change of potential energy, a result that is made possible by the use of the Gmk = rkr start ; k in Equation (35) instead of an alternative such as Gmk = r 2 k . Thus, Equation (46) can be written as  \n1 2 dmk v 2 k GLYPH<0> v 2 start ; k GLYPH<14> t GLYPH<0> Gmkdmk GLYPH<14> t 1 rk GLYPH<0> 1 r start ; k ! = GLYPH<0> Ak ˆ vk ( Pk GLYPH<0> 1 GLYPH<0> Pk ) + 4 GLYPH<25> ˆ vk rk ! ( Q visc ; k GLYPH<0> 1 GLYPH<0> Q visc ; k ) : (48)  \nSimilarly, multiplying equation (37) by dmk also yields an equation with units of luminosity:  \n( ek GLYPH<0> e start ; k ) GLYPH<14> t dmk = GLYPH<0> ( Lk GLYPH<0> Lk + 1) GLYPH<0> Pk ( Ak ˆ vk GLYPH<0> Ak + 1 ˆ vk + 1) + 4 GLYPH<25> Q visc ; k ˆ vk rk GLYPH<0> ˆ vk + 1 rk + 1 ! + ( GLYPH<15> nuc ; k GLYPH<0> GLYPH<15>GLYPH<23>; k + GLYPH<15> extra ; k ) dmk ; (49)  \nAdding Equations (48) and (49) and summing over the grid index k gives  \nX k GLYPH<0> Gmkdmk GLYPH<14> t 1 rk GLYPH<0> 1 r start ; k ! + 1 2 dmk v 2 k GLYPH<0> v 2 start ; k GLYPH<14> t + ( ek GLYPH<0> e start ; k ) GLYPH<14> t dmk X GLYPH<0> ( Lk GLYPH<0> Lk + 1) (50)  \n= k + PkAk + 1ˆ vk + 1 GLYPH<0> Pk GLYPH<0> 1 Ak ˆ vk + 4 GLYPH<25> Q visc ; k GLYPH<0> 1 ˆ vk rk ! GLYPH<0> 4 GLYPH<25> Q visc ; k ˆ vk + 1 rk + 1 ! + ( GLYPH<15> nuc ; k GLYPH<0> GLYPH<15>GLYPH<23>; k + GLYPH<15> extra ; k ) dmk :  \nThe sum over the pressure terms is  \nX k ( PkAk + 1 ˆ vk + 1 GLYPH<0> Pk GLYPH<0> 1 Ak ˆ vk ) = GLYPH<0> GLYPH<2> ( PA ˆ v ) surface GLYPH<0> ( PA ˆ v ) center GLYPH<3> = GLYPH<0> ( L acoustic ; surface GLYPH<0> L acoustic ; center) ; (51)  \nwhich cancels term by term except at the boundaries. We define L acoustic ; surface as the work done by the model on the atmosphere at the surface and L acoustic ; center as the work done on the model at the center, for example, by an artificial piston. The sum over the artificial viscosity terms leads to a similar expression, but because Q visc vanishes at the surface and the center, the sum equals zero. That is, the energy added by artificial viscous heating in the energy equation exactly balances the loss of kinetic energy by artificial viscous acceleration in the momentum equation.  \nThe terms on the left hand side of Equation (50) are the difference in the total energy between the start and end of a step divided by the length of the step, in other words, the average rate of change of the total energy of the model. Therefore Equation (50) can be written as  \n( E final GLYPH<0> E initial) =GLYPH<14> t = GLYPH<0> ( L surface GLYPH<0> L center) GLYPH<0> ( L acoustic ; surface GLYPH<0> L acoustic ; center) + X k GLYPH<0> GLYPH<15> nuc ; k GLYPH<0> GLYPH<15>GLYPH<23>; k + GLYPH<15> extra ; k GLYPH<1> dmk : (52)  \nThis equation embodies conservation of energy in MESAstar : the rate of change of total energy equals the combined energy sources and sinks. This demonstrates that in the given form, the algebraic equations intrinsically conserve energy in the sense that failure to get energy balance can only arise from the residual numerical errors that are inherent in using approximate solutions to the equations. This in turn means that to control energy balance errors, we can focus on reducing residuals either by changes in the Newton solver or by timestep reductions.  \n4.6. Controlling the Accuracy of Energy Conservation  \nThe Newton solver considers both the sizes of incremental changes to the variables and the sizes of residual errors for the equations. For the energy equation, the residual used by the solver is defined to be the timestep GLYPH<14> t times the di GLYPH<11> erence between the left and right sides of Equation (37) divided by e start ; k ; in other words, the residual is the error as a fraction of the specific energy at the start of the step. By adjusting tolerances for the average and maximum size of residuals, we force the Newton solver to take extra iterations to reduce the residuals which will in turn reduce the total error in energy conservation.  \nA second and related way to control energy conservation errors is to use the average and maximum energy residuals to adjust the timesteps. For example, if the maximum magnitude for an energy residual exceeds a specified hard limit, then the proposed solution is rejected and the step is retried with a smaller timestep. If the maximum is smaller than the hard limit but exceeds another specified limit, then there is no forced retry for this step, but the next timestep is reduced by the ratio of the limit divided by the maximum magnitude. If the maximum is smaller than both limits, then other factors determine the next timestep. Later in this section, we will show that these approaches in combination with the intrinsic conservation form of the equations yield a solution for a shock in an envelope that evolves with reasonably large timesteps while conserving energy to a high degree of accuracy.', '4.7. Limiting Acceleration of Convective Velocity': 'When using hydrodynamics, we often require timesteps that are so small that we need to limit the increase in convective velocities as calculated in the standard instantaneous mixing length theory (MLT) so that they do not assume unphysically large accelerations. If convection velocities are allowed to adjust instantaneously, then our methods for artificially creating shocks will fail since however rapidly we inject energy over a limited region, convection will be able to transport the energy away. To be able to simulate shocks we need to have a way to limit convection velocity acceleration.  \nThe primary scheme we use for this is derived from Arnett (1969) and Wood (1974). The MLT implementation in MESA has been extended to take as additional arguments the timestep and the previous convection velocity at the same mass location ( vc ; prev). It calculates a provisional convection velocity ( vc 0) using the standard instantaneous MLT, then defines a convective timescale ( GLYPH<28> MLT) as the local pressure scale height ( H ) divided by the sum of the provisional plus previous velocities. If GLYPH<14> t is less than GLYPH<28> MLT, then the next convection velocity ( vc ) is only incremented from the previous one by the di GLYPH<11> erence of the provisional minus the previous velocities times the ratio of the timestep divided by the time scale  \nvc = vc ; prev + min 1 ; GLYPH<14> t GLYPH<28> MLT ! ( vc 0 GLYPH<0> vc ; prev) ; (53)  \nGLYPH<28> MLT = H ( vc 0 + vc ; prev) : (54)  \nwhere  \nAs an alternative scheme for limiting convection acceleration, we also allow the maximum rate of change of convection velocity to be set as a fraction, g GLYPH<18> , of the local gravitational acceleration. If vc 0 > vc ; prev, then  \nvc = min( vc 0 ; vc ; prev + GLYPH<14> t g GLYPH<18> g ) ; (55)  \nThe final vc is used to recalculate the convection e GLYPH<14> ciency, which is used to calculate the MLT temperature gradient.  \nThese methods for limiting the acceleration of convective velocities reduce the energy transport rate as well as the rate of compositional mixing. Both schemes seem to give at least qualitatively reasonable results and avoid the problems of unphysically large accelerations that are possible with standard instantaneous MLT. Hopefully this ad hoc solution will soon be replaced by a quantitatively correct formulation.', '4.8. Surface Boundary Conditions': 'MESA provides a variety of options for surface boundary conditions (see, e.g., section 5.3 of Paper I), and several more have been added for use with hydrodynamics. The simplest allow specification of a particular value for the surface pressure, the surface temperature, or the T e GLYPH<11> if the surface is not at the photosphere. In the case of a given fixed surface pressure, the corresponding surface temperature is set using the surface luminosity and radius based on the usual black body relation. For the second case, where the surface temperature is fixed, the surface pressure is set to the corresponding radiation pressure. For both of these, if the surface is not at the photosphere, T e GLYPH<11> is set using the Eddington T -GLYPH<28> relation. Finally, for specified T e GLYPH<11> when the surface is not at the photosphere, the corresponding surface temperature is also derived using the Eddington T -GLYPH<28> relation, and the surface pressure is set to the radiation pressure for that temperature.  \nFor computations involving shocks at the surface, there is an option to use boundary conditions that specify a vanishing gradient for compression at the surface and a temperature corresponding to black body radiation. The outermost cells ( k = 1 ; 2) satisfy the equation  \n1 GLYPH<26> 1 GLYPH<0> 1 GLYPH<26> start ; 1 = 1 GLYPH<26> 2 GLYPH<0> 1 GLYPH<26> start ; 2 ; (56)  \nwhich represents the vanishing of the surface compression gradient.  \nFinally, for computations involving interior shocks but low velocities at the surface, there is an option to use the surface pressure from the selected atmosphere prescription with the momentum equation relating the surface velocity to the surface pressure gradient. This form for the surface boundary conditions is used in the shocked massive star example in Section 6 and in the following envelope shock test.', '4.9. Shock Test': 'To test the implementation, we shock the extended envelope of a 6.93 M GLYPH<12> asymptotic giant branch (AGB) star evolved from a 7 M GLYPH<12> main-sequence star without rotation and an initial metallicity of 0.001. This case is chosen because of the uniform properties of the extended envelope (i.e., small density range, smooth density, and uniform composition). Our interest is to study the propagation of the shock, the properties of the shocked material, and the magnitude of energy conservation errors. In Section 6, we present results that mimic core-collapse supernovae.  \nExplosion simulations with MESA start from a converged model. The core is excised by removing inner shells of the model and setting new inner boundary conditions for mass, radius, velocity, and luminosity. For the current test, we remove the center just above the helium core at a mass of 2.40 M GLYPH<12> which corresponds to an inner radius of 27.2 R GLYPH<12> . The stellar surface lies at a radius of 282.7 R GLYPH<12> . During the following  \n3  \n3  \nFigure 16. Multi-epoch snapshots for the hydrodynamical simulation of a 10 49 erg shock in the envelope of a 6 : 93 M GLYPH<12> AGB star. We show the density (top row), temperature (middle row), and velocity (bottom row), versus Lagrangian mass (left column) and radius (right column). In each panel, the solid line refers to the MESA results and the dashed line to the V1D results.  \n<!-- image -->  \nevolution, the excised region is treated as a point mass and is linked to the above layers by the inner boundary conditions which can be changed at each timestep to simulate various core events. The model grid was adjusted at each step to give higher resolution in the vicinity of the shock. The total number of cells stayed at about 1000, with cell masses dropping to about 10 GLYPH<0> 4 of the total in approximately 100 cells around the shock.', 'thermal bomb:': "- GLYPH<15> The first option changes the inner boundary conditions for velocity and radius to mimic a piston. A corecollapse supernova can be simulated by moving the inner radius inwards (collapse) at a free-fall speed and then violently outwards (bounce and explosion). The parameterization for the piston-driven explosion is the same as that described in Woosley & Weaver (1995), and includes the infall piston time, the final inward piston radius, the initial outward piston speed, and the final  \nIn MESA , the artificial explosion that creates a shock can be produced in three ways: a piston, a luminosity flash, or a  \npiston radius.  \n- GLYPH<15> The second option increases the inner boundary luminosity over a specified time in order to deliver the desired total energy. In this approach, the inner boundary radius is fixed at all times and becomes a zero-flux inner boundary once the flash is over. In this approach, the inner boundary radius is fixed (zero velocity) and we inject the energy within the first zone of the domain.\n- GLYPH<15> The third option deposits energy at a constant rate during a specified time and in a region bounded by two specified Lagrangian-mass coordinates. As in the second option, the inner boundary radius is fixed at all times.  \nThe di GLYPH<11> erences amongst these three options can alter the properties of the shocked envelope.  \nOnce the stellar core has been excised, the remaining envelope has a smooth density profile, resembling a power law whose exponent is GLYPH<0> 1 at depth and decreases outwards to become about GLYPH<0> 10 at the photosphere (top row panels of Figure 16). Because convective accelerations are limited, the energy deposited increases the internal energy within the innermost 0.05 M GLYPH<12> of the grid. The pressure build-up leads to the sudden expansion of the innermost layers and the formation of a mildly supersonic shock (Mach number GLYPH<25> 2). The shock propagates at a velocity in excess of 1000 km s GLYPH<0> 1 initially, but slows to a few 100 km s GLYPH<0> 1 by the time it reaches the stellar surface after 3 GLYPH<2> 10 5 s. The density contrast across this somewhat weak shock is GLYPH<25> 6. For a strong shock, one expects a density jump of 4 for an ideal gas with an adiabatic index of 5 / 3 and a value of 7 for a radiation-dominated gas ( GLYPH<13> = 4 / 3).  \nTo benchmark MESA for these shock tests, we have used the explicit radiation-hydrodynamics code V1D (Livne 1993; Dessart et al. 2010a,b). Options 1 and 3 are implemented in V1D . For the present envelope shock test, and subsequently for the explosion tests, we use option 3 in both codes. We initiate the explosion by depositing a total of 10 49 erg at a constant rate over 10 s between the Lagrangian mass coordinates of 2.40 and 2.45 M GLYPH<12> . This energy deposited is well in excess of the initial binding energy, which is approximately GLYPH<0> 2 GLYPH<2> 10 47 erg. Once the energy injection is over, we save a model which is then used as initial conditions for a shock evolution simulation.  \nThis simulation is analogous to a shock-tube test. However, in the stellar context (realistic stellar envelope, realistic equation of state, spherical expansion), there is no analytical solution for comparison. We thus run the same simulation with the code V1D and include the results in Figure 16. The results agree at multiple times spanning the progression of the shock towards the stellar surface (the times used for comparison are the same to within 1% and the grid resolution is comparable). The sharpness of the shock in the two simulations di GLYPH<11> ers with time and location. In V1D , the artificial viscosity has a physical spread of two grid zones, irrespective of radius, while in this MESA run, the spread is set to 0.1% of the local radius (see Section 4.2).  \nSince the explosion is started as a thermal bomb, the bulk of the energy is initially internal, see Figure 17. As the material expands and accelerates, the kinetic energy increases, mirroring the decrease in internal energy (essentially no energy is used to unbind the envelope). At the time of shock emergence, the internal and kinetic energies are comparable.  \nFigure 17. Top: Evolution of internal energy E int, gravitational energy E grav, kinetic energy E kin, and their sum E tot for the envelope shock test simulated with MESA and V1D . Bottom: Log of cumulative relative error CRE (Equation 58) of the total energy E tot. We neglect sources (nuclear burning) and sinks (radiation losses), which are negligible.  \n<!-- image -->  \nFigure 18. Top: Normalized values of velocity, artificial acceleration g visc, artificial viscous heating GLYPH<15> visc, and Mach number in the vicinity of the shock. The dashed vertical line marks where the Mach number is unity. Bottom: Dependency of the shock morphology on changes in the viscosity parameter l 2. The dots shown for the model with l 2 = 0.004 denote the location of the MESA grid points at that time. For all these tests, we deposit an energy of 10 50 erg at a constant power over 1 s.  \n<!-- image -->  \n/circledot  \nIn the present case, we can preserve good accuracy while  \nstill allowing time steps an order of magnitude greater than the Courant time. 13 The error in energy conservation at timestep i is  \nGLYPH<14> E error ; i = Ei GLYPH<0> Ei GLYPH<0> 1 GLYPH<0> E sources ; i (57)  \nwhere E sources ; i is the right hand side of Equation (52) multiplied by GLYPH<14> t . The cumulative relative error in energy at a timestep n is  \nCRE ( tn ) = 1 En n X i = 1 GLYPH<14> E error ; i : (58)  \nIn the test case, after about 15,100 timesteps when the shock reaches 6.6 M GLYPH<12> , the cumulative relative error has grown to about GLYPH<0> 1 : 4 GLYPH<2> 10 GLYPH<0> 6 , corresponding a roughly linear growth rate of about GLYPH<0> 1 GLYPH<2> 10 GLYPH<0> 10 per step (bottom panel in Figure 17). Note that at this stage of evolution the shock is nearing the outer edge of the envelope but has not actually broken out through the surface. Issues of shock break out are beyond the scope of the current implementation. Using the parameters selected for the test, the energy conservation with V1D is not as good as with MESA (the jumps in cumulative error correspond to times when the limit on the time step are loosened); comparable accuracy can be obtained with V1D by reducing the explicit time step well below the Courant limit.  \nFinally, to illustrate the e GLYPH<11> ects of artificial viscosity, we vary the quadratic term l 2 that controls the spread of the shock in response to compression (see Equation 33), with the explosion energy increased to 10 50 erg in order to produce a stronger shock, and otherwise the same parameters and initial conditions. The top panel of Figure 18 shows the artificial acceleration ( g visc) and energy ( GLYPH<15> visc) terms that enter the momentum and the energy equations for l 2 = 0.001. The acceleration term is positive ahead of the shock, causing a preacceleration of the unshocked material, and negative behind the shock causing a deceleration of the post-shock material. The energy corresponding to those changes in momentum is balanced by the extra term for artificial viscous heating in the energy equation ( GLYPH<15> visc). The lower panel of Figure 18 shows the expected increase in the width of the shock as we raise the parameter l 2. For the model with l 2 = 0.004, dots locate grid cells. Note that with the smallest value ( l 2 = 0.001), the velocity is showing small oscillations ('ringing') behind the shock indicating that we have reached a practical lower bound for the shock spread given the other parameter choices and the nature of the specific problem.", '5. ADVANCED BURNING': 'For the advanced stages of stellar burning, we show here that more accurate summations yield more e GLYPH<14> cient time integrations. This development allows MESA to use large insitu reaction networks. It o GLYPH<11> ers an improvement by providing a single solution methodology that avoids the challenges of stitching together di GLYPH<11> erent solution methods such as nuclear statistical equilibrium (NSE) or co-processing a reaction network. We discuss this development and apply it to the evolution of an X-ray burst on a neutron star. In Section 6 we discuss the pre-supernova progenitors and combine the new capability for advanced burning with the implicit hydrodynamics module to discuss the explosion of core-collapse supernovae.  \n13 The Courant time, equal to the minimum sound crossing time through a grid zone. In this envelope test, it is of the order of 10 s initially, increasing progressively to 40 s prior to shock emergence.  \nThe equations that describe the continuum limit of reacting nuclei are  \n˙ Yi = X j ci GLYPH<21> jYj + X j ; k ci j cj j ! j ck j ! GLYPH<21> j ; kYjYk + X j ; k ; l ci j cj j ! j ck j ! j cl j ! GLYPH<21> j ; k ; lYjYkYl ; (59)  \nwhere Yi is the abundance of isotope i , GLYPH<21> is a reaction rate, and the three sums are over reactions which produce or destroy a nucleus of species i with 1, 2, and 3 reacting nuclei, respectively (e.g., Meyer et al. 1998; Hix & Meyer 2006; Guidry et al. 2013; Longland et al. 2014). The positive or negative stoichiometric coe GLYPH<14> cients ci account for the numbers of nuclei created or destroyed in a reaction. The factorials in the denominators avoid double counting of identical particles.  \nFigure 19. Evolution of the composition for a one-zone burn at constant T = 9.6 GLYPH<2> 10 9 K and GLYPH<26> = 6.0 GLYPH<2> 10 9 g cm GLYPH<0> 3 for 10 6 s starting with a pure 28 Si composition. The calculation uses the mesa 204.net isotope listing (see Table 2), the most abundant isotopes are drawn with thick lines, and several isotopes are labeled. The initial composition is quickly erased as NSE for Y e GLYPH<25> 0.5 is established by GLYPH<25> 10 GLYPH<0> 8 s. Several orders of magnitude in time pass before weak reactions drive a second period of rearrangement. By GLYPH<25> 10 s a second NSE quiescent period with Y e GLYPH<25> 0.403 is established.  \n<!-- image -->  \nFigure 19 shows the evolution of the mass fractions for a MESA one-zone burn at constant T = 9.6 GLYPH<2> 10 9 Kand GLYPH<26> = 6.0 GLYPH<2> 10 9 g cm GLYPH<0> 3 for 10 6 s starting with a pure 28 Si composition. The 204 isotope network, mesa 204.net , used in the calculation is listed in Table 2, and includes the isotopes identified in Heger et al. (2001) as important for Y e in core-collapse models. The thermonuclear reaction rates are from JINA reaclib version V2.0 2013-04-02 (Cyburt et al. 2010). Implementation of reaction rates and associated quantities are described in Paper I and Paper II.  \nThe thermodynamic conditions used in Figure 19 are representative of the central regions of massive stars during the advanced stages of evolution. At such temperatures the initial composition of pure 28 Si undergoes a rapid readjustment. The timescale for the initial Y e GLYPH<25> 0.5 composition to relax to an NSE composition is roughly GLYPH<28> nse GLYPH<25> GLYPH<26> 1 = 5 exp(179 : 7 = T 9 GLYPH<0> 40 : 5) s = 3 GLYPH<2> 10 GLYPH<0> 8 s (Khokhlov 1991; Calder et al. 2007), commensurate with the first burning phase in Figure 19. Between GLYPH<25> 10 GLYPH<0> 8 s and GLYPH<25> 10 GLYPH<0> 4 s the isotopes 4 He and 54 Fe dominate the Y e GLYPH<25> 0.5 NSE composition.  \nTable 2 204 Isotope Network Listing.Table 3 Final Y e for Figure 19  \nTable 3 shows the sensitivity of the final Y e in this calculation to the number of isotopes in the network. Each successively larger network encompasses the previous smaller network and was crafted to yield approximately the same final Y e value as given by the largest network. The 204 isotope network used in Figure 19 is in the regime where larger networks give the same final Y e to 3 significant figures.  \nSince T and GLYPH<26> are constant, only changes to Y e can change the abundances. A second period of intense rearrangement begins at GLYPH<25> 10 GLYPH<0> 4 s and ends at GLYPH<25> 10 s. This activity is driven primarily by p ( e GLYPH<0> ; GLYPH<23> ) n and n ( e + ; ¯ GLYPH<23> ) p and other weak reactions that change Y e. Beyond GLYPH<25> 10 s the isotopes 48 Ca, 49 Ca, and 51 Sc dominate the Y e GLYPH<25> 0.403 NSE composition.', '5.1. More Accurate Summations Yield More E GLYPH<14> cient Integrations': "Wenowtest di GLYPH<11> erent summation methods for Equation (59) and demonstrate that improved accuracy of the summations reduces the number of time steps with a commensurate reduction in the execution time GLYPH<0> while producing the same answers to within the specified integration accuracy.  \nWhen the summations in Equation (59) are accumulated in IEEE 64-bit arithmetic (16 significant figures, real*8 precision in Fortran on most architectures; more specifically binary64 with round to nearest and round ties to even) the integration in Figure 19 takes 3062 time steps using a variableorder Bader-Deuflhard integrator with a specified accuracy of GLYPH<28> int = 10 GLYPH<0> 4 and a scaling value y scale = 10 GLYPH<0> 3 . The specified accuracy GLYPH<28> int limits the maximum error over one time step for any isotope. Other potential, but less demanding, choices for the meaning of GLYPH<28> int include limiting the average or rootmean-square error over one time step for all isotopes. When  \nan abundance is greater than y scale a relative error is calculated, while for abundances smaller than y scale, the absolute error is calculated (e.g., Press et al. 1992). In essence, only abundances greater than y scale can exert control on the size of the time step.  \nAt temperatures larger than GLYPH<25> 5 GLYPH<2> 10 9 K, integrating Equation (59) can be challenging as terms in the summation usually become large and opposite in sign. As shown above, the classic symptom during an integration under these thermodynamic conditions is the integrator taking an excessive number of very small time steps in order to satisfy the specified integration accuracy criteria. The traditional workaround to this numerical problem is abandoning a network integration at elevated temperatures and deploying equilibrium solution methods. This switching of methods raises its own numerical issues when used within the larger context of multi-dimensional simulations or stellar evolution models (see Section 5.2).  \nWhen the summations are accumulated in IEEE 128-bit arithmetic (32 significant figures, real*16 precision in Fortran on most architectures), the same integration takes only 55 time steps, a factor GLYPH<25> 50 improvement in the number of time steps, and a factor of GLYPH<25> 30 less execution time. Both calculations returned the same answers to within the specified integration error tolerances. For tighter integration tolerances of GLYPH<28> int = 10 GLYPH<0> 6 and y scale = 10 GLYPH<0> 5 , the evolution with summations in IEEE 64-bit arithmetic takes 10,081 time steps while the evolution with summations in IEEE 128-bit arithmetic takes 88 time steps. This is a factor of GLYPH<25> 100 improvement in the number of time steps, a factor of GLYPH<25> 150 in execution time, with both calculations again producing the same abundances to within the specified integration error tolerances. Both sets of integration tolerances are practical, everyday usage tolerances; they are not extreme cases of hypothetical interest only. Using low-order Rosenbrock and first-order Euler integrators also showed similar improvements in the number of time steps when the summations were performed in IEEE 128-bit arithmetic instead of IEEE 64-bit arithmetic. We achieve a reduction in the number of timesteps and execution times regardless of the number of isotopes, choice of integrator, integration tolerances, or linear algebra solver. This improvement in e GLYPH<14> -ciency is fundamentally driven by a reduction in the numerical noise of the function being integrated.  \nUnless precautions are taken the summation of large sets of numbers can be very inaccurate due to the accumulation of rounding errors. Methods for accurate summation within the bounds of a given arithmetic remain an active field of research (e.g., Demmel & Hida 2003; McNamee 2004; Ogita et al. 2005; Rump et al. 2008; Graillat & M'enissier-Morain 2012; Collange et al. 2014). These summation discrepancies also worsen on heterogeneous architectures GLYPH<0> such as clusters with NVIDIA GPUs or Xeon Phi accelerators GLYPH<0> which combine programming environments that may obey various floating-point models and o GLYPH<11> er di GLYPH<11> erent precision results.  \nThe summations in Equation (59) for the neutron, proton, and GLYPH<11> -particle abundances are especially prone to inaccuracies because every isotope in a network reacts with these three particles. We report on the summation errors for these three isotopes. Each term in the summations of Equation (59) is calculated using IEEE 64-bit arithmetic and then copied into a IEEE 128-bit variable using the Fortran promotion rules. Each IEEE 128-bit term is then imported into the MP (Brent 1978) and MPf90 (Bailey 1995) multiple precision packages. All other aspects of the integration were executed in IEEE 64-  \nTable 4 Results of Summation Experiments.  \n- a Relative to the 100 digit sum by the MP and MPf90 multiple precision packages.\n- b For a Bader-Deuflard integrator in IEEE 64-bit arithmetic.\n- c For a single thread on one 2.7 GHz Intel Xeon E5 core with the Intel 15.0.1 Fortran compiler, and relative to the execution time for the integration with 128-bit summations with terms in the order given.  \nbit arithmetic. The summations are then accumulated with  \n- GLYPH<15> IEEE 64-bit: terms in the order as given.\n- GLYPH<15> IEEE 64-bit: terms sorted by their absolute value in ascending order and the Kahan (1965) algorithm, which reduces the numerical error in summation by retaining a separate variable to accumulate the errors.\n- GLYPH<15> IEEE 64-bit: terms sorted by their absolute value in ascending order.\n- GLYPH<15> IEEE 128-bit: terms in the order as given.\n- GLYPH<15> MP and MPf90 100 digits: terms sorted by their absolute value in ascending order.\n- GLYPH<15> IEEE 128-bit: terms sorted by their absolute value in ascending order.  \nThere are many summation methods and alternative multiple precision packages we did not deploy in these studies (e.g., Knuth 1997; Higham 2002; Li et al. 2002; Muller et al. 2010; Collange et al. 2014).  \nTable 4 summarizes these summation experiments. We confirm that 100 digits are su GLYPH<14> cient to prevent errors in our multiple precision sums. Column 4 gives the minimum number of correct digits in a summation. There are three time periods in the evolution of Figure 19 where summations performed in IEEE 64-bit arithmetic greatly increase the number of time steps taken by the integration. One is during the first rearrangement into the NSE state ending around 10 GLYPH<0> 8 s, another is during the second rearrangement around 10 GLYPH<0> 1 s and the third time period is when the abundances do not change much ( ˙ Y GLYPH<25> 0) and reaction rates reach equilibrium. It is during these equilibrium periods where a summation in IEEE 64-bit arithmetic with the terms summed in the order given may yield only 6 accurate digits. As a result, 3062 time steps are needed to complete the integration (e.g., row 2 of Table 4). It is important to note that this strategy and choice of arithmetic is commonly used by nuclear reaction networks (e.g., Timmes 1999) GLYPH<0> and the most inaccurate choice. Row 4 of Table 4 is an important case, sorted plus Kahan summation, because it demonstrates that a marginal improvement in the accuracy of the summation (8 minimum correct digits) has a major reduction on the number of time steps (1174 time steps) and execution time (a factor of GLYPH<25> 2.5 smaller). This establishes the general trend that improved accuracy of the summations reduces the number of time steps with a commensurate reduction in the execution  \ntime GLYPH<0> while producing the same answers to within the specified integration accuracy.  \nThe left hand side of Equation (59) for ˙ Yi is a IEEE 64-bit array to be filled with one of the summations. Setting ˙ Yi equal to one of the IEEE 128-bit summations (at least 22 digits of accuracy) or the multiple precision package summations gives the most e GLYPH<14> cient integration (55 time steps, rows 5 and 6 in Table 4) because the Fortran precision demotion rules assure ˙ Yi is accurate to the limit of IEEE 64-bit arithmetic. The next best strategy, but a distant second, is setting ˙ Yi equal to the sorted, Kahan summation. The worst case is setting ˙ Yi equal to the 64-bit arithmetic sum with the terms in the order they appear GLYPH<0> which is a common approach (e.g., Timmes 1999).  \n4  \nFigure 20. Number of accurate digits in 64-bit and 128-bit summations for ˙ Y ( 4 He) as measured by the 100 digit sum calculated by the multiple precision packages. The x-axis gives the time step number for the integration done with IEEE 128-bit summations. The number of accurate digits in ˙ Y (p) and ˙ Y (n) are within a few digits of ˙ Y ( 4 He).  \n<!-- image -->  \nFigure 20 shows the number of correct digits in 64-bit and 128-bit summations for ˙ Y ( 4 He) with the terms accumulated in the order they are given. The number of correct digits is measured against the 100 digit sum calculated by the multiple precision packages MP and MPf90 . The choices for the inte-  \nrator and integration tolerances are the same as in Figure 19. Figure 20 shows that the minimum number of accurate digits is usually within a few digits of the limit of IEEE 64-bit arithmetic, but degrades to 6 digits (see row 1 of Table 4) during a time period of intense isotope rearrangement. These relatively large inaccurate summations cause the right hand side of Equation (59) to be poorly defined in IEEE 64-bit arithmetic. As a direct result, the integration of Equation (59) with IEEE 64-bit summations takes 3062 time steps to complete. Sorting the terms in the sums in ascending order and using the Kahan summation algorithm results in an accurate digit pattern that is very similar except the number of accurate digits is improved by one or two (see row 3 of Table 4). As a direct result of the improved accuracy of the summations, the number of timesteps is reduced from 3062 to 1141.  \nWe found no notable improvements by increasing the accuracy of the summations for the Jacobian matrix used by the sti GLYPH<11> ordinary di GLYPH<11> erential equation integrators. For this problem, it is evidently more important to better define the function GLYPH<0> the right hand side of Equation (59) GLYPH<0> than the Jacobian matrix holding the derivatives of the function.  \nFor the IEEE 128-bit sum relative to the 100 digit sum in Figure 20, the minimum number of accurate digits is usually near the limit of IEEE 128-bit arithmetic, but degrades to 21 digits (see row 5 of Table 4) during the second period of intense isotope rearrangement. Relative to the IEEE 64-bit summations the number of correct digits is improved by at least 15, consistent with our conversion of the IEEE 64-bit terms in the sum to IEEE 128-bit using the Fortran promotion rules. As a direct result of the improved accuracy of the summations the integration that used the IEEE 128-bit summation took only 55 timesteps to complete the evolution.  \nBased on these experiments we currently choose to improve the accuracy of the summations by converting the terms in the order they appear from IEEE 64-bit to IEEE 128-bit using the Fortran promotion rules and adding the terms in IEEE 128-bit arithmetic. IEEE 128-bit precision is presently almost always implemented in software by a variety of techniques (e.g., double-double methods), since direct hardware support for IEEE 128-bit precision is presently rare. However, Table 4 shows the reduction in the number of time steps from accumulating the sums in IEEE 128-bit arithmetic far exceeds the extra computational cost per addition.", '5.2. A Uniform Solution Method for Nuclear Burning in Stellar Evolution': 'At high temperatures, the traditional workaround for the numerical problem of inaccurate summations in IEEE 64-bit arithmetic is to forgo using a reaction network integration to evolve the abundances and nuclear energy generation rate and to replace it with equilibrium solution techniques. An example of such an equilibrium calculation is NSE, where a root-find for the neutron and proton chemical potentials is performed. Once these two chemical potentials are known, all the abundances can be determined from nuclear Saha equations (e.g., Cli GLYPH<11> ord & Tayler 1965; Hartmann et al. 1985; Meyer et al. 1998; Nadyozhin & Yudin 2004; Seitenzahl et al. 2008; Odrzywolek 2012). Equilibrium solution methods by themselves are e GLYPH<14> cient, robust, and inexpensive.  \nHowever, combining reaction networks and equilibrium solution methods creates its own numerical issues, especially when the temperature and density are spatial and time dependent. For example, the temperature of a cell may start relatively low, move into quasi-static equilibrium (QSE) range  \nabove 3 GLYPH<2> 10 9 K, and then move into NSE range above 5 GLYPH<2> 10 9 K. Ad-hoc decision trees must be created for switching between a network integration, QSE solutions, and NSE solutions. These switches can introduce unphysical discontinuities in the abundances either from one timestep to the next or in the abundance spatial profiles from one cell to the next.  \nFinally, equilibrium methods determine the composition at a fixed electron fraction Y e. It then becomes necessary to solve an ordinary di GLYPH<11> erential equation for ˙ Y e based on weak reaction rates in order to advance the abundance solution with a time varying Y e (McLaughlin et al. 1996; Townsley et al. 2009; Arcones et al. 2010, also see Section 8). Switching between integration and equilibrium methods mid-stream is a liability, not a positive asset.  \nFurthermore, cells near the transition between a network integration and an equilibrium method can be unstable in the sense that the equilibrium solution can evolve a cell to lower temperatures pushing the cell into using a network integration, while the solution from the network integration can evolve the cell towards higher temperatures evolving the cell back towards using the equilibrium solution. Moreover, the reaction network used for the time integration is di GLYPH<11> erent (usually smaller) than the isotope listing used for the equilibrium methods. This necessitates crafting a delicate mapping between two abundance vectors, which may also introduce unphysical discontinuities. In addition, care must be taken to assure the reaction rate screening corrections used in the time integration are properly taken into account in the equilibrium solution method, otherwise a fundamental incompatibility exists between the abundance vectors.  \nThe need for traditional workarounds forced by limited accuracy of the summation is now avoided. The summation experiments in Section 5.1 demonstrate that network integration can be robust and e GLYPH<14> cient, even at very high temperatures, when the accuracy of the summations is improved. We stress this is not just a solution to issues of limited accuracy. It also o GLYPH<11> ers an improvement in MESA by providing a single solution methodology, network integration, that avoids the challenges of stitching together di GLYPH<11> erent solution methods.', '5.3. X-ray Burst Models and Adaptive Nets': "The new capabilities described above allow MESAstar to use large in-situ reaction networks (i.e. fully coupled to the stellar evolution rather than uncoupled co-processing). A demonstration is Type 1 X-ray bursts, a class of objects with unstable nuclear burning on the surface of a neutron star (NS). These bursts are sensitive functions of accretion rate (Chen et al. 1997), accretion composition (Galloway et al. 2006), the spatial distribution of burning on the surface of the NS (Bildsten 1995), the type of burning that occurs between bursts (Galloway et al. 2008) as well as possibly other conditions, for instance 'superbursts' where carbon, rather than H / He, burns (Cumming & Bildsten 2001). Here we focus on a simplified model of constant accretion rate, where the burning occurs over the whole surface of the NS. GS 1826-24 (Tanaka 1989), also known as the 'clocked burster' (Ubertini et al. 1999), provides an example of such a system due to its regular Type 1 X-ray bursts.  \nAs material is accreted at the surface of a NS it is compressed and heats the underlying material. The accreted hydrogen (from a low mass main sequence (MS) star (Chen et al. 1997)) burns via the hot CNO cycle. However, with high enough accretion rates the hydrogen will be accreted faster than the hot CNO cycle, which is limited by the GLYPH<12> -decay  \nTable 5 Recurrence times of X-ray bursts.  \nFigure 21. Kippenhahn plot during two X-ray bursts for the rp 305 net with the solar metallicity accretion model. The x-axis values are times relative to the peak of each burst, note the non-linearity of the scale. The y-axis values are the column depth and the color coding shows the temperature of the NS envelope. The dashed contours show the extent of the convective regions.  \n<!-- image -->  \ntimescale (of order minutes), can process the material. The accreted helium ignites unstably in a hydrogen rich environment, allowing rapid proton ( rp ) captures onto seed nuclei (Wallace & Woosley 1981). This process forms nuclei along the proton drip line up to and beyond the iron group (Schatz et al. 1999), peaking at 107 Te, when GLYPH<11> -decays prevent heavier elements from being formed (Schatz et al. 2001; Fisker et al. 2008). Once the burst begins, convection will commence, mixing the freshly burnt material with the ashes of previous burning episodes (Weinberg et al. 2006).  \nWemodel the NS envelope using inner boundary conditions for mass and radius of M c = 1 : 4 M GLYPH<12> and R c = 11 : 2 km (Heger et al. 2007), implying a gravitational redshift of 1 + z = 1 : 26. The base of the envelope is composed of an inert layer that does not undergo reactions. The luminosity at the base of the envelope is set to L = 1 : 6 GLYPH<2> 10 34 ergs s GLYPH<0> 1 (Woosley et al. 2004). We base our nuclear networks on the 304 species rp.net net-  \nGS 1826-24 has been studied by the Rossi X-Ray Timing Explorer (RXTE) over several years (Galloway et al. 2004, 2008). The bursts showed a decrease in the recurrence time between bursts, from 4.1 hr in 2000 to 3.56 hr in 2002, though during each observational epoch the bursts were consistent with each other. Based on the ratio of the burst energy to the persistent flux, it is assumed that the bursts are powered by hydrogen burning of solar metallicity material.  \nFigure 22. The folded burst profiles for the di GLYPH<11> erent nuclear networks as compared to GS 1826-24 for an accretion rate of 3 GLYPH<2> 10 GLYPH<0> 9 M GLYPH<12> yr GLYPH<0> 1 with 2% metals with a solar composition, Three rp network models are shown and one of the adaptive net models. The insert shows a zoom in of the first 30 s during the burst.  \n<!-- image -->  \nFigure 23. Folded burst light curve for the rp 305 net, with a solar metallicity accretion composition, shown for ˙ M = 2 : 4 GLYPH<2> 10 GLYPH<0> 9 M GLYPH<12> yr GLYPH<0> 1 and ˙ M = 3 GLYPH<2> 10 GLYPH<0> 9 M GLYPH<12> yr GLYPH<0> 1 , normalized to the peak flux measured for GS 182624. The insert shows the first 30 s of the burst.  \n<!-- image -->  \nwork of Fisker et al. (2008), which includes proton rich isotopes up to 107 Te. Isotopes above 66 Zn, which is the peak isotope in the mesa 204.net , are included due to the proton captures possible on highZ isotopes during the peak of the burst (Fisker et al. 2006). We also include the e GLYPH<11> ects of rotational mixing by setting a minimum amount of mixing in the NS envelope. This mixing, while having a physical motivation (Piro & Bildsten 2007; Keek et al. 2009), is there primarily to improve the convergence of MESA models by smoothing out the compositional gradients that form in the ashes of previous bursts. We include the post-Newtonian correction to correct the local gravity in each cell for GR e GLYPH<11> ects. During the burst we allow the accretion to continue.  \nFigure 21 shows the temperature profile during two X-ray bursts, for the rp 305 net, accreting solar metallicity material at a rate of 3 : 0 GLYPH<2> 10 GLYPH<0> 9 M GLYPH<12> yr GLYPH<0> 1 . At t 0 GLYPH<25> GLYPH<0> 10 s the envelope ignites material and drives the formation of the first convection zone. This zone expands outwards in the envelope mixing the ashes from the burning at the base of the envelope outwards to lower pressures (Weinberg et al. 2006). As the burst decays the convection zone recedes outwards and by t 0 GLYPH<25> 150 s the envelope returns to its pre-burst temperature profile.  \nOur results are compared to the RXTE observations of GS 1826-24 over bursts 9-20. Time resolved spectra were binned during the bursts' rise time and decay (Galloway et al. 2008; Zamfir et al. 2012). Data output by MESA is not GR time corrected, thus we set the burst times to be t 0 = t (1 + z), and average multiple bursts to produce a scaled light-curve.  \nWe test three reaction networks, rp 53 , rp 153 and rp 305 , each a modified form of that in Fisker et al. (2008). Table 5 shows that increasing the number of isotopes in the reaction network increases the recurrence time and that all (for ˙ M = 3 : 0 GLYPH<2> 10 GLYPH<0> 9 M GLYPH<12> yr GLYPH<0> 1 ) have recurrence times GLYPH<25> 1 hr less than that of GS 1826-24.  \nTo achieve a better match to the GS 1826-24 recurrence time (see Table 5), we reduce the accretion rate to ˙ M = 2 : 4 GLYPH<2> 10 GLYPH<0> 9 M GLYPH<12> yr GLYPH<0> 1 . However, Figure 23 shows that the light curve comparisons are not as good as for the higher ˙ M .  \nFigure 22 and its insert show the folded light curves for each of the three rp reaction networks plus the GS 1826-24 observations. The rise time is sensitive to the net, with the largest net matching the observed slow rise. The observed decay profile is also best matched by rp 305 . Burst to burst variations of the models decrease with increasing net size and can be further reduced by increasing the temporal resolution of the models. However, increasing the size of the net reduces the variation without having to increase the temporal resolution and also highlights the impact of MESA capability to include large nuclear networks.  \nGS 1826-24 was also modeled by Heger et al. (2007) with accretion of hydrogen, helium and 2% 14 N. For comparison, we run a model with this same composition with the rp 305 net. Table 5 shows that the recurrence time decreases slightly when accreting 2% 14 N rather than 2% metals. The model with metal accretion is in better agreement with both the light curve rise and decay.  \nWe now explore adaptive nets (Woosley et al. 2004), where we allow MESA to determine which isotopes (and reactions) are necessary by assessing the available reaction pathways for the most abundant isotopes. The network is constructed by first finding those isotopes with an abundance above a threshold, X keep, and then introducing those isotopes which are connected by adding or removing protons, neutrons, or GLYPH<11>  \nparticles. That determination is made via the additional parameters X n (i.e. neutron reactions) and X p (i.e. proton and GLYPH<11> reactions) potentially re-adding isotopes removed with the initiating X keep threshold.  \nAccreting solar composition material at ˙ M = 3 : 0 GLYPH<2> 10 GLYPH<0> 9 M GLYPH<12> yr GLYPH<0> 1 we follow the model to the second burst, finding a recurrence time of 3 : 1 hrs, comparable to that from the rp 305 net (Table 5). The adaptive net has a better rise time profile than the rp nets, while the rp 305 net has a better fit to the decay. This gives us confidence that the rp 305 net includes all relevant isotopes which drive the X-ray burst and thus is a useful approximation. For suitable values for the sensitivity of the adaptive net, the net limits itself to GLYPH<25> 400 isotopes between bursts, which increases to GLYPH<25> 600 isotopes during the burst. Variations of a factor 100 in the threshold parameters only change the isotope count by at most 50 isotopes and do not a GLYPH<11> ect the final results.", '6. CORE-COLLAPSE SUPERNOVAE': 'The capability of using large, in-situ reaction networks without the need for equilibrium or co-processing techniques was described in Section 5 and applied to X-ray burst models. Weextend our demonstration of this capability by first considering pre-supernova models. We then combine the advanced burning development with the implicit treatment of shocks discussed in Section 4 to core-collapse supernovae models.  \nFigure 24. Evolution of T c and GLYPH<26> c in solar metallicity, non-rotating Mi = 15 and 30 M GLYPH<12> pre-supernova models. The curves are calculated using an in-situ 204 isotope reaction network. Locations of the core carbon, neon, oxygen, and silicon ignition are labeled, as is the scaling relation T 3 c / GLYPH<26> c, and the E F = k B T GLYPH<25> 4 electron degeneracy curve. Regions dominated by electronpositron pairs, photodisintegration, and rapid electron capture are shaded and labeled.  \n<!-- image -->', '6.1. Pre-Supernova Evolution without QSE or NSE': "Figure 24 shows the T c GLYPH<0> GLYPH<26> c evolution of Mi = 15 and 30 M GLYPH<12> models from the onset of carbon burning until ironcore collapse. These non-rotating, solar metallicity models used the 204 isotope reaction network described in Section 5 and MESAstar 's 'Dutch' mass loss prescription with GLYPH<17> = 0.8. These models have GLYPH<25> 2200 cells on the main-sequence, GLYPH<25> 3500 cells as the star becomes a red supergiant, and GLYPH<25> 2300 cells at the onset of core collapse. At core collapse the final  \nFigure 27 shows the mass fraction profiles of the ten most abundant isotopes within the iron core at the onset of core collapse for the Mi = 15 M GLYPH<12> model evolved with the 204 isotope network. Each isotope shown dominates the NSE composition at some location within the iron core, although we stress that no NSE or QSE approximation was used; the same 204 isotope reaction network was used throughout the entire model from the pre-main-sequence to the onset of core collapse.  \n<!-- image -->  \n/circledot  \nFigure 25. Radial velocity and Y e profiles at the onset of core collapse for the Mi = 15 M GLYPH<12> model. Dashed curves show the results using a 22 isotope network and solid curves show the results using a 204 isotope network. Both models are evolved from the pre main-sequence to the onset of core collapse with their respective reaction network. The vertical gray lines mark the mass of the iron core as defined by the Y e jump.  \n<!-- image -->  \n/circledot  \nFigure 26. Thermodynamic profiles at the onset of core collapse for the Mi = 15 M GLYPH<12> model. Dashed curves show the results using a 22 isotope network and solid curves show the results using a 204 isotope network. Both models are evolved from the pre main-sequence to the onset of core collapse with their respective reaction network. The vertical gray lines mark the mass of the iron core as defined by the Y e jump.  \nmasses are Mf = 13.0 and 15.2 M GLYPH<12> . The curves fall below the T 3 c / GLYPH<26> c scaling relation because the core becomes partially electron degenerate, as indicated by tracks crossing the Fermi energy E F = ( k B T ) GLYPH<25> 4 curve. Evolution towards lower density at nearly constant temperature signals ignition of a nuclear fuel.  \nFigure 25 shows the radial velocity and Y e profiles at the onset of core collapse for the Mi = 15 M GLYPH<12> model. Dashed curves show the results using a 22 isotope network and solid curves show the results using a 204 isotope network. Both models are evolved from the pre main-sequence to the onset of core collapse with their respective reaction network. The vertical gray lines mark the mass of the iron core as defined by the Y e jump, which is m GLYPH<25> 1.43 M GLYPH<12> for the 204 isotope model and m GLYPH<25> 1.59 M GLYPH<12> for the 22 isotope model. The infall speed has reached GLYPH<25> 1000 km s GLYPH<0> 1 just inside these iron core locations.  \n<!-- image -->  \n/circledot  \nFigure 27. Mass fraction profiles of the ten most abundant isotopes within the iron core at the onset of core collapse for the Mi = 15 M GLYPH<12> model evolved with the 204 isotope network. The entire iron core is in NSE and the mass fractions adapt to the changing temperature, density and Y e.  \nFigure 26 shows the thermodynamic profiles at the onset of core collapse for the Mi = 15 M GLYPH<12> model. Dashed curves again show the results using a 22 isotope network and solid curves show the results using a 204 isotope network. The vertical gray lines again mark the mass of the iron core as defined by the Y e jump in Figure 25. The impact of these di GLYPH<11> erences remains to be explored.  \nThe most abundant isotopes in an NSE distribution generally have an individual Y e that is within a small range of the local Y e. A small spread usually exists due to nuclear structure e GLYPH<11> ects. For example, the dominant isotopes at the center in Figure 27 are 49 Sc and 48 Ca. These isotopes have individual Y e of 0.429 and 0.417, respectively; commensurate with the central Y e GLYPH<25> 0.428 shown in Figure 26. The dominant isotope changes as the NSE distribution adapts to the rapidly decreasing density profile and increasing Y e profile. All the isotopes in the iron core eventually become part of the compact remnant after the explosion. However, the thermodynamic and composition profiles near the mass cut depend on the profiles interior to the mass cut.", '6.2. Core-collapse Supernova Explosions': 'The envelope shock tests described in Section 4.9 show that the hydrodynamic solver in MESA meets the basic requirements for shock propagation in a star. The AGB star model was selected because of the well behaved conditions of its envelope - a density structure that is smooth and monotonically declining towards the stellar surface, and a uniform composition.  \nHere we explore the more challenging conditions associated with a strong shock born deep in the stellar interior of a massive star. We study the dynamics of such a supernova  \n3  \nFigure 28. Multi-epoch snapshots of the hydrodynamical simulation from energy injection mimicking core-collapse supernova. The initial model is a 15 M GLYPH<12> star at solar metallicity, evolved with mass loss but no rotation, and employing a nuclear network of 22 isotopes. We show the density (top row), temperature (middle row), and velocity (bottom row), versus Lagrangian mass (left column) and radius (right column). In each panel, the solid line shows the MESA results and the dashed line the V1D results.  \n<!-- image -->  \n/circledot  \n/circledot  \nshock and the explosive nucleosynthesis that takes place in the wake of the shock during the first second. The yields from explosive nucleosynthesis depend on both the energy and the power (characteristic energy deposition timescale), while the dynamics of the shock are primarily dependent on the total energy deposited.  \nThe starting conditions for the explosion simulations are the two 15 M GLYPH<12> pre-supernova models discussed above; one for the approximate 22 isotope network and one for the 204 isotope network.  \n6.2.1. Explosion Dynamics  \nSince the focus here is on the dynamics rather than nucleosynthesis, we expedite the MESA simulation by using the 22 isotope network for both the pre-supernova and the explosion phases. Before triggering the explosion, we remove the iron core of the red-supergiant star by placing the inner boundary of the grid at a Lagrangian mass of 1.75 M GLYPH<12> . We trigger the explosion by depositing 1.52 GLYPH<2> 10 51 erg at a constant rate during 1 s. The artificial viscosity is raised during the energy deposition phase (i.e., l 2 = 0.01), when the shock is at small radii, and lowered in the subsequent evolution until shock breakout (i.e., l 2 = 0.003). Since the binding energy of the envelope to be shocked is GLYPH<0> 3.2 GLYPH<2> 10 50 erg at the time we trigger the explosion, this choice of energy deposition yields a total energy at the end of the deposition phase of 1.2 GLYPH<2> 10 51 erg. This is generally considered a standard value for a core-collapse supernova.  \nThe radii of the innermost shells are initially very small since they lie at the outer edge of the iron core. Consequently, they su GLYPH<11> er considerable cooling from expansion. Figure 28 shows a drop in temperature from a few 10 9 K down to GLYPH<24> 10 4 K at GLYPH<24> 1 day. In addition, the supernova shock splits into a reverse / forward shock structure when it encounters the density drop at the transition between the He-rich core and the H-rich envelope. The reverse shock is the new feature, absent in the envelope shock test, that causes a significant deceleration of He-core material. The conversion of kinetic energy into internal energy causes this inner material to heat up, erasing the cooling e GLYPH<11> ect from expansion. The innermost layers, which travel the slowest, will be shocked last. These innermost zones can evolve to temperatures GLYPH<24> 10 4 K. It is in these innermost regions at late times that the di GLYPH<11> erences between MESA and V1D are the largest. The o GLYPH<11> set occurs in a region of relatively high density ( GLYPH<24> 10 GLYPH<0> 7 g cm GLYPH<0> 3 ) and low temperature ( GLYPH<24> 10 4 K). The o GLYPH<11> set in temperature between the MESA and V1D simulations at late times stems from a di GLYPH<11> erence in the equation of state for metal-rich regions. MESA accounts for ionization through the OPAL equation of state table for metallicities z < 0.04. For higher metal abundances where OPAL tables are unavailable, MESA currently assumes full ionization while V1D solves for the ionization state of the gas. Note that such density / temperature regimes are normally not encountered in stellar interiors. For other quantities and / or locations / times, the agreement between MESA and V1D is excellent.  \nFigure 28 shows that the development of the explosion is analogous to the tests using the (low-density) envelope of an AGB star in Section 4.9, but with significant quantitative differences. Here, the shock born at the edge of the iron core first travels through the dense CO-rich core. At the outer edge of the He-rich shell, the shock traverses a steep density gradient to enter the low-density H-rich envelope. Hence, the shock crosses regions with densities ranging from GLYPH<24> 10 6 g cm GLYPH<0> 3 at the edge of the iron core down to 10 GLYPH<0> 10 g cm GLYPH<0> 3 at the stellar surface.  \nWe also note that in the MESA simulation, two small spikes appear in the temperature and density profiles at < 2.5 M GLYPH<12> at GLYPH<25> 1000 s after the energy deposition phase. This feature is absent in V1D because V1D uses a much larger viscous spread when the shock is in the helium-rich core ( R < R GLYPH<12> ). One can reduce or eliminate such spikes by increasing the viscous spread, although this may visibly smear the shock when it crosses the H-rich envelope - the current choice seems a suitable compromise.  \nIn contrast to the envelope shock test, this supernova explosion configuration raises the temperature by a factor of about ten. Consequently, because P rad = P gas / T 3 =GLYPH<26> , the post-shock material becomes completely radiation dominated ( P rad GLYPH<29> P gas). If we neglect the binding energy and the kinetic energy of the post-shock material, the post-shock energy is of the order of the explosion energy. We indeed find a good correspondence between the post-shock temperature computed by MESA and the temperature obtained from ( E 0 = aV ) 1 = 4 (where a is the radiation constant, E 0 is a fitting parameter, typically of the order of the explosion energy, and V = 4 GLYPH<25> R 3 sh = 3 is the volume within the shock radius R sh). As expected, we also find that the shock accelerates (decelerates) in regions where GLYPH<26> sh R 3 sh decreases (increases) outward.', '6.2.2. Explosive Nucleosynthesis': "Here we compare the shock nucleosynthesis results from the two independent codes, MESA and V1D . The same initial 204 isotope pre-supernova model was the starting point. Our first test case is a strong explosion triggered by injecting 1.57 GLYPH<2> 10 51 erg for 0.05 s and within 0.02 M GLYPH<12> of the mass cut, which is positioned at the outer edge of the iron core at 1.5 M GLYPH<12> . The exact choice of explosion energy, deposition time scale, and mass cut is not strictly relevant.  \nOverall the agreement between MESA and V1D on this strong explosion is very good. The small di GLYPH<11> erences between MESA and V1D in Figures 29 and 30 are due to the di GLYPH<11> erence in mapping procedures.  \nFigures 29 and 30 compare the mass fraction profiles of MESA with a 22 isotope network, MESA with a 204 isotope network, and V1D with a 54 isotope network. The first comparison at 0.0 s shows the impact of mapping from the presupernova 204 isotope network to the networks used in the shock nucleosynthesis test. The next comparison at 0.05 s is at the end of the energy deposition phase. The final comparison at 42.7 s is after explosive nucleosynthesis has completed. In all cases, the silicon-rich and oxygen-rich shells are strongly influenced by the explosion; the former primarily for the production of 56 Ni and the latter primarily for the production of 28 Si and 32 S. The 56 Ni yields at 42.7 s are 0.092 M GLYPH<12> for V1D , 0.087 M GLYPH<12> for MESA with 22 isotopes, and 0.096 M GLYPH<12> for MESA with 204 isotopes.  \nMESA follows rules for mapping isotopes from one network to another network: If an isotope present in the old network is also present in the new network, then the abundance from the old network is copied to the abundance for the new network. Isotopes in the new network that are not present in the old network are initially given a mass fraction of zero. MESA then separately renormalizes classes of isotopes to have the same total mass fraction in the new network as in the old network. The classes are neutrons, hydrogen, helium, carbon, nitrogen, oxygen, and other metals. This procedure guarantees that the sum of the mass fractions in a given class will be the same in the new network as in the old network. V1D 's mapping procedure for isotopes is the following. For any isotope present in the V1D network but absent in the MESA input the mass fraction is set to the solar metallicity value. When an isotope included in the MESA input is absent in the V1D network, this isotope is left out in the V1D simulation. After completing the mappings, the resulting composition is renormalized so that the sum of the mass fractions is unity.  \nOur second test case is a lower power explosion. We inject 1.326 GLYPH<2> 10 51 erg in 1.0 s over 0.05 M GLYPH<12> of the mass cut, which  \nlog X  \nlog X  \nFigure 31. Composition profiles of the eight most abundant isotopes at 1 s in the inner 0.4 M GLYPH<12> of the ejecta for MESA (solid) and V1D (dashed) simulations of a 10 51 erg explosion in the 15 M GLYPH<12> model. The Si-rich and O-rich shells have been influenced by explosive nucleosynthesis, the former primarily for the production of 56 Ni and the latter primarily for the production of 28 Si and 32 S.  \n<!-- image -->  \n/circledot  \nFigure 29. Nucleosynthesis profiles of selected isotopes for the 1.57 GLYPH<2> 10 51 erg energy deposition test case at 0.0 s (top) and 0.05 s (bottom). The dashed lines show the MESA results with a 22 isotope network, the solid lines show the MESA results with a 204 isotope network, and the long dashed lines show the V1D results with a 54 isotope network.Figure 30. Same as Figure 29 but at a time after all nucleosynthesis has completed.  \n<!-- image -->  \nis also positioned at the outer edge of the iron core at 1.5 M GLYPH<12> . The total energy after the deposition phase is 10 51 erg.  \n<!-- image -->  \n/circledot  \nIn Figure 31, we show the composition profiles for the eight most abundant isotopes in the inner GLYPH<25> 0.4 M GLYPH<12> at the end of the energy-deposition phase (i.e., at 1 s). The correspondence between MESA and V1D is again very good. The 56 Ni mass fraction approaches unity - it would reach unity if we appreciably increased the power (see Figure 29 for example). Some 58 Ni is produced in the same region, while 54 Fe is synthesized in the layers immediately above. The 56 Ni yields at 1.0 s are 0.0041 M GLYPH<12> for V1D , and 0.011 M GLYPH<12> for MESA with 204 isotopes.  \nThis work shows that the power of the explosion has a a significant impact on the abundance profiles. In the high power explosion, the yield of 56 Ni is GLYPH<25> 10 times larger and the 4 He is several orders of magnitude more abundant. The nucleosynthesis of the low power explosion is completed at end of deposition phase at 1.0 s, while nucleosynthesis in the high power case continues for GLYPH<25> 30 s. This sensitivity suggests potentially observable signatures between low and high power explosions. In addition, the explosive nucleosynthesis that takes place in core-collapse supernovae is sensitive to the way the explosion is triggered. With the approach we use (fixed power during the energy deposition phase), we find that increasing the explosion energy (at a given power), the power (at a given explosion energy), or both alters the amount of mass burnt. Moving the mass cut deeper into denser layers considerably enhances the amount of burnt material but this material may fall back rather than be ejected. Moving the mass cut further out into lower-density regions may completely quench the production of 56 Ni, in favor, for example, of 28 Si. It is thus important to keep in mind that the piston or thermal explosion trigger is artificial and that the yields from explosive nucleosynthesis bear significant uncertainties.", '7. IMPROVED TREATMENT OF MASS ACCRETION': "Adding mass to a star requires a way to accurately and e GLYPH<14> -ciently compute the thermal state of the freshly accreted material in the outermost layers. This is simplified by a hierarchy of timescales. For accretion at ˙ M , there are two important timescales at a given location, m , the thermal time GLYPH<28> th ' ( M GLYPH<0> m ) CPT = L , where CP is the specific heat at constant  \nP , and the time to accrete this same layer, GLYPH<28> acc ' ( M GLYPH<0> m ) = ˙ M . Near the surface, L GLYPH<29> CPT ˙ M , implying that GLYPH<28> th GLYPH<28> GLYPH<28> acc, so that these layers have ample time to relax to the thermal equilibrium configuration fixed by L (Nomoto & Sugimoto 1977; Nomoto 1982; Townsley & Bildsten 2004).  \nThe timescale hierarchy GLYPH<28> th GLYPH<28> GLYPH<28> acc implies that the outer regions evolve nearly homologously in the fractional mass coordinate q = m = M (Sugimoto & Nomoto 1975). Hence, the thermal profile (e.g. the run of T with P or GLYPH<26> ) of the outer layer is nearly constant in time even as fluid elements are compressed to higher pressures and have T ( m ) increase. More formally, T ( q ) varies slowly in time near the surface, where (1 GLYPH<0> q ) GLYPH<28> 1. This motivates reformulating the Lagrangian based form of  \nIn cases where L arises solely from compression of material, such as a very rapidly accreting star of high ˙ M or an old, cold accreting WD, then L GLYPH<24> CPT b ˙ M , where T b is the temperature at the degenerate / nondegenerate transition in a WD (Townsley & Bildsten 2004) or of the core in a normal star. Even in these cases, the outer layers have T GLYPH<28> Tb , allowing the inequality GLYPH<28> th GLYPH<28> GLYPH<28> acc to hold. This also implies that the thermal state of the arriving material is unimportant, allowing us to safely use the approximation that material arrives with the same entropy as the photosphere, since material relaxes toward this on the very short GLYPH<28> th at the photosphere. 14 Even when ˙ M varies on short timescales, using an averaged accretion rate is a good approximation for computing the evolution of the interior layers due to their long GLYPH<28> th (Piro et al. 2005; Townsley & Gansicke 2009).  \nGLYPH<15> grav = GLYPH<0> T Ds Dt GLYPH<17> GLYPH<0> T @ s @ t ! m (60)  \nthat is needed in the energy equation, @ L =@ m = GLYPH<15> grav + GLYPH<15> nuc GLYPH<0> GLYPH<15>GLYPH<23> (Sugimoto 1970; Sugimoto & Nomoto 1975; Paper I), to a version in the coordinate q ,  \nGLYPH<15> grav = GLYPH<0> T @ s @ t ! q + T d ln M ( t ) dt @ s @ ln q ! t : (61)  \nSugimoto & Nomoto (1975), and later works based on it, denote the second term on the right hand side the 'homologous' term. Physically it is the local loss of entropy in the fluid element as it is compressed to higher pressure. They label the first term on the right hand side the 'non-homologous' term. It arises from the much slower departure of the outer layers from simple homologous evolution on a timescale M = ˙ M .  \nIn Paper II, following the work of Townsley & Bildsten (2004), only the homologous term in Equation (61) was included in GLYPH<15> grav in and near regions of newly accreted material. We also described the huge advantage of such an implementation, as it allows the mass added per timestep to be much larger than the smallest cell mass near the surface, while maintaining accurate thermal profiles at low pressures. However, leaving out the non-homologous term can create a discontinuity in GLYPH<15> grav at the location where the standard Lagrangian derivative, Equation (60), begins to be used. We now de-  \nMESAstar includes the ability to have an inner inert core of mass Mc . In this situation q = ( m GLYPH<0> Mc ) = ( M GLYPH<0> Mc ) rather than the more typical m = M . For simplicity here, we will use Mc = 0. Approximate homology holds in either case for 1 GLYPH<0> q GLYPH<28> 1.  \n14 Apossible exception to this case is rapidly accreting pre-main sequence stars where the accretion shock is so optically thick that the material's entropy remains high and is advected inward (Palla & Stahler 1990).  \nscribe the improvement we have made to MESAstar so that it now includes both the homologous and non-homologous terms. Hence, the two forms of Ds = Dt are physically equivalent and there is no longer any discontinuity.", '7.1. Lagrangian and Homologous Regions': 'The independent coordinates used for writing the timedependent structure of the star are m and t , and for fluid elements deep within the star at both timesteps, the conventional form of Equation (60) is adequate. One numeric subtlety of accretion is that the derivative at constant m cannot be evaluated for material that is not present in the star at the beginning of the timestep. However, the simplification available when GLYPH<28> th GLYPH<28> GLYPH<28> acc, manifest in Equation (61), enables GLYPH<15> grav to be evaluated in the outer regions.  \nWhen T and GLYPH<26> are used as independent variables, we write  \nGLYPH<15> grav = GLYPH<0> CPT " (1 GLYPH<0> r ad GLYPH<31> T ) @ ln T @ t ! m GLYPH<0> r ad GLYPH<31>GLYPH<26> @ ln GLYPH<26> @ t ! m # (62)  \nin the interior of the star, as in Paper I, and near the surface we choose to write, using Equation (61),  \nGLYPH<15> grav = GLYPH<15> grav ; nh + GLYPH<15> grav ; h ; (63)  \nwhere  \nGLYPH<15> grav ; nh = GLYPH<0> CPT 2 6 6 6 6 4 (1 GLYPH<0> r ad GLYPH<31> T ) @ ln T @ t ! q GLYPH<0> r ad GLYPH<31>GLYPH<26> @ ln GLYPH<26> @ t ! q 3 7 7 7 7 5 (64)  \nGLYPH<15> grav ; h = CPTGm ˙ M 4 GLYPH<25> r 4 P ( r ad GLYPH<0> r T ) : (65)  \nand  \nHere r T = d ln T = d ln P is the T -P profile in the star, and we have used the thermodynamic derivatives r ad = ( @ ln T =@ ln P ) s , GLYPH<31> T = ( @ P =@ T ) GLYPH<26> , and GLYPH<31>GLYPH<26> = ( @ P =@GLYPH<26> ) T . There is also a transition region where a weighted combination of these forms is used, with weights varying linearly in m .  \n<!-- image -->  \nFigure 32. Illustration of the MESAstar mesh in both m and q for a timestep in which mass is added to the star over the time interval t 1 to t 2. Vertical and slanted lines indicate cell boundaries. Cell size is exaggerated; there can be many cells in the newly added material. Three regions are chosen in the process of expanding the mesh for the added material, an inner Lagrangian region, an outer homologous region, and a transition region.  \n<!-- image -->  \nPlacement of the transition is related to the mesh. The MESAstar mesh structure is unchanged from that discussed in Paper I and Paper II. An illustration of the mesh regions and an indication of the behavior of cell boundaries during accretion is shown in Figure 32. In the case of mass loss, analogous operations are performed; here we focus only on mass gain.  \nTime derivatives appearing in the equations for physical evolution are then estimated using first order di GLYPH<11> erences. In the Lagrangian mesh region, the finite di GLYPH<11> erence form of ( @=@ t ) m involves a simple same-cell di GLYPH<11> erence. Similarly, in the homologous mesh region, the finite di GLYPH<11> erence form of ( @=@ t ) q involves a same-cell di GLYPH<11> erence. In most cases, by design, these same-cell di GLYPH<11> erences are for values whose changes, e.g. GLYPH<14> ln T , are directly available from the iterative solution of the new structure, allowing us to avoid the numerical problems inherent in subtracting two almost identical numbers. In the transition region both m and q coordinates of cells have been modified, so we cannot do a same-cell difference for either ( @=@ t ) m or ( @=@ t ) q . Instead, we interpolate values from the model at the start of the step to corresponding locations in m or q at the end of the step.  \nA mass GLYPH<14> M = ˙ M ( t 2 GLYPH<0> t 1) is added to the star from time t 1 to t 2. Sizes are exaggerated for clarity; there are generally many zones in each region, possibly hundreds in the newly added material. The diagram is shown in both mass coordinate, m , and homology coordinate, q . Before each timestep, MESAstar adjusts the initial mesh resolution by splitting or merging cells based on local gradient conditions, producing an adjusted resolution mesh, which we show at t 1. No simulation time elapses during that process. When constructing the mesh that will represent the star at t 2, MESAstar divides the new mesh into three regions: an inner Lagrangian region, in which the m boundaries of each cell are preserved during the timestep, a transition region, and an outer, homologous region, in which the q boundaries of each cell are the same across the timestep. The result of this operation is an expanded mass domain shown at t 2 in Figure 32.  \nA smooth and accurate value for GLYPH<15> grav in the transition region is important. To ensure this, the location of the transition region is selected to reduce the di GLYPH<11> erences between the constant m and constant q forms of the time derivative and maintain accurate finite di GLYPH<11> erencing. As a simple mechanism to control these, we limit, in units of cell size, the o GLYPH<11> set in the interpolation used to translate locations from the beginning to the end of the timestep (Miles et al. 2015). Using the cell size implicitly takes advantage of the limits imposed by mesh controls on the maximum possible magnitude of cell-to-cell changes in key variables, including the variables of interest for GLYPH<15> grav time derivatives.', '7.2. Testing': 'In order to demonstrate that the interface between the outer homologous region and the inner Lagrangian region provides a smooth profile that is independent of timestep size, we have repeated the test shown in Paper II section 5.3. This test involves accretion of solar composition material onto a WD at 10 GLYPH<0> 10 M GLYPH<12> yr GLYPH<0> 1 . We use the same starting model as in Paper II, which was produced by accreting hydrogen-rich material through several hydrogen shell flashes on a 0 : 6 M GLYPH<12> WD with an initial core temperature of about 10 7 K. As accretion proceeds, the total accumulated accreted mass, M acc, increases up to a maximum which causes the hydrogen flash and  \nnova runaway, M ign.  \nFigure 33. Profiles spanning the Lagrangian-homologous grid transition region in a 0 : 6 M GLYPH<12> WDwith a core temperature of 10 7 K undergoing accretion of solar composition material at a rate of 10 GLYPH<0> 10 M GLYPH<12> yr GLYPH<0> 1 . At the time shown M acc = M ign = 0 : 2. Shown for comparison are the results of the treatment discussed here for a 5000 yr timestep and a 1000 yr timestep and the treatment discussed in Paper II. For the current treatment the transition from homologous to Lagrangian is indicated by two open circles at either end of the region. The location GLYPH<14> M in from the surface, the base of material newly accreted this step, is indicated by the triangle.  \n<!-- image -->  \n-  \nProfiles near the transitions region are shown in Figure 33 at the time when M acc = M ign = 0 : 2, which had the most severe discontinuity in Paper II. The first panel shows r T and the second panel shows 4 GLYPH<25> r 2 GLYPH<26> HP GLYPH<2> GLYPH<15> grav, where HP is the pressure scale height. This is the amount of energy being released due to the T Ds = Dt term in the energy equation within a scale height, and has units of luminosity. We have chosen a timestep of 5 ; 000 yr, which places the homology-Lagrangian transition region in a similar place to the location of the discontinuity in Paper II.  \nThe orange curve in Figure 33 was computed using MESA version r4664, as used in Paper II. This displays the discontinuity due to using GLYPH<15> grav from Equation (62) and Equation (63) with only the homologous term and a transition point a factor of 5 deeper in pressure than GLYPH<14> M . The black curve shows the same simulation with the same timestep for the treatment discussed here, in which GLYPH<15> grav ; nh is included and the transition region, indicated with solid dots at either end, is placed as described in Section 7.1. We see that, away from the transitions, GLYPH<15> grav is unchanged from the values found in Paper II, in which only the homologous term was used in the exterior. We also show the result for a timestep of 1,000 yr is indistinguishable on this plot; the GLYPH<15> grav profile di GLYPH<11> ers by a fraction of a percent at the edges of the transition region, and less elsewhere. In the current treatment the profiles are independent of timestep  \nglyph[arrowbothv]  \nglyph[arrowbothv]  \nsize.', '8. WEAK REACTIONS': "An erratum for this paper was published as Paxton et al. (2016). The errors reported there have been corrected in this version of the paper.  \nThese grids include the thermodynamic conditions where the electrons are degenerate and relativistic, which are realized for example in massive white dwarfs and cores of intermediate mass stars. Under these conditions, the rates of electron-capture and beta-decay reactions are su GLYPH<14> ciently sensitive to density and temperature that they can change by tens of orders of magnitude between adjacent points in these tables. Linear or cubic interpolation cannot accurately reproduce the value of the rate between the tabulated points. glyph[arrowbothv]  \nThe rates module provides weak reaction rates for hundreds of isotopes. By default, when atoms are fully ionized, these rates are based (in order of precedence) on the tabulations of Langanke & Mart'ınez-Pinedo (2000), Oda et al. (1994), and Fuller et al. (1985). These tables span a wide range of density and temperature, 1 GLYPH<20> log( GLYPH<26> Y e = g cm GLYPH<0> 3 ) GLYPH<20> 11 and 7 GLYPH<20> log( T = K) . 10 : 5, but are relatively coarse, with 11 points in the GLYPH<26> Y e dimension ( GLYPH<1> log GLYPH<26> Y e = 1) and 12 points in the T dimension ( GLYPH<1> log T GLYPH<25> 0 : 25). glyph[arrowbothv]  \nThe di GLYPH<14> culty of interpolating in coarse rate tabulations was discussed by Fuller et al. (1985), who proposed a physically-motivated interpolation scheme, hereafter referred to as Fuller, Fowler, and Newman (FFN) interpolation. Their procedure assumes the rate has the form given by a single transition between the parent and daughter nuclear ground states. However, the true rate may be dominated by allowed transitions to or from excited states in the parent or daughter nucleus. This is almost always the case when the ground state to ground state transition is highly forbidden. The specific transition that dominates the rate may change over the range of thermodynamic conditions covered by the table. The FFN interpolation method does not account for these complications.  \nIn recent years, a number of authors have discussed the importance of well-sampled weak rates in capturing the influence of these processes on stellar evolution. This can be achieved by generating denser tables for the specific reactions of interest or by using analytic approximations to the rates (e.g., Toki et al. 2013; Mart'ınez-Pinedo et al. 2014). We now present a capability by which MESA can calculate weak reaction rates on-the-fly from input nuclear data. This removes the potential for interpolation artifacts. It also enables easy experimentation in cases where the input nuclear data may not be well-measured. We begin with an overview of how we calculate these weak rates and illustrate their utility and a few applications.  \nFigure 34 compares the results of the interpolation methods described in the preceding paragraphs with the on-thefly approach that we have implemented in MESA and will be described here. It shows the electron-capture rate on 24 Mg and beta-decay rate of 24 Na at fixed temperature. Linear interpolation of these coarse tables fails to reproduce the rapid variation in the rate. The FFN interpolation method produces curves with characteristic shapes more similar to the true rate, but because the Q -value is that of the ground state to ground state transition and not that of the transition that dominates the rate, the density dependence is not correct in detail.", '8.1. Calculation of Weak Rates': 'Figure 34. The top panel (bottom panel) shows the rate of electron capture on 24 Mg (beta decay of 24 Na) as a function of density at a fixed temperature of log( T = K) = 8 : 6. The Oda et al. (1994) tabulated points are shown as black dots. The dotted line shows the result of using linear interpolation between the tabulated points. The dashed line shows the result of using the physicallymotivated interpolation method suggested by Fuller et al. (1985). The solid line shows the rate calculated using the on-the-fly rate calculation capability of MESA documented in this section. Slight di GLYPH<11> erences between points and the line are due to di GLYPH<11> erences in the input nuclear data.  \n<!-- image -->  \nConsider two nuclei A GLYPH<17> ( Z ; N ) and B GLYPH<17> ( Z GLYPH<0> 1 ; N + 1) that have two states connected by an electron-capture transition  \nA + e GLYPH<0> ! B + GLYPH<23> e (66)  \nand beta-decay transition  \nB ! A + e GLYPH<0> + ¯ GLYPH<23> e : (67)  \nThe energy di GLYPH<11> erence between the ground states can be written as  \nQ g = ( ( MA GLYPH<0> MB ) c 2 for electron capture ; ( MB GLYPH<0> MA ) c 2 for beta decay : (68)  \nwhere MA and MB are the nuclear rest masses of the ground states. The total energy di GLYPH<11> erence between any two states can be written as  \nQij = Q g + Ei GLYPH<0> Ef ; (69)  \nwhere Ei and Ef are the energies of the initial and final states measured relative to the ground state. For the transitions that we consider here, Q g < 0 and Qij < 0 for electron capture and Q g > 0 and Qij > 0 for beta decay.  \nThe total rate of the process (electron capture or beta decay) is the sum of the individual transition rates from the i -th parent  \nIn this section, we use J to represent the nuclear spin. We work in the allowed approximation, which neglects all total lepton angular momentum ( L = 0). This restricts us to Fermi transitions, where the total lepton spin is S = 0, and therefore the initial and final nuclear spins are equal ( Ji = Jf ), and Gamow-Teller transitions, where S = 1, and therefore Ji = Jf ; Jf GLYPH<6> 1 (excluding Ji = Jf = 0). In both cases, there is no parity change: GLYPH<25> i GLYPH<25> f = + 1 (e.g., Commins 1973).  \nstate to the j -th daughter state, GLYPH<21> i j , weighted by the occupation probability of the i -th parent state, pi .  \nGLYPH<21> = X i pi X j ln 2 ( f t ) i j GLYPH<8> ( GLYPH<22> e ; T ; Qij ) ; (70)  \nwhere ( f t ) is the comparative half-life and can be either measured experimentally or calculated from theoretical weakinteraction nuclear matrix elements. The i -sum is over all parent states and the j -sum is over all daughter states. The occupation probability is  \npi = (2 Ji + 1) exp ( GLYPH<0> GLYPH<12> Ei ) P k (2 Jk + 1) exp( GLYPH<0> GLYPH<12> Ek ) ; (71)  \nwhere we define GLYPH<12> = ( k B T ) GLYPH<0> 1 . The quantity GLYPH<8> is a phase space factor which depends on the electron chemical potential GLYPH<22> e (including the electron rest mass), on the temperature T , and the energy di GLYPH<11> erence Qij . The value of GLYPH<8> for electron capture is  \nGLYPH<8> ec = exp( GLYPH<25>GLYPH<11> Z ) ( m e c 2 ) 5 Z 1 GLYPH<0> Qij E 2 e ( E e + Qij ) 2 1 + exp[ GLYPH<12> ( E e GLYPH<0> GLYPH<22> e)] dE e ; (72)  \nwhere GLYPH<11> is the fine structure constant. For beta decay it is  \nGLYPH<8> GLYPH<12> = exp( GLYPH<25>GLYPH<11> Z ) ( m e c 2 ) 5 Z Qij m e c 2 E 2 e ( E e GLYPH<0> Qij ) 2 1 + exp[ GLYPH<0> GLYPH<12> ( E e GLYPH<0> GLYPH<22> e)] dE e : (73)  \nSimilarly, the total rate of energy loss via neutrinos is  \n"GLYPH<23> = X i pi X j m e c 2 ln 2 ( f t ) i j GLYPH<9> ( GLYPH<22> e ; T ; Qij ) ; (74)  \nThe value of GLYPH<9> for electron capture is  \nGLYPH<9> ec = exp( GLYPH<25>GLYPH<11> Z ) ( m e c 2 ) 6 Z 1 GLYPH<0> Qij E 2 e ( E e + Qij ) 3 1 + exp[ GLYPH<12> ( E e GLYPH<0> GLYPH<22> e)] dE e ; (75)  \nand for beta decay it is  \nGLYPH<9> GLYPH<12> = exp( GLYPH<25>GLYPH<11> Z ) ( m e c 2 ) 6 Z Qij m e c 2 E 2 e ( E e GLYPH<0> Qij ) 3 1 + exp[ GLYPH<0> GLYPH<12> ( E e GLYPH<0> GLYPH<22> e)] dE e : (76)  \nIn order to implement these equations in MESA , we rewrite the integrals in terms of Fermi-Dirac integrals, following appendix A of Schwab et al. (2015). MESA implements fast quadrature routines to evaluate integrals of this form. Each time a weak rate is needed, it is calculated on-the-fly. We discuss the computational cost of this procedure in Section 8.3.  \nAssuming thermal equilibrium, the energy released by weak reactions depends only on total reaction rate, total neutrino loss rate, energy di GLYPH<11> erence between the nuclei, and the electron chemical potential. Therefore, the total specific heating rate from a reaction is  \nGLYPH<15> ec = nA GLYPH<26> h ( Q g + GLYPH<22> e) GLYPH<21> ec GLYPH<0> "GLYPH<23>; ec i ; (77)  \nGLYPH<15>GLYPH<12> = nB GLYPH<26> h ( Q g GLYPH<0> GLYPH<22> e) GLYPH<21>GLYPH<12> GLYPH<0> "GLYPH<23>;GLYPH<12> i ; (78)  \nwhere nA and nB are the number densities of the species undergoing electron capture and beta decay, respectively, and GLYPH<26> is the total mass density.  \nTherefore, given a list of nuclear levels and the ( f t )-values for the transitions between them, MESA can calculate the rates  \nFigure 35. The electron capture (solid lines) and beta decay (dashed lines) rates of the 23 Na23 Ne Urca pair as calculated by MESA , using the on-the-fly methods described in this section. The value of log( T = K) is shown next to each electron capture line; the beta decay line of matching color is at the same temperature. The rates vary rapidly, with both temperature and density, near the threshold density, which is roughly in the center of the plot.  \n<!-- image -->  \nof electron capture and beta decay and the corresponding energy generation rates. Typically only a few low-lying states and the transitions between them are needed. As an example, Figure 35 shows the rates for the 23 Na23 Ne Urca pair.', '8.1.1. Coulomb Corrections': 'In a dense plasma, the electrostatic interactions of the ions and electrons introduce corrections to the weak rates relative to those which assume a Fermi gas of electrons and an ideal gas of ions. Our treatment of these e GLYPH<11> ects, which is presented in appendix B of Schwab et al. (2015), is similar to appendix A of Juodagalvis et al. (2010).  \nThe electron density relevant to the reaction rate is not the average electron density, but rather the electron density at the position of the nucleus. This correction is accounted for as a shift in the value of the electron chemical potential that enters the phase space factor, GLYPH<22> 0 e = GLYPH<22> e + Vs . Values of Vs have been calculated by Itoh et al. (2002). This correction does not enter the energy generation rates because it has not changed the energy cost to add or remove an electron.  \nSince electron capture and beta decay change the ion charge, the Coulomb interaction energy changes the energy di GLYPH<11> erence between the parent and daughter nuclear states. To calculate this shift, we use the excess ion chemical potential GLYPH<22> ex from Potekhin et al. (2009). We incorporate this e GLYPH<11> ect by shifting the value of Qij , as defined in Equation (69), by an amount GLYPH<1> E = GLYPH<22> ex ; parent GLYPH<0> GLYPH<22> ex ; daughter. This shift, Q 0 i j = Qij +GLYPH<1> E , enters the calculation of the phase space factors and the energy generation rates.', '8.2. Applications': 'When GLYPH<22> e . j Q j , only the few electrons in the tail of the Fermi-Dirac distribution have su GLYPH<14> cient energy to overcome the energy gap and capture on A to form B . Thus, the rate of electron capture is small compared to the rate of beta decay, and so isotope B is favored in the equilibrium. When GLYPH<22> e & j Q j , there are only a few unoccupied states available to accept the energetic electron from the beta decay. This final state blocking means the rate of beta decay is small compared to  \nthe rate of electron capture, and so isotope A will be favored in the equilibrium.  \nThe shift in this equilibrium can have profound consequences when it occurs in stellar interiors. It modifies the composition, reduces the electron fraction, and alters the thermal state of the plasma. We now discuss two applications of our on-the-fly treatment of the weak rates: the Urca process and accretion-induced collapse.', '8.2.1. Urca Process': 'When the ground state to ground state transition is allowed (odd nuclei), the rates of electron capture and beta decay are both significant when GLYPH<22> e GLYPH<25> j Q g j . Since each reaction produces a neutrino which free-streams out of the star, this can lead to significant cooling. With a total number density of an Urca species nU = nA + nB , assuming the abundances are given by the detailed balance condition nA GLYPH<21> ec + nB GLYPH<21>GLYPH<12> = 0, the volumetric neutrino cooling rate will be nUC , where  \nC = "GLYPH<23>; ec GLYPH<21>GLYPH<12> + "GLYPH<23>;GLYPH<12> GLYPH<21> ec GLYPH<21>GLYPH<12> + GLYPH<21> ec : (79)  \nIn the limit k B T GLYPH<28> j Q j , the maximum value of the Urca cooling rate at a given temperature has a simple form (e.g., Tsuruta &Cameron 1970)  \nC max = 7 GLYPH<25> 4 ln 2 60 m e c 2 ( f t ) GLYPH<12> + ( f t )ec ! k B T m e c 2 ! 4 Q m e c 2 ! 2 exp( GLYPH<25>GLYPH<11> Z ) : (80)  \nThus, when the Urca process is important, well-resolved weak rates are necessary to correctly capture the temperature evolution of the core (Toki et al. 2013; Jones et al. 2013). Jones et al. (2013) used MESA r3709 along with a denser table described in Toki et al. (2013) to do their work. The Toki et al. (2013) table is not publicly available, so to reproduce the results of Jones et al. (2014) we save a model of an 8.8 M GLYPH<12> star at log( GLYPH<26> c = g cm GLYPH<0> 3 ) = 8 : 95 from our run with MESA (version r3709) using the Jones et al. (2014) inlists. We then load this model into a newer MESA version (r7503) that has access to the on-the-fly weak rates and evolve this model using a network with only the Urca process reactions. During this phase other nuclear reactions are not important to the central evolution.  \nWell-sampled rates such as those shown in Figure 35 are necessary to reproduce the correct Urca cooling rates. We illustrate this in Figure 36, which shows C max for the 23 Na23 Ne Urca process for temperatures 10 8 - 10 9 K. The circles show the results using the on-the-fly treatment described in this paper; the squares show the results using the coarse tables of Oda et al. (1994). The dashed line shows the cooling rate expected from Equation (80) which is in excellent agreement with the results of the on-the-fly method. The Urca cooling rates calculated from interpolating in coarse tables severely underestimate the true cooling rate when k B T GLYPH<28> j Q j .  \nFigure 37 shows the central temperature and density of the core. The solid lines show the evolution using the on-the-fly rates, the dashed lines show the results when interpolating in coarse tables. The drops in temperature at log( GLYPH<26> c = g cm GLYPH<0> 3 ) GLYPH<25> 9 : 1 and GLYPH<25> 9 : 25 correspond to cooling from the 25 Mg25 Na and 23 Na23 Ne Urca pairs, respectively. The corresponding shifts in composition can be clearly seen in the lower panel. These results demonstrate the importance of densely-sampled weak rates to the evolution of the core.  \nFigure 36. The e GLYPH<11> ect of the interpolation method on the Urca process cooling rates. The circles show the maximum value of C (Equation 79) calculated using the on-the-fly methods discussed in this section; the squares show the results using the coarse tables of Oda et al. (1994). Interpolation in these coarse tables severely underestimates the Urca cooling rates at low temperatures. The dashed line shows the expected value of the cooling rate given by Equation (80).  \n<!-- image -->  \nFigure 37. The top panel shows the evolution of T c and GLYPH<26> c in an 8.8 M GLYPH<12> star. The bottom panel shows the central 25 Mg and 23 Na mass fractions. The solid lines show the evolution using the on-the-fly rates, the dashed lines show the results when interpolating in coarse tables. The locations of the changes in mass fraction match the locations of cooling in the top panel. This demonstrates the importance of densely-sampled weak rates to the evolution of the core.  \n<!-- image -->  \nWhen the ground state to ground state transition is forbidden (even nuclei), the first transition to become significant is typically an allowed transition into an excited state. In these cases, the beta decays from the daughter ground state are blocked and decays from daughter excited states are strongly suppressed by the Boltzmann factor. Therefore significant cooling via the Urca process does not occur. Instead, since the captures are preferentially to an excited state, significant heating occurs via gamma-ray emission as level populations relax to a thermal distribution.  \nTwo important capture chains occur in oxygen-neonmagnesium (ONeMg) cores: 24 Mg ! 24 Na ! 24 Ne and 20 Ne ! 20 F ! 20 O. For these sequences of captures, the excess electron energy is thermalized. These are the key reactions in electron-capture supernovae and the accretion-  \nFigure 38. The evolution of a cold ONeMg WD towards AIC. The top panel shows the evolution of GLYPH<26> c and T c. The bottom shows the central 24 Mg and 20 Ne mass fractions. The solid lines show the evolution using the on-the-fly rates; the dashed lines show the results when interpolating in coarse tables.  \n<!-- image -->  \ninduced collapse (AIC) of ONeMg white dwarfs (e.g., Miyaji et al. 1980). As the degenerate core approaches the Chandrasekhar mass, the electron captures remove the pressure support and heat the plasma. Figure 38 shows the evolution of a cold ONeMg WD ( X O = 0 : 5, X Ne = 0 : 45, X Mg = 0 : 05) accreting at ˙ M = 10 GLYPH<0> 6 M GLYPH<12> yr GLYPH<0> 1 . The solid lines show the evolution using the on-the-fly rates described in this section; the dashed lines show the results when interpolating in coarse tables. When using the coarse tables, the electron captures on 24 Mg do not occur until approximately a factor of two larger density. At this greater density, the energy deposition from each capture is higher and this leads to a large temperature change due to the A = 24 captures alone. In contrast, the on-the-fly rates show the behavior demonstrated in previous studies of this evolution that did not use sparse tables (e.g., Miyaji & Nomoto 1987): the A = 24 captures heat the plasma and accelerate the contraction; the A = 20 captures, due to the higher 20 Ne abundance and a higher energy release per capture, cause a thermal runaway and the formation of an oxygen deflagration (Schwab et al. 2015).', '8.3. Guidelines': "MESA provides the nuclear data used in the calculation of the reactions specifically discussed in this section (Tilley et al. 1998; Firestone 2007a,b, 2009; Shamsuzzoha Basunia 2011; Mart'ınez-Pinedo et al. 2014). To consider additional reactions, a list of nuclear levels and ( f t )-values must be specified.  \nThe expressions in Section 8.1 assume degenerate, relativistic electrons. As GLYPH<22> e increases, additional transitions to higher energy states of the daughter nuclei and must be included. At higher temperatures, excited states of the parent nucleus will begin to be thermally populated and captures or decays from those states and must be included. At temperatures and densities where the composition approaches NSE, these methods are particularly inappropriate, as it is necessary to consider large pools of isotopes (Juodagalvis et al. 2010).", '9. CHEMICAL DIFFUSION': "MESA 's early implementation of microscopic element di GLYPH<11> usion incorporated the approach used by Thoul et al. (1994) in their seminal work on understanding the sedimentation of helium in the solar interior. The fundamental starting point for  \nthis treatment of di GLYPH<11> usion is the Boltzmann equation with the assumption of binary collisions where the particle's mean free path is much larger than the average particle spacing. This formalism, encoded in the Burgers equations (Burgers 1969), assumes that ions interact with an e GLYPH<11> ective potential that governs isolated interactions between only two particles at a time. For more strongly coupled plasmas, as GLYPH<0> GLYPH<25> e 2 = ( GLYPH<21> ion k B T ) exceeds unity (where GLYPH<21> ion = (3 = 4 GLYPH<25> n ion) 1 = 3 is the mean inter-ion spacing, and n ion is the total ion number density), it is no longer clear that this assumption remains valid. Later updates to MESA incorporated the work of Hu et al. (2011) on radiative levitation and incorporated the resistance coe GLYPH<14> cients calculated by Paquette et al. (1986) for approaches to the denser plasma regime as GLYPH<0> ! 1.  \nHere we describe MESA 's current implementation of chemical di GLYPH<11> usion and then discuss the path forward for di GLYPH<11> usion implementations in the GLYPH<0> > 1 regime, needed for accurate studies of di GLYPH<11> usion in the interiors of white dwarfs or surfaces of neutron stars. Recent theoretical work in this strongly coupled regime (Baalrud & Daligault 2013, 2014; Beznogov & Yakovlev 2014) provides support for a future update of MESA .", '9.1. Current Methods in MESAstar': 'We now describe the formalism and assumptions underlying the approach to di GLYPH<11> usion currently present in MESA . This is followed by a discussion of the framework for numerical implementation of this formalism provided by Thoul et al. (1994) and key modifications present in the current version of the MESA di GLYPH<11> usion routine.', '9.1.1. Burgers Equations and the Low Density Limit': 'The Burgers equations for di GLYPH<11> usion in an ionized plasma are derived using the Boltzmann equation for the distribution function Fs ( x ; GLYPH<24> ; t ) for particles of type s  \n@ Fs @ t + X i GLYPH<24> i @ Fs @ xi + X i fsi ms @ Fs @GLYPH<24> i = dFs dt ! collision ; (81)  \nwhere xi are the components of the position vector, GLYPH<24> i are the components of the velocity vector, fsi are components of the forces on particles of type s , and ms is the mass for those particles. Throughout this section, the indices s and t refer to particle species, while i and j are used to index other quantities such as spatial components of vectors.  \nBurgers adopts the 13-moment approximation due to Grad (1949) as a closure scheme for taking moments of the Boltzmann equation. Burgers also assumes an approximately Maxwellian distribution function  \nFs = ns GLYPH<25> 3 = 2 a 3 s exp GLYPH<0> c 2 s a 2 s ! (1 + GLYPH<30> s ) ; (82)  \nwhere as = (2 k B T = ms ) 1 = 2 , csi = GLYPH<24> i GLYPH<0> usi represents the components of the deviation of the velocity from the mean flow velocity u s of the species, and  \nGLYPH<30> s = X i ; j Bsijcsicsj + X i Csi c 2 s GLYPH<0> 5 2 a 2 s ! csi (83)  \nis the small deviation ( GLYPH<30> s GLYPH<28> 1) from the Maxwellian distribution. The coe GLYPH<14> cients Bsij and Csi are defined such that the distribution function has a total of 13 free parameters corresponding to the 13 moments of the closure scheme (see Burgers 1969).  \nBurgers derives the collision integrals ( S ( l ) st ) and crosssections ( GLYPH<6> ( l j ) st ) that result from taking moments of the right hand side of the Boltzmann equation  \nS ( l ) st = 2 GLYPH<25> Z 1 0 (1 GLYPH<0> cos l GLYPH<31> st ) b db ; (84)  \nGLYPH<6> ( l j ) st = 4 GLYPH<25> GLYPH<25> 3 = 2 Z 1 0 dv exp GLYPH<0> v 2 GLYPH<11> 2 st ! v 2 j + 3 GLYPH<11> 2 j + 4 st S ( l ) st ; (85)  \nwhere GLYPH<11> 2 st = 2 k B T =GLYPH<22> st , GLYPH<22> st = msmt = ( ms + mt ), v represents the relative velocity of colliding particles, and the angle of deviation GLYPH<31> st is a function of both v and the impact parameter b that depends on the physics of the two-particle interaction between colliding particles in the gas. Burgers then defines the dimensionless coe GLYPH<14> cients zst , z 0 st , z 00 st , and z 000 st , along with resistance coe GLYPH<14> cients ( Kst ) in terms of the collision integrals:  \nKst = Kts = 2 3 nsnt GLYPH<22> st GLYPH<11> st GLYPH<6> (11) st ; GLYPH<6> (12) st = GLYPH<6> (11) st = 5 2 (1 GLYPH<0> zst ) ; GLYPH<6> (13) st = GLYPH<6> (11) st = 25 4 GLYPH<0> 25 2 zst + 5 2 z 0 st ; GLYPH<6> (22) st = GLYPH<6> (11) st = z 00 st ; GLYPH<6> (23) st = GLYPH<6> (11) st = z 000 st : (86)  \nIn the \'single-fluid picture\' the di GLYPH<11> usion velocities are defined with reference to the mean velocity of the gas as a whole ( u ), rather than with respect to the mean species velocity ( u s ):  \nusi = 1 ns Z d GLYPH<24> GLYPH<24> i Fs ; u = 1 GLYPH<26> X s GLYPH<26> s u s ; w s = u s GLYPH<0> u : (87)  \nBurgers defines residual heat flow vectors  \nrsi = ms 2 nsk B T Z d 3 GLYPH<24> ( GLYPH<24> i GLYPH<0> ui ) j GLYPH<24> GLYPH<0> u j 2 Fs ! GLYPH<0> 5 2 wsi : (88)  \nAs shown in section 18 of Burgers (1969) if we assume j w s j GLYPH<28> as and the absence of magnetic fields, the basic equations of di GLYPH<11> usion are  \nr ps GLYPH<0> GLYPH<26> s g GLYPH<0> GLYPH<26> es E = X t , s Kst ( w t GLYPH<0> w s ) + X t , s Kstzst mt r s GLYPH<0> ms r t ms + mt ; (89)  \n5 2 nsk B r T = GLYPH<0> 2 5 Kssz 00 ss r s GLYPH<0> 5 2 X t , s Kstzst mt ms + mt ( w t GLYPH<0> w s ) GLYPH<0> X t , s Kst " 3 m 2 s + m 2 t z 0 st ( ms + mt ) 2 + 4 5 msmt ( ms + mt ) 2 z 00 st # r s + X t , s Kst msmt ( ms + mt ) 2 3 + z 0 st GLYPH<0> 4 5 z 00 st ! r t ; (90)  \nwhere E is the quasi-static electric field and GLYPH<26> es is the average charge density of species s . These equations are still general, with the form of the resistance coe GLYPH<14> cients not yet fully specified. The physics of the particular types of interactions within ideal gases is fully contained in the coe GLYPH<14> cients Kst , zst , z 0 st , z 00 st , and z 000 st .  \nFor ionized gases, the resistance coe GLYPH<14> cients require evaluation of collision integrals that diverge for a pure Coulomb  \npotential. However, since the two-particle interaction potential is only truly applicable on short length scales, an integration cuto GLYPH<11> or screened potential is commonly adopted. Burgers chooses to calculate resistance coe GLYPH<14> cients using a pure Coulomb potential truncated at the Debye radius  \nR D = 0 B B B B B @ 4 GLYPH<25> X s nsZ 2 s e 2 k B T 1 C C C C C A GLYPH<0> 1 = 2 ; (91)  \nwhich is assumed to be much larger than the inter-ion spacing. Indeed, for a plasma of one species, R D =GLYPH<21> ion = (3 GLYPH<0> ) GLYPH<0> 1 = 2 . Applying this form of interaction to the collision integrals, the l = 1 integrals defined in Equation (84) can be evaluated (Baalrud & Daligault 2014)  \nS (1) st = 2 GLYPH<25> R 2 D GLYPH<11> 4 st GLYPH<3> 2 st v 4 ln 2 6 6 6 6 6 4 1 + GLYPH<3> 2 st v GLYPH<11> st ! 4 3 7 7 7 7 7 5 ; (92)  \nwhere GLYPH<3> st = GLYPH<22> st GLYPH<11> 2 st R D = ( ZsZte 2 ). In order to perform the integral in Equation (85), Burgers notes that the dependence of S ( l ) st on v inside the logarithmic term is weak, so that we can replace v 2 there with its average value h v 2 i = 3 k B T =GLYPH<22> st . Assuming a very dilute plasma, so that GLYPH<3> 2 st h v 2 i 2 =GLYPH<11> 4 st GLYPH<29> 1, Burgers then writes  \nS (1) st GLYPH<25> 4 GLYPH<25> R 2 D GLYPH<11> 4 st GLYPH<3> 2 st v 4 ln 3 k B TR D ZsZte 2 ! ; (93)  \nand the final result for the resistance coe GLYPH<14> cients follows as  \nKst GLYPH<25> 16 p GLYPH<25> 3 nsntZ 2 s Z 2 t e 4 GLYPH<22> st GLYPH<11> 3 st ln 3 k B TR D ZsZte 2 ! ; (94)  \nzst = 3 5 ; z 0 st = 13 10 ; z 00 st = 2 ; z 000 st = 4 : (95)  \nWith these coe GLYPH<14> cients now fully specified, Burgers di GLYPH<11> usion equations along with constraints such as charge neutrality and current neutrality form a closed set of equations, which can be solved for w s , r s , E , and g from the input of a stellar profile.', "9.1.2. MESA 's Implementation of Thoul et al.'s Approach": "The di GLYPH<11> usion routine originally implemented in MESA was based on the work of Thoul et al. (1994). They start with the Burgers equations, written in a compact notation following Noerdlinger (1977, 1978) that is equivalent to Equations (89) and (90) in one dimension. However, the approach of Thoul et al. (1994) di GLYPH<11> ers from Burgers' original treatment in one important respect: the resistance coe GLYPH<14> cients are based on a modified result for the collision integrals. They follow Equation (95) for the various zst coe GLYPH<14> cients, which uses a pure Coulomb potential with a cuto GLYPH<11> at the Debye length, but the Kst coe GLYPH<14> cients were derived from an alternative fitting of the Coulomb logarithms introduced by Iben & MacDonald (1985). For these coe GLYPH<14> cients, they define GLYPH<21> = max( R D ; GLYPH<21> ion), and use  \nKst = 16 p GLYPH<25> 3 nsntZ 2 s Z 2 t e 2 GLYPH<22> st GLYPH<11> 3 st GLYPH<2> 1 : 6249 2 ln 2 6 6 6 6 6 4 1 + 0 : 18769 4 k B T GLYPH<21> ZsZte 2 ! 1 : 2 3 7 7 7 7 7 5 : (96)  \nThis expression is a fit to the numerical results of Fontaine & Michaud (1979), motivated by white dwarf conditions where  \nBurgers' approximations for dealing with Equation (92) are not valid ( GLYPH<0> > 1). Since this fit focuses on the strong coupling regime, and di GLYPH<11> ers from Equation (94), these results can be incorrect in the limit of a dilute plasma as we discuss later. Nevertheless, Thoul et al. (1994) elected to use Equation (96) under all conditions, since it provides an approximately correct solution in a convenient closed form.  \nUsing Equations (89) and (90) along with the constraints of current neutrality ( P s GLYPH<26> esws = 0) and local mass conservation ( P s GLYPH<26> sws = 0), Thoul et al. (1994) express an entire closed system of equations in a dimensionless matrix form suitable for numerical evaluation:  \np K 0 GLYPH<11> i d ln p dr + GLYPH<23> i d ln T dr + S X j = 1 j , e GLYPH<13> i j d ln Cj dr ! = 2 S + 2 X j = 1 GLYPH<1> i jWj ; (97)  \nwhere S is the total number of species in the gas (including electrons) and Cj = nj = n e is the concentration of the j th species. Consult Thoul et al. (1994) for definitions of K 0, GLYPH<11> i , GLYPH<23> i , GLYPH<13> i j , and GLYPH<1> i j . The definition of Wj is  \nWj = 8 > > > > > > < > > > > > > : wj for j = 1 : : : S ; rj for j = S + 1 : : : 2 S ; K GLYPH<0> 1 0 n e eE for j = 2 S + 1 ; K GLYPH<0> 1 0 n e m p g for j = 2 S + 2 : (98)  \nThis is the vector containing the unknown quantities solved for after specifying K 0, GLYPH<11> i , GLYPH<23> i , GLYPH<13> i j , and GLYPH<1> i j . The routine provided by Thoul et al. (1994) inverts Equation (97) for one term in the left hand side at a time so as to find the 'generalized di GLYPH<11> usion coe GLYPH<14> cients,' which can be used to construct di GLYPH<11> usion velocities or contributions from pressure, temperature, or concentrations individually.", '9.1.3. Modified Coe GLYPH<14> cients and Radiative Levitation as Implemented by Hu et al.': "Hu et al. (2011) extend the methods of Thoul et al. (1994) by introducing some key modifications. First, they include an extra force term due to radiative levitation, so that Equation (89) becomes  \ndps dr + GLYPH<26> s ( g GLYPH<0> g rad ; s ) GLYPH<0> ns ¯ ZseE = X t , s Kst ( wt GLYPH<0> ws ) + X t , s Kstzst mtrs GLYPH<0> msrt ms + mt ; (99)  \nwhere g rad ; s refers to the radiative acceleration on species s . ¯ Zs is the average charge of species s , allowing an account of partial ionization so that ns ¯ Zse = GLYPH<26> es . They do not modify Equation (90). 15  \nVst ( r ) = ¯ Zs ¯ Zte 2 exp( GLYPH<0> r =GLYPH<21> ) r ; (100)  \nIn contrast to Thoul's original routine, Hu et al. (2011) use the resistance coe GLYPH<14> cients from Paquette et al. (1986), which were generated based on substantial improvements to Fontaine & Michaud (1979). In evaluating the collision integrals, Paquette et al. (1986) use a screened Coulomb potential of the form  \n15 As written in equation (3) of Hu et al. (2011), their expression has two errors in the first term on the right hand side of the first line: the sign is wrong, and it is missing resistance coe GLYPH<14> cients Kij . Since neither of these errors propagates into later sections of the paper, it appears that both are simply typos, and otherwise their expression matches Equation (90) exactly.  \nwhere, once again, GLYPH<21> = max( R D ; GLYPH<21> ion). As we note below, this choice of GLYPH<21> makes a substantial di GLYPH<11> erence in strongly coupled plasmas, where the Debye radius no longer corresponds to a distance at which other nearby charged particles can significantly screen the Coulomb field. After setting up the algebra for a matrix solution very similar to that of Thoul et al. (1994), Hu et al. (2011) solve for the vector Wj (as defined in Equation 98) appearing in the equation  \np K 0 GLYPH<0> GLYPH<11> imig rad ; i k B T + GLYPH<11> i d ln p dr + GLYPH<23> i d ln T dr + S X j = 1 j , e GLYPH<13> i j d ln Cj dr ! = 2 S + 2 X j = 1 GLYPH<1> i jWj : (101)  \nMany of the quantities appearing in this equation are defined di GLYPH<11> erently than in Thoul et al. (1994); see Hu et al. (2011) for details. We can also solve this equation directly for the vector Wj to obtain  \nW 2 S + 1 W 2 S + 2 = K GLYPH<0> 1 0 n e eE K GLYPH<0> 1 0 n e m p g = eE m p g ; (102)  \nthe strength of the electric field relative to gravity.", '9.2. Analytic Expression for the Electric Field': "In some simple cases, Burgers equations can be solved to yield an analytic expression for the electric field, providing a useful test for MESA . Starting directly with his di GLYPH<11> usion equations, Burgers (1969) arrives at the following expressions for a pure plasma of electrons along with one species of ions (charge Ze ):  \nr p e + n e e E ZK 0 = w + 3 5 r e ; (103)  \n5 2 n e k B r T ZK 0 = GLYPH<0> 3 2 w GLYPH<0> 2 5 K ee z 00 ee K ie + z 0 ie ! r e ; (104)  \nwhere w = w i GLYPH<0> w e. For a plasma with only one ion species in di GLYPH<11> usion equilibrium, the constraints of current neutrality and local mass conservation give w = 0. In the case of a pure hydrogen plasma, p = 2 p e, and in hydrostatic equilibrium r p e = r p = 2 = GLYPH<26> g = 2. Hence, we can solve the above set of equations to find  \ne E = GLYPH<0> 1 2 m p g GLYPH<0> 3 2 2 5 K ee z 00 ee K ie + z 0 ie ! GLYPH<0> 1 k B r T : (105)  \nThe coe GLYPH<14> cient for the temperature gradient term depends directly on the nature of the resistance coe GLYPH<14> cients in the Burgers formalism, so di GLYPH<11> erent models of Coulomb collisions in ionized plasma will lead to di GLYPH<11> erent results for the electric field.  \nAs a slight generalization of Equation (105) in one dimension, we write  \neE m p g = 1 2 GLYPH<0> GLYPH<11> e k B m p g dT dr : (106)  \nIf we calculate the coe GLYPH<14> cient GLYPH<11> e using the Burgers' formalism with Equations (95) and (96), we find  \nGLYPH<11> e = 3 2 2 5 K ee z 00 ee K ie + z 0 ie ! GLYPH<0> 1 = 0 : 804 (107)  \nA comparable analytic expression for the electric field is provided by Roussel-Dupr'e (1981), who applies a BoltzmannFokker-Planck approach to finding di GLYPH<11> usion coe GLYPH<14> cients for trace elements in hydrogen plasma. His treatment of di GLYPH<11> usion is more precise than the Burgers' formalism, but has the limitation of only being applicable in the case of nearly pure hydrogen with a di GLYPH<11> using trace element. His result for the electric field matches the form of Equation (106) with the coe GLYPH<14> cient GLYPH<11> e = 0 : 703. This provides another useful point of comparison in the specific case of nearly pure hydrogen plasmas. Below we use this analytic expression as a test of the updated resistance coe GLYPH<14> cients employed by Hu et al. (2011).", '9.3. Results and Comparisons': 'We have constructed several simple MESA test cases in order to illustrate the e GLYPH<11> ects of radiative levitation and di GLYPH<11> erent resistance coe GLYPH<14> cients. Where possible, we compare MESA output to corresponding analytic expressions.', '9.3.1. Electric Fields': "By default, MESA uses the resistance coe GLYPH<14> cients provided by Paquette et al. (1986), but it can also use the resistance coe GLYPH<14> cients defined by Iben & MacDonald (1985), given here in Equation (96). In the case of a pure hydrogen star, the coe GLYPH<14> -cients given in Equation (96) lead directly to Equation (107), so these coe GLYPH<14> cients are especially useful in performing simple comparisons of MESA output to a corresponding analytic expression. Due to the complicated numerical methods used to obtain the resistance coe GLYPH<14> cients of Paquette et al. (1986), it is not possible to write down a directly corresponding closed form analytic expression for the electric field, but results based on these more precise calculations compare favorably to those of Roussel-Dupr'e (1981) in the case of a pure hydrogen plasma. Starting with the MESA test suite, we constructed a solar mass pure hydrogen star, and we ran just long enough to turn on the di GLYPH<11> usion routine and gather output for electric and gravitational fields. For such a star, we can compare MESA results for the electric field directly to the analytic expression given in Equation (106), with GLYPH<11> e = 0 : 804 in the solution of Burgers (1969) and GLYPH<11> e = 0 : 703 for Roussel-Dupr'e (1981).  \nThe Sun provides another interesting test case for comparing the e GLYPH<11> ects of using di GLYPH<11> erent resistance coe GLYPH<14> cients. An example solar model from the MESA test suite was run with different choices of the resistance coe GLYPH<14> cients. Figure 40 shows a slight di GLYPH<11> erence between the electric field strengths relative to gravity given by the Paquette et al. (1986) coe GLYPH<14> cients and those by Iben & MacDonald (1985).  \nFigure 39 plots the result of Equation (106) for both values of GLYPH<11> e , along with the results from the di GLYPH<11> usion routine (Equation 102) for each type of resistance coe GLYPH<14> cients available in MESA . As expected, the curve calculated from the MESA diffusion routine output using the resistance coe GLYPH<14> cients of Iben &MacDonald (1985) closely matches the analytic expression with GLYPH<11> e = 0 : 804 as calculated by Burgers (1969) using his similar coe GLYPH<14> cients. When using the more detailed numerical calculations for the resistance coe GLYPH<14> cients provided by Paquette et al. (1986), the di GLYPH<11> usion routine output closely resembles the more precise analytic calculation given by Roussel-Dupr'e (1981).", '9.3.2. Gravitational Fields': 'The MESA di GLYPH<11> usion routine treats both the electric field and local gravitational acceleration as unknown quantities. MESA  \nFigure 40. Comparison of electric field strengths relative to gravity using di GLYPH<11> erent resistance coe GLYPH<14> cients in a solar model.  \n<!-- image -->  \n/circledot  \nFigure 39. Comparison of electric field strengths relative to gravity in a pure hydrogen star ( M = 1 : 0 M GLYPH<12> , T e GLYPH<11> = 5 : 74 GLYPH<2> 10 3 K, L = 0 : 576 L GLYPH<12> ) with nuclear burning artificially suppressed in the MESA routine to avoid any helium contamination. Solid lines represent the analytic expression given by Equation (106) for two di GLYPH<11> erent values of the coe GLYPH<14> cient GLYPH<11> e . Dashed lines represent output from the MESA di GLYPH<11> usion routine as described in Equation (102), with the only di GLYPH<11> erence being the resistance coe GLYPH<14> cients used to solve the Burgers equations.  \n<!-- image -->  \n/circledot  \nrecords the quantity W 2 S + 2 (Equation 98), used to calculate the gravitational acceleration from the di GLYPH<11> usion routine:  \ng di GLYPH<11> = K 0 W 2 S + 2 n e m p : (108)  \nThis expression for g di GLYPH<11> is independent of the simpler expression for local gravitational acceleration g gauss = Gm = r 2 . Figure 41 compares g gauss and g di GLYPH<11> for a typical profile found using the example solar model from the MESA test suite. In Figure 41 a profile from a star of larger mass ( M = 1 : 5 M GLYPH<12> ) shows disagreement between the gravity outputs in the convective core. The Burgers formalism assumes heat transfer that is correlated with temperature gradients through the residual heat flow vectors defined in Equation (88). This assumption breaks down when most of energy is transported by convection; however, the e GLYPH<11> ects of di GLYPH<11> usion in this region are  \ncompletely overwhelmed by convective mixing and are therefore inconsequential.  \nFigure 43. Abundance profile of a CO WD ( M = 0.611 M GLYPH<12> , T e GLYPH<11> = 5.16 GLYPH<2> 10 3 K, L = 9 : 29 GLYPH<2> 10 GLYPH<0> 5 L GLYPH<12> ) after 4 Gyr of WD evolution. The region left of the blue, dashed line is the interior of the WD, where GLYPH<0> GLYPH<21> 1. Left of the red, dashed line GLYPH<0> GLYPH<21> 50, and di GLYPH<11> usion has been turned o GLYPH<11> for this region. The electrons are an ideal gas to the right of the black dot-dash line.  \n<!-- image -->  \n/circledot  \nFigure 41. Comparison of gravitational fields obtained from g di GLYPH<11> and g gauss in two MESA test suite cases. The two lines representing the Sun (age = 4 : 57 Gyr) show good agreement, while the two lines representing a 1 : 5 M GLYPH<12> star disagree in regions with large convective flux where di GLYPH<11> usion is inconsequential.', '9.3.3. Radiative Levitation': "MESA 's implementation of radiative levitation is based on Hu et al. (2011). Figure 42 shows an abundance profile of a subdwarf B star model produced by MESA , where radiative levitation is responsible for the presence of 56 Fe, 58 Ni, and other metals near the surface (as also seen in figure 3 of Hu et al. 2011).", '9.3.4. White Dwarf Sedimentation': "In a cooling WD, di GLYPH<11> usion governs sedimentation over long timescales. The assumptions behind the formalism of the Burgers equations do not hold under white dwarf conditions:  \n- GLYPH<15> The Burgers equations assume all particle species satisfy an ideal gas equation of state. In the context of a degenerate WD both electrons and ions violate this assumption.\n- GLYPH<15> The very dense, strongly coupled ( GLYPH<0> > 1) conditions of a WD call into question the validity of the two-particle scattering picture used to calculate the ion resistance coe GLYPH<14> cients.  \nNevertheless, for lack of a better option, previous studies have relied on the Burgers equations with the coe GLYPH<14> cients of Paquette et al. (1986). For example, see C'orsico et al. (2002).  \nFigure 43 shows an abundance profile produced by MESA for a CO WD after 4 Gyr of evolution, where di GLYPH<11> usion governs sedimentation in the outer layers. The vertical lines in Figure 43 mark the outer boundaries of regions where the two concerns listed above become significant. Nearly all of the WD resides inside at least one of these regimes, and much  \n<!-- image -->  \n(  \n/circledot  \n)  \nFigure 42. Abundance profile of a subdwarf B star model ( M = 0 : 462 M GLYPH<12> , T e GLYPH<11> = 2 : 67 GLYPH<2> 10 4 K, L = 1 : 12 L GLYPH<12> , age = 5 Myr) showing the e GLYPH<11> ects of radiative levitation with a layer of 56 Fe = 58 Ni at the surface.  \nof the interesting di GLYPH<11> usion sedimentation occurs inside regions that are both significantly coupled and highly degenerate. Thus, improvements to the treatment of di GLYPH<11> usion are clearly necessary before we are able to describe di GLYPH<11> usion in WDsadequately. This MESA run turns o GLYPH<11> di GLYPH<11> usion for GLYPH<0> GLYPH<21> 50, where we expect strong coupling to substantially modify the underlying equations.  \n<!-- image -->  \n(  \n/circledot  \n)", '9.4. Expanding the Domain of Validity and Next Steps': "The validity of the Boltzmann approach becomes questionable as GLYPH<0> > 1 and the ions become a liquid. Bildsten & Hall (2001) estimated the di GLYPH<11> usion coe GLYPH<14> cient in this liquid regime by using the Stokes-Einstein relation. However, for a broadbased code such as MESA , we need to implement di GLYPH<11> usion into the GLYPH<0> > 1 regime in a manner that allows for a smooth transition between coupling regimes.  \nCan these approximations be extrapolated to the strongly coupled regime of GLYPH<0> > 1? Baalrud & Daligault (2013) provide a method for numerically calculating resistance coe GLYPH<14> -cients using a hypernetted chain (HNC) approximation from e GLYPH<11> ective potentials. Figure 44 compares their HNC results (diamonds) to their Molecular Dynamics (MD) simulations of a one-component plasma (OCP, circles) for the self-di GLYPH<11> usion coe GLYPH<14> cient D GLYPH<3> , defined by  \nPaquette et al. (1986) successfully described di GLYPH<11> usion in a regime of intermediate coupling through the use of screened potentials, which are a way to account for the collective nature of interactions in a dense plasma. Though there is no rigorous reason to expect that a formalism based on the two-particle scattering picture should work well as GLYPH<0> ! 1, their comparison to simulations verified that this description of di GLYPH<11> usion is very accurate for GLYPH<0> . 1.  \nD GLYPH<3> = D GLYPH<21> 2 ion ! p ; (109)  \nwhere ! p is the plasma frequency and D = 2 D (2) ss (the factor of 2 in this definition ensures that if we redefine species s in terms of two subspecies s 1 and s 2, then D = D (2) s 1 s 2 ). The general expression for the interdi GLYPH<11> usion coe GLYPH<14> cient is  \nD (2) st = nsnt ns + nt k B T Kst (1 GLYPH<0> GLYPH<1> ) ; (110)  \nwhere the 1 GLYPH<0> GLYPH<1> term in the denominator accounts for a second order correction that can be defined using  \nGLYPH<1> = (2 GLYPH<6> (12) st GLYPH<0> 5 GLYPH<6> (11) st ) 2 = GLYPH<6> (11) st 55 GLYPH<6> (11) st GLYPH<0> 20 GLYPH<6> (12) st + 4 GLYPH<6> (13) st + 8 GLYPH<6> (22) st : (111)  \nFor reference, we also include a direct fit of Daligault & Murillo (2005) to the MD data of Ranganathan et al. (2003), given by  \nD GLYPH<3> = 0 : 0028 + 0 : 00525 173 GLYPH<0> GLYPH<0> 1 ! 1 : 154 : (112)  \nThe agreement between the HNC and MD simulations shows that the HNC does a better job of accounting for correlation physics in strongly coupled plasmas than a simple screened Coulomb potential and allows for a surprising (and still physically unexplained) extension of the Burgers formalism into the strongly coupled regime. This recent work allows us to go into the large GLYPH<0> limit with the Burgers formalism, but the question remains as to how we obtain di GLYPH<11> usion coe GLYPH<14> cients in a reliable manner.  \nThe self-di GLYPH<11> usion coe GLYPH<14> cients from the two options in MESA are shown in Figure 44 and correlate with the MD data better than expected for the high GLYPH<0> regime. In particular, the agreement is much better than that shown in figure 2 of Baalrud & Daligault (2013) for either 'cuto GLYPH<11> ' or 'screened' Coulomb methods. The reason for this agreement is that both MESA implementations use the inter-ion spacing rather than the Debye length once GLYPH<0> > 1 = 3, which yields favorable scalings in  \nFigure 44. Compilation of the self-di GLYPH<11> usion coe GLYPH<14> cients obtained from different methods. 'MD Data' and 'HNC' points are taken from Baalrud & Daligault (2013). The solid black line is the result of the MESA calculation using the coe GLYPH<14> cients of Paquette et al. (1986). The dashed green line is the result of the calculation using the resistance coe GLYPH<14> cients from the original routine of Thoul et al. (1994) based on the fit to the Coulomb logarithm found in Iben & MacDonald (1985), given here in Equation (96). The dashed purple line represents the fit to MD data given here in Equation (112).  \n<!-- image -->  \nthe high GLYPH<0> limit. Iben & MacDonald (1985) constructed their fitting formula based on a few numerical results for GLYPH<0> > 1. Paquette et al. (1986) also showed that their formalism can be extended to GLYPH<0> > 1 as long as the inter-ion spacing is used rather than the Debye radius for the screening length.  \nThough MESA does not yet provide the capability of implementing resistance coe GLYPH<14> cients based on the HNC method, we hope to accomplish this in the near future by means of a table similar to that provided for the coe GLYPH<14> cients of Paquette et al. (1986). For a more thorough discussion of these methods and the likely path of application to mixtures, consult Beznogov &Yakovlev (2014). We will also need to correctly account for the electron degeneracy and the non-ideal equation of state for the ions, both of which modify the electrostatic field needed to correctly determine the forces that drive di GLYPH<11> usion.", '10. SOFTWARE INFRASTRUCTURE': "Here we describe a number of changes to MESA that have occurred since Paper II and are of potential interest to users of MESA or developers of similar software.  \nPreviously, test cases in the MESA test suite accepted different results so long as they were within a certain tolerance, an appropriate choice for testing scientific results where the physical uncertainties are much greater than the numerical ones. However, we found that this made detecting and tracking bugs across platforms di GLYPH<14> cult. For the purposes of code testing, it is much better to insist that any inconsistency is a problem, no matter how small.  \nMESA can be compiled with either the GNU or Intel Fortran compilers, runs on multiple operating systems (Windows, OS X, and Linux), and can use di GLYPH<11> erent numbers of OpenMP threads. It is necessary to regularly test that the code is performing correctly across the di GLYPH<11> erent combinations of compiler, OS, thread count. To this end, developers and engaged users run the MESA test suite on a wide range of systems before each release.  \nMotivated by this challenge, MESA now provides bit-for-bit consistency for all results across all the supported platforms.  \nIt is essential to emphasize that the goal of this achievement is to enable better testing. It allows users to exactly reproduce the results of others, independent of platform di GLYPH<11> erences, which is especially useful to developers attempting to reproduce bugs. The achievement of bit-for-bit consistency is not a claim that the results of MESA calculations are physically accurate or numerically converged to any specific degree.  \nThis bit-for-bit consistency was achieved via the following choices:  \n- GLYPH<15> Using parallel algorithms that give identical results independent of number of threads or order of thread execution. MESA 's linear algebra solver is based on BCYCLIC (Hirshman et al. 2010). It sub-divides the work between threads based on the the size of the matrix rather than on the number of threads available. It is also necessary to avoid OpenMP reduction clauses, which provide no guarantees on ordering of operations.\n- GLYPH<15> Specifying compiler flags that forbid the compiler from making any optimization that can a GLYPH<11> ect floating point precision (e.g., forbid re-association and fast math operations). Most optimizations are still allowed.\n- GLYPH<15> Using an I / O library that does precise conversion from binary to ASCII for double precision numbers.\n- GLYPH<15> Using a math library that gives consistent results for operations such as log , exp , sin , cos , pow . MESA uses CRLIBM 16 in round towards zero mode. The choice to use a math library that gives exact results is not because 16 digit accuracy from the math routines in necessary. Rather, we want consistent results across supported platforms and this is the best way to achieve this consistency.\n- GLYPH<15> Replacing integer power expressions (i.e., x**3 ) by repeated multiplications (i.e., x*x*x ). Di GLYPH<11> erent compilers implement integer powers di GLYPH<11> erently, giving di GLYPH<11> erent results.  \nHaving achieved bit-for-bit identical results, we can test files for exact equality. This applies both to the module-bymodule tests that run at installation time and the case by case tests in the star and binary test suites. These test cases compare the final model from the test run to a saved result from a previous MESA version. If they are not exactly the same, the test fails. The test is also restarted from an intermediate state to confirm that runs which are stopped and restarted yield exactly the same results as those that are not.  \nSince Paper II the main MESA website 18 has undergone significant revision, making it easier for new users to get started  \nWhile MESAstar is parallelized via OpenMP, the install process has historically been serial. MESA contains approximately 1000 Fortran files and so the ability to compile more than one file simultaneously has the potential to provide significant reductions in the time needed to install MESA . Recently the compilation step has been parallelized, enabled by the automated dependency generation tool makedepf90 , 17 allowing multiple instances of the Fortran compiler to be invoked simultaneously. This is of particular utility for developers who may recompile MESA frequently.  \nwith MESA . This restructuring has also made it easier for the developers to keep material up-to-date as MESA evolves. One of the most important improvements is that the files that document the default value of each MESA option use the Markdown 19 markup language. This allows documentation web pages to be generated automatically for each MESA release.  \nImprovements have also been made to the distribution of MESA . Previously, MESA was available only by checking out the source code using the Subversion 20 version control system. Now, every release version of MESA (including past releases) is available for download as a ZIP archive. This is simpler and saves bandwidth and disk space. It has quickly become the preferred way to install MESA with the ZIP file of the current release being downloaded tens of times per week.", '11. SUMMARY AND CONCLUSIONS': "We have explained and, where possible, verified or validated, major new capabilities and improvements implemented in MESA since the publication of Paper I and Paper II. These advancements include interacting binary systems (Section 2), implicit hydrodynamics and shocks (Section 4), in-situ usage of large reaction networks, especially for X-ray bursts and core-collapse supernova progenitors (Section 5), and the explosion of massive stars (Section 6). These new capabilities will allow for extended exploration of core collapse progenitors and the sensitivity of shock nucleosynthesis to their explosion mechanism. The full coupling of MESA to the GYRE non-adiabatic pulsation instrument (Section 3) has already revealed the richness of the instability strips for massive stars and enables the continued growth of astero-seismology across the HR diagram. Progress in the treatment of mass accretion (Section 7) and weak reaction rates (Section 8) will improve studies of their impact on stellar evolution. We also discuss the domain of validity for particle di GLYPH<11> usion within MESA and describe a path forward for extending di GLYPH<11> usion into the regime relevant to WD sedimentation (Section 9). We also describe significant improvements to the infrastructure of MESA (Section 10). MESAstar input files and related materials for all the figures are available at http: // mesastar.org.  \nAs a community software instrument for stellar astrophysics new directions for MESA will be driven by: features useful to the MESA user community, advances in the physics modules, algorithmic developments, and architectural evolution. Potential examples include a treatment of ionization in  \nThese hitherto unpublished advancements have already enabled a number of studies in interacting binary systems (Wolf et al. 2013; Pavlovskii & Ivanova 2015; Vos et al. 2015) and stellar pulsations (P'apics et al. 2014; Stello et al. 2014; Quinn et al. 2015; Cunha et al. 2015), and led to the discovery of new features in the thermal runaway during the evolution of ONeMg cores towards AIC (Schwab et al. 2015). It also enabled the first three dimensional simulations of the final minutes of iron core growth in a massive star up to and including the point of core gravitational instability and collapse (Couch et al. 2015). In addition, these enhanced capabilities have allowed for applications of MESAstar that were not initially envisioned, such as the treatment of Magneto-Rotational Instability in stars (Wheeler et al. 2015), e GLYPH<11> ects of axions on nucleosynthesis in massive stars (Aoyama & Suzuki 2015), and particle physics beyond the Standard Model (Curtin & Tsai 2014).  \nthe equation of state for an arbitrary composition across an expanded region in the GLYPH<26> -T plane, non-linear pulsations, Monte Carlo based thermonuclear reaction rates, modules for subsonic flame propagation, ports to additional architectures, and a web-interface to MESA for education.  \nIt is a pleasure to thank Nilou Afsari, Dave Arnett, Warrick Ball, Ed Brown, Jieun Choi, Joergen Christensen-Dalsgaard, Andrew Cumming, Sebastien Deheuvels, Aaron Dotter, Chris Fryer, Duncan Galloway, Pascale Garaud, Alfred Gautschy, Samuel Jones, Max Katz, Eli Livne, Marcin Mackiewicz, Chris Mankovich, Casey Meakin, Broxton Miles, Ehsan Moravveji, Kevin Moore, Eliot Quataert, Jeremy Sakstein, Richard Stancli GLYPH<11> e, Willie Strickland, Anne Thoul, and Joris Vos. We also thank the participants of the 2013 and 2014 MESA Summer Schools for their willingness to experiment with new capabilities.  \nThis project was supported by NSF under the SI 2 program grants (ACI-1339581, ACI-1339600, ACI-1339606) and NASA under the TCAN program grants (NNX14AB53G, NNX14AB55G, NNX12AC72G). The work at UC Santa Barbara was also supported by the NSF under grants PHY 1125915, AST 11-09174, AST 12-05574. J.S. is supported by the NSF Graduate Research Fellowship Program under grant DGE 11-06400. L.D. acknowledges financial support by the 'Agence Nationale de la Recherche' under grant ANR-2011Blanc-SIMI-BS56-0007. D.T. acknowledges support under HST-GO-12870.14-A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. R.H.D.T. acknowledges resources provided by of the University of Wisconsin-Madison Advanced Computing Initiative. F.X.T. acknowledges support from the Simons Foundation.", 'REFERENCES': "- Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Springer Science + Business Media B.V) Aoyama, S., & Suzuki, T. K. 2015, ArXiv e-prints Arcones, A., Mart'ınez-Pinedo, G., Roberts, L. F., & Woosley, S. E. 2010, A&A, 522, A25 Arnett, W. D. 1969, Ap&SS, 5, 180 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411 Baalrud, S. D., & Daligault, J. 2013, Physical Review Letters, 110, 235001 -. 2014, Physics of Plasmas, 21, 055707 Baglin, A., Auvergne, M., Barge, P., et al. 2009, in IAU Symposium, Vol. 253, IAU Symposium, ed. F. Pont, D. Sasselov, & M. J. Holman, 71-81 Bailey, D. H. 1995, ACM Trans. Math. Softw., 21, 379 Ball, W. H., & Gizon, L. 2014, A&A, 568, A123 Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55 Bedding, T. R., Huber, D., Stello, D., et al. 2010, ApJL, 713, L176 Bedding, T. R., Mosser, B., Huber, D., et al. 2011, Nature, 471, 608 Beznogov, M. V., & Yakovlev, D. G. 2014, Phys. Rev. E, 90, 033102 Bildsten, L. 1995, ApJ, 438, 852 Bildsten, L., & Hall, D. M. 2001, ApJL, 549, L219 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 Brent, R. P. 1978, ACM Trans. Math. Softw., 4 Burgers, J. M. 1969, Flow Equations for Composite Gases (Academic Press) Calder, A. C., Townsley, D. M., Seitenzahl, I. R., et al. 2007, ApJ, 656, 313 Cantiello, M., & Langer, N. 2010, A&A, 521, A9 Cantiello, M., Yoon, S.-C., Langer, N., & Livio, M. 2007, A&A, 465, L29 Chaplin, W. J., & Miglio, A. 2013, ARA&A, 51, 353 Chen, W., Shrader, C. R., & Livio, M. 1997, ApJ, 491, 312 Christensen-Dalsgaard, J. 2008, Ap&SS, 316, 113 Christensen-Dalsgaard, J., & Thompson, M. J. 2011, in IAU Symposium, Vol. 271, IAU Symposium, ed. N. H. Brummell, A. S. Brun, M. S. Miesch, & Y. Ponty, 32-61 Cli GLYPH<11> ord, F. E., & Tayler, R. J. 1965, MmRAS, 69, 21 Collange, S., Defour, D., Graillat, S., & Iakymchuk, R. 2014, Full-Speed Deterministic Bit-Accurate Parallel Floating-Point Summation on Multiand Many-Core Architectures, Tech. rep., https: // hal.archives-ouvertes.fr / hal-00949355v1\n- Commins, E. 1973, Weak Interactions (McGraw-Hill)\n- Couch, S. M., Chatzopoulos, E., Arnett, W. D., & Timmes, F. X. 2015, ArXiv e-prints\n- C'orsico, A. H., Benvenuto, O. G., Althaus, L. G., & Serenelli, A. M. 2002, MNRAS, 332, 392\n- Cox, A. N., Morgan, S. M., Rogers, F. J., & Iglesias, C. A. 1992, ApJ, 393, 272\n- Cunha, M. S., Stello, D., Avelino, P. P., Christensen-Dalsgaard, J., & Townsend, R. H. D. 2015, ArXiv e-prints\n- Creevey, O. L., Th'evenin, F., Berio, P., et al. 2015, A&A, 575, A26 Cumming, A., & Bildsten, L. 2001, ApJL, 559, L127  \nCurtin, D., & Tsai, Y. 2014, Journal of High Energy Physics, 11, 136  \nDaligault, J., & Murillo, M. S. 2005, Phys. Rev. E, 71, 036408  \nCyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240  \nDaszy'nska-Daszkiewicz, J., Ostrowski, J., & Pamyatnykh, A. A. 2013,  \n- de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764  \nMNRAS, 432, 3153  \n- De Ridder, J., Barban, C., Baudin, F., et al. 2009, Nature, 459, 398 Deheuvels, S., Bruntt, H., Michel, E., et al. 2010, A&A, 515, A87 Demmel, J., & Hida, Y. 2003, SIAM J. Sci. Comput., 25, 1214 Dessart, L., Livne, E., & Waldman, R. 2010a, MNRAS, 408, 827\n- Detmers, R. G., Langer, N., Podsiadlowski, P., & Izzard, R. G. 2008, A&A, 484, 831 Donati, J.-F., & Landstreet, J. D. 2009, ARA&A, 47, 333\n- -. 2010b, MNRAS, 405, 2113\n- Dorfi, E. A. 1998, in Saas-Fee Advanced Course 27: Computational\n- Dziembowski, W. A., Moskalik, P., & Pamyatnykh, A. A. 1993, MNRAS, 265, 588\n- Methods for Astrophysical Fluid Flow., ed. O. Steiner & A. Gautschy, 263 Dupret, M.-A., Grigahc'ene, A., Garrido, R., Gabriel, M., & Scuflaire, R. 2004, A&A, 414, L17\n- Dziembowski, W. A., & Pamiatnykh, A. A. 1993, MNRAS, 262, 204 Eggleton, P. P. 1971, MNRAS, 151, 351 -. 1983, A&A, 268, 368  \n-. 2007b, Nuclear Data Sheets, 108, 2319  \n- Firestone, R. B. 2007a, Nuclear Data Sheets, 108, 1\n- -. 2009, Nuclear Data Sheets, 110, 1691\n- Fisker, J. L., Schatz, H., & Thielemann, F.-K. 2008, ApJS, 174, 261  \nFisker, J. L., Gorres, J., Wiescher, M., & Davids, B. 2006, ApJ, 650, 332  \n- Fontaine, G., & Michaud, G. 1979, in IAU Colloq. 53: White Dwarfs and Variable Degenerate Stars, ed. H. M. van Horn & V. Weidemann, 192-196 Fraley, G. S. 1968, Ap&SS, 2, 96  \nGalloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty,  \n- Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1 Galloway, D. K., Cumming, A., Kuulkers, E., et al. 2004, ApJ, 601, 466  \nD. 2008, ApJS, 179, 360  \n- Gansicke, B. T., Dillon, M., Southworth, J., et al. 2009, MNRAS, 397, 2170 Gehrels, N., Ramirez-Ruiz, E., & Fox, D. 2009, Annual Review of Astronomy and Astrophysics, 47, 567\n- Galloway, D. K., Psaltis, D., Muno, M. P., & Chakrabarty, D. 2006, ApJ, 639, 1033\n- Godart, M., Noels, A., Dupret, M.-A., & Lebreton, Y. 2009, MNRAS, 396, 1833\n- Grad, H. 1949, Communications on Pure and Applied Mathematics, 2, 331 Graillat, S., & M'enissier-Morain, V. 2012, Information and Computation,\n- Goldstein, J., & Townsend, R. H. D. 2015, in preparation\n- 216, 57 , special Issue: 8th Conference on Real Numbers and Computers Grigahc'ene, A., Uytterhoeven, K., Antoci, V., et al. 2010, Astronomische Nachrichten, 331, 989\n- Hartmann, D., Woosley, S. E., & El Eid, M. F. 1985, ApJ, 297, 837\n- Grott, M., Chernigovski, S., & Glatzel, W. 2005, MNRAS, 360, 1532 Guidry, M. W., Billings, J. J., & Hix, W. R. 2013, Computational Science and Discovery, 6, 015003\n- Heger, A., Cumming, A., Galloway, D. K., & Woosley, S. E. 2007, ApJL, 671, L141\n- Heger, A., Woosley, S. E., Mart'ınez-Pinedo, G., & Langanke, K. 2001, ApJ, 560, 307\n- Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368\n- Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 Herschel, W. 1802, Royal Society of London Philosophical Transactions Series I, 92, 477\n- Hirshman, S., Perumalla, K., Lynch, V., & Sanchez, R. 2010, Journal of Computational Physics, 229, 6392\n- Higham, N. J. 2002, Accuracy and Stability of Numerical Algorithms, 2nd edn. (Philadelphia, PA, USA: Society for Industrial and Applied Mathematics)\n- Hix, W. R., & Meyer, B. S. 2006, Nuclear Physics A, 777, 188 Howell, S. B., Nelson, L. A., & Rappaport, S. 2001, ApJ, 550, 897 Hu, H., Tout, C. A., Glebbeek, E., & Dupret, M.-A. 2011, MNRAS, 418, 195 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 Hut, P. 1981, ApJ, 99, 126\n- Itoh, N., Tomizawa, N., Tamamura, M., Wanajo, S., & Nozawa, S. 2002, ApJ, 579, 380\n- Iben, Jr., I., & MacDonald, J. 1985, ApJ, 296, 540\n- Jones, S., Hirschi, R., & Nomoto, K. 2014, ApJ, 797, 83 Jones, S., Hirschi, R., Nomoto, K., et al. 2013, ApJ, 772, 150\n- Juodagalvis, A., Langanke, K., Hix, W. R., Mart'ınez-Pinedo, G., & Sampaio, J. M. 2010, Nuclear Physics A, 848, 454 Kahan, W. 1965, Commun. ACM, 8, 40 Keek, L., Langer, N., & in't Zand, J. J. M. 2009, A&A, 502, 871 Khokhlov, A. M. 1991, A&A, 245, 114 Kippenhahn, R., Ruschenplatt, G., & Thomas, H.-C. 1980, A&A, 91, 175 Kippenhahn, R., & Weigert, A. 1967, Zeitschrift fur Astrophysik, 66, 251 Kjeldsen, H., Bedding, T. R., & Christensen-Dalsgaard, J. 2008, ApJL, 683, L175\n- Kobulnicky, H. A., Kiminki, D. C., Lundquist, M. J., et al. 2014, ApJS, 213, 34\n- Knuth, D. E. 1997, The Art of Computer Programming, Volume 2 (3rd Ed.): Seminumerical Algorithms (Boston, MA, USA: Addison-Wesley Longman Publishing Co., Inc.)  \nKolb, U., & Ritter, H. 1990, A&A, 236, 385  \n- Langer, N., Wellstein, S., & Petrovic, J. 2003, Proceedings of IAU Symposium, 212, 275\n- Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2006, Science, 314, 97 Langanke, K., & Mart'ınez-Pinedo, G. 2000, Nuclear Physics A, 673, 481 Langer, N. 2012, ARA&A, 50, 107\n- Law, N. M., Kulkarni, S. R., Dekany, R. G., et al. 2009, PASP, 121, 1395 Li, X. S., Demmel, J. W., Bailey, D. H., et al. 2002, ACM Trans. Math. Softw., 28, 152  \nLongland, R., Martin, D., & Jos'e, J. 2014, A&A, 563, A67  \n- Lin, J., Rappaport, S., Podsiadlowski, P., et al. 2011, ApJ, 732, 70 Livne, E. 1993, ApJ, 412, 634  \nLubow, S. H., & Shu, F. H. 1975, ApJ, 198, 385  \nMadhusudhan, N., Justham, S., Nelson, L., et al. 2006, ApJ, 640, 918  \nMacFadyen, A. I., & Woosley, S. E. 1999, ApJ, 524, 262  \nMaeder, A., & Meynet, G. 2000, ARA&A, 38, 143  \nMart'ınez-Pinedo, G., Lam, Y. H., Langanke, K., Zegers, R. G. T., &  \n- Marsh, T. R., Nelemans, G., & Steeghs, D. 2004, MNRAS, 350, 113  \nSullivan, C. 2014, Phys. Rev. C, 89, 045806  \n- Meyer, B. S., Krishnan, T. D., & Clayton, D. D. 1998, ApJ, 498, 808 Meynet, G., Ekstrom, S., Georgy, C., Chiappini, C., & Maeder, A. 2010, Evolution of Massive Stars along the Cosmic History (Wiley-VCH Verlag GmbH), 97-126\n- McLaughlin, G. C., Fuller, G. M., & Wilson, J. R. 1996, ApJ, 472, 440 McNamee, J. M. 2004, SIGSAM Bull., 38, 1 Mestel, L. 1968, MNRAS, 138, 359  \nMiles, B. J., Townsley, D. M., & Paxton, B. 2015, in preparation Miyaji, S., & Nomoto, K. 1987, ApJ, 318, 307  \nMuller, J.-M., Brisebarre, N., de Dinechin, F., et al. 2010, Handbook of  \n- Miyaji, S., Nomoto, K., Yokoi, K., & Sugimoto, D. 1980, PASJ, 32, 303  \nFloating-Point Arithmetic (Birkhauser Boston), 572, ACM G.1.0; G.1.2;  \n- Nadyozhin, D. K., & Yudin, A. V. 2004, Astronomy Letters, 30, 634 Nelder, J. A., & Mead, R. 1965, Comput. J, 7, 308 Noerdlinger, P. D. 1977, A&A, 57, 407  \nG.4; B.2.0; B.2.4; F.2.1., ISBN 978-0-8176-4704-9  \n-. 1978, ApJS, 36, 259  \n- Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457\n- Nomoto, K. 1982, ApJ, 253, 798  \nNomoto, K., Saio, H., Kato, M., & Hachisu, I. 2007, ApJ, 663, 1269  \n- Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231  \nNomoto, K., & Sugimoto, D. 1977, PASJ, 29, 765  \nOdrzywolek, A. 2012, Atomic Data and Nuclear Data Tables, 98, 852 Ofek, E. O., Sullivan, M., Shaviv, N. J., et al. 2014, ApJ, 789, 104 Ogita, T., Rump, S. M., & Oishi, S. 2005, SIAM J. Sci. Comput, 26, 2005 Oluseyi, H. M., Becker, A. C., Culliton, C., et al. 2012, AJ, 144, 9 Palla, F., & Stahler, S. W. 1990, ApJL, 360, L47 Pamyatnykh, A. A. 1999, Acta Astr., 49, 119  \nP'apics, P. I., Moravveji, E., Aerts, C., et al. 2014, A&A, 570, A8 Paquette, C., Pelletier, C., Fontaine, G., & Michaud, G. 1986, ApJS, 61, 177 Pavlovskii, K., & Ivanova, N. 2015, MNRAS, 449, 4415 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208 Paxton, B., Marchant, P., Schwab, J., et al. 2016, ApJS, 223, 18 Peters, P. C. 1964, Physical Review, 136, 1224 Petrovic, J., Langer, N., & van der Hucht, K. A. 2005, A&A, 435, 1013 Piro, A. L., Arras, P., & Bildsten, L. 2005, ApJ, 628, 401 Piro, A. L., & Bildsten, L. 2007, ApJ, 663, 1252 Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964 Potekhin, A. Y., Chabrier, G., & Rogers, F. J. 2009, Phys. Rev. E, 79, 016411 Powell, M. J. D. 2004, The NEWUOA software for unconstrained optimization without derivatives, Tech. Rep. NA05, Department of  \n- Papadopoulos, A., D'Andrea, C. B., Sullivan, M., et al. 2015, MNRAS, 449, 1215  \nApplied Mathematics and Theoretical Physics, University of Cambridge Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing, Vol. 2nd ed. (Cambridge: University Press)  \n- Quinn, S. N., White, T. R., Latham, D. W., et al. 2015, ApJ, 803, 49 Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1 Ranganathan, S., Johnson, R., & Woodward, C. 2003, Physics and Chemistry of Liquids, 41, 123\n- Rappaport, S., Joss, P. C., & Verbunt, F. 1983, ApJ, 275, 713\n- Rappaport, S., Joss, P. C., & Webbink, R. F. 1982, ApJ, 254, 616 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003 Ritter, H. 1988, A&A, 202, 93\n- Rump, S. M., Ogita, T., & Oishi, S. 2008, SIAM J. Sci. Comput., 31, 189 Sacco, G. G., Je GLYPH<11> ries, R. D., Randich, S., et al. 2015, A&A, 574, L7 Saio, H. 2011, MNRAS, 412, 1814\n- Roussel-Dupr'e, R. 1981, ApJ, 243, 329\n- Saio, H., Kuschnig, R., Gautschy, A., et al. 2006, ApJ, 650, 1111 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 444, 337 Schatz, H., Bildsten, L., Cumming, A., & Wiescher, M. 1999, ApJ, 524, 1014\n- Schwab, J., Quataert, E., & Bildsten, L. 2015, ArXiv e-prints\n- Schatz, H., Aprahamian, A., Barnard, V., et al. 2001, Physical Review Letters, 86, 3471\n- Seitenzahl, I. R., Timmes, F. X., Marin-Lafl'eche, A., et al. 2008, ApJL, 685, L129  \nSkumanich, A. 1972, ApJ, 171, 565  \n- Shamsuzzoha Basunia, M. 2011, Nuclear Data Sheets, 112, 1875 Shen, K. J., & Bildsten, L. 2007, ApJ, 660, 1444  \nSmartt, S. J. 2009, ARA&A, 47, 63  \n- Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, A&A, 327, 620\n- Smith, N. 2014, Annual Review of Astronomy and Astrophysics, 52, 487 Smolec, R., & Moskalik, P. 2008, Acta Astron., 58, 193\n- Stancli GLYPH<11> e, R. J., & Eldridge, J. J. 2009, MNRAS, 396, 1699\n- Sugimoto, D. 1970, ApJ, 159, 619  \nStello, D., Compton, D. L., Bedding, T. R., et al. 2014, ApJL, 788, L10  \n- Sugimoto, D., & Nomoto, K. 1975, PASJ, 27, 197\n- Tauris, T. M., & van den Heuvel, E. P. J. 2006, Formation and evolution of compact stellar X-ray sources (Cambridge University Press), 623-665 Thoul, A. A., Bahcall, J. N., & Loeb, A. 1994, ApJ, 421, 828\n- Tanaka, Y. 1989, in ESA Special Publication, Vol. 296, Two Topics in X-Ray Astronomy, Volume 1: X Ray Binaries. Volume 2: AGN and the X Ray Background, ed. J. Hunt & B. Battrick, 3-13\n- Tilley, D. R., Cheves, C. M., Kelley, J. H., Raman, S., & Weller, H. R. 1998, Nuclear Physics A, 636, 249\n- Toki, H., Suzuki, T., Nomoto, K., Jones, S., & Hirschi, R. 2013, Phys. Rev. C, 88, 015806\n- Timmes, F. X. 1999, ApJS, 124, 241\n- Torres, G., Andersen, J., & Gim'enez, A. 2010, A&A Rev., 18, 67 Tout, C. A., Aarseth, S. J., Pols, O. R., & Eggleton, P. P. 1997, MNRAS, 291, 732\n- Townsley, D. M., & Bildsten, L. 2004, ApJ, 600, 390\n- Townsend, R. H. D., & Teitler, S. A. 2013, MNRAS, 435, 3406\n- Townsley, D. M., & Gansicke, B. T. 2009, ApJ, 693, 1007\n- Tsuruta, S., & Cameron, A. G. W. 1970, Ap&SS, 7, 374\n- Townsley, D. M., Jackson, A. P., Calder, A. C., et al. 2009, ApJ, 701, 1582\n- Ubertini, P., Bazzano, A., Cocchi, M., et al. 1999, ApJL, 514, L27  \nVerbunt, F., & Zwaan, C. 1981, A&A, 100, L7  \n- Ulrich, R. K., & Burger, H. L. 1976, ApJ, 206, 509\n- Vos, J., Østensen, R. H., Marchant, P., & Van Winckel, H. 2015, A&A, 579, A49  \nWeinberg, N. N., Bildsten, L., & Schatz, H. 2006, ApJ, 639, 1018  \n- Vreeswijk, P. M., Savaglio, S., Gal-Yam, A., et al. 2014, ApJ, 797, 24 Wallace, R. K., & Woosley, S. E. 1981, ApJS, 45, 389 Weaver, T. A., Zimmerman, G. B., & Woosley, S. E. 1978, ApJ, 225, 1021  \nWeisberg, J. M., & Taylor, J. H. 2005, ASP Conference Series, 328, 25  \n- Wheeler, J. C., Kagan, D., & Chatzopoulos, E. 2015, ApJ, 799, 85  \nWellstein, S., Langer, N., & Braun, H. 2001, A&A, 369, 939  \n- Whyte, C. A., & Eggleton, P. P. 1980, MNRAS, 190, 801\n- Woosley, S., & Bloom, J. 2006, Annual Review of Astronomy and Astrophysics, 44, 507\n- Wolf, W. M., Bildsten, L., Brooks, J., & Paxton, B. 2013, ApJ, 777, 136 Wood, P. R. 1974, ApJ, 190, 609\n- Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015 Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181\n- Zamfir, M., Cumming, A., & Galloway, D. K. 2012, ApJ, 749, 69 Zdravkov, T., & Pamyatnykh, A. A. 2008, Journal of Physics Conference Series, 118, 012079\n- Woosley, S. E., Heger, A., Cumming, A., et al. 2004, ApJS, 151, 75 Yuan, F., Lidman, C., Davis, T. M., et al. 2015, ArXiv e-prints Zahn, J.-P. 1977, A&A, 57, 383"}