We update the capabilities of the software instrument Modules for Experiments in Stellar Astrophysics MESA and enhance its ease of use and availability. Our new approach to locating convective boundaries is consistent with the physics of convection and yields reliable values of the convectivecore mass during both hydrogen and heliumburning phases. Stars with Mlt 8 MSUBSUB become white dwarfs and cool to the point where the electrons are degenerate and the ions are strongly coupled a realm now available to study with MESA due to improved treatments of element diffusion latent heat release and blending of equations of state. Studies of the final fates of massive stars are extended in MESA by our addition of an approximate Riemann solver that captures shocks and conserves energy to high accuracy during dynamic epochs. We also introduce a 1D capability for modeling the effects of RayleighTaylor instabilities that in combination with the coupling to a public version of the STELLA radiation transfer instrument creates new avenues for exploring Type II supernova properties. These capabilities are exhibited with exploratory models of pairinstability supernovae pulsational pairinstability supernovae and the formation of stellarmass black holes. The applicability of MESA is now widened by the capability to import multidimensional hydrodynamic models into MESA. We close by introducing software modules for handling floating point exceptions and stellar model optimization as well as four new software tools  MESAWeb MESADocker pyMESA and mesastar.org  to enhance MESAs education and research impact.
['convection', 'diffusion', 'hydrodynamics', 'methods: numerical', 'stars: evolution', 'supernovae: general', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena']
Modules for Experiments in Stellar Astrophysics MESA Convective Boundaries Element Diffusion and Massive Star Explosions
{'MODULES FOR EXPERIMENTS IN STELLAR ASTROPHYSICS (MESA): CONVECTIVE BOUNDARIES, ELEMENT DIFFUSION, AND MASSIVE STAR EXPLOSIONS': "Bill Paxton, 1 Josiah Schwab, 2, ∗ Evan B. Bauer, 3 Lars Bildsten, 1, 3 Sergei Blinnikov, 4, 5, 6 Paul Duffell, 7 R. Farmer, 8, 9 Jared A. Goldberg, 3 Pablo Marchant, 10 Elena Sorokina, 4, 5 Anne Thoul, 11 Richard H. D. Townsend, 12 and F. X. Timmes 9  \n- 1 Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA\n- 2 Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA\n- 3 Department of Physics, University of California, Santa Barbara, CA 93106, USA\n- 4 M.V. Lomonosov Moscow State University, Sternberg Astronomical Institute, Universitetsky pr., 13, Moscow, 119234, Russia\n- 5 NIC Kurchatov Institute - Institute for Theoretical and Experimental Physics (ITEP), 117218 Moscow, Russia\n- 6 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo, Kashiwa, Chiba 277-8583, Japan\n- 7 Department of Astronomy and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720, USA\n- 8 Anton Pannenkoek Institute for Astronomy, University of Amsterdam, NL-1090 GE Amsterdam, the Netherlands\n- 9 School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA\n- 10 Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA\n- 11 Space sciences, Technologies and Astrophysics Research (STAR) Institute, Universit'e de Li'ege, All'ee du 6 Ao ˆ u t 19C, Bat. B5C, 4000 Li'ege, Belgium\n- 12 Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA", 'ABSTRACT': "We update the capabilities of the software instrument Modules for Experiments in Stellar Astrophysics ( MESA ) and enhance its ease of use and availability. Our new approach to locating convective boundaries is consistent with the physics of convection, and yields reliable values of the convective core mass during both hydrogen and helium burning phases. Stars with M < 8 M glyph[circledot] become white dwarfs and cool to the point where the electrons are degenerate and the ions are strongly coupled, a realm now available to study with MESA due to improved treatments of element diffusion, latent heat release, and blending of equations of state. Studies of the final fates of massive stars are extended in MESA by our addition of an approximate Riemann solver that captures shocks and conserves energy to high accuracy during dynamic epochs. We also introduce a 1D capability for modeling the effects of Rayleigh-Taylor instabilities that, in combination with the coupling to a public version of the STELLA radiation transfer instrument, creates new avenues for exploring Type II supernovae properties. These capabilities are exhibited with exploratory models of pair-instability supernova, pulsational pair-instability supernova, and the formation of stellar mass black holes. The applicability of MESA is now widened by the capability of importing multi-dimensional hydrodynamic models into MESA . We close by introducing software modules for handling floating point exceptions and stellar model optimization, and four new software tools -MESA-Web , MESA -Docker, pyMESA , and mesastar.org -to enhance MESA 's education and research impact.  \nKeywords: convection - diffusion - hydrodynamics - methods: numerical - stars: evolution supernovae: general  \n∗  \nHubble Fellow", '1. INTRODUCTION': "Over the next decade multi-messenger astronomy will probe the rich stellar astrophysics of transient phenomena in the sky, including gravitational waves from the mergers of neutron stars and black holes, light curves and spectra from core-collapse supernovae, and the oscillation modes of stars. On the observational side of this new era, the Laser Interferometer GravitationalWave Observatory has demonstrated the existence of binary stellar-mass black hole systems (Abbott et al. 2016a,b,c, 2017a,b) and continues to monitor the sky with broadband detectors for gravitational waves from compact binary inspirals and asymmetrical exploding massive stars (Fryer et al. 2002; Gossan et al. 2016; Abbott et al. 2016d,e,f, 2017a). The Gaia Data Release 1,  \ncontaining about one billion stars, begins the process of converting the spectrophotometric measurements to distances, proper motions, luminosities, effective temperatures, surface gravities, and elemental compositions (Gaia Collaboration et al. 2016a,b). This stellar census will provide the observational data to tackle a range of questions related to the origin, structure, and evolutionary history of stars in the Milky Way (Creevey et al. 2015; Sacco et al. 2015; Lindegren et al. 2016; van Leeuwen et al. 2017). The Neutron star Interior Composition Explorer mission, delivered to the International Space Station in June 2017, will provide rotationresolved spectroscopy of the thermal and non-thermal emissions of neutron stars in the soft X-ray band with over 15 million seconds of exposures (Gendreau et al. 2012; Arzoumanian et al. 2014; Gendreau et al. 2016) to open a new window into the interior structure and dynamics that underlie neutron stars (e.g., Ozel et al. 2016; Miller 2016). With first light at Palomar Observatory in 2017, the Zwicky Transient Facility (Kulkarni 2016) will scan more than 3750 deg 2 hr -1 to a depth of about 20 mag to discover young supernovae less than 24 hours after explosion each night, hunt for electromagnetic counterparts of gravitational-wave events (Ghosh et al. 2017), and search for rare and exotic transients. Repeated imaging of the Northern sky, including the Galactic Plane, will produce a photometric variability catalog with nearly 300 observations each year (Laher et al. 2017) for detailed studies of variable stars and binary systems. From its unique high earth orbit, the Transiting Exoplanet Survey Satellite aims to survey about 200,000 nearby G, K and M type stars with apparent magnitudes brighter than about 12 mag with a 1 minute cadence across a 400 deg 2 area of the sky (Ricker et al. 2016; Sullivan et al. 2015, 2017) to open a new era on stellar variability. The Large Synoptic Survey Telescope will image the entire Southern Hemisphere deeply in multiple optical colors every week with a 3.5 deg 2 , three billion pixel digital camera (LSST Science Collaboration et al. 2017) to open new perspectives on transient objects such as tidal disruption events (Bade et al. 1996; Stern et al. 2004; Arcavi et al. 2014; Komossa 2015) and interacting close binary systems (Oluseyi et al. 2012; Korol et al. 2017). The Jiangmen Underground Neutrino Observatory will usher in a new generation of multipurpose neutrino detectors (Li 2014; Brugi'ere 2017) designed in part to open a new avenue on neutrinos from pre-supernova massive stars (e.g., Odrzywolek 2009; Misch & Fuller 2016; Patton et al. 2017a,b) and core-collapse supernova explosions (e.g., Hirata et al. 1987; Janka 2017).  \nThis ongoing explosion of activity in multi-messenger stellar astronomy powers theoretical and computational developments, in particular the evolution of the community software instrument Modules for Experiments in Stellar Astrophysics ( MESA ) for research and education. We introduce MESA in Paxton et al. (2011, Paper I) and significantly expand its range of capabilities in Paxton et al. (2013, Paper II) and Paxton et al. (2015, Paper III). These prior papers, as well as this one, are 'instrument' papers that describe the capabilities and limitations of MESA while also comparing to other available numerical or analytic results. This paper describes the major new advances to MESA for modeling convective boundaries, element diffusion, implicit shock hydrodynamics, massive star explosions and light curves, pulsational pair-instability supernovae, and black hole formation. We do not fully explore these results and their implications here. The scientific potential of these new capabilities will be unlocked in future work via the efforts of the MESA user community.  \nThe convective regions of stars remain a rich site of fascinating challenges including the interplay between mixing, composition gradients, and element diffusion. A convection region transports energy through the vertical exchange of matter. The location where the radial velocity of the bulk motions goes to zero is a natural way to define the edge of a convection region (Vitense 1953; Bohm-Vitense 1958). It is necessary to ensure that convective boundaries are properly positioned (e.g., Eggleton 1972; Gabriel et al. 2014), because their exact placement can have a strong influence on the evolution of the stellar model (Salaris & Cassisi 2017). An important new addition to MESA is an improved treatment of convective boundaries, allowing them to evolve toward a state where the radiative gradient equals the adiabatic gradient on the convective side of the boundary. As a consequence, the Schwarzchild and Ledoux criteria now give the same position for convective boundaries.  \nGradients can drive changes in the composition profile of a star. For example, if gradients occur in the concentrations of chemical elements, then diffusion tends to smooth out the differences. Temperature gradients can push heavier species towards regions of higher temperature, while pressure gradients can propel heavier species to diffuse towards regions of higher pressure (Thoul et al. 1994; Hansen et al. 2004; Kippenhahn et al. 2012; Michaud et al. 2015). Treatments of diffusion typically assume that all diffusing species are ideal gases (e.g., Burgers 1969; Thoul et al. 1994). For white dwarf interiors and neutron star envelopes, degenerate electrons violate this assumption (Deloye & Bildsten 2002; Chang et al. 2010). In addition, strong Coulomb cou-  \npling in plasmas requires modifications to the binary scattering formalism for calculating cross-sections used to obtain diffusion coefficients (Paquette et al. 1986a; Stanton & Murillo 2016; Daligault et al. 2016; Shaffer et al. 2017). MESA 's extensions of element diffusion for degenerate and strongly coupled plasmas open a pathway into the regime relevant to sedimentation in the interiors of white dwarfs (Iben & MacDonald 1985; Iben et al. 1992; Koester 2009; Hollands et al. 2017) and the surfaces of neutron stars (Chang & Bildsten 2003, 2004; Beznogov et al. 2016).  \nMassive ( M glyph[greaterorsimilar] 8 M glyph[circledot] ) stars explode when energy from the collapse of their core to a compact object emerges as an outgoing shock wave into the outer parts of the star. The outward propagation of this shock wave generates Rayleigh-Taylor instabilities that can mix material behind the shock front (Chevalier 1976; Chevalier & Klein 1978; Weaver & Woosley 1980; Benz & Thielemann 1990; Herant & Benz 1991; Hammer et al. 2010; Wongwathanarat et al. 2015; Utrobin et al. 2017). The resulting light curves of Type II supernovae can be sub-divided into multiple classes but we focus here on Type IIP supernovae (e.g., Smartt 2009a, 2015; Smith et al. 2016). Our improvements to MESA - implicit shock capturing hydrodynamics, Rayleigh-Taylor instability modeling in 1D (Duffell 2016), and radiative transfer using the public version of the STELLA instrument (Blinnikov & Sorokina 2004; Baklanov et al. 2005; Blinnikov et al. 2006) open up new avenues for researching the diverse set of Type II supernovae.  \nPair-instability leads to a partial collapse, which in turn causes runaway thermonuclear burning in the carbon-oxygen core (Fowler & Hoyle 1964; Rakavy & Shaviv 1967; Barkat et al. 1967; Rakavy et al. 1967; Fraley 1968). A wide variety of outcomes is possible depending on the star's mass and rotation. A single energetic burst from nuclear burning can disrupt the entire star without leaving a black hole remnant behind to produce a pair-instability supernova (Ober et al. 1983; Fryer et al. 2001; Scannapieco et al. 2005; Kasen et al. 2011; Chatzopoulos et al. 2013). Alternatively, a series of bursts can trigger a cyclic pattern of nuclear burning, expansion and contraction, leading to a pulsational pair-instability supernova that leaves a black hole remnant (Barkat et al. 1967; Woosley et al. 2007a; Chatzopoulos & Wheeler 2012; Woosley 2017; Limongi 2017). Many of these variations can now be explored in MESA , as can lower mass progenitors that do not pulse before collapse to a black hole.  \nMESA is a community-driven software instrument for stellar astrophysics. New directions will be motivated by features useful to the MESA user community, ad-  \nvances in the physics modules, algorithmic developments, and architectural evolution. Potential examples for expanding MESA 's scientific, computational, and educational capabilities include seamlessly leveraging many-core architectures, an improved treatment of the equation of state, Jupyter/Python notebooks for education, and continued integration with software instruments useful to the astronomy and astrophysics community. Examples include ADIPLS (Christensen-Dalsgaard 2008; Christensen-Dalsgaard & Thompson 2011), GYRE (Townsend & Teitler 2013), and STELLA (Blinnikov et al. 1998; Blinnikov & Sorokina 2004; Baklanov et al. 2005; Blinnikov et al. 2006).  \nThe paper is organized as follows. Section 2 introduces a new treatment of convective boundaries. In Section 3 we present an implementation of element diffusion that accounts for electron degeneracy and strongly coupled interactions. Section 4 describes the Riemann solver for shock capturing in MESA 's new implicit hydrodynamics solver, and Section 5 presents a model for approximating the 3D effects of the Rayleigh-Taylor instability. In Section 6 we introduce the coupling of MESA and an implementation of the STELLA radiative transfer instrument to explore the modeling of Type IIP supernova light curves from post-explosion to post-plateau. In Section 7 we show advances to model pair-instability supernova, pulsation pair-instability supernova, and black hole formation. Section 8 discusses energy accounting in stellar evolution.  \nAppendix A discusses improvements to estimating a model's absolute magnitude in a chosen color filter, Appendix B offers guidance on importing multidimensional models into MESA , and Appendix C details the implementation of element diffusion in MESA . Appendix D introduces two new software modules for handling floating point exceptions and stellar model optimization, and four new software tools for education and research: MESA -Web, MESA -Docker, pyMESA , and mesastar.org .  \nImportant symbols are defined in Table 1. Acronyms used are denoted in Table 2. We denote components of MESA , such as modules and routines, in typewriter font e.g., colors .  \nTable 1 . Important Variables. Single character symbols are listed first, symbols with modifiers are listed second. Some symbols may be further subscripted by c (indicating a central quantity), by a cell index k , or by an index that runs over species ( i , j , s , or t ).  \nTable 1 (continued)Table 2 . Acronyms used in this paper.  \nTable 2 (continued)", '2. CONVECTIVE BOUNDARIES': "Gabriel et al. (2014) discuss the correct positioning of convective boundaries in stellar evolution models. Following earlier work (e.g., Roxburgh 1978) they argue that a convective boundary should be defined as the point where the convective velocity vanishes. Within local mixing-length theory (MLT), this condition is equivalent to the requirement ∇ rad = ∇ ad , where ∇ rad and ∇ ad are the radiative and adiabatic temperature gradients, respectively. Critically, this equality must be satisfied on the convective side of the boundary, because the MLT convective velocity is only well defined there. Moreover, because the fluid on the convective side is presumed to be well-mixed, the Ledoux temperature gradient ∇ L = ∇ ad + B (Equation 11 of Paper II) can play no part in setting the location of the boundary.  \nIf the chemical composition is continuous across the convective boundary, then so too are ∇ rad and ∇ ad , and requiring ∇ rad = ∇ ad on the convective side of the boundary results in the same equality on the radiative side. However, a composition discontinuity produces a jump in density and opacity, and in turn a discontinuity in ∇ rad and ∇ ad . Hence, it is generally the case that ∇ rad = ∇ ad on the radiative side of the boundary.  \nglyph[negationslash]  \nIn numerical codes based on discrete grids, the nuance of the foregoing discussion is often overlooked in favor of a simple approach for locating convective boundaries based on sign changes in the discriminant y = ∇ rad -∇ ad (or y = ∇ rad - ∇ L , if the Ledoux stability criterion is used). This approach works well when the chemical composition remains continuous, but is problematic when the composition - and hence y - is discontinuous at the boundary; it typically leads to configurations where ∇ rad > ∇ ad on the convective side, which is unphysical and ultimately retards the growth of the convective region. Previous versions of MESA have taken this approach; the outcome is evident in Figure 15 of Paper II, which shows the convective core mass as a function of age during the He-burning evolution of a 3 M glyph[circledot] star. In the model with no overshoot and the  \nSchwarzschild stability criterion, the core grows only modestly in mass before reaching a plateau. Inspection of the model confirms that ∇ rad > ∇ ad on the convective side of the core boundary, signifying that the core growth is being impeded.  \nGabriel et al. (2014) highlight a further issue with this simple sign-change approach, whereby the location of a convective boundary is not uniquely determined but rather depends on the mixing history near the boundary. We have confirmed this issue is present in MESA when using the sign-change approach. This manifests itself as a lack of convergence in some models (e.g., the 3 M glyph[circledot] He-burning example) when the resolution is increased and/or the timestep shortened.  \nTo resolve these issues we implement a new 'predictive mixing' scheme in MESA . It is inspired both by the 'maximal overshoot' scheme introduced by Constantino et al. (2015), and by the procedure described by Bossini et al. (2015). In the new scheme, the extent of a convection region is allowed to expand at each time step until the boundaries reach the point where ∇ rad = ∇ ad on their convective side. We describe the new scheme in detail in the following section and then present results obtained with this scheme in four scenarios: a growing convective core in a low-mass star on the main sequence (MS), a retreating convective core in a high-mass star on the MS, growing He-burning cores in intermediateand low-mass stars, and a surface convective region in a low-mass star on the MS. In all cases, we assume an initial He mass fraction Y = 0 . 28, an initial metal mass fraction Z = 0 . 02, and we neglect rotation and mass loss.", '2.1. Predictive Mixing': "The MESA predictive mixing scheme initially proceeds in the same manner as the simple sign-change approach, by finding the cells where y > 0 on one face (convective) and y < 0 on the other face (radiative). For each of these candidate boundary cells, the algorithm considers how y would change if the cell were completely mixed with the rest of the adjoining convection region. This prediction involves re-evaluating opacities, densities and other data throughout the mixed region, under the assumption that the composition is completely uniform. If y would become positive on both faces of the candidate boundary cell, then the adjacent cell in the radiative region becomes the new candidate boundary cell and a new round of predictive mixing begins. The process continues iteratively until the candidate cell after the predictive mixing still has a negative y on the radiative face. The code reverts to the previous candidate, identifies it as the final convective boundary cell, recalculates con-  \nvective diffusivities and convective velocities using MLT, and writes these into the model for use in the composition solver (see Paper I, Section 6.2). No abundances are directly modified in the model during the predictive iterations. Below, we demonstrate that this algorithm leads MESA to a solution of the stellar structure equations in which 0 ≤ y glyph[lessmuch] ∇ ad on the convective side of each boundary cell.  \nThe physical justification for our predictive mixing scheme traces back to a narrative advanced by Castellani et al. (1971). Focusing on He core-burning, these authors argue that any gentle mixing outside the core boundary irreversibly alters the composition there, and the resulting increase in opacity raises the local ∇ rad from sub-adiabatic to super-adiabatic. The outcome is a 'self-driving mechanism for the extension of the convective region', which continues until ∇ rad = ∇ ad on the convective side of the core boundary. While Castellani et al. (1971) invoked overshoot as the source of the mixing outside the boundary, Michaud et al. (2007) show that element diffusion can serve equally well and leads to the same outcome.  \nFor MS stars with growing convective cores, the extension of the core boundary cannot be driven in exactly the same way as the He-burning case, because helium has a lower opacity than hydrogen. However, gentle mixing outside the core boundary erases any composition gradients there, and it is the loss of these gradientsand their accompanying stabilizing effect-that drives the extension of the convective region until ∇ rad = ∇ ad on the convective side of the core boundary.  \nThe predictive mixing scheme doesn't specify the nature of the gentle mixing beyond convective boundaries, instead focusing on its effects. Tied in with this agnosticism is the presumption that the mixing-driven expansion of convective boundaries is so rapid that it can be approximated as instantaneous. This is likely a reasonable approach during core H and He burning; Castellani et al. (1971) argue that the growth of the core boundary in the latter case should proceed on a timescale which is much shorter than the burning lifetime. However, there may be circumstances where the finite timescale for boundary growth cannot be neglected.  \nBecause uniform composition is assumed during the predictive mixing iterations, there is no functional distinction between the Schwarzschild and Ledoux criteria when evaluating the discriminant y . However, the preliminary search for sign changes in y , before any predictions are made, does take into account composition gradients when the Ledoux criterion is used. As a result, the initial candidate boundary cells can differ between the two criteria. In many cases this difference is  \nFigure 1. The mass coordinate m bdy of the convective-core boundary plotted as a function of MS age, for the 1 . 5 M glyph[circledot] stellar model discussed in Section 2.2. Different line styles/colors show the separate runs described in the text.  \n<!-- image -->  \nunimportant, with the final location of the boundaries being insensitive to the choice of criterion. The one exception is when a region with ∇ ad < ∇ rad < ∇ L is bounded on both sides by radiative regions; then, it will be completely overlooked during a preliminary search with the Ledoux criterion. As we shall demonstrate later, such scenarios arise in our calculations outside convective cores during MS evolution.", '2.2. Evolution of a Growing Convective Core on the Main Sequence': 'We evolve a 1 . 5 M glyph[circledot] star from the zero age main sequence (ZAMS) to the terminal age main sequence (TAMS) using the predictive mixing scheme at the convective core boundary; this is the same mass and evolutionary stage considered in Section 5.1.2 of Gabriel et al. (2014). Figure 1 plots the mass coordinate of the convective-core boundary as a function of MS age, showing results from separate runs using Schwarzschild and Ledoux criteria, and additional runs with the incremental inclusion of semi-convection (in just the Ledoux case) and then element diffusion (in both cases). The semiconvection is modeled using the Langer et al. (1985) scheme with an efficiency parameter α sc = 0 . 1 (see Paper II for a complete description of the semi-convection implementation in MESA ). For comparison, the figure also shows the outcome of using the Ledoux criterion but no predictive mixing; in contrast to the other cases which broadly agree with one another, the core growth is inhibited and the H-burning lifetime correspondingly truncated.  \nNote that in Figure 13 of Paper II the results obtained with the Ledoux criterion show a shrinking convective  \nFigure 2. Profiles of ∇ rad , ∇ ad , ∇ L , and X as a function of mass coordinate, in the inner part of the 1 . 5 M glyph[circledot] stellar model at X c = 0 . 42. The panels show the separate runs described in the text. Gray shading indicates regions undergoing convection. Unless otherwise indicated, all models used predictive mixing.  \n<!-- image -->  \ncore; this behavior was due to a separate problem arising from over-smoothing of the composition gradient (see Moore & Garaud 2016) and has since been rectified in MESA . For completeness, we include this case in Figure 1.  \nFigure 2 plots the profiles of ∇ rad , ∇ ad , ∇ L and X , in the inner part of the 1 . 5 M glyph[circledot] star nearing the halfway point of its MS evolution (a core H mass fraction X c = 0 . 42). In the upper row, the left panel il-  \ntes the run with the Ledoux criterion plus predictive mixing (the dotted curve in Figure 1), while the right panel shows the run with the Ledoux criterion but without predictive mixing (the black curve, ibid.). Clearly, without predictive mixing ∇ rad remains significantly larger than ∇ ad on the convective side of the boundary, which as discussed previously is physically inconsistent. When using predictive mixing, however,  \nthe profiles satisfy ∇ rad = ∇ ad on the convective side, and closely match those seen in the left panel in Figure 6 of Gabriel et al. (2014). The small bump in ∇ rad just above the boundary is Schwarzschild unstable but Ledoux stable.  \nThe middle panels of Figure 2 show the runs with the Ledoux criterion and predictive mixing, and the incremental addition of semi-convection (left) and then element diffusion (right). Inside the core boundary, the profiles are almost identical to those shown in the upper-left panel but just outside the boundary, the semiconvection converts the composition discontinuity into a steep gradient and flattens the bump in ∇ rad to a neutral, ∇ rad = ∇ ad profile. Element diffusion further softens the abundance profile, as shown in the middle-right panel. Note that element diffusion has only a small effect on the location of the convective boundary; this is barely noticeable in Figure 2, but a slight extension of the boundary can be seen in Figure 1 toward the later part of the MS, for the two cases including diffusion.  \nThe lower panels of Figure 2 show the runs using the Schwarzschild criterion and predictive mixing, without (left) and with (right) element diffusion. In the left panel, the abundance profile shows a chaotic staircaselike profile, due to mixing by transient convective shells that appear and disappear from one timestep to the next (two of these shells can be seen in the figure). The shells do not appear in the Ledoux plots (middle and upper panels) because the region outside the core is stabilized in its entirety by the abundance gradient: ∇ rad < ∇ L . This serves as a good illustration of the earlier discussion (Section 2.1) of how the Schwarzschild and Ledoux criteria can sometimes lead to different outcomes. It is important to note, however, that the location of the core boundary is the same in all cases with predictive mixing; the differences only appear in the inhomogeneous region beyond the boundary which arises from slow H burning outside the core.  \nThe lower-right panel of Figure 2 shows that adding element diffusion removes the abundance discontinuities, replacing them with a smooth gradient. The resulting profiles appear almost identical to the Ledoux case shown in the middle-right panel of the figure (and compare also the curves with diffusion in Figure 1).', '2.3. Evolution of a Retreating Convective Core on the Main Sequence': "We now evolve a 16 M glyph[circledot] star from ZAMS to TAMS using the new predictive mixing scheme at the convective core boundary; this is the same mass and evolutionary stage considered in Section 5.1.1 of Gabriel et al. (2014). Figure 3 plots the mass of the convective core  \nFigure 3. The mass coordinate m bdy of the convectivecore boundary as a function of MS age, for the 16 M glyph[circledot] stellar model discussed in Section 2.3. Different line styles/colors show the separate runs described in the text.  \n<!-- image -->  \nas a function of MS age, showing results from separate runs using the Schwarzschild and Ledoux criteria, and with and without predictive mixing. The agreement between these four cases is very close. However, as was the case in the preceding section, there are differences outside the convective core. These can be seen in Figure 4, which plots the profiles of ∇ rad , ∇ ad , ∇ L and X near the end of the star's MS evolution ( X c = 0 . 15), for the two runs with predictive mixing.  \nEven though both runs exhibit the same core structure, with ∇ rad = ∇ ad at the convective side of the core boundary, the inhomogeneous region left behind by the retreating core is very different. The H abundance obtained with the Schwarzschild criterion shows the same staircase-like profile seen in the lower-left panel of Figure 2, again due to mixing by transient convective shells. These shells are not present when the Ledoux criterion is used, with the exception of a persistent solitary shell at the top of the inhomogeneous region (corresponding to where the core boundary was located at the ZAMS); the behavior of this shell is discussed by Gabriel et al. (2014, their Section 5.5.1; and compare also against their Figure 4). Between the shell and the core boundary, the abundance profile from the Ledoux run remains relatively smooth. The different abundance profiles in the two runs will have a direct influence on the Brunt-Vaisala frequency profile, and therefore on the oscillation frequencies of the stellar model.  \n2.4. Evolution of the Convective Core during Core He Burning  \n<!-- image -->  \n/circledot  \n/circledot  \nFigure 4. Profiles of ∇ rad , ∇ ad , ∇ L , and X as a function of mass coordinate, in the inner part of the 16 M glyph[circledot] stellar model at X c ≈ 0 . 15. The panels show the separate runs described in the text. Gray shading indicates regions undergoing convection. Both models use predictive mixing at the convective core boundary.Figure 5. The mass coordinate m bdy of the convectivecore boundary plotted as a function of CHeB age, for the 1 M glyph[circledot] stellar model discussed in Section 2.4. Different line styles/colors show the separate runs described in the text.  \n<!-- image -->  \nAs reviewed by Salaris & Cassisi (2017), the modeling of mixing in low- and intermediate-mass stars during core He burning (CHeB) is particularly challenging. Correct treatment of convective boundaries is complicated by the fact that the ∇ rad profile within the core convection region develops a local minimum at some point during CHeB evolution (see the middle and lower panels of Figure 6). This is a consequence of the complex behavior of the physical quantities (opacity, temperature, density, etc.) involved in the expression for ∇ rad . With further outward propagation of the convective boundary, the mixing of fresh He into the core can lower the radiative gradient throughout the core to  \nsuch an extent that ∇ rad = ∇ ad at the local minimum of ∇ rad . When this happens, the part of the convection region interior to the minimum becomes decoupled from the part exterior to the minimum: the convection region has split. This phenomenon was first discussed by Eggleton (1972), and a variety of ad-hoc approaches have been proposed to follow the subsequent evolution, mostly focused around the narrative that the exterior part undergoes partial mixing with the adjacent radiative region until it reaches convective neutrality (see, e.g., Castellani et al. 1985, and references therein).  \nAnother problem appears near the end of CHeB. At that point, even small amounts of He added to the core (which is almost totally depleted in He) will enhance the rate of energy production and thus the luminosity, resulting in an increase in ∇ rad . This increase leads to a sudden growth in the core boundary and a 'breathing pulse'. The He is then quickly burned in the core and the star re-adjusts itself. The existence of these breathing pulses remains controversial and it is still unclear whether they are numerical or physical (Caputo et al. 1989; Cassisi et al. 2003; Farmer et al. 2016; Constantino et al. 2017). All of these problems are clearly described and illustrated in Salaris & Cassisi (2017).  \nTo manage these complexities, the predictive mixing scheme must be modified. When a convection region splits, it is no longer meaningful to re-evaluate y using opacities and other data calculated on the assumption of uniform composition throughout (Section 2.1), because the radiative region appearing at the split point prevents the free exchange of material between the adjacent convection regions. Although in principle we could resort to the partial mixing mentioned above, in practice it is  \nFigure 6. Profiles of ∇ rad , ∇ ad , and Y as a function of mass coordinate, for the 1 M glyph[circledot] stellar model. The panels correspond to different stages during CHeB: Y c = 0 . 9 (upper), Y c = 0 . 6 (middle), and Y c = 0 . 3 (lower). Gray shading indicates regions undergoing convection.  \n<!-- image -->  \n/circledot  \nnot clear how this might be implemented within a diffusive mixing framework. Constantino et al. (2015) have developed an overshoot-like prescription which appears useful for mimicking the convective neutrality achieved by partial mixing (see their Section 2.3.3), but it in-  \nFigure 7. The mass coordinate m bdy of the convectivecore boundary plotted as a function of CHeB age, for the 3 M glyph[circledot] stellar model discussed in Section 2.4. Different line styles/colors show the separate runs described in the text.  \n<!-- image -->  \nvolves a number of unconstrained parameters. Therefore, on grounds of simplicity and pragmatism - and recognizing that better approaches may become apparent in the future - we modify the predictive scheme to prevent it from causing a convection region to split in the first place. This involves a new control parameter, predictive\\_superad\\_thresh ; if during the predictive mixing iterations the super-adiabaticity ∇ rad / ∇ ad -1 drops below this threshold anywhere in the mixed region, then the code backs off the mixing by one cell and updates the model convective diffusivities and convective velocities in the usual manner.  \nFurther functionality, controlled by a new parameter predictive\\_avoid\\_reversal , also helps to prevent splitting and breathing pulses. When this parameter is set to the name of a MESA isotope, then the code monitors how the predictive mixing alters the abundance evolution of that isotope in the convection region. If it would cause this evolution to reverse (i.e., switch from decreasing to increasing, or vice-versa), then the code backs off the mixing by one cell and updates the model as before. Thus, for instance, setting this parameter to 'he4' during CHeB ensures that the predictive mixing scheme does not cause the core He abundance to increase across a timestep.  \nTo illustrate the preceding discussion, we evolve a 1M glyph[circledot] star through CHeB; this is the same mass considered by Constantino et al. (2015). Figure 5 plots the mass of the convective core as a function of CHeB age (defined as the time elapsed since the central Y drops below 0.98), showing results from separate runs with and without predictive mixing, and  \nFigure 8. Evolution of the 3 M glyph[circledot] stellar model in the HR diagram, from ZAMS to the beginning of CHeB (left), and throughout the CHeB phase (right). Separate tracks show the different cases considered in the text; in the left panel, the Schwarzschild track without predictive mixing lies beneath the tracks with predictive mixing.  \n<!-- image -->  \nusing the Schwarzschild and Ledoux criteria. For the cases with predictive mixing, we adopt a value of 0.005 for the predictive\\_superad\\_thresh parameter, and set predictive\\_avoid\\_reversal to 'he4' to prevent any reversal in the core He abundance. Figure 5 also shows the results from an additional Ledoux/predictive mixing run where we allow the core to split by not setting the predictive\\_avoid\\_reversal and predictive\\_superad\\_thresh controls.  \nFigure 5 shows that without predictive mixing the core is prevented from growing, and the CHeB lifetime significantly curtailed, irrespective of whether the Schwarzschild or Ledoux criteria are used (see also Figure 15 of Paper II). With predictive mixing but no splitting allowed, however, the core grows steadily until He is exhausted, and no breathing pulses are seen. There is almost no difference between the Schwarzschild and Ledoux cases. When the core is allowed to split, the evolution is much noisier. Starting at an age ≈ 25 Myr, the core undergoes episodes of splitting and rejoining that repeat on a short timescale. Toward the end of the evolution, as the core helium abundance becomes very small, the timescale between successive splittings becomes longer, until the core finally splits without rejoining. The overall CHeB lifetime of the model is shortened by ≈ 6 Myr relative to the cases where splitting is avoided.  \nFigure 6 plots the profiles of ∇ rad , ∇ ad and Y for the 1 M glyph[circledot] star at three points during its CHeB evolution, corresponding to core helium mass fractions Y c = 0 . 9, 0.6 and 0.3. The profiles are all from the run with the Ledoux criterion and predictive mixing. In the upper panel, a local minimum in ∇ rad has yet to de-  \nand the core boundary satisfies the ∇ rad = ∇ ad equality on its convective side. In the middle and lower panels, the local minimum in ∇ rad can clearly be seen; in these cases, the predictive mixing has extended the convection region as far as possible without pushing the minimum ∇ rad below the threshold set by the predictive\\_superad\\_thresh control. MESA treats the region between the ∇ rad minimum and the convective boundary as fully convective. On the convective side of this boundary ∇ rad > ∇ ad , which is physically inconsistent but cannot be remedied with predictive mixing alone: any further extension of the boundary would cause the convection region to split. As discussed above, fixing this inconsistency requires some way of modeling the partial mixing expected to occur in the part of the convection region between the ∇ rad minimum and the boundary.  \nThe abundance profiles plotted in Figure 6 show a sharp transition between the He-depleted core and the He-rich radiative region above. Although not shown, the carbon and oxygen abundance profiles exhibit corresponding jumps at the core boundary. Similar results are obtained by Constantino et al. (2015) with their 'maximal overshoot' scheme (cf. their Figure 2); and those authors also find a core mass evolution during CHeB that closely resembles the outcome from predictive mixing (cf. their Figure 8 and our Figure 5). These similarities are not coincidental; although the predictive mixing and maximal overshoot schemes have different narratives and implementations, both have the effect of growing the core boundary during CHeB to the greatest extent permitted without causing the convection region to split. The larger cores that result from this growth  \nappear to provide a better match to Kepler asteroseismic period spacings, when compared with other mixing schemes that produce smaller cores (Constantino et al. 2015); and with certain assumptions about post-CHeB evolution, the larger cores can also provide a satisfactory fit to observational cluster counts (Constantino et al. 2016).  \nTo explore whether the predictive mixing performs equally well for a higher-mass star that has not passed through the He flash, we also evolve a 3M glyph[circledot] star through CHeB; this is the same mass and evolutionary stage considered in Figure 15 of Paper II. Figure 7 plots the mass of the convective core as a function of CHeB age, showing results from separate runs with and without predictive mixing, and using the Schwarzschild and Ledoux criteria. For the cases with predictive mixing, we again adopt a value of 0.005 for the predictive\\_superad\\_thresh parameter, and set predictive\\_avoid\\_reversal to 'he4' to prevent any reversal in the core He abundance. As before, we find that the predictive mixing allows the core to grow steadily; and that the Schwarzschild and Ledoux criteria give essentially the same outcome.  \nAs a visual summary of how predictive mixing influences a star's evolution, Figure 8 plots evolutionary tracks of the 3 M glyph[circledot] model in the Hertzsprung-Russell (HR) diagram, for the same combinations of mixing and stability criteria considered in Figure 7. The left panel focuses on the MS and red giant branch phases, and the right panel on the CHeB phase. In the left panel, the case with the Ledoux criterion but without predictive mixing stands out from the other three as having a slightly reduced luminosity. This behavior arises because the boundary of the hydrogen-burning convective core is incorrectly positioned during the early MS evolution, retarding the growth of the core (the same effect can be seen for the 1 . 5 M glyph[circledot] model in the upper panels of Figure 2). During the subsequent CHeB phase, all four tracks are similar until slightly after the luminosity minimum, when the helium-burning convective core starts to grow; this growth is retarded in both cases without predictive mixing, leading to reduced luminosities and the shorter CHeB lifetimes seen in Figure 7. For the cases with predictive mixing, there is no difference between the Schwarzschild and Ledoux, either on the MS or after.", '2.5. Evolution of the Bottom of the Surface Convective Region in a Low-Mass Star': 'We now evolve a 1 M glyph[circledot] star from ZAMS to TAMS, using the predictive mixing scheme to position the lower boundary of the convective envelope. We include ele-  \nFigure 9. The mass coordinate m bdy of the lower boundary of the envelope convection region plotted as a function of MS age, for the 1 M glyph[circledot] stellar model discussed in Section 2.5. Different line styles/colors show the separate runs described in the text.  \n<!-- image -->  \nment diffusion in these calculations; when it is excluded, the composition remains completely uniform throughout the stellar envelope, and predictive mixing makes no difference whatsoever to the evolution. Figure 9 plots the mass coordinate of the convective boundary as a function of MS age, showing results from separate runs with and without predictive mixing, and using the Schwarzschild and Ledoux criteria.  \nThe four runs are in agreement until an age ≈ 6.5 Gyr; after this point, the downward growth of the region boundary is slower in the run that does not include predictive mixing with the Ledoux criterion. Figure 10 plots the profiles of ∇ rad , ∇ ad , ∇ L and X , in the outer part of the 1 M glyph[circledot] star at an age 8 . 40 Gyr. The left panel illustrates the run with the Ledoux criterion plus predictive mixing, while the right panel shows the run with the Ledoux criterion but without predictive mixing. The former shows that ∇ rad = ∇ ad on the convective (upper) side of the convective boundary, while the latter has ∇ rad > ∇ ad consistent with the boundary growth being retarded.', '2.6. Effect of Timesteps and Mesh Size': 'We now demonstrate how limiting the maximum timestep δt max (set by the max\\_years\\_for\\_timestep control) and changing the mesh resolution parameter ∆ (set by the mesh\\_delta\\_coeff control; see Section B.4 of Paper II for further details) influences the results presented in the previous sections.  \nFirst we consider the effects of changing timestep and resolution on the position of the convective envelope boundary in the 1 M glyph[circledot] model considered in Section 2.5,  \nFigure 10. Profiles of ∇ rad , ∇ ad , ∇ L , and X as a function of mass coordinate, in the outer envelope of the 1 M glyph[circledot] stellar model at an age 8.40 Gyr. The panels show the separate runs described in the text. Gray shading indicates regions undergoing convection.  \n<!-- image -->  \nFigure 11. The mass coordinate m bdy of the lower boundary of the envelope convection region plotted as a function of MS age, for the 1 M glyph[circledot] stellar model. Different line styles/colors show the separate runs with alternative timestep ( δt max ) and mesh resolution (∆) choices. The choices adopted in Section 2.5 are marked with an asterisk [*].  \n<!-- image -->  \nfocusing specifically on the case with the Ledoux criterion and predictive mixing. The results presented previously in Figure 9 are calculated using δt max = 5Myr and ∆ = 0 . 5. Figure 11 demonstrates that halving either δt max or ∆ has little effect on these results, confirming that the calculations are converged. Such settings need to be applied when a converged result is desired from MESA for this calculation.  \nFigures 12 and 13 repeat this exercise for the position of the core convection boundary in the 1 . 5 M glyph[circledot] MS  \nFigure 12. The mass coordinate m bdy of the convectivecore boundary plotted as a function of MS age, for the 1 . 5 M glyph[circledot] stellar model. Different line styles/colors show the separate runs with alternative timestep ( δt max ) and mesh resolution (∆) choices. The choices adopted in Section 2.2 are marked with an asterisk [*].  \n<!-- image -->  \nmodel and 1 M glyph[circledot] CHeB model, respectively. The results presented previously are clearly converged, and this exercise clarifies the MESA settings that should be used for this calculation.', '3. ELEMENT DIFFUSION': 'Section 9 of Paper III describes in detail the old implementation of element diffusion in MESA . Section 9.3.4 points out limitations to those methods, namely: (1) electron degeneracy was not properly accounted for in the diffusion equations, and (2) strong Coulomb interaction introduced theoretical uncertainties for the diffusion  \nFigure 13. The mass coordinate m bdy of the convectivecore boundary plotted as a function of CHeB age, for the 1 M glyph[circledot] stellar model. Different line styles/colors show the separate runs with alternative timestep ( δt max ) and mesh resolution (∆) choices. The choices adopted in Section 2.4 are marked with an asterisk [*].  \n<!-- image -->  \ncoefficients. These two issues are especially important when modeling diffusion in WDs. Here we describe the impact of degeneracy and present new methods to incorporate its effects. We also discuss recent updates to diffusion coefficients and potential approaches for further improvements.', '3.1. Degeneracy and the Approach in Paper III': "The approach to diffusion presented in Section 9 of Paper III assumes all particles obey the ideal gas law. Electron degeneracy pressure can significantly modify the EOS and violate this assumption.  \nFor a plasma species s (i.e., electrons and ions) with partial pressure P s , mass density ρ s , charge density ρ es , number density n s , and temperature T , the Burgers (1969) equations for diffusion are  \nglyph[negationslash]  \nd P s d r + ρ s g -ρ es E = ∑ t = s K st ( w t -w s ) + ∑ t = s K st z st m t r s -m s r t m s + m t , (1)  \nglyph[negationslash]  \nglyph[negationslash]  \nglyph[negationslash]  \nglyph[negationslash]  \n5 2 n s k B d T d r = -2 5 K ss z '' ss r s -5 2 ∑ t = s K st z st m t m s + m t ( w t -w s ) -∑ t = s K st [ 3 m 2 s + m 2 t z ' st ( m s + m t ) 2 + 4 5 m s m t ( m s + m t ) 2 z '' st ] r s + ∑ t = s K st m s m t ( m s + m t ) 2 ( 3 + z ' st -4 5 z '' st ) r t . (2)  \nThe resistance coefficients K st , z st , z ' st , and z '' st are defined in Equation (86) of Paper III. With S representing the total number of plasma species, we must solve for 2 S +2 unknowns: S diffusion velocities ( w s ), S heat flow vectors ( r s ), the electric field ( E ), and the gravitational acceleration ( g ). The Burgers equations above for each species provide 2 S equations, so we can close the system with two additional constraints, which are no net flow of mass or electric current due to diffusion,  \n∑ s ρ s w s = 0 , (3)  \nThis gives a total of 2 S +2 equations.  \nglyph[negationslash]  \n∑ s ρ es w s = 0 . (4)  \nWhen electrons are degenerate, Equation (1) is difficult to apply since d P e / d r no longer takes a simple analytic form. Moreover, the temperature term appearing on the left hand side of Equation (2) clandestinely assumes an ideal gas law. Burgers (1969) defines the temperature for each species as T s ≡ P s /n s k B and assumes thermal equilibrium between all species so that T ≡ T s . The quantities P s and n s are defined in terms of moments of a Maxwellian distribution function, but the Fermi-Dirac distribution function for electrons no longer reduces to a Maxwellian form when they are degenerate, and hence T e = P e /n e k B . If the electrons remain in thermal equilibrium with their surroundings while failing to satisfy an ideal-gas relation for their temperature, the Burgers treatment assigns an incorrect temperature to degenerate electrons for the d T/ d r term in Equation (2).  \nFurthermore, the approach to diffusion described in Paper III follows Thoul et al. (1994) in rearranging and rescaling all equations into one matrix system with units convenient for solving numerically,  \nglyph[negationslash]  \nP K 0 ( α i dln P d r + ν i dln T d r + S ∑ j =1 j =e γ ij dln C j d r ) = 2 S +2 ∑ j =1 ∆ ij W j . (5)  \nThe sum on the left hand side skips the electron index because C e ≡ 1 by construction, and so we save  \nresources by not evaluating its gradient unnecessarily. Here, indices i = 1 , 2 , . . . S encode the S equations given by Equation (1), indices i = S +1 , S +2 , . . . , 2 S encode the S equations given by Equation (2), and indices i = 2 S + 1 , 2 S + 2 encode the 2 constraints of no current or mass flux. For definitions of the various coefficients and matrices in Equation (5), consult Paper III and Thoul et al. (1994). We repeat a few particularly relevant definitions here. First, let C s = n s /n e denote the species concentration, where n e is the electron number density. Second, define the total concentration as C = ∑ s C s . Then the quantity α i appearing in Equation (5) above is defined as  \nThe term α i dln P/ d r in Equation (5) is meant to capture contributions of the driving terms d P s / d r in Equation (1). But this correspondence only holds if the ratio of the partial pressure P s for species s to the total pressure P is given by  \nα i =    C i /C i = 1 , 2 , . . . S, 0 i = S +1 , . . . 2 S +2 . (6)  \nP s P = C s C = n s /n e ∑ t n t /n e = n s n , (7)  \nwhere n is the total number density. This holds as long as all pressures are ideal-gas. However, once electron degeneracy modifies the equation of state, P does not scale linearly with n , and so Equation (7) fails for all species in the plasma. This means the α i term no longer accurately represents the information in the Burgers equations for the diffusion velocity of any species.  \nMoreover, the prefactor P/K 0 in Equation (5) also assumes ideal gas for each species. The quantity K 0 = 1 . 144 × 10 -40 ( T/ 10 7 K) -3 / 2 n 2 e simply scales out some of the information common to all diffusion coefficients in the units used for Equation (5). Thoul et al. (1994) assume an ideal gas to simplify the prefactor in Equation (5) to  \nP K 0 = 2 . 00 ( T/ 10 7 K) 5 / 2 ( ρ/ 100 g cm -3 ) ( ∑ s C s )( ∑ s A s C s ) , (8)  \nwhere A s is the mass of species s in atomic mass units. This scaling was propagated into the MESA diffusion routine described in Paper III. Since ideal gas pressure can be significantly smaller than total pressure when electrons are degenerate, this prefactor for Equation (5) is systematically too small for degenerate plasmas. This can result in diffusion velocities that are many orders of magnitude smaller than obtained by a proper solution.  \nFigure 14. The gravitational acceleration reported by the diffusion routine described in Paper III compared with g Gauss = Gm/r 2 for a 0 . 6 M glyph[circledot] MESA WD model.  \n<!-- image -->  \n-  \n/circledot  \nWe can verify that there are problems in the degenerate regime by looking at the local gravitational acceleration g diff , which is solved for simultaneously with the diffusion velocities in the diffusion routine described by Paper III. MESA also reports the gravitational acceleration independent of the diffusion routine, g Gauss = Gm/r 2 . For a MESA WD model, layers below the surface quickly become degenerate, and the difference between g diff and g Gauss is significant (Figure 14). This reflects the fact that the solutions given by the diffusion routine scale with a pressure that is far too small in the interior.", '3.2. New Methods': 'We now describe new methods that have been introduced to avoid the limitations discussed in Section 3.1.', '3.2.1. Recasting the Burgers Equations': 'The problems with Equation (5) demonstrated in Figure 14 can be circumvented by solving the Burgers equations directly as presented in Equations (1) and (2). When avoiding the rescaling of the Burgers equations that was originally adopted from Thoul et al. (1994), no limitations on the form of total pressure are present.  \nTo that end, we recast the diffusion solver into the form given in Appendix C. This form closely follows the general approach presented by Thoul et al. (1994) for arranging the full set of equations into a single matrix equation, but enters the Burgers equations into that matrix structure without rescaling any quantities. We therefore avoid making any additional ideal-gas assumptions beyond those already present in the Burgers equations.', '3.2.2. Resolving the Degeneracy Problem': 'Electron degeneracy makes it difficult to evaluate the term d P s / d r in Equation (1) in the case of electrons, but it is possible to form a closed set of diffusion equations that makes no explicit reference to this equation for the electrons. Even in many applications involving WDs, each ion species can be treated as approximately ideal, and hence Equation (1) remains useful for ions. We are then left with just two problematic equations out of the system of 2 S + 2 equations: Equations (1) and (2) for the electrons.  \nFor the S -1 species of ions in the system, we can write S -1 Equations (1) in the form  \nglyph[negationslash]  \nglyph[negationslash]  \nn s k B T dln T d r + n s k B T dln n s d r + n s A s m p g -n s ¯ Z s q e E = ∑ t = s K st ( w t -w s ) + ∑ t = s K st z st A t r s -A s r t A s + A t , (9)  \nwhere ¯ Z s is the average charge of species s obtained using Paquette et al. (1986b). Taking this together with S Equations (2) and the two constraints on current and mass flux, we have a total of 2 S + 1 equations. If we drop g as an unknown and treat it as a fixed input to the diffusion routine in MESA using g = Gm/r 2 , we are left with 2 S +1 unknowns. This gives a closed system of diffusion equations with no explicit reference to the problematic Equation (1) for electrons. This is the form of diffusion equations described in Appendix C.  \nThe thermal diffusion terms (those including d T/ d r in Equation 2) still contain ideal-gas assumptions as described in Section 3.1. Fortunately, in WD cores where strong electron degeneracy occurs, electron conduction leads to efficient thermal transport, resulting in small temperature gradients. With d T/ d r glyph[lessmuch] T/H , where H = P/ρg is the local scale height, the heat flow vectors (representing kinetic energy carried along a temperature gradient by diffusing particles) become negligible: r s glyph[lessmuch] w t for all w t . Thus for WD interiors the system of diffusion equations can be simplified by dropping the S heat flow terms, removing the need for the S Equations (2). Indeed, according to Iben & MacDonald (1985) and Paquette et al. (1986b), thermal diffusion leads only to small corrections to the diffusion velocities for degenerate WD interiors.  \nTherefore, following Iben & MacDonald (1985), we provide options for neglecting thermal diffusion in electron degenerate regions, setting r s = 0 and dropping Equation (2) for each species. Equation (9) then simplifies to the following S -1 equations that no longer  \ndepend on r s for the ions:  \n1 n s ∑ t K st ( w t -w s ) + ¯ Z s q e E = A s m p g + k B T dln T d r + k B T dln n s d r , (10)  \nwhich matches Equation (10) from Iben & MacDonald (1985). Together with the 2 constraints, this leaves a simplified set of S +1 equations for S +1 unknowns: S diffusion velocities w s and the electric field E .  \nThermal diffusion terms tend to enhance gravitational settling velocities (Iben et al. 1992). This can be seen in Figure 15 for a 1 . 25 M glyph[circledot] star on the MS, where the solvers that include thermal diffusion speed the sedimentation of 16 O away from the surface relative to the solver that neglects thermal diffusion. MESA also provides options for smoothly transitioning between diffusion velocities obtained with and without thermal diffusion (averaging between the two solutions in a blending region as a function of electron degeneracy parameter). By default, this transition region occurs when the electron chemical potential is near µ e ∼ k B T , but it is left to the user to decide on an appropriate range of electron degeneracy over which thermal diffusion should be shut off, if at all. The effect of blending between solvers with and without thermal diffusion is to suppress the thermal enhancements to diffusion velocities, smoothly pushing the enhancements to zero as electrons reach a degeneracy threshold. The implementation for the simplified set of diffusion Equations (10) and the smooth turn-off of thermal diffusion terms as a function of degeneracy are described in Appendix C.  \nIn order to confirm that we recover the correct behavior on the MS, we compare results obtained with different diffusion routines for a 1 . 25 M glyph[circledot] star in Figure 15. Here the results based on Thoul et al. (1994) are valid, since no significant departures from ideal-gas behavior are present near the surface. The results obtained with the new scheme are in agreement.', '3.2.3. Diffusive Equilibrium': 'Paper II and Paper III show abundance profiles for WDs that have reached diffusive equilibrium in their outer layers. Figure 23 of Paper II compares the diffusive tails of H and He to an analytic expression from Althaus et al. (2003) and finds good agreement. However, Althaus et al. (2003) note that their analytic expression for diffusive equilibrium follows Arcoragi & Fontaine (1980) in assuming an ideal gas, and the equilibrium abundance profiles from their evolutionary models deviate from the analytic expression due to the inclusion of electron degeneracy. Similarly, the He layer of the WD  \nFigure 15. Surface 16 O mass fraction of a 1 . 25 M glyph[circledot] star over its MS lifetime. It first decreases as diffusion causes sedimentation. Then it increases after the small surface convection zone begins to grow, catching the receding 16 O and mixing it back toward the surface.  \n<!-- image -->  \nmodel shown in Figure 43 of Paper III is partially degenerate, and hence the driving forces for diffusion should be modified in this region.  \nFor a fully-ionized isothermal ideal gas the electric field that serves as one of the driving forces for diffusion in Equation (9) takes the form q e E = [ A/ ( Z +1)] m p g . In contrast, in the limit of strong electron degeneracy, the electric field approaches q e E = ( A/Z ) m p g . When He is the background material, the electric-to-gravitational force ratio q e E/m p g increases from 4 / 3 to 2. In this limit, any trace isotopes with A/Z = 2 see no net sedimentation force ( Zq e E -Am p g = 0), while H with A/Z = 1 sees a significant upward sedimentation force ( Zq e E -Am p g > 0). This extra buoyant force on H in a degenerate He background pushes the diffusive tail further toward the surface relative to the ideal-gas case, as shown in Figure 16. With the proper handling of electron degeneracy described in Section 3.2, our MESA models now agree with the time-dependent diffusion models shown in Figure 18 of Althaus et al. (2003).', '3.2.4. Radiative Levitation': 'Radiative levitation is included as an optional extra term. The Burgers equations are modified with a extra forcing term by taking ρ s g → ρ s ( g -g rad ,s ), as shown in Equation (99) of Paper III. Our implementation continues to follow Hu et al. (2011) but no longer employs their matrix structure for the Burgers equations; details of how the g rad ,s terms are handled with the updated diffusion schemes can be found in Appendix C.  \nFigure 16. Abundance profiles in 0 . 6 M glyph[circledot] MESA WD models at T eff = 5 , 000 K after evolving for 4 Gyr to approach diffusive equilibrium in the outer layers. The old equations assume an ideal gas; the new equations include the effects of electron degeneracy.  \n<!-- image -->  \n-  \n/circledot', '3.3. Updated Diffusion Coefficients': 'The Paquette et al. (1986a) diffusion coefficients have served as the standard for stellar diffusion problems. The scattering cross-sections for these coefficients are calculated using a screened Coulomb potential  \nV 12 ( r ) = ¯ Z 1 ¯ Z 2 q 2 e r exp( -λ/r ) , (11)  \nwith the screening length chosen as λ = max( λ D , ¯ a i ), where λ D is the Debye length, ¯ a i = (3 / 4 πn i ) 1 / 3 is the average interionic distance, and n i is the ion density. This choice is a crude but effective way to handle the strongly coupled regime; as shown in Paper III, this yields reasonable agreement with diffusion coefficients calculated from molecular dynamics.  \nStanton & Murillo (2016) provide updated calculations of collision integrals for screened Coulomb interactions and suggest improvements to the treatment of screening length. They provide fitting functions and tables that can be used with any choice of screening length. In MESA we follow their suggested screening prescription. The electron screening length is given by a Thomas-  \nFermi approximation that accounts for non-relativistic degeneracy:  \nλ e =   4 πq 2 e n e √ ( k B T ) 2 + ( 2 3 E F ) 2   -1 / 2 , (12)  \nλ i = ( 4 π ¯ Z 2 i q 2 e n i k B T ) -1 / 2 . (13)  \nwhere E F = glyph[planckover2pi1] 2 (3 π 2 n e ) 2 / 3 / 2 m e is the electron Fermi energy. The direct inclusion of degeneracy increases λ e . The ion screening lengths are the Debye lengths for each species,  \nTo prevent ions from screening below the inter-ionic spacing, Stanton & Murillo (2016) introduce an approximate ion-sphere for each species a i ≡ (3 ¯ Z i / 4 πn e ) 1 / 3 , and define an ion-sphere coupling parameter  \nΓ i ≡ ( ¯ Z i q e ) 2 a i k B T . (14)  \nTheir net effective screening length is then  \nλ eff ≡ [ 1 λ 2 e + ∑ i 1 λ 2 i ( 1 1 + 3Γ i ) ] -1 / 2 . (15)  \nThis construction enforces a minimum on the screening length at approximately the ion-sphere radius a i for each species, similar to the strict minimum at ¯ a i set by Paquette et al. (1986a). Stanton & Murillo (2016) point out that this adjustment to the ion screening length is physically motivated by the ion pair distribution functions in a strongly coupled plasma, where the occupation probability within the ion-sphere radius is negligible, and hence no ions are present to provide screening beneath that cutoff. The proper handling of degeneracy in the electron screening length makes it unnecessary to impose any particular minimum there, so there is no longer any ad hoc appeal to a universal minimum screening length.  \nFor repulsive Coulomb potentials of the form in Equation (11), Stanton & Murillo (2016) provide fits and tables of collision integrals and coefficients that we now use to calculate the resistance coefficients K st for inclusion in the Burgers equations in MESA . They do not provide fits for attractive potentials, and Paquette et al. (1986a) note that interactions with these potentials behave significantly differently from those with repulsive potentials when screened. Hence, MESA continues to use the Paquette et al. (1986a) coefficients for electron-ion terms, and adopts Stanton & Murillo (2016) for all ion-ion coefficients. In any case, it is evident from Equation (94) in Paper III that the resistance coefficients approximately  \nfollow K st ∝ µ 1 / 2 st , where µ st is the reduced mass of particles s and t ; so, electron-ion resistance coefficients are generally negligible compared to the ion-ion terms.  \nThe calculations of Paquette et al. (1986a) overestimate the electron-ion resistance coefficients in the case where electrons are degenerate. This is because diffusion and resistance coefficients are generally calculated assuming that the velocity distributions of all particles are Maxwellian, and the coefficients roughly scale as K st ∝ v -2 s v -2 t . When the electrons become degenerate, their characteristic kinetic energies are of order E F glyph[greatermuch] k B T , and so their velocity distribution skews toward larger velocities. This results in smaller resistance coefficients K st , overestimating the impact of electronion drag. However, the overestimate results in coefficients that remain negligible compared to ion-ion terms, and no attempt is made to correct it in MESA .  \nFor repulsive potentials, the coefficients from Stanton & Murillo (2016) generally agree with those of Paquette et al. (1986a) to within a few percent. In strongly coupled WD interiors the Stanton & Murillo (2016) coefficients lead to ∼ 10% shorter diffusion timescales due to a screening length that is allowed to be somewhat smaller than the minimum value imposed by Paquette et al. (1986a): λ eff < ¯ a i . Future prospects for further improvements to diffusion coefficients include the recent progress on effective potential methods from Daligault et al. (2016) and Shaffer et al. (2017).', '3.4. Diffusion-Induced Flashes on He WDs': "Diffusion-induced H shell flashes on low-mass ( M glyph[lessorsimilar] 0 . 4 M glyph[circledot] ) He WDs are known to alter their cooling times (Althaus & Benvenuto 2000; Althaus et al. 2001) and seismic properties (Althaus et al. 2013). Istrate et al. (2016a,b) use MESA to model this process, generating tables of cooling timescales and comparing MESA models with those of Althaus et al. (2013).  \nFigure 17 shows an exploration of the H shell flash domain for a large grid of Z = 0 . 02 MESA models over a range of He-core and H-envelope masses. Here the envelope mass is defined as the total mass of H-rich material ( X > 0 . 01) at the surface at the beginning of the WD cooling track. Lines show the minimum envelope masses for which H shell flashes occur given various diffusion prescriptions.  \nFor a given core mass, there is a range of envelope masses that exhibit shell flashes only if diffusion is included, but this range depends on the diffusion prescription. The two lower lines for models including diffusion in Figure 17 differ only in the handling of electron degeneracy in the diffusion scheme. This illustrates the importance of properly handling degeneracy as described in  \nFigure 18 shows the delay in WD cooling from introducing glyph[epsilon1] 22 into 0 . 6 M glyph[circledot] models. These models turn off diffusion for Γ > 175, so glyph[epsilon1] 22 is only active in material for which crystallization has not yet occurred. The timedelays shown in Figure 18 are in good agreement with those shown by Deloye & Bildsten (2002) and Garc'ıaBerro et al. (2008) for comparable cases.  \n<!-- image -->  \n/circledot  \nFigure 17. Minimum envelope mass M env for which a H shell flash occurs on a He WD core mass M core for Z = 0 . 02 MESA models with and without diffusion. The regime for the phenomenon of diffusion-induced flashes lies between the boundaries for models with and without diffusion.  \nSection 3.2, since the diffusion-induced flashes are typically ignited by CNO burning in the diffusive tail of H that reaches into the partially degenerate He layers. WDs in this mass range often experience cycles of many H flashes, depleting H incrementally until insufficient H remains to ignite another flash. The disagreement between diffusion prescriptions on the minimum envelope mass for flashes is therefore significant, as this will determine the total number of flashes and final H mass that sets the ultimate cooling timescale for an object.  \nTo explore the full range of parameters presented in Figure 17, our WD models were built by artificially stripping the H envelope down to a specific mass coordinate above the He core of a 1 . 0 M glyph[circledot] model ascending the RGB. For a discussion of MESA models including protoWD formation and the resulting H envelope masses, see Istrate et al. (2016b).", '3.5. Heating from 22 Ne Settling': "glyph[negationslash]  \nIn the strongly degenerate limit, q e E/m p g ≈ 2 for C/O WD cores. For an isotope where A/Z = 2, the electric and gravitational fields result in a net force that drives diffusion. For 22 Ne in cooling WD interiors, this force is F = Zq e E -Am p g ≈ -2 m p g , causing 22 Ne to sediment toward the center and deposit energy as it moves deeper into the gravitational potential (Bildsten & Hall 2001; Deloye & Bildsten 2002; Garc'ıa-Berro et al. 2008, 2010). This heating can prolong the WD cooling timescale, especially at late times when the WD is very dim and radiates away the energy slowly. This effect may be especially important for explaining WD  \n<!-- image -->  \n/circledot  \nFigure 18. Extra cooling time required to reach a given luminosity for 0 . 6 M glyph[circledot] WD models including heating from 22 Ne settling, relative to models neglecting this heating. For comparison, we also show a result from Garc'ıa-Berro et al. (2008) for a 0 . 6 M glyph[circledot] WD with an Oxygen-dominated core composition. Figure 57 shows the same quantity including other physical processes such as crystallization for the same Z = 0 . 02 WD model shown here.  \nluminosity functions in old and metal-rich open clusters such as NGC 6791, where abundant 22 Ne is available in WD interiors to provide heating.  \nMESA now offers an option to include this heating term in the energy equation (see Section 8.7) when diffusion is enabled. The specific rate at which energy is deposited is  \nglyph[epsilon1] 22 = | F | v 22 ( Am p ) /X 22 = (22 m p g -10 q e E ) X 22 v 22 22 m p . (16)  \nThe 22 Ne diffusion velocity ( v 22 ) and electric field are calculated in the diffusion routine and then used to evaluate the above heating term. Note that the updates to diffusion described in Section 3.2 are essential for correctly calculating both the diffusion velocity and magnitude of the driving force in the degenerate interior of the WD.", '4. IMPLICIT HYDRODYNAMICS': 'In Paper III we describe implicit shock-capturing hydrodynamics capabilities based on the use of an artificial viscosity. We now add an option for using an approximate Riemann solver, the HLLC (Harten-Lax-van Leer-  \nContact) solver introduced by Toro et al. (1994). (See also Batten et al. 1997 for an early implicit implementation of HLLC.) The HLLC method provides improved shock capturing and energy conservation by avoiding the need for artificial viscosity. However, the methods presented in Paper III are still included in MESA so that users may continue to apply them.', '4.1. Implementation of HLLC': "Accurate shock-capturing methods evaluate the flux of hydrodynamical conserved quantities by extrapolating the solution on each interface between zones over the course of the timestep. The different methods for projecting the solution into the future are known as different 'Riemann solvers'. HLLC is designed to accurately capture the evolution of contact discontinuities. When implemented on a Lagrangian grid, HLLC is able to evolve purely advective flows without any contact smearing (Cheng & Shu 2007; Duffell & MacFadyen 2011; Cheng et al. 2012; Cheng & Shu 2014).  \nPaper I and Paper III discussed the evolution of a velocity variable v , defined at cell faces. When using HLLC, MESA instead evolves a cell-centered velocity u .  \nWe solve a Riemann problem at the cell face with index k . The cell to the left (toward the center) is the cell with index k ; the cell to the right (toward the surface) is the cell with index k -1. The cell face radius is r k . The mass contained within an individual cell is dm k . The mass enclosed from the center of the star to the cell face is m k . For convenience, we define the face area as A k = 4 πr 2 k . Thermodynamic variables (e.g., P k , ρ k ) are defined at cell centers by mass. Figure 19 shows the layout of cells.  \nFigure 19. Cell and face variables relevant for hydrodynamics in MESA when using HLLC.  \n<!-- image -->  \nThe calculation begins by making estimates for the density and velocity at the left and right of the face. Explicit codes sometimes use multipoint polynomial interpolation based on values in neighboring cells to improve the reconstruction of the values at the face. However for an implicit code such as MESA , that would introduce dependencies in the partial derivatives for the Jacobian that would violate the necessary block tridiagonal structure (see Appendix B in Paper II). To avoid this problem, we use the cell center density and velocity alone to estimate the values at the edges of that cell. The variables for the left and right values are named relative to the edge rather than the cell, that is  \nρ L = ρ k ρ R = ρ k -1 u L = u k u R = u k -1 . (17)  \nThis choice limits the solution to be first-order accurate in space.  \nUsing an approach similar to Kappeli & Mishra (2014), we reconstruct the pressure at the faces assuming hydrostatic equilibrium. The pressure derivative implied by hydrostatic equilibrium at the face is  \n( d P d m ) hse = -Gm k A k r 2 k , (18)  \nand we reconstruct the pressure to the left and right of the face  \nP L = P k + dm k 2 ( d P d m ) hse , P R = P k -1 -dm k -1 2 ( d P d m ) hse . (19)  \nThis choice improves the timescale over which hydrostatic equilibrium can be maintained when using HLLC, and facilitates the process of switching from a hydrostatic model to one in which a velocity variable is evolved.  \nThe 1D Lagrangian context makes the implementation of HLLC straightforward. In a Lagrangian code like MESA there is no mass flux across cell faces. In hydrodynamics, there is no mass flux across a contact discontinuity. HLLC includes the contact wave, so we simply associate the contact wave with the cell face. 1 As given by Toro (2009), the HLLC estimate of the contact wave speed is  \nS ∗ = u R ρ R ( S R -u R ) + u L ρ L ( u L -S L ) + ( P L -P R ) ρ R ( S R -u R ) + ρ L ( u L -S L ) , (20)  \n1 In Section 5, where we consider the effects of mass diffusion, we will need to slightly revise this association.  \nand hence u face = S ∗ . Likewise, the pressure at the cell face is the pressure at the contact wave, P face = P ∗ , where  \nP ∗ = 1 2 [( ρ R ( u R -S R )( u R -S ∗ ) + P R ) + ( ρ L ( u L -S L )( u L -S ∗ ) + P L )] , (21)  \nand S L and S R are the fastest wavespeeds moving to the left and right, respectively. To evaluate these, we assume the simple and most conservative bounds on the signal velocities,  \nS L = min( u L -c s , L , u R -c s , R ) , S R = max( u L + c s , L , u R + c s , R ) , (22)  \nwhere c s , L and c s , R are the sound speeds on the left and right sides of the cell boundary, respectively. Having obtained values for u face and P face , we now formulate the versions of the equations used when HLLC is enabled.  \nIn the Lagrangian picture, the cell boundaries move with the fluid velocity, such that the net fluxes for mass, momentum, and energy from cell k to cell k -1 are extremely simple (Cheng & Shu 2014) and given by  \nF ρ,k = 0 , F p,k = A k P face ,k , F e,k = A k P face ,k u face ,k + L k . (23)  \nThe L k term in the energy flux does not come from the solution of the Riemann problem, but from the fact that MESA also evolves a luminosity variable that reflects the radiative or convective transport of energy.  \nThe finite volume form of the mass conservation equation remains the same as that given in Paper I,  \nln r k = 1 3 ln [ r 3 k +1 + 3 4 π dm k ρ k ] . (24)  \nHowever, the equation for the radius (Equation 29 in Paper III) has changed. The new equation for the radius is  \nr k = r start ,k + u face ,k δt , (25)  \nwhere δt is the timestep. For numerical precision, we re-write this as  \nexp(ln r k -ln r start ,k ) -1 ︸ ︷︷ ︸ expm = u face ,k δt r start ,k , (26)  \nwhere this recasting allows use of crlibm (de Dinechin et al. 2007, see also Paper III) function expm to evaluate the function exp( x ) -1 to machine precision (as indicated by the underbrace).  \nThe local radial momentum equation for cell k is  \nu k -u start ,k δt = -1 2 ( Gm k r 2 k + Gm k +1 r 2 k +1 ) + P k dm k ( A k -A k +1 ) + 1 dm k ( F p,k +1 -F p,k ) . (27)  \nOn the right hand side, the first term is gravitational, the second is a geometric source term that arises from putting the equation in conservation-law form, and the final term is the momentum flux in and out of the cell found by HLLC.  \nThe local total energy conservation equation for cell k is  \ne k -e start ,k + 1 2 ( u 2 k -u 2 start ,k ) -Gm C ( 1 r C -1 r C , start ) = δt [ 1 dm k ( F e,k +1 -F e,k ) + glyph[epsilon1] nuc ,k -glyph[epsilon1] ν,k + glyph[epsilon1] extra ,k ] . (28)  \n(See Section 8.3 for a discussion of how this energy accounting is related to that typically used in stellar evolution calculations.) We define the cell center quantities m C and r C to be r and m at the center of mass of the cell. The terms on the left split the local total energy into internal, kinetic, and potential components. The right side gives the energy in and out of the cell and the energy sources and sinks in the cell. Energy loss from neutrinos due to nuclear reactions is already subtracted from the nuclear burning term, glyph[epsilon1] nuc , so only the neutrino energy loss rate from thermal processes, glyph[epsilon1] ν , is explicitly accounted for in Equation (28). Other processes are accounted for via glyph[epsilon1] extra .  \nAs in Paper I, the temperature differences of interior cells T k are set by energy transport across boundaries,  \nT k -1 -T k = dm k [ ∇ T,k ( d P d m ) hse T k P k ] , (29)  \nwhere ∇ T,k is provided by MESA module mlt (see Section 5.1 in Paper I) and the overbars indicate quantities at the cell faces (see Figure 19). This equation relates temperatures of neighboring cells; the actual temperature in each cell is then fixed by a surface boundary condition.  \nMESA 's HLLC includes the effects of rotation in the shellular approximation (see Paper II, Section 6.1) and can also include a post-Newtonian correction to the gravitational force. (For an example application to neutron star envelopes, see Paper III, Section 5.3). These capabilities require modifications to the momentum equation. In both cases, they can be treated as a  \nrescaling of the local gravitational constant G → fG . In the case of rotation, the rescaling factor is f P (Paper II, Equation 23). In the post-Newtonian case, it is (1 -2 Gm/ ( rc 2 )) -1 / 2 . Therefore, when either of these is used with the hydrodynamics capabilities described in this section, the rescaling is applied to the G in pressure reconstruction (Equation 18) and separately to each G (for cell k and k +1) in the momentum equation (Equation 27).", '4.2. Mesh Refinement': "During a typical stellar evolution run, MESA controls its meshing using 'mesh functions' that limit the maximum allowed change of various quantities between adjacent cells (see Section 6.5 in Paper I and Section B.4 in Paper II). With HLLC, the criteria to split or merge cells are written solely in terms of the radial coordinate in order to simplify the adjustments to the mesh in response to large changes in density before and after a shock. Cells split when they decompress enough that their radial extent becomes too large, and they merge with a neighbor when they compress enough that their radial extent becomes too small.  \nThe refinement criteria can use either linear ( x = r ) or logarithmic ( x =ln r ) radius. The user selects a target number of cells, N target . MESA translates this into a target cell size, dx target =( x surface -x center ) /N target . A cell is considered too large if dx k /dx target > f long and a cell is considered too small if dx target /dx k > f short . 2 The refinement then proceeds iteratively. At each iteration, MESA selects the smallest and largest cells. If the largest cell is too large, it is split. If the smallest cell is too small, it is merged unless the magnitude of the difference between its velocity and that of either neighbor is a significant fraction of the local sound speed: this prevents merging in the immediate vicinity of the shock where there are sharp jumps in velocity. The refinement proceeds up to some maximum number of iterations, though in practice the procedure typically stops before then because no more cells satisfy the criteria to be split or merged.  \nA cell merges with its smaller neighbor, unless they have a different most-abundant chemical species, in which case the cell merges with the other neighbor instead. When a cell is split, differences in quantities such as density and chemical abundances between the two child cells are determined by interpolation from the neighboring cells.", '4.3. Time Resolution': 'Since the hydrodynamics equations are being solved implicitly, MESA is not subject to the Courant-FriedrichsLewy (CFL) timestep condition for numerical stability. The size of the MESA timestep is instead limited by the restrictions on the allowed changes in the structure of the star. The usual timestep controls continue to apply.  \nWhile numerical stability does not require the restrictive CFL timestep condition, the choice of timestep does affect the accuracy of the solution. A CFL-like limit is often also applied because it can be a convenient additional way to restrict timesteps along with the other options. Such a restriction allows for well-converged solutions. The timestep can be limited by the requirement that  \nwhere f t is a user parameter. Unlike in an explicit code where a similar minimum must be evaluated over all cells, in MESA the minimum is taken only over cells for which  \nmin { dr k | u k | + c s , k } < f t × δt , (30)  \nmax {| u k -u k +1 | , | u k -u k -1 |} c s , max < f u , (31)  \nwhere c s , max is the maximum evaluated over nearby cells, and f u is a user parameter. 3 This means that only regions near the shock front limit the timestep. The option for additional limitations on where this condition is evaluated (e.g., in mass) are provided.', '4.4. Hydrodynamic Test Problems': 'In order to test the HLLC implementation, comparisons are now made to problems with known solutions.', '4.4.1. Sedov Blast Wave': 'In the Sedov blast wave problem, an energy E blast is deposited at the origin at time zero in a domain with a non-uniform density profile ρ = ρ 0 r -ω , where ρ 0 and ω are constants. We assume an ideal-gas EOS with a constant adiabatic index γ , that is P = ( γ -1) ρe .  \nGeneration of numerical Sedov solutions is discussed in Kamm & Timmes (2007). A constant initial density profile, ω =0, is frequently used in verification tests (e.g., Gehmeyr & Mihalas 1994; Fryxell et al. 2000). Although a power-law initial density profile is more challenging for verification studies, we explore such a profile because density gradients are prevalent in astrophysics.  \nFigure 20. Sedov blast wave density, pressure and velocity profiles at the labeled times. The analytic (black) and MESA (colored) curves show solutions for a shock propagating down a ρ = ρ 0 r -1 density profile with an adiabatic index γ =7/5. Deviations from the analytic solutions are glyph[lessorsimilar] 2%.  \n<!-- image -->  \nTo model a shock propagating down a linear density gradient in spherical geometry we set ω =1, γ =7/5, ρ 0 =1gcm -3 , and P 0 =0ergcm -3 in the analytic solution; while we set P 0 =10 -6 erg cm -3 in MESA as a stable numerical approximation to zero pressure. The initial blast energy, E blast =1.464 erg, is determined by choosing that r shock =1 cm at t =1 s and then calculating the Sedov energy integral. Figure 20 shows the evolution of the density, pressure, and velocity. As the shock propagates outward from the origin these quantities monotonically decrease as mass is swept up by the shock. The spherical Sedov problem admits a similarity solution. Figure 21 demonstrates that MESA maintains the analytic self-similar profiles at different times.', '4.4.2. Noh Problem': "Noh (1987) describes a standard verification problem that tests the ability to transform kinetic energy into internal energy, and the ability to follow supersonic flows. A sphere of cold gas with an ideal-gas EOS and constant adiabatic index γ , that is P = ( γ -1) ρe , is initialized  \nFigure 21. Sedov blast wave self-similarity of the analytic (black curves) and MESA (colored symbols) solutions. Scaled velocity v/v shock , pressure P/P shock , and density ρ/ρ shock profiles for a shock propagating down a ρ = ρ 0 r -1 density profile at the three different times are overlaid. Symbols for each epoch mark cell locations. Deviations from the analytic self-similar solutions are glyph[lessorsimilar] 2%.  \n<!-- image -->  \n~  \nFigure 22. Analytic (colored solid) and MESA (black dashed) solutions for the density, pressure, and velocity at t =0.3s in the Noh problem. Disagreements near the center are due to wall-heating as discussed in the text.  \n<!-- image -->  \nwith a uniform, radially inward speed of 1 cm s -1 . A shock forms at the origin and propagates outward as the gas stagnates. For an initial gas density ρ 0 =1gcm -3 , the analytic solution in spherical geometry for γ =5/3 predicts a density 64 g cm -3 in the stagnated gas.  \nFigure 22 shows the analytic and MESA profiles for the density, pressure, and material speed at t =0.3 s. Most implementations, including MESA 's, produce anomalous 'wall-heating' near the origin (although see Gehmeyr  \net al. 1997). As the shock forms at the origin the momentum equation tries to establish the correct pressure level. However, numerical dissipation generates additional entropy. The density near the origin drops below the correct value to compensate for the excess internal energy (e.g., Noh 1987; Rider 2000). Thus, the density profile is altered near the origin while the pressure profile remains at the correct constant level in the post-shock region.  \nFigure 22 shows the analytic solution and MESA solution for a fixed timestep of δt =5 × 10 -6 s and 10,000 cells. Deviations from the analytic solutions are glyph[lessorsimilar] 1%, except for the density near the origin and the shock front. A convergence exercise with different fixed timesteps and spatial resolutions suggests spatial resolution is relatively more important in the MESA solutions than temporal resolution for the Noh problem.", '4.4.3. Supernova Shock': 'The problem of a supernova (SN) shock moving through a stellar envelope has been extensively studied. For a radiation-dominated strong shock, a simple analytic expression for the shock velocity is provided by Matzner & McKee (1999),  \nv sh = α ( E m -M center ) 1 / 2 [ m -M center ρr 3 ] 0 . 19 , (32)  \nwhere we adopt α = 0 . 736 as suggested by Tan et al. (2001). The explosion has an energy E . The mass that enters into this expression is the mass entrained by the shock and so differs from the Lagrangian mass coordinate ( m ) by the mass of the remnant ( M center ). Since the material in the shocked envelope has an adiabatic index of 4 / 3, the Matzner & McKee (1999) prediction for the post-shock velocity is v MM = 6 v sh / 7.  \nMESA defines the shock location to be the outermost point where the fluid Mach number exceeds 1, as measured in the rest frame of the star. Since the primary application of these capabilities are blast waves propagating into approximately static stellar envelopes, this shock detection criterion suffices. Figure 23 compares the velocity in a MESA model (the 19 M glyph[circledot] model of SN1999em; see Section 6) with v MM . We show explosions with two different energies, E = 0 . 9 × 10 51 erg and E = 2 . 7 × 10 51 erg. Both cases have M center = 1 . 5 M glyph[circledot] . Typical differences are at the few percent level.', '4.4.4. Weak Shock Propagation': "We now explore weak shocks with Mach numbers M = 1 . 2 -2 . 2 propagating outward in the hot stellar envelope of a classical nova progenitor. The model is a 0 . 8 M glyph[circledot] WD. The H/He envelope extends from  \nFigure 23. Comparison of MESA with the analytic results of Matzner & McKee (1999) and Tan et al. (2001). This started as a 19 M glyph[circledot] model at ZAMS; at explosion it is 17.79 M glyph[circledot] with M center = 1 . 5 M glyph[circledot] . The upper panel shows an explosion with E = 0 . 9 × 10 51 erg; the lower panel shows E = 2 . 7 × 10 51 erg. The gray dashed curve shows the analytic prediction for the post-shock fluid velocity given the density profile of the initial model. The solid curves show velocity profiles from the MESA calculation at specific times. The unfilled diamonds indicate where on the dashed curve the two should be compared.  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \nr = 7 . 1 × 10 8 cm to r = 7 . 85 × 10 8 cm with densities ρ = 10 -100 g cm -3 and temperatures T ≈ 10 7 K.  \nAfter excising the core we run the model with HLLC enabled for 100 s, corresponding to ≈ 50 sound crossing times in the outer envelope, to allow the envelope to settle. Afterwards, the envelope has a total energy of -9 . 16 × 10 45 erg, with -9 . 38 × 10 45 erg in potential energy, 2 . 22 × 10 44 erg in thermal energy, and a negligible kinetic energy 1 . 2 × 10 29 erg. We turn off convective energy transport to study the properties of weak shocks. To create weak shocks, we inject 0 . 5% -5% of the total thermal energy into a single cell with mass dm = 1 . 6 × 10 25 g at r = 7 . 3475 × 10 8 cm over 10 -4 s. Fig-  \nFigure 24. Velocity at different times in the envelope after 6 × 10 42 ergs has been injected. One shock front travels upward and grows as it enters the outer atmosphere, and another pulse travels downward and reflects off of the inner boundary. The thin gray lines denote the region of study.  \n<!-- image -->  \nhows the resulting upward and downward propagating shocks. We restrict our region of study to the region where the upward and downward shocks are wellseparated, in the radius range of r = 7 . 4 -7 . 65 × 10 8 cm (denoted with thin gray lines in Figure 24). We do not study the properties of the downward shock and its artificial reflection from the 'floor' of our model.  \nWe define u 0 = d r peak / d t as the shock velocity, where r peak is defined as the radial location with maximum fluid velocity. We compare the properties of the shock to analytic expectations for cases where γ is identical in the pre- and post-shock material. Pre-shocked quantities carry a 0 subscript, and shocked quantities carry a 1, and we use the sound speeds, c s , and pressures, P , on either side of the discontinuity. Following Zel'dovich & Raizer (1967), in the rest frame of the shock front, the pre-shock gas travels into the shock front at velocity  \n( u 0 c s , 0 ) 2 = ( γ -1) + ( γ +1) P 1 /P 0 2 γ . (33)  \nThe post-shock velocity u 1 has magnitude | u 1 | = | u 0 -u peak | , where u peak denotes the fluid velocity u at r peak . The analytic expression is  \n( u 1 c s , 1 ) 2 = ( γ -1) + ( γ +1) P 0 /P 1 2 γ . (34)  \nLocal shocked quantities are evaluated at the cell with the maximum Lagrangian fluid velocity, while preshocked quantities are evaluated at the cell in the initial MESA profile (before the shock has propagated) with the same mass coordinate as the shock front when it reaches r peak . The thin black lines in the upper and lower panels  \nFigure 25. Comparison of the MESA calculation (colored lines) with analytic expressions (thin black lines) for ( u 0 /c s , 0 ) 2 (upper) and ( u 1 /c s , 1 ) 2 (lower) for different energies injected.  \n<!-- image -->  \nof Figure 25 are the right-hand side of Equations (33) and (34), respectively, for shocks produced by different amounts of injected energy. Colored lines show the lefthand side of each equation as calculated from the MESA model.  \nWe now compare the temperatures T 0 and T 1 of the pre- and post-shock gas. We expect  \nT 1 T 0 = P 1 P 0 [ ( γ -1)( P 1 /P 0 ) + ( γ +1) ( γ -1)( P 1 /P 0 ) + ( γ -1) ] . (35)  \nThe thin black lines in the upper panel of Figure 26 show the right hand side of Equation (36) and the solid colored lines correspond to quantities calculated by MESA . The colored dotted lines in the top plot show the temperature change for an adiabatic compression,  \n( T 1 T 0 ) S =const = ( P 1 P 0 ) (1 -1 /γ ) , (36)  \nmaking it clear that for the weakest shocks, the temperature jump is that expected from an adiabatic compression. However, for stronger shocks, the temperature is higher due to the entropy increase associated with the shock. For a gas with specific heat capacity c V , this entropy jump is  \nS 1 -S 0 = c V ln { P 1 P 0 [ ( γ -1)( P 1 /P 0 ) + ( γ +1) ( γ -1)( P 1 /P 0 ) + ( γ -1) ] γ } , (37)  \nFigure 26. Comparison of the MESA calculation (colored lines) to the expectations of shock theory (thin black lines) for the temperature increase (upper) and entropy increase (lower) for different energies injected. Colored dotted lines in the upper plot indicate the temperature change for a purely adiabatic compression.  \n<!-- image -->  \nshown by the thin black lines in the bottom plot of Figure 26. The colored lines correspond to quantities calculated by MESA . The agreement is excellent for large entropy jumps, but becomes noisy at lower injection energies because extracting small changes from the background is then challenging. For the weakest shock, the entropy changes are orders of magnitude smaller than the background entropy 1 . 2 -1 . 6 × 10 9 erg g -1 K -1 in the region of interest.", '5. RAYLEIGH-TAYLOR INSTABILITIES': 'The outward moving shock in a core collapse SN explosion encounters multiple composition boundaries. Across these boundaries the density gradient is steep, especially at the H/He boundary. Post-shock, these regions become unstable to the Rayleigh-Taylor instability (RTI). Early analytics and 2D simulations (Chevalier 1976; Chevalier & Klein 1978; Weaver & Woosley 1980; Benz & Thielemann 1990; Herant & Benz 1991) and modern 3D calculations (Hammer et al. 2010; Wongwathanarat et al. 2015; Utrobin et al. 2017) show significant impact on the density, velocity, and composition structure of the ejecta.  \nIt has been known for decades that the resulting compositional mixing can significantly alter the photometry of the SN. This effect has been roughly included in 1D modeling of Type IIP light curves resulting from explo-  \nsions deep within a red supergiant (Eastman & Pinto 1993; Utrobin 2007; Dessart & Hillier 2010, 2011). The mass densities and energy densities are also smeared out by the mixing from the RTI (see Bersten et al. 2011 for an early discussion raising this concern). In their recent modeling of the Type IIP SN 1999em, Utrobin et al. (2017) capture the impact of the RTI using a 3D model pre-breakout and connect to observable SN properties with a 1D post-breakout radiation calculation.  \nTo approximate the 3D effects of the RTI, we implement a scheme by Duffell (2016) that modifies the 1D spherical hydrodynamics equations. This scheme has been recently applied to the specific case of core collapse SN by P. Duffell et al. (2017, in preparation) and is now implemented in MESA for use along with the HLLC scheme. In this section, we describe the MESA implementation and compare to 3D calculations of Wongwathanarat et al. (2015). The use of the resulting RTI-mixed ejecta for SN lightcurves and velocities will be discussed in Section 6.', '5.1. Implementation of Duffell RTI': 'The Duffell (2016) scheme evolves an additional scalar quantity α R representing the relative strength of turbulent fluctuations. 4 The evolution equation for α R is an advection-diffusion equation with source terms. In Eulerian form, this is  \n∂ ∂t ( ρα R ) + 1 r 2 ∂ ∂r [ r 2 ( ρα R u -η R ∂ ∂r ( ρα R ) )] = S + α + S -α , (38)  \nwhere  \nS + α = ( A R + B R α R ) √ max ( 0 , -∂P ∂r ∂ρ ∂r ) , S -α = -D R α R ρc s r -1 , η R = C R α R c s r . (39)  \nThe source and sink terms S + α and S -α represent growth and decay of the turbulence, respectively. These terms along with a diffusion coefficient η R are determined via scaling arguments. The dimensionless coefficients in front of these quantities (growth coefficients A R , B R , diffusion coefficient C R , and decay coefficient D R ) are determined by calibrating a suite of 1D models against 3D hydrodynamics simulations. The original model of Duffell (2016) calibrates against 2D simulations; see P. Duffell et al. (2017, in preparation), for the recalibration of these constants to 3D simulations. The  \nvalues of the constants found by that 3D calibration are A R = 10 -3 , B R = 2 . 5, C R = 0 . 2, and D R = 2 . 0. In MESA , these constants are adjustable so that the user may explore the effect of varying them. For example, we show later the effect of D R = 2 . 0 , 3 . 0 , and4 . 0 on the mass fractions in massive star SN models at shock breakout.  \nAdditionally, a diffusive term (with diffusivity η R ) appears in each of the mass, momentum, and energy equations. For the sake of exploration in MESA , we allow each diffusivity to be scaled by an independent factor. With the diffusive term, the mass flux becomes (cf. Equation 23)  \nF ρ = A ρ ( u -˙ r ) -A η R ∂ρ ∂r , (40)  \nand the choice u = ˙ r (i.e., u face = S ∗ ) no longer causes this quantity to vanish. If no correction were applied, MESA would no longer preserve the mass coordinates of zone faces. In order to preserve the Lagrangian nature of the equations, we allow for an additional velocity between the cell face and the fluid. The advective flux introduced by the relative motion of the face will then exactly cancel this diffusive flux, restoring the Lagrangian nature of the scheme. Assuming ˙ r = u + δu , the no mass flux condition can be rewritten as  \nδu = η R 1 ρ ∂ρ ∂r = η R ∂ρ ∂m ( 1 ρ dm dr ) . (41)  \nThe term in parentheses is equal to A . In the finite volume form evolved by MESA , evaluating this condition at the cell face gives  \nδu = A η R ( ρ L -ρ R ) dm . (42)  \nTherefore, we modify the HLLC equation u face = S ∗ to  \nu face = S ∗ + δu , (43)  \nand proceed as in Section 4.1 (see Equation 20 and surrounding discussion). Usually | δu | glyph[lessmuch] S ∗ , so in practice this is a small modification and the HLLC scheme still works well.  \nFor a scalar quantity f , the flux is the sum of the diffusive flux plus the advective flux ( A ρfδu ) created by the velocity shift δu , that is  \nF f = A η R ∂ρ ∂r f ︸ ︷︷ ︸ advective -A η R ∂ ( ρf ) ∂r ︸ ︷︷ ︸ diffusive = -A η R ρ ∂f ∂r . (44)  \nRewriting the spatial derivative as a mass derivative gives  \nF f = -η R ( A ρ ) 2 ∂f ∂m . (45)  \nTo evaluate the fluxes for a cell k , we define  \nσ R ,k = η R ,k ( A k ρ k ) 2 1 dm k , (46)  \nwhere  \nη R ,k = C R α R ,k c s ,k r k . (47)  \nThe fluxes across faces are  \nF ρ,k = 0 , F p,k = A k P face ,k -σ R ,k ( u k -1 -u k ) , F e,k = A k P face ,k u face ,k + L k -σ R ,k ( e k -1 -e k ) . (48)  \nThe finite volume version of Equation (38) evolved by MESA is  \nα R -α R , start = δt [ 1 dm k ( F α,k +1 -F α,k ) +( A R + B R α R ) [ 1 ρ k √ max ( 0 , -∂P ∂r ∂ρ ∂r ) ] start -D R α R ,k ( 2 c s , k r k + r k +1 ) start ] , (49)  \nwhere  \nF α,k = -σ R ,k ( α R ,k -1 -α R ,k ) . (50)  \nWe evaluate the product of the P and ρ spatial derivatives as  \n∂P ∂r ∂ρ ∂r = ( P k -P k -1 dr k )( ρ k -ρ k -1 dr k ) , (51)  \nwhere  \ndr k = dm k 4 πr 2 C ρ k , (52)  \nwhich is numerically preferable to a subtraction of radial coordinates. At sharp jumps in density and pressure, these source terms can diverge, and therefore options to smooth ∂P ∂r ∂ρ ∂r are available, though they are off by default. In practice, smoothing does not appear to be necessary in cases we have explored, as HLLC typically smears out these sharp jumps over several cells in the model at the shock, and RTI mixing then smooths out the jumps more post-shock.', '5.2. Comparing a Munich 3D Model to MESA with Duffell RTI': 'We now develop a MESA analog to a specific 3D simulation of Wongwathanarat et al. (2015). This provides a comparison of the predictions from the MESA implementation of the RTI mixing described in the previous subsection (which we refer to as Duffell RTI) with those obtained in a 3D simulation. The Wongwathanarat et al.  \n(2015) progenitor model we use, L15-1-cw, has a mass of 15 M glyph[circledot] based on Limongi et al. (2000). We refer to this as the Munich L15 model. Note, as made clear in Wongwathanarat et al. (2015), most prior studies simulating RTI in SN envelopes disregard early-time asymmetries, relying on explosions that are initiated assuming spherical symmetry. Since those explosion asymmetries appear to have significant consequences, it is important to start from a 3D model like L15-1-cw when evaluating the use of MESA for SNe.  \nTo compare with the Munich L15 model, we construct a MESA starting model with similar parameters. Future studies of a variety of 3D models will be necessary to assess the impact on our 1D results of a variety of 3D asymmetries in the initial explosion. Just as Duffell RTI allows 1D simulations to capture many of the effects of the 3D RTI, it may be possible to extend 1D codes in the future to include relevant effects of explosion asymmetries in a self-consistent manner rather than by expediencies such as we describe below for initializing the 56 Ni abundance.  \nWe now compare the 3D shell averages of Wongwathanarat et al. (2015) to MESA with Duffell RTI enabled. The left panel of Figure 27 shows the resulting abundances when the shock is at 4 . 8 M glyph[circledot] , with the thin lines from MESA with D R = 3 and the thick lines the 3D averages from the Munich L15 model. For H, He, and O, the MESA lines agree with the Munich model. If nothing is done to take into account the initial asymmetry of the explosion, the radial extent of the 56 Ni in the Munich model far exceeds what can be achieved in MESA by Duffell RTI mixing. Hence, at this moment in the model evolution, we use the Munich L15 results to fix the extent of the MESA distribution of 56 Ni. Later mixing in the MESA run is done by Duffell RTI. The right panel of Figure 27 shows the comparison with the Munich model just before shock breakout. For this case, we show three simulations with D R = 2 , 3 and 4.  \nIn Figure 28 we show the MESA profiles of density (upper panel) and velocity (lower panel) at the moment when the shock is at 14 . 7 M glyph[circledot] . The solid lines are with Duffell RTI enabled, while the dotted lines are with it turned off. As shown by Wongwathanarat et al. (2015), Utrobin et al. (2017), and P. Duffell et al. (2017, in preparation) the operation of the RTI removes the unphysical density feature produced in 1D simulations without it. Such features can be seen in Figures 2 of Eastman et al. (1994) and Dessart & Hillier (2011) and in the dotted black line in the upper panel of Figure 28. Duffell RTI also alters the velocity structure of the material near the H/He boundary, as we discuss more in Section 6.6. The thick gray lines in both plots show the  \n1D shell averages of the 3D Munich L15 model. The fainter gray lines show the density and velocity profiles for a variety of angles in the Munich model. The asymmetries of the shock in the Munich model lead to its location varying between mass coordinates 10 . 5 M glyph[circledot] and 14 . 5 M glyph[circledot] . This variation with angle leads to 1D shell averages that do not show a sharp shock feature, but instead have more rounded shapes. Since the 1D MESA results have the shock at a single mass coordinate, they are similar to Munich profiles at a particular angle. This difference must be considered when comparing results from MESA to shell averages from the Munich model. It also shows that the time of shock breakout, which is well-defined in the 1D model, varies with angle in the 3D model.', '6. LIGHT CURVES AND VELOCITY EVOLUTION OF CORE COLLAPSE SUPERNOVAE': "We now present MESA modeling of the ejecta evolution triggered by core collapse in massive stars (roughly M > 8 M glyph[circledot] ). The new MESA capabilities enable selfconsistent calculations of photometric evolution of core collapse supernovae (SNe) using the STELLA code (Blinnikov et al. 1998; Blinnikov & Sorokina 2004; Baklanov et al. 2005; Blinnikov et al. 2006). A public version of STELLA is now included with the MESA distribution, and the interface from MESA to STELLA has been customized for ease of use. 5  \nOur main emphasis in this section is on the commonly observed Type IIP 'plateau' SNe that originate from energy deposited deep in the core of a M ≈ 8 -20 M glyph[circledot] red supergiant (Smartt 2009b). We also exhibit how these new capabilities enable simulations of core collapse events that occur after the star has lost the majority of its outer hydrogen envelope, the Type IIb and Ib SNe.  \nThe new capabilities we present are provided by a powerful combination of MESA and STELLA . Post core collapse evolution proceeds in two distinct phases. First we use MESA to evolve models from a few seconds after the central explosion triggered by core collapse to a time just before the outgoing shock reaches the stellar surface. These calculations make use of HLLC (Section 4) and Duffell RTI (Section 5). Subsequently, we use STELLA to evolve models through shock breakout and beyond the end of the plateau, generating light curves and velocities of the material at the photosphere and above.  \nSimulations using 3D models from the core collapse event to shock breakout are computationally expensive  \nFigure 28. Comparison of density (upper panel) and velocity (lower panel) profiles. The solid black line shows the MESA model using the Duffell RTI capabilities documented here with D R = 3; the dotted black line shows the same model run without the effects of the RTI. The thick gray lines show the 1D shell averages of the Munich model, while the fainter gray lines show the Munich model densities and velocities at different angles.  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \nFigure 27. Comparison of abundances in MESA models (thin lines) with 3D shell averages from the Munich L15 model (thick lines). This is a comparison of analogous models at similar times, so the goal is to illustrate qualitative agreement. Left panel: for D R = 3, the time when the forward shock is at 4.8 M glyph[circledot] and the reverse shock at about 4.0 M glyph[circledot] . The MESA 56 Ni curve is the result of artificially inserting 56 Ni in the model at this time. Right panel: the time the forward shock is at 14.7 M glyph[circledot] , so near breakout, and for D R = 2 , 3 and 4.  \n<!-- image -->  \n/circledot  \nbut now feasible (Wongwathanarat et al. 2015; Utrobin et al. 2017), and it will be a significant contribution to have more of them available in the future. To explore  \nthe subsequent glyph[greaterorsimilar] 100 days of photometric and spectroscopic evolution, 1D approximations are common. The new capabilities with MESA and STELLA also use a 1D approximation for both the pre- and post-breakout evolution. This provides a less computationally costly alternative for initial exploration of the parameter space for potential progenitors prior to or instead of doing a more realistic but more computationally costly 3D simulation. The pair MESA and STELLA can produce useful results in a few hours running on a modern multicore desktop workstation (see Section 6.7), while the 3D prebreakout evolution and post-breakout spectral analysis can take weeks running on a supercomputer. MESA and STELLA are not a replacement for the more computationally expensive codes but will be useful in conjunction with them.  \nThroughout this section, we present models that we characterize as 'similar to' observed SNe. We list the properties and parameters of these models in Table 3. As we discuss in Section 6.8, where we describe the procedure by which these models were generated, they are not 'best-fit' models. Rather, they simply serve as illustrative cases of these new capabilities.", '6.1. From Core Collapse to Near Breakout with MESA': 'Models of massive stars can be evolved in MESA up to the onset of the rapid infall of the iron core (see Paper I, Paper II, Paper III). However, MESA cannot model the core collapse event itself. Hence, to transition from the onset of core infall to the explosion phase, we rely on a a variety of approximate procedures (Paper III).  \nFor the current efforts, our approach is as follows. We remove the center section of the model at the location where the entropy per baryon is 4 k B , excising the portion of the model that will have collapsed to form a proto-neutron star. This corresponds approximately to the iron core, typically at about 1 . 5 M glyph[circledot] . We allow the model to continue infall until its inner boundary (IB) reaches 200 km, near the location of the stalled shock (thanks to H.-T. Janka, private communication, for suggesting this scheme). After the first few seconds, we account for further fallback by removing negative velocity material at the IB. We are not seeking a numerical model of realistic fallback since that depends on 3D details of the explosion that are beyond what MESA can simulate.  \nThe stellar explosion is induced by injecting energy in a thin layer of approximately 0 . 01 M glyph[circledot] at the IB for 5 ms, at a rate sufficient to raise the total energy of the model to a user-specified value. In the subsequent evolution, nuclear reactions are allowed to change abundances but not to generate energy. This choice is suitable because we are not seeking accurate nucleosynthetic yields. The explosion energy spent to photodisintegrate the core to a mix of protons, neutrons, and alpha particles is soon after roughly repaid by energy released as those particles recombine to form products such as 56 Ni. Getting an accurate accounting of the energy balance of that complex process is beyond the scope of this paper and is not attempted in the following examples. Our choice to exclude nuclear energy generation can be seen as a simplifying assumption that the cost of photodisintegration is balanced by the return from later recombination. For users wishing to refine this, any excess change in energy from nuclear reactions can be included in the specification of the post-explosion total energy of the model.  \nThe conservation of total energy throughout the run is estimated by summing the per-step errors from postexplosion to near breakout. At each timestep, we compare the actual change in total energy between initial and final models for the step, to the change expected from surface luminosity and neutrino losses over the duration of the step. The runs for the models reported below typically show relative cumulative errors in conservation of total energy of less than 1%, with most of that error happening in the first few minutes post explosion when the shock is most extreme. For later stages, the cumulative relative error is orders of magnitude smaller.  \nThe post-explosion evolution of the MESA model is determined by the shock traversal through the star and the resulting Duffell RTI. Figure 29 illustrates the difference between models with and without the effects of the RTI by showing density and pressure profiles. They are shown when the forward shock is about halfway through  \nthe star and when the reverse shock originating at the H/He boundary has reached ≈ 4 M glyph[circledot] on its way to the center. The reverse shock is primarily responsible for the large RTI effects evident in the plots. The online animated figure shows the time evolution of these and many other quantities of interest from seconds after explosion to near shock breakout.', '6.2. From Near Breakout through the Plateau: Introducing STELLA': 'To follow the evolution of the model through shock breakout and beyond we use a multigroup (i.e., frequencydependent) radiation hydrodynamics code. 6 When the shock is near breakout, we hand the MESA model off to STELLA in an appropriate form which involves interpolating to the desired grid and optionally adding circumstellar material (CSM) according to user specifications. With that done, MESA is finished, and STELLA takes over (see Section 6.3 for a discussion of how we select when to hand off).  \nSTELLA (Blinnikov et al. 1998; Blinnikov & Sorokina 2004; Baklanov et al. 2005; Blinnikov et al. 2006) is able to model SN evolution at early times, before the expansion is homologous. It can also handle shock breakout and interaction with circumstellar material outside the conventional stellar photosphere. STELLA is an implicitly differenced hydrodynamics code that incorporates multigroup radiative transfer. The time-dependent equations are solved implicitly for the angular moments of intensity averaged over fixed frequency bands. STELLA takes about the same amount of time for near-breakout to post-plateau evolution as MESA takes to simulate from explosion to near-breakout: about an hour on current workstations.  \nSTELLA solves the radiative transfer equations in the intensity momentum approximation in each frequency bin. We use from 40 to 200 frequency groups, enough to produce bolometric luminosities and broad-band colors, but not sufficient to produce spectra. Better broad-band light curves can be produced with the larger number of frequency groups, but 40 is sufficient for a bolometric lightcurve and gives faster runtimes since each group must be represented by a variable and an equation at each zone. The opacity is computed based on over 153,000 spectral lines from Kurucz & Bell (1995) and Verner et al. (1996). The expansion opacity formalism from Eastman & Pinto (1993) is used for line opacities taking high velocity gradients into account. The opac-  \nFigure 29. Density (left y-axis, orange curves) and pressure (right y-axis, blue curves) for MESA models with Duffell RTI (left panel) and without any RTI effects (right panel) at a time when the forward shock is approximately halfway through the star. The online animated Figure shows the time evolution of these and other quantities for each case.  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \nity also includes photoionization, free-free absorption, and electron scattering. LTE is assumed in the plasma, which allows the use of the Boltzmann-Saha distribution for ionization and level populations. STELLA does not include a nuclear reaction network except for the radioactive decay chain initiated from 56 Ni. For calculating the overall opacity, the code uses 16 species: H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, a sum of stable Fe and radioactive 56 Co, and stable Ni and radioactive 56 Ni. Energy from nickel and cobalt radioactive decay is deposited as positrons and gamma-rays and is treated in a one-group transport approximation according to Swartz et al. (1995).  \nSTELLA solves the conservation equations for mass, momentum, and total energy on a Lagrangian comoving grid. It employs artificial viscosity based on the standard von Neumann artificial viscous pressure used for stabilizing solutions (Von Neumann & Richtmyer 1950) and a cold artificial viscosity used to smear shocks (Blinnikov et al. 1998; Moriya et al. 2013). The coupled equations of radiation hydrodynamics are solved through an implicit high-order predictor-corrector procedure based on the methods of Gear (1971) and Brayton et al. (1972); see Blinnikov & Panov (1996) and Stabrowski (1997) for details.  \nWe explore the sensitivity of bolometric light curves ( L bol ) reported by STELLA to the number of frequency bins, spatial zoning, and error tolerances. The result of our sensitivity study is that 40 frequency bins, 300 spatial zones, and an error tolerance 0 . 001 for the GearBrayton method typically give a converged model. In our experience using MESA and STELLA for Type IIP SNe, we have not found cases that require different values for number of frequency bins and error tolerance. Some  \ncases may need a larger number of zones in order to minimize numerical artifacts producing spurious oscillations in the light curve. This problem can often be fixed by a relatively small increase in the number of zones; this is shown for a case similar to SN 2012A in Figure 30.  \nFigure 30. Effect of number of STELLA zones on a model (12A) similar to SN 2012A. The inset zooms in on the region of the plateau that exhibits numerical oscillations. A small increase in the number of zones significantly reduces this artifact and doubling the number of zones almost completely removes it.  \n<!-- image -->', '6.3. Handing Off from MESA to STELLA': "A time must be chosen to hand off the MESA model to STELLA . This choice is driven by a compromise between two considerations. First, RTI modeling ceases once STELLA is running even though the effects of RTI may not be complete at that time. Therefore, one wants  \nFigure 31. Bolometric lightcurves for Type IIP SNe obtained where the transition between MESA and STELLA has been done at different times relative to shock breakout. Performing the handoff just before breakout (blue curve, day 0) is the recommended choice.  \n<!-- image -->  \nthe model to remain in MESA as long as possible. But second, STELLA more accurately handles shock breakout and the outermost layers, especially if any matter is placed above the photosphere or if significant radiation is free-streaming from just below the photosphere prior to shock breakout. Moreover, the sophisticated multigroup radiation transfer of STELLA will do a much better job than (gray) MESA can at later times post breakout. Hence for longer-term lightcurve evolution, this motivates the default choice to perform this handoff just before breakout.  \nIn order to illustrate the effects of this choice, Figure 31 shows bolometric light curves for cases where the handoffs are done at different times. Note that this plot shows MESA being forced to run post-breakout even though that is not recommended. The deviation of the light curves for later handoff are primarily the result of STELLA doing a better job because of its multigroup radiation transfer rather than any late-stage RTI effects being captured by MESA that are missed by STELLA . That is because, for this case, the H envelope is of normal thickness and the reverse shock from the H/He boundary has time to reach the center, completing essentially all of the RTI effects before breakout.  \nIn the runs presented in the remainder of this section, we choose to do the MESA -toSTELLA handoff shortly before breakout, as determined by the outgoing shock front reaching a location 0 . 11 M glyph[circledot] below the surface of the model (this location is a user-defined parameter). Again, we note that in some cases the reverse shock is still far from the center at this moment, and not all of  \nthe RTI mixing has completed. In particular, this is true for models with a partially stripped envelope (see Section 6.9). For now, this remains a caveat for the user; a solution would be to have the post-breakout radiation hydrodynamics code include a treatment of the effects of RTI. When presenting the results, we define t = 0 as the time of shock breakout-which we identify using the peak of the bolometric luminosity-and not the (earlier) time of the MESA -toSTELLA handoff.  \nBecause of STELLA 's treatment of radiation hydrodynamics, we have not had to take the progress of the model toward homologous expansion into consideration in selecting a time to hand off from MESA . However this is a consideration for doing a handoff to radiative transfer codes that assume homology. More accurate spectral and lightcurve modeling with full radiative transfer, such as EDDINGTON (Eastman & Pinto 1993), SEDONA (Kasen et al. 2006), and CMFGEN (Dessart & Hillier 2010), assume homologous expansion in their current applications to SNe, and this should be considered when deciding the time to hand off from another simulation. Indeed, Eastman et al. (1994) and Dessart & Hillier (2011) discuss this challenge, especially for the inner-most material that has not reached a homologous stage and can still have a reverse shock running through it. Approximations made in mapping to a thereafter homologous code can impact the late-time photospheric velocity evolution and the nebular line width predictions associated with the innermost ejecta.  \nIn contrast, STELLA does not assume homologous expansion, so early handoffs are fine; it can handle the effects of remaining pressure gradients as the model moves toward homologous expansion. This is important, as the time it takes to reach homology in these models can be quite long. Figure 32 shows velocity evolution results for a model similar to SN 1999em (see discussion in Section 6.8). Homologous expansion would imply that v/r is flat, whereas a 20% variation from simple homology is evident at 20 days. An additional way in which homology can be violated long after shock breakout is from 56 Ni decay, especially in Type Ia SNe (Woosley et al. 2007b). As is evident in Figure 32, the much smaller mass fractions of 56 Ni in Type IIP SNe do not cause such a problem. The contrasting light curves with and without 56 Ni are shown in Figure 33, exhibiting the prolonging of the plateau due to radioactive decay (Kasen & Woosley 2009; Sukhbold et al. 2016a).", '6.4. Connecting to Observations: Photospheric Properties from STELLA': 'To set the stage for the rest of this section, we describe a particular model in detail. Figure 34 shows the  \nFigure 32. Profile of v/r throughout the model, normalized to the values at the surface, at a number of epochs. In homologous expansion, the profiles would be constant and equal to unity. The models are still significantly non-homologous at 20 days. The bottom panel shows the 56 Ni mass fraction for this model at MESA -toSTELLA handoff (near shock breakout). The extended 56 Ni profile (total nickel mass ≈ 0 . 04 M glyph[circledot] ) indicates that it has been mixed by RTI effects out to near the surface. However, as illustrated by the dotted lines in the upper panel, the decay energy does not have a significant effect on the approach to homology.  \n<!-- image -->  \n/circledot  \nevolution of a model we construct to be similar to the Type IIP SN 1999em (99em 19 in Table 3). The quantities shown are those generated during the STELLA phase of the evolution. Panel (a) is the bolometric luminosity, while panel (b) shows velocity at the location of the photosphere (where τ Ros = 2 / 3) and panels (c) and (d) show the mass and radius coordinate of this location. This illustrates the familiar result that the photosphere only reaches the deeper parts of the ejecta after about day 50. The radiation and gas temperatures at the photosphere are shown in panel (e), as is an effective temperature defined by the bolometric luminosity leaving the photosphere. Panel (f) shows the optical depth to the IB, highlighting that the radiative diffusion approximation is excellent (since τ IB glyph[greatermuch] 1) until day 120, at which point the plateau ends and the IB temperature (panel g) approaches that of the photosphere. (Curves showing photospheric quantities stop once τ IB < 3.) Meanwhile, the photospheric radius (panel d) stays remarkably constant throughout the plateau.  \nOur emphasis is on bolometric luminosities, where 40 STELLA frequency bins is adequate. However, broad-  \nFigure 33. Lightcurves and velocities for a model similar to SN 1999em with ≈ 0 . 04 M glyph[circledot] 56 Ni and without 56 Ni. The main effect of the radioactive decay is to prolong the plateau.  \n<!-- image -->  \nband light curves are also reported by STELLA . Figure 35 shows how the STELLA colors change as one goes from 40 to 200 frequency bins in a model approximately matching the bolometric luminosity of SN 1999em (99em 19 in Table 3). This reflects that a given band is spanned by only a small number of frequency bins. The nonpublic research version of STELLA can opt to use many more frequency bins to address under-resolution issues. There are no current plans to include that capability in MESA . We also show what a blackbody would predict using the MESA colors module (see Appendix A). This makes it clear that the line-blanketing in the U band is well handled by STELLA . We do not include colors in our subsequent discussions, but we expect they may be useful to users who have access to observations in one or two bands, but not enough data to produce a bolometric light curve from observations.', '6.5. Connecting to Observations: Fe II Line Velocities': "It is important to be able to interpret the ejecta velocities measured by observers, which are often inferred from the absorption minimum in the Fe II 5169 ˚ A line. Modeling these absorption features requires more detailed radiative transfer than available in STELLA . However, rather than assume the photospheric velocity reported by STELLA is identical to the Fe II 5169 ˚ A line, we have added the capability of finding the location (and hence the velocity) of material above the photosphere where the Sobolev optical depth in the Fe II 5169 ˚ A line  \nFigure 34. Photosphere and IB properties of a model (99em 19 in Table 3), similar to the Type IIP SN 1999em, as a function of time. From top-to-bottom, the figure shows the bolometric luminosity, the velocity, mass coordinate, and radius coordinate of the photosphere, three temperatures (gas, radiation, and 'effective') at the photosphere, the optical depth to the IB, and the IB temperature. The gray line in panel (c) shows the Lagrangian mass coordinate of the IB.  \n<!-- image -->  \nis a specified value. 7 This will prove to be most important after day 30 or so, when the photosphere has started to move inward in mass coordinate into ejecta with a shallow density profile.  \nThe strength of a line in a homologously expanding atmosphere is quantified by the Sobolev optical depth (Sobolev 1960; Castor 1970; Mihalas 1978; Kasen et al.  \nFigure 35. Comparison of model 99em 19 with the multicolor lightcurve of SN 1999em, showing colors from STELLA and blackbody colors from MESA . Circles indicate observational data. This demonstrates the effect of the number of STELLA frequency bins on the predicted colors.  \n<!-- image -->  \n2006), which for the Fe II line at any position is  \nτ Sob = πq 2 e m e c n Fe η i ft exp λ 0 , (53)  \nwhere λ 0 = 5169 ˚ A is the line center wavelength for the Fe II line, f = 0 . 023 is its oscillator strength, n Fe is the number density of iron atoms, and t exp is the time since breakout. The quantity η i is the fraction of iron atoms that are in the lower level of the transition of interest and depends on the properties of the gas. D. Kasen (2017, private communication) provided an η i ( ρ, T ) table for post-processing to produce the Fe II line velocities, calculated under the assumption of LTE and covering log( ρ/ g cm -3 ) = -16 to -8 and log( T/ K) = 3 . 3 to 4 . 3.  \nWe use Equation (53) after the STELLA run to provide the velocity of material that satisfies a chosen value of τ Sob . This yields a velocity that can be compared to the measured Fe II line velocities. Figure 36 shows the resulting comparisons for various choices of τ Sob for a model similar to the Type IIP SN 2012A found solely by matching the bolometric luminosity (upper panel). The lower panel displays the Fe II 5169 ˚ A data and the velocities derived from the photosphere and for a range of values of τ Sob . At early times, there is little difference between the photospheric velocity and that of the Fe line. However, as the photosphere moves deeper into the ejecta, the two velocities substantially diverge. The velocity inferred from the Sobolev argument gives a much better match to observations than the photospheric velocity. Motivated by this compari-  \nwe choose τ Sob = 1 for our later plots, a parameter that the user is free to adjust.  \nFigure 36. Comparison of different definitions of velocity for a model (12A) similar to SN 2012A. The upper panel shows the data and the model for the bolometric luminosity. The lower panel shows the velocity of a few different locations depending on the Sobolev optical depth in the Fe II line.  \n<!-- image -->", '6.6. The Impact of Pre-breakout RTI Mixing': "We have previously outlined the inclusion of a method for RTI mixing in MESA (the Duffell scheme; Section 5), the use of MESA to evolve models pre-breakout (Section 6.1), the use of STELLA to evolve models postbreakout (Section 6.2) and described how to connect the models to observations (Sections 6.4 and 6.5).  \nIn this way, MESA plus STELLA allows users to explore the impact of RTI mixing on Type IIP light curves and velocities. Prior work in this direction (Eastman et al. 1994; Utrobin 2007; Dessart & Hillier 2010, 2011; Morozova et al. 2015) focused on the impact of compositional mixing, often with averaging approaches to achieve various levels of mixing. Only the recent work of Utrobin et al. (2017) incorporated compositional mixing from a 3D model and also included the modified density and velocity structures, also seen in the 1D RTI mixing (P. Duffell et al. 2017, in preparation).  \nFigure 37 shows the lightcurve and velocities of model 99em 19. The luminosity without RTI mixing has a dis-  \ntinctive rise just before the plateau as shown by Eastman et al. (1994) and Utrobin (2007). As RTI causes many associated changes in composition, density, velocity, and energy density for the innermost material, we cannot specifically identify the immediate cause of the lengthening of the plateau phase when RTI is incorporated without further experiments. These are now possible using MESA and STELLA but are beyond the scope of this paper. The lower panel shows the photospheric and Fe II line velocities with and without RTI mixing. The most evident change is at the end of the plateau, when the material that was near the H/He boundary in the red supergiant is approaching the SN photosphere. That material is strongly affected by RTI mixing as shown in Figure 28 and discussed in P. Duffell et al. (2017, in preparation).  \nFigure 37. Effect of Duffell RTI on lightcurves (upper panel) and velocities (lower panel) of a Type IIP SN model (99em 19) similar to SN 1999em. The inset shows the time near day 120 where the altered density structure causes a significant difference in the Fe II line velocities.  \n<!-- image -->  \nTo enable exploration of the impact of various components of the RTI mixing, we explicitly allow for the diffusion coefficients for density, momentum, energy, and composition to be scaled by independent constant factors relative to the value η R given in Equation (47). We show in Figure 38 the impact of varying the coefficient in the internal energy flux in Equation (48), which we refer to as η R ,e . These plots show the energy density and density of the ejecta just before shock breakout in one of our models (99em 19) in Table 3). The blue line  \nFigure 38. Structural effect of the RTI energy diffusion coefficient η R ,e . Other parameters are the same as model 99em 19. The profiles are shown just before shock breakout. The effects on the lightcurve are shown in Figure 40.  \n<!-- image -->  \n/circledot  \nis for the fiducial value, whereas the red line is for an extreme increase of a factor of 100. The only locations that are sensitive to these changes are the innermost mass coordinates where RTI was most active, the same regions where the lightcurve and velocities seem to be sensitive to changes related to RTI mixing. The variable η R ,e was found to be a useful 'knob' to vary for modeling of specific SNe.", '6.7. Exploring the Explosive Landscape': "A strength of the new MESA plus STELLA capabilities is their ease of use. This enables detailed quantitative studies of large numbers of core collapse SNe. The open source nature of MESA , the inclusion of STELLA in the MESA distribution, and the repository of examples contained within the MESAstar test suite allow a user to obtain models that can be compared directly to observations. Indeed, with minimal manual intervention, a user can take a star from the pre-MS to a SN light curve within a few hours of computer runtime. To emphasize this point, we describe here how this might be done (Section 6.7.1). To demonstrate how parameter choices affect lightcurves, we show a large sample of variations of a standard case for 'high-middle-low' settings of some of the main parameters (Section 6.7.2). In Section 6.8 we will exhibit a few specific models created to be roughly similar to known Type IIP SNe. The potential is clear for an extensive database of such SNe models created  \nFigure 39. Effect of single-parameter variations to the progenitor std 16. The upper three panels vary initial properties of the star; the lower three vary modeling assumptions during evolution to core-collapse.  \n<!-- image -->  \nFigure 40. Effect of single-parameter variations associated with the SN explosion. Unvaried parameters have the values of the model listed in the upper right of that panel.  \n<!-- image -->  \nusing MESA and STELLA ; its actualization is beyond the scope of this paper.", '6.7.1. Generating Models with MESA plus STELLA': "The first step in generating a core collapse SN lightcurve is to use MESA to make a pre-SN stellar model that is undergoing core collapse. The test case example make pre ccsn can serve as a useful template. As part of the required inlists, the user must select values for the main variables: initial mass ( M ZAMS ), initial metallicity ( Z ), initial rotation (( v/v c ) ZAMS ), overshooting parameter ( f ov ), wind scaling factor ( η wind ), and the mixing length for MLT in the H envelope ( α MLT , H ). Of course, the user may tune other MESA parameters of interest. The run from pre-MS to Fe core infall runs automatically given these parameters and, depending on the case, takes roughly an hour on a modern multicore desktop workstation. Users interested in details of pre-SN models may require settings that lead to significant additional runtime (e.g., Farmer et al. 2016; Renzo et al. 2017).  \nThe second step loads the model at core infall into MESA , emulates the core collapse explosion by excising the core and injecting energy and Ni (as described in Section 6.1) and evolves until near shock breakout. The test case example ccsn IIp can serve as a useful template. Again, the user must set the value of the various 'knobs' controlling the properties of the explosion such as the total energy E and the 56 Ni mass M Ni . Early ( t < 20 days) lightcurves of core collapse SNe are better fit when large amounts of CSM are placed outside the conventional photosphere (Morozova et al. 2016; Dessart et al. 2017; Morozova et al. 2017a,b). We provide an option to include CSM. We also provide the option for 'boxcar' smoothing of the model abundances before the handoff from MESA to STELLA (Kasen & Woosley 2009; Dessart et al. 2013; Morozova et al. 2015). The end result of this step is a model suitable for input into STELLA , so one must also indicate the number of STELLA zones to be used. This MESA phase from after explosion to near breakout typically takes about 30 minutes on a modern multicore desktop workstation.  \nThe final step uses the results produced in the previous step as input to STELLA and evolves the model through shock break-out to the post-plateau phase. A script to execute STELLA is provided. This stage takes about an additional 30 minutes on a modern multicore desktop workstation for typical cases. When STELLA finishes, a post-processing step produces data for comparison to observational results.  \nFigure 39 exhibits the std 16 model lightcurves as progenitor parameters are varied. Many variations behave as expected from previous analytical and numerical scalings (Popov 1993; Kasen & Woosley 2009; Sukhbold et al. 2016a). For example, the decrease in the plateau duration with lower ZAMS masses or higher mass loss (increased η wind ) is as expected. The increase in plateau luminosity with decreasing α MLT , H is because those stars with lower α MLT , H have a larger stellar radius at time of explosion. However, other variations in these figures are not as easily diagnosed.  \nFigure 40 exhibits model lightcurves as explosion parameters are varied. Again, many cases lead to the expected outcomes, such as the increase in the plateau luminosity with increasing explosion energy and the increased duration of the plateau with increasing nickel mass. The changes caused by varying the RTI parameters are slight for the compositional mixing and boxcars, though, as we discussed in Section 6.6, modifying the diffusion of energy density during RTI does impact the shape at the end of the plateau. The impact of the CSM is similar to that shown by Morozova et al. (2017a) and Dessart et al. (2017).  \nWith experience in the effects of varying the parameters (knobs) shown in Figures 39 and 40, it is sometimes possible to get a rough match between model and observations after a dozen or so attempts. That is about the amount of effort we undertook to get the models similar to various observed SNe presented in Section 6.8. Of course the effects of the various knobs do not combine in any simple manner, so it can be a nontrivial challenge to find a combination that gives a good match for both velocities and lightcurve. Our experience suggests that it is a good strategy to match velocities before lightcurves since there are few ways available to shift velocities and many ways to change lightcurves. It is important to include velocities in judging potential matches because of the multiple degeneracies, as will be seen below where we show two models similar to SN 1999em with quite different ejecta masses and explosion energies. Even when using both velocities and lightcurves, it remains a challenge to find a unique 'best' match.", '6.8. Applications to a Few Type IIP SNe': "To show examples of what can be accomplished with these new capabilities, we have modeled four Type IIP SNe: 1999em, 2005cs, 2009N and 2012A. These cover a range of luminosities, plateau durations and nickel masses and have readily available data (Pejcha & Prieto 2015a,b) for bolometric luminosities and Fe II ve-  \nlocities. 8 We follow the steps described in Section 6.7, iterating to reach the matches shown. The models are not intended to demonstrate the best matches that can be achieved using MESA and STELLA . An investment of more effort could produce better matches, but is beyond the scope of this paper. The parameters we choose are shown in Table 3.  \nFigure 41. Effect of CSM on a model (99em 19) similar to SN 1999em. The influence is particularly apparent in the early time lightcurves and velocities.  \n<!-- image -->  \nWe note a few general insights gained from our modeling. We found that the radii of red supergiant models from MESA were too large for these Type IIP SNe models unless we set α MLT , H = 3. All models benefited at early times by having some CSM present. Figure 41 shows how the early 1999em model predictions change as CSM is added to the value shown in Table 3. The luminosity at early times is a far better match, as are the earliest velocity data. As expected, by day 50 and beyond there is no impact of the CSM on the model predictions. Comparisons of how the luminosity collapsed at the end of the plateau drove us to prefer an enhancement in η R ,e in several cases.  \nTo exhibit some of the possible degeneracies, we constructed two distinct models for 1999em. As shown in Figure 42, they are both reasonable models for the bolometric luminosity and Fe II velocities. However,  \nTable 3. Key properties and parameters associated with the core-collapse SN models. The column 'case' identifies the model. The 'progenitor parameters' sub-table lists input parameters used during the MESA evolution of the models to core infall: initial mass ( M ZAMS ), initial metallicity ( Z ), initial rotation (( v/v c ) ZAMS ), overshooting parameter ( f ov ), wind scaling factor ( η wind ), and the mixing length for MLT in the H envelope ( α MLT , H ). The 'stellar properties at the time of explosion' sub-table lists physical quantities evaluated in the MESA model at the time the Fe core begins to infall: mass ( M exp ), effective temperature ( T eff ), radius ( R exp ), luminosity ( L exp ), mass of the He core ( M He ), and initial mass of the Fe core that will be excised ( M c , i ). The 'explosion properties and parameters' sub-table lists input parameters like the total energy after explosion E exp and the 56 Ni mass M Ni as well as properties of the model including the final core mass after fallback ( M c , f ) and the total ejecta mass ( M ej ). This sub-table also lists input parameters used in the MESA plus STELLA modeling such as the RTI parameter ( η R ,e ) and the number of boxcar smoothing passes ('boxcar'). Parameters controlling the extent of the CSM are also needed; for a wind profile this includes the wind duration ( t CSM ), mass loss rate ( ˙ M CSM ), and velocity ( v CSM ). Many properties are omitted for the stripped case because this is an ensemble of models with a range of envelope stripping (see Section 6.9).stellar properties at time of explosion", 'explosion parameters and properties': "Figure 42. Two models similar to SN 1999em (99em 16 and 99em 19) with significantly different ejecta masses and total energies.  \n<!-- image -->  \nFigure 43. Model (05cs) similar to the low luminosity example SN 2005cs.  \n<!-- image -->  \ntheir ejected masses and radii differ significantly, one has 12 . 9 M glyph[circledot] and 770 R glyph[circledot] , whereas the other has 15 . 9 M glyph[circledot] and 600R glyph[circledot] . Utrobin (2007) gave an ejected mass of 19 . 0 ± 1 . 2 M glyph[circledot] , a radius of 500 R glyph[circledot] , and an explosion energy of 1 . 3 × 10 51 erg. Bersten et al. (2011) gave an ejecta mass of 17 . 6 M glyph[circledot] , radius of 800 R glyph[circledot] , and explosion energy  \nFigure 44. Model (09N) similar to SN 2009N. Further experiments might produce a model with a better match to the drop at the end of the plateau. Alternatively, this model might be a useful start when looking for a match to an observed light curve with a slow decline from the plateau.  \n<!-- image -->  \nof 1 . 25 × 10 51 erg. Utrobin et al. (2017) model this event with a 3D simulation from explosion to shock breakout, similar to the Munich L15 model we discuss in Section 5, but with an explosion energy of about 0 . 5 × 10 51 erg. For comparison, the MESA models for 1999em have total energies after explosion of 0 . 60 × 10 51 erg for the case with 12 . 9 M glyph[circledot] ejected mass, and 0 . 78 × 10 51 erg for the case with 16 . 3 M glyph[circledot] ejected mass.  \nWe previously showed 2012A in Figure 36. Our model had an ejected mass of 10 . 1 M glyph[circledot] to compare with 7 . 8 M glyph[circledot] from Morozova et al. (2017b), 12 . 5 M glyph[circledot] from Tomasella et al. (2013) and 13 . 1 ± 0 . 7 M glyph[circledot] from Utrobin & Chugai (2015). Tomasella et al. (2013) also reported a progenitor luminosity of log( L/ L glyph[circledot] ) = 4 . 73 ± 0 . 13, just a bit fainter than our model's value. Figure 43 shows our model for 2005cs. Our model has an ejected mass of 10 . 4 M glyph[circledot] , slightly higher than the 9 . 5 M glyph[circledot] reported by Spiro et al. (2014) and the 7 . 8 M glyph[circledot] reported by Morozova et al. (2017b). Figure 44 shows our model for 2009N, which has an ejected mass of 9 . 9 M glyph[circledot] , whereas Morozova et al. (2017b) found 9 . 3 M glyph[circledot] and Tak'ats et al. (2014) found 11 . 5 M glyph[circledot] .", '6.9. Partially Stripped Core Collapse SNe': "There is a well-defined class of core collapse SNe where either much (Type IIb) or nearly all (Type Ib and Ic) of the H envelope was lost prior to the core collapse event.  \nFigure 45. Models of partially stripped SNe. These M ZAMS = 17M glyph[circledot] models have a range of H envelope masses giving rise to a range of plateau durations. The upper panel shows the bolometric luminosity while the lower panel shows the velocity. All models have the same total energy post-explosion of 0 . 65 × 10 51 erg and a 56 Ni mass of 0 . 037 M glyph[circledot] . Other model parameters are indicated in Table 3 (case 'stripped'). The lowest mass case has about 0 . 1 M glyph[circledot] of the H envelope remaining, similar to that of a Type IIb SN (Ergon et al. 2015).  \n<!-- image -->  \nDessart et al. (2015) performed detailed radiative transfer models for a large set of progenitors from binary evolution, while Morozova et al. (2015) carried out diffusive calculations with varying amounts of mass loss. Yoon et al. (2017) explored MESA models constructed from binary transfer scenarios and applied them to a set of well observed Type IIb events. We have not yet been able to deal successfully with Ic models because of numerical problems related to the extreme ejecta velocities that occur at shock breakout. However, it is possible to do both IIb and Ib models as shown here.  \nIn Figure 45, we show the MESA plus STELLA predictions for luminosities and photospheric velocities for a range of models with varying amounts of mass stripped from a 17 M glyph[circledot] ZAMS model, ranging from the entire initial H envelope still remaining down to only 0 . 1 M glyph[circledot] of the H envelope left at the time of explosion. Similar to Figure 7 of Morozova et al. (2015), the plateau period becomes shorter as the residual H shell mass declines. Our smallest mass model has an H envelope mass comparable to typical models of Type IIb SNe and generates a light curve comparable to observed Type IIb SNe (Ergon et al. 2015). Figure 46 shows the interior properties  \nof these same models near the moment of shock breakout. For models which have been stripped, the reverse shock has not reached the IB at the time the forward shock reaches the surface. Since RTI mixing does not occur in STELLA , these models would incompletely include the effects of the RTI.  \nFigure 46. Models of M ZAMS = 17M glyph[circledot] stars (case 'stripped') that have experienced a range of stripping. The density, velocity and pressure profiles are shown at the time of handoff from MESA to STELLA , very close to shock breakout. The gray band shows the range of locations of the H/He boundary at the time of explosion.  \n<!-- image -->  \n/circledot  \nCao et al. (2013) discovered the fully stripped Type Ib SNe iPTF13bvn in the nearby spiral galaxy NGC 5806 with the intermediate Palomar Transient Factory (Law et al. 2009). This is one of only a few stripped SNe with a progenitor detection. Using data from Cao et al. (2013) and Fremling et al. (2014), we show in Figure 47 our model that approximately matches the iPTF13bvn light curve. The model is derived from an 11 M glyph[circledot] ZAMS model and has a remaining mass of only 3 . 4 M glyph[circledot] at the time of explosion with total energy after explosion of 0 . 95 × 10 51 ergs and a 56 Ni mass of 0 . 11 M glyph[circledot] distributed throughout the remaining star (ejecta mass 1.8 M glyph[circledot] ). Fremling et al. (2014) also modeled this lightcurve, finding the total energy to be 0 . 85 +0 . 5 -0 . 4 × 10 51 ergs, with a a 56 Ni mass of 0 . 049 +0 . 02 -0 . 012 M glyph[circledot] , and total ejecta mass of  \n1 . 94 +0 . 50 -0 . 58 M glyph[circledot] . Our parameters are similar, falling within the range of the quoted uncertainties, except for the 56 Ni mass.Figure 47. Comparison of model (13bvn) similar to the stripped Type Ib SN iPTF13bvn with the observed bolometric lightcurve.  \n<!-- image -->", '7. BLACK HOLE FORMATION': "Compact objects are a natural product of the evolution of massive stars. A broad consensus on which massive stars produce black holes (BH) has not yet been reached (Timmes et al. 1996; Fryer & Kalogera 2001; Heger et al. 2003; Eldridge & Tout 2004; Zhang et al. 2008; Ugliano et al. 2012; Clausen et al. 2015; Sukhbold et al. 2016b; Muller et al. 2016; Limongi 2017).  \nThe lack of consensus is due to a variety of differences in the modeling, including stellar wind treatments during the pre-supernova stage (Renzo et al. 2017); shellular rotation prescriptions (e.g., Limongi 2017); sensitivity to the initial metallicity (e.g., O'Connor & Ott 2011), number of isotopes in the reaction network (Farmer et al. 2016), adopted values of critical reaction rates (deBoer et al. 2017; Fields et al. 2017), and ignition of core carbon burning (Farmer et al. 2015; Cristini et al. 2017; Petermann et al. 2017); variations from spatial and temporal resolution (Farmer et al. 2016); convection during core-collapse (e.g., Couch et al. 2015); and effects from binary partners (e.g., Marchant et al. 2016; Batta et al. 2017). In addition, current estimates of the neutron star and BH initial mass function chiefly rely on parameterized explosion models and not on first principles calculations.  \nThis section explores MESA models that can produce BHs. First, we consider M ZAMS ≤ 60M glyph[circledot] models that  \ncan form a BH without encountering dynamical instability due to e + e -pair production. Second, we survey M ZAMS ≥ 60M glyph[circledot] models that encounter dynamical instability, either entering the Γ 1 ≤ 4/3 regime once to produce a pair-instability supernova (PISN) (Fowler & Hoyle 1964; Rakavy et al. 1967; Rakavy & Shaviv 1967; Barkat et al. 1967; Fraley 1968; Ober et al. 1983; Fryer et al. 2001; Scannapieco et al. 2005; Kasen et al. 2011; Chatzopoulos et al. 2013), or multiple times to produce a pulsational pair-instability supernova (PPISN) and a BH remnant (Barkat et al. 1967; Woosley et al. 2007a; Chatzopoulos & Wheeler 2012; Woosley 2017).  \nFigure 48. He core mass (upper panel) and mass location where the gravitational binding energy is equal to 10 48 erg (lower panel) for ZAMS masses of 30, 45, and 60 M glyph[circledot] . Three stellar wind treatments, two wind scaling factors, and two rotation rates are shown for each ZAMS mass. The variation, illustrated by the tan band, induced by these modelling choices increases with ZAMS masses. Also shown are models from the literature (Renzo et al. 2017; Farmer et al. 2016; Sukhbold & Woosley 2014; Limongi & Chieffi 2003, 2006; Woosley et al. 2002), although each adopts different modelling choices and definitions of the He core mass.  \n<!-- image -->  \n<!-- image -->  \n/circledot", '7.1. Progenitors that Do Not Pulse': 'The upper panel of Figure 48 shows the He core mass from Z =0.02 M ZAMS =30, 45, and 60 M glyph[circledot] models. The lower panel shows the mass location where the binding energy is 10 48 erg. Each ZAMS mass uses an exponential convective overshoot parameter f ov =0.004 applied at all convective boundaries, a mixing length α MLT =1.5, MLT++ enabled (see Paper II), and is run to the onset of core-collapse (infall velocity ≥ 1000 km s -1 ). We illustrate the variation in the He core mass and mass location where the binding energy is 10 48 erg from the effects of rotation, wind strength, and the wind schemes of Nieuwenhuijzen & de Jager (1990), van Loon et al. (2005), and de Jager et al. (1988).  \nTo estimate a BH mass from the structure at core collapse, we use the mass location where the binding energy integrated from the surface exceeds 10 48 erg. This is motivated by neutrinos removing ≈ 10 53 erg during core-collapse, reducing the gravitational mass of the core by ≈ 0.3 M glyph[circledot] . The outer part of the star responds to the sudden decrease in the gravitational field by driving a sound wave that steepens into a shock that unbinds some of the outer envelope (Coughlin et al. 2017). Mass with binding energy glyph[lessorsimilar] 10 47 erg is likely to be ejected (Nadezhin 1980; Lovegrove & Woosley 2013) while mass that is not ejected will likely become part of the BH. Figure 48 suggests that BH masses estimated in this simple way can be significantly larger than the final He core mass, and more sensitive to the assumed model parameters. For example, there is wide variation in the expected BH mass for the 60 M glyph[circledot] progenitor depending on choice of wind scheme and scaling factor, whereas modest rotation has a smaller effect.', '7.2. Pulsational Pair-instability Supernovae': 'Stars with M ZAMS glyph[greaterorsimilar] 60 M glyph[circledot] are expected to become dynamically unstable before core O depletion as e + e -pair production leads to regions where the adiabatic index Γ 1 ≤ 4/3 (Fowler & Hoyle 1964; Rakavy & Shaviv 1967). The ensuing collapse results in explosive O burning, with a variety of possible outcomes. Stars can produce PISNe where the energy injected from explosive O burning completely unbinds the star without leaving a compact remnant. Alternatively, stars can undergo a cyclic pattern of entering the pair instability region, contracting, burning, and expanding, leading to PPISNe.  \nIndividual pulses in a PPISN can remove a large fraction of the mass of the star at velocities of several thousand kms -1 , with the remaining material settling down into hydrostatic equilibrium at a lower central temperature than before the pulse. The star then contracts as it loses energy due to radiation and neutrino emission un-  \nFigure 49. Central temperature and density for models with metallicity Z = 0 . 001 undergoing core collapse (50 M glyph[circledot] ), PPISN (90M glyph[circledot] ), and PISN (140 M glyph[circledot] ).  \n<!-- image -->  \ntil it undergoes an additional pulse or collapses to form a BH. Depending on its initial mass, the time between pulses varies from a fraction of a year to millennia, with the outer ejected layers expanding to very low densities and becoming optically thin.  \nMESA currently cannot simultaneously follow the long term evolution of the bound core and the ejecta, making it necessary to remove the unbound layers from the stellar model. To do this we model individual pulses using both the Riemann solver hydrodynamics (Section 4), as well as the 1D treatment of the Rayleigh-Taylor instability (Section 5), until the star is approximately in hydrostatic equilibrium. We then relax a new stellar model using the methods described in Appendix B, such that it has the same mass, entropy, and composition profiles as the layers that remained bound in the hydrodynamical model. This model is then evolved assuming hydrostatic equilibrium until the onset of another pulse or the final core-collapse to a BH.  \nAs an example, we compute models at a metallicity Z = 0 . 001, using similar parameters as in Section 7.1. The van Loon scheme is used for low-temperature winds with a scaling factor of 0 . 4. In addition, convection is modeled as a time-dependent process by limiting changes in convective velocities as in Arnett (1969) and Wood (1974).  \nFigure 49 shows the evolution in the ρ c -T c plane during late burning stages for a 90 M glyph[circledot] model undergoing a PPISN, a 50 M glyph[circledot] model experiencing iron-core collapse, and a 140 M glyph[circledot] model producing a PISN. Although the center of the 90 M glyph[circledot] star does not evolve into the region where Γ 1 < 4 / 3, the outer layers of the CO core do. Coupled with enhanced neutrino losses from pair-  \nFigure 50. Late-time evolution of a M ZAMS = 90M glyph[circledot] star with a metallicity Z = 0 . 001 undergoing a PPISN. (upper panel) Evolution of the central temperature showing a zoomed-in region covering 2.2 days which contains the first two pulses, as well as an additional zoom-in covering 1.8 hours which shows the first pulse and its ring-down into hydrostatic equilibrium. (lower panel) Total mass of the star below the escape velocity, He core mass and CO core mass during pulsations. The online animated Figure shows the time evolution of these quantities and the interior structure of the star.  \n<!-- image -->  \nnihilation, this causes the star to collapse and undergo four distinct pulses before finally collapsing into a BH. At the onset of the first pulse, the star has a mass of 87 . 1 M glyph[circledot] , with a He core of 45 . 6 M glyph[circledot] and a CO core of 41 . 1 M glyph[circledot] . As shown in Figure 50 the first two pulses happen within two days of each other and they remove the entire H envelope. The remaining two pulses remove almost the entire He envelope, resulting in a final mass of 41 . 2 M glyph[circledot] when the star collapses into a BH.  \nFigure 51 shows key masses on a grid encompassing ZAMS masses for which PPISNe occur under our model assumptions. Our PPISN progenitors have He core masses in the range of 28 M glyph[circledot] -67 M glyph[circledot] , and no BHs with masses above 50 M glyph[circledot] are formed. These results are in broad agreement with Woosley (2017). However, Figure 51 shows that the range of ZAMS masses that result in a PPISN is significantly different to the one computed by Woosley (2017). This can be attributed to a different choice of input physics such as core overshooting, as well as a different initial metallicity.', '8. ENERGY ACCOUNTING IN STELLAR EVOLUTION': 'Figure 51. Total mass, He core mass, and CO core mass at carbon depletion for single stars at a metallicity Z = 0 . 001. The 50M glyph[circledot] model undergoes iron-core collapse, while the 140 M glyph[circledot] model experiences complete disruption through a PISN. All other models experience PISNe. For comparison, the models with no mass loss from Woosley (2017) at a metallicity Z = 0 . 0016 are also shown.  \n<!-- image -->  \n/circledot  \nPaper I describes the stellar structure equations and their implementation in MESA . In order to provide physically and numerically accurate solutions of these equations, it is often necessary to evaluate them in different ways depending on the details of the star being simulated. In particular, there are a number of different ways to formulate and evaluate the equations solved by MESA that encode local and global energy conservation. The goal of this section is to clarify the available options, discuss when and why they are used, and describe how various forms of energy are tracked and accounted for in stellar evolution. While in places this section reads like a tutorial, it is in fact the first time we have presented a detailed description of a complex and critical aspect of how MESA works, information that is important for intelligent use of this software tool.  \nIn Section 8.1 we describe the fundamental equations we are solving, and in Section 8.2 we describe choices associated with their numerical implementation. In Section 8.3 we describe the connection between the form of the energy equation typically used in stellar evolution calculations and the version used when the hydrodynamics options discussed in Section 4 are enabled. In Section 8.4 we clarify how the energy associated with ionization is included in MESA . In Section 8.5 we describe the numerical approach necessary to ensure that the latent heat associated with crystallization in a white dwarf (WD) is included in MESA . In Section 8.6 we discuss the difficulties introduced by the necessity to blend between  \ndifferent equations of state (EOS) as the thermodynamic conditions in the stellar interior change, and how MESA minimizes artifacts associated with these blends. In Section 8.7 we discuss the energy associated with gravitational settling.', '8.1. Fundamental Equations': "In the stellar structure equations (e.g., Cox & Giuli 1968; Kippenhahn et al. 2012), energy conservation is typically formulated by considering the energy flow in and out of a fluid parcel. In this Lagrangian picture, to understand how the energy of a fluid parcel is changing, we account for the specific (i.e., per unit mass) rate of energy injection into the parcel, glyph[epsilon1] , and the specific rate of energy flow through the boundaries ( ∂L/∂m ; L ( m ) is the luminosity profile and m the Lagrangian mass coordinate). The specific heating rate (D q/ D t ) for the parcel must then satisfy  \nD q D t = glyph[epsilon1] -∂L ∂m , (54)  \nwhere D / D t is the Lagrangian time derivative. Except in the case of hydrodynamics described in Section 4 (where a total energy equation is solved; see Section 8.3), the basic equation to be solved is always some form of Equation (54). By tradition, the negative of the lefthand side of Equation (54) is called glyph[epsilon1] grav .  \nThermodynamics relates the heating of material to the changes in its properties. The first law of thermodynamics states that the total heat added δQ for a parcel is  \nδQ ≡ d E + P d V, (55)  \nwhere E is the internal energy, P is the pressure, and V is the volume. Let N i be the number of particles of species i in the parcel. Then expanding E in terms of the independent thermodynamic basis variables ( S, V, N i ) yields the following thermodynamic identity:  \nd E + P d V = T d S + ∑ i µ i d N i , (56)  \nwhere S is the entropy, and T is the temperature. The sum runs over all species present, and  \nµ i ≡ ( ∂E ∂N i ) S,V (57)  \nis the chemical potential for species i .  \nThe number abundance of every species is defined with reference to the total number of baryons N B as Y i ≡ N i /N B . Denoting Avogadro's number by N A , the atomic mass unit is m amu = 1g /N A . The specific  \n(i.e., per unit mass) form of Equation (56) is then given by multiplying by the invariant N A /N B to find  \nδq ≡ d e + P d ( 1 ρ ) = T d s + ∑ i ( ∂e ∂Y i ) s,ρ d Y i . (58)  \nThe total baryonic mass density is ρ , so that 1 /ρ is the specific volume, and e and s are specific energy and entropy respectively. Local thermodynamic equilibrium (LTE) determines a unique solution for the ionization state of each isotope. Thus, composition is completely specified by a set of number abundances { Y i } for all nuclear isotopes.  \nEquation (58) is relativistically correct when the rest mass is included in the energy and the chemical potentials. Therefore, in principle, changes in nuclear rest masses due to nuclear reactions could be accounted for via this equation. However, in MESA , the energetic effects associated with composition changes due to nuclear reactions are not included in glyph[epsilon1] grav . Instead, these important terms are accounted for via glyph[epsilon1] nuc (the specific energy generation rate of nuclear reactions) which is evaluated separately and included as part of the local source term glyph[epsilon1] in Equation (54) (see Paper I).  \nIt is often convenient to specify compositions in terms of the baryonic mass fractions { X i } via the relation X i = A i Y i , where A i is the mass number for isotope i . Since rest-mass changes due to nuclear reactions are handled separately from glyph[epsilon1] grav , ρ and { X i } can be treated as independent basis variables without introducing any ambiguity into the chemical potential term in Equation (58). Some EOS options express the composition dependence in terms of aggregate quantities; examples include hydrogen abundance X , helium abundance Y , metallicity Z , average mass number ¯ A , and average atomic number ¯ Z .  \nThe value for glyph[epsilon1] grav can be computed beginning from either the left or right hand side of the equals sign in Equation (58). Usually, some form of the left hand side is used, but in Section 8.5 we will describe a case where it is more convenient to use the right hand side.", '8.2. Implementation': "Basic variables are those quantities directly calculated by MESAstar 's solver. Examples include velocity, radius, and the thermodynamic variables. MESA offers options for selecting ( ρ, T, { X i } ) or ( P gas , T, { X i } ) as the thermodynamic variables. The EOS routines calculate other thermodynamic quantities as a function of the chosen variables, e.g., e = e ( ρ, T, { X i } ). MESA solves the stellar structure equations implicitly, thus it is possible to approximate total time derivatives of any quantity calculated in the stellar model simply by differencing its  \nvalue at the start and end of a timestep. Therefore, one way to evaluate glyph[epsilon1] grav would be to directly calculate the time derivatives in Equation (58). Two possible versions of glyph[epsilon1] grav would then be  \n-glyph[epsilon1] grav = T D s D t + ∑ i ∂e ∂Y i D Y i D t , (59)  \nand  \nWhile simple to construct, the finite differences necessary to calculate these equations are often numerically problematic.  \n-glyph[epsilon1] grav = D e D t + P D D t ( 1 ρ ) . (60)  \nTo see the potential numerical issues, consider the implementation of Equation (54) using Equation (60) in cell k with mass dm k over a timestep δt . The derivative of a quantity D y/ D t is typically constructed as a finite difference of y over the timestep, so after integrating over the mass of zone k , we have  \n0 = ( glyph[epsilon1] k -e k, end -e k, start δt -P k 1 /ρ k, end -1 /ρ k, start δt ) dm k -( L k -L k +1 ) . (61)  \nThe implicit solver scheme in MESA attempts to reduce the residual from evaluating the right hand side of this equation below some tolerance.  \nWhile the implicit scheme in MESA may sometimes find acceptable results for an equation such as Equation (61), finite numerical precision can result in troublesome behavior for the time derivatives involving subtractions. In particular, over a small timestep where the change in e k or ρ k is small compared to the overall magnitude of these quantities, floating point arithmetic can suffer significant loss of precision. When energy scales arising from these types of finite difference derivatives are comparable to glyph[epsilon1] k , the implicit solver may be unable to converge to an acceptable solution.  \nTo avoid these problems, the equations can be cast in terms of derivatives that are not evaluated using subtractions. Such derivatives are available only for the basic variables, since the Jacobian matrix for an evolution step satisfying the equations of stellar structure in MESA is written in terms of the basic variables and their derivatives (see Paper I, Section 6.2 and Paper II, Section B.2 and Figure 47). For MESA , ρ and T serve as default variables.  \nModifying Equation (60) to take advantage of ρ as a basic variable yields  \n-glyph[epsilon1] grav = D e D t -P ρ Dln ρ D t , (62)  \nbut the change in e is still evaluated using subtraction. Another related form, obtained by application of mass continuity, is  \n-glyph[epsilon1] grav = D e D t + P ∂ ∂m ( v A ) , (63)  \nwhere v is the cell velocity and A is the area of the cell face. This is the form used in the artificial viscosity based hydrodynamics options described in Paper III.  \nExpanding the total derivative of energy and thus eliminating the subtraction motivates the following alternative forms. Expanding e in terms of its dependence on the basic variables ρ and T and dropping the dependence on composition gives  \n-glyph[epsilon1] grav = c V T Dln T D t + [ ρ ( ∂e ∂ρ ) T -P ρ ] Dln ρ D t , (64)  \nwhere c V ≡ ( ∂q/∂T ) ρ = ( ∂e/∂T ) ρ . One can also choose to expand e in terms of its dependence on P and T (dropping composition dependence) and then convert to a form given in terms of ρ instead of P to obtain  \n-glyph[epsilon1] grav = c P T [ (1 -∇ ad χ T ) Dln T D t -∇ ad χ ρ Dln ρ D t ] , (65)  \nwhere c P ≡ ( ∂q/∂T ) P and ∇ ad ≡ ( ∂ ln T/∂ ln P ) s . The derivation for this expression in terms of P and T is given in Chapter 4 of Kippenhahn et al. (2012), from which it is straightforward to obtain Equation (65) using χ T ≡ ( ∂ ln P/∂ ln T ) ρ and χ ρ ≡ ( ∂ ln P/∂ ln ρ ) T .  \nSince ρ and T are basic variables, the time derivatives appearing in Equations (64) or (65) involve no subtractions. Hence, solving Equation (54) with glyph[epsilon1] grav as defined by those two equations will not be susceptible to the same losses of numerical precision as other forms, at the cost of dropping the composition terms. Similarly, Equation (4.47) of Kippenhahn et al. (2012) will yield the same stability when P and T are used as basic variables. When P gas and T are selected as basic variables, the identification P = P gas + aT 4 / 3 allows writing  \n-glyph[epsilon1] grav = c P T [( 1 -4 ∇ ad P rad P ) Dln T D t -∇ ad P gas P Dln P gas D t ] . (66)  \nSection 4.5 in Kippenhahn et al. (2012) also shows how this local energy treatment of glyph[epsilon1] grav results in global energy conservation, including total gravitational potential energy from which the name glyph[epsilon1] grav is derived.  \nThe superior numerical stability of Equations (64)(66) comes at the cost of using derivative quantities such as c V and χ ρ . The Jacobian matrix of an implicit method thus requires the partial derivatives of c V  \nTable 4. Summary of glyph[epsilon1] grav Options  \nand χ ρ . An EOS must therefore be capable of returning the state functions P , e , and s along with their first derivatives (e.g., c V and χ ρ ) and second derivatives (e.g., ∂c V /∂T ).  \nAs noted above, Equations (64)-(66) drop the composition terms, which is justifiable if the derivatives ( ∂e/∂X i )(D X i / D t ) are negligible for each X i . Dropping composition terms is often justified in stellar evolution scenarios where timescales for these changes are very slow or their associated energies are negligible, such as MS burning where energy release from nuclear burning dominates any small changes in internal energy due to composition evolution over a single step (Kippenhahn et al. 1965; Garc'ıa-Berro et al. 2008). Making this assumption, MESA also offers an option for calculating glyph[epsilon1] grav in terms of a simplified form of Equation (59):  \n-glyph[epsilon1] grav = T D s D t , (67)  \nwhich drops composition dependence to offer an expression that is more convenient to evaluate.  \nHowever, even after composition dependence related to nuclear burning is accounted for with a separate glyph[epsilon1] nuc term as discussed in Section 8.1, other processes that change abundances (e.g., mixing) may be important. In cases where dropping these terms is not justifiable, it may be necessary to add a compensating local source term glyph[epsilon1] in Equation (54).  \nIn summary, MESA currently offers options for solving Equation (54) with glyph[epsilon1] grav defined in any of the ways given in Equations (60)-(67). Figure 52 schematically summarizes the relationships between these equations and Table 4 shows the inlist commands necessary for invoking each of these options. Usually, the superior numerical stability gained by using Equation (65) is to be preferred, and hence it is the MESA default, but users should be aware of the possibility that other forms may be necessary to capture important physics. One such case for Equation (67) is described in Section 8.5. An-  \nother is the artificial viscosity-based implicit hydrodynamics described in Paper III (see Section 4, Equation 41), where choosing Equation (63) helps ensure intrinsic energy conservation.", '8.3. Relationship to the Riemann Solver-Based Hydrodynamics Implementation': 'When using the Riemann solver-based hydrodynamics capabilities described in Section 4, MESA does not cast the stellar structure equations in terms of local heating as in Equation (54). Instead, it combines Equation (54) with the constraint of fluid momentum conservation to form a local total energy equation.  \nWe begin with the mass continuity equation,  \nD ρ D t = -ρ r 2 ∂ ∂r ( r 2 u ) , (68)  \nand the momentum equation,  \nD u D t = -1 ρ ∂P ∂r -∂ Φ ∂r , (69)  \nwritten in Lagrangian form and assuming spherical symmetry. The variable u is the radial velocity and Φ is the gravitational potential. The Lagrangian derivative operator is D / D t = ∂/∂t + u∂/∂r .  \nMultiplying Equation (69) by u gives  \nD D t ( 1 2 u 2 ) = -u ρ ∂P ∂r -u ∂ Φ ∂r . (70)  \nThe gravitational potential does not explicitly depend on time ( ∂ Φ /∂t = 0), so DΦ / D t = u∂ Φ /∂r . This implies  \nUsing Equations (54) and (60) we have  \nD D t ( 1 2 u 2 +Φ ) = -u ρ ∂P ∂r . (71)  \nD e D t -P ρ 2 D ρ D t = glyph[epsilon1] -∂L ∂m . (72)  \nAdding Equations (71) and (72) gives  \nD D t ( e + 1 2 u 2 +Φ ) = P ρ 2 D ρ D t -u ρ ∂P ∂r + glyph[epsilon1] -∂L ∂m . (73)  \nUsing mass continuity (Equation 68) this becomes  \nD D t ( e + 1 2 u 2 +Φ ) = -1 ρr 2 ∂ ∂r ( Pur 2 ) + glyph[epsilon1] -∂L ∂m . (74)  \nIn spherical coordinates  \n1 ρr 2 ∂ ∂r ( r 2 f ) = ∂ ( A f ) ∂m , (75)  \nwhere A = 4 πr 2 . Thus we arrive at the equation that MESA solves,  \nD D t ( e + 1 2 u 2 +Φ ) = glyph[epsilon1] -∂ ∂m ( L + P A u ) . (76)  \nFigure 53. The value of glyph[epsilon1] grav in the pure-H pre-MS model, evaluated over a region near the stellar surface that includes an ionization zone (where glyph[epsilon1] ion = 0). The solid colored lines indicate the individual energy terms. Their sum (dotted black line) agrees with the value calculated by MESA (solid gray line).  \n<!-- image -->  \n[  \n]  \nFigure 52. Schematic showing the relationships in Equations (59)-(65).  \n<!-- image -->', '8.4. Ionization': "The internal energy reported the EOS should include the energy associated with ionization 9 and molecular dissociation. The assumption of LTE specifies the ionization state given ( ρ, T, { X i } ). Since MESA does not regard a change in ionization as a change in composition, it is not necessary to include separate composition derivatives in glyph[epsilon1] grav in order to account for the energetic effects of changes in ionization state.  \nTo demonstrate a specific scenario where MESA accounts for ionization energy, we evolve a 1 M glyph[circledot] pre-MS model composed of pure H. We compare quantities calculated by MESA with other, simpler estimates. We calculate the thermal energy assuming a monatomic ideal gas,  \ne thermal = 3 N A k B T 2 µ . (77)  \nWe calculate the ionization energy for pure H as  \ne ion = (1 -f H ) N A E H + N A E H2 2 , (78)  \nwhere we assume the ionization fraction of H is given by the Saha equation. The variable f H represents the neutral fraction of H. The H ionization energy is E H = 13 . 6 eV and Equation (78) also includes the dissociation energy of molecular H ( E H2 = 4 . 52 eV) assuming that no H is in the molecular state.  \n<!-- image -->  \n-  \nglyph[negationslash]  \n/circledot  \nDuring the evolution, we record glyph[epsilon1] grav calculated by MESA using Equation (65). We also evaluate the quantity  \nD D t ( e thermal ) ︸ ︷︷ ︸ glyph[epsilon1] thermal + D D t ( e ion ) ︸ ︷︷ ︸ glyph[epsilon1] ion + P D D t ( 1 ρ ) ︸ ︷︷ ︸ glyph[epsilon1] PdV  \nthat separates out the thermal and ionization energy. In Figure 53 we compare these two approaches, making it clear that all three terms in the above expression play an important role. Their sum agrees with the MESA glyph[epsilon1] grav , indicating that each of these terms is accounted for in the MESA calculation.  \nFigure 54. Specific internal energy at a fixed Lagrangian coordinate in the pre-MS model (upper panel). The solid colored curves indicate the individual energy terms. The internal energy reported by MESA (solid gray curve) exceeds the thermal energy because of the ionization energy. The lower panel shows the neutral fraction of the H.  \n<!-- image -->  \nFigure 54 shows the history of the material at the Lagrangian coordinate ( M -m ) / M glyph[circledot] = 10 -5 . We plot e reported by the MESA EOS along with e thermal and e ion (calculated in same manner as above). At this location, the specific internal energy is dominated by the ionization energy. The lower panel of this figure shows the neutral fraction of H; towards the left of the plot, the H is fully neutral. In this region the ionization energy plateaus at the dissociation energy of molecular H (see Equation 78).  \nFor a star in hydrostatic equilibrium, the virial theorem states that  \n-1 2 ∫ M 0 Gm r dm + ∫ M 0 3 P 2 ρ dm = 0 . (79)  \nThe right term's integrand, 3 P/ (2 ρ ), is the specific thermal energy of an ideal monatomic ideal gas. Figure 55 shows the total internal energy and gravitational potential energy reported by MESA for the pure-H pre-MS model. On the same scale we show half the total potential energy plus the internal energy. This quantity is not zero; rather, by the virial theorem, it should sum to the non-thermal and non-ideal internal energy (e.g., the ionization energy). This value, recorded from the MESA model, agrees well with our estimate of the ionization energy. Also note that at early times the total energy of the  \nFigure 55. Total potential and internal energy in the pureH pre-MS model. The sum of half the total potential energy plus the internal energy (solid gray curve), which by the virial theorem should be the non-thermal internal energy, agrees well with our estimate of the ionization energy (dashed black curve). The deviation at glyph[greaterorsimilar] 10 5 yr is caused by non-ideal gas effects.  \n<!-- image -->  \n-  \nFigure 56. Specific (red) and cumulative (black) total energy (IE + PE) in the envelope of an AGB model ( M = 1 . 0 M glyph[circledot] , L = 4 . 97 × 10 3 L glyph[circledot] , T eff = 2 , 920 K, R = 276R glyph[circledot] ). This energy is positive in the envelope due to the inclusion of ionization energy in the internal energy reported by MESA .  \n<!-- image -->  \n/circledot  \nstar (internal + potential, not shown) is positive. The phenomenon of positive total energy when ionization energy is included also occurs for envelopes of stars on the asymptotic giant branch (AGB; Paczy'nski & Zi'oglyph[suppress]lkowski 1968). Figure 56 shows the total energy in the envelope of a 1 . 0 M glyph[circledot] MESA model on the AGB. This confirms that the ionization energy is included when MESA reports the total energy of a model.", '8.5. Latent Heat': 'Figure 57. Cooling for a 0 . 6 M glyph[circledot] C/O WD ( M H = 2 . 7 × 10 -5 M glyph[circledot] , M He = 1 . 5 × 10 -2 M glyph[circledot] ) with different treatments of the latent heat of crystallization. The default treatment smoothly injects the latent heat over the range 150 ≤ Γ ≤ 175. The right panel shows differences in cooling time (relative to the default shown in the left panel) required to reach a given luminosity for other treatments.  \n<!-- image -->  \n<!-- image -->  \n/circledot  \nPaper II discusses the inclusion of the latent heat of crystallization for long term WD cooling. Crystallization is a first-order phase transition that manifests in the PC EOS (Potekhin & Chabrier 2010) as an entropy discontinuity at a plasma coupling parameter of Γ = 175, and can be captured in stellar evolution with glyph[epsilon1] grav in the form of Equation (67). Since the publication of Paper II, controls have been added to MESA to allow smoothing out the injection of latent heat in glyph[epsilon1] grav over a user-specified range of Γ. By default, the range for crystallization is softened to 150 ≤ Γ ≤ 175 to avoid numerical difficulties with sudden energy injection associated with a sharp transition at Γ = 175. The controls allow for tightening this range for more precise timing on the occurrence of crystallization if necessary. Figure 57 shows the small impact on cooling time for a 0 . 6 M glyph[circledot] WDfrom spreading the latent heat over this broader range of Γ relative to a tighter phase transition for 174 ≤ Γ ≤ 176.  \nThe spreading of the phase transition is accomplished by calculating both the liquid and solid solutions within the PC EOS and linearly blending the entropy s and internal energy e over the specified range of Γ. With glyph[epsilon1] grav expressed in the form of Equation (67), the energy of the phase transition is captured as fluid elements smoothly traverse from liquid-phase to solid-phase. By default MESA automatically switches to using glyph[epsilon1] grav in the form of Equation (67) for Γ > 150. This choice ensures the capture of latent heat release.  \nTheoretical and observational work has suggested that crystallization in C/O mixtures may occur at higher Γ than the classical one component plasma value of Γ = 175 (Horowitz et al. 2007; Winget et al. 2009; Medin & Cumming 2010; Althaus et al. 2012). Our updated crystallization controls allow for investigating the effect  \non stellar evolution of crystallization at Γ ≈ 240. Figure 57 shows the potential effects on WD cooling times of varying the Γ for crystallization. Because the heating from crystallization is released very late in the WD evolution, its effects on cooling times are on the order of Gyr, and variations in crystallization treatment can lead to changes that are a significant fraction of this timescale.  \nThe composition terms in Equation (59) that were dropped to form Equation (67) are negligible as long as there is no mixing in the crystallization region. Phase separation may violate this assumption and require a modified treatment, but we do not consider this process here. Detailed phase-diagrams for crystallization and the possible associated phase-separation effects are not currently supported in MESA , so our investigation here is limited to the effects of crystallization as a function of a fixed Γ range.', '8.6. EOS Blending': 'As shown in Figure 1 of Paper I, MESA employs a patchwork of several EOSs to provide coverage of a maximal amount of ρ -T space. When blending from one EOS region into another, care is required to avoid introducing spurious contributions into glyph[epsilon1] grav . At high density, MESA blends from the Helmholtz EOS (HELM, Timmes & Swesty 2000) for Γ < 10 to the PC EOS (Potekhin & Chabrier 2010) for Γ > 20 by default. This default has been changed from the original default of 40 ≤ Γ ≤ 80 given in Paper I due to the optimal agreement between relevant quantities shown in Figure 59, as explained below. Overall, the two EOSs agree well on thermodynamic quantities in the blending region ( ∼ 1% for e and s ), but Figure 58 shows that the absolute magnitude of  \nthe disagreement can still be large enough to influence glyph[epsilon1] grav for a cooling WD when glyph[epsilon1] grav is expressed in the form of Equations (59)-(63).  \nThe left panel of Figure 58 indicates that typically the internal energy difference is ∆ e ∼ 10 15 erg g -1 , while c P T ∼ 10 14 -10 15 erg g -1 in the region of the blend. As a WD model cools, most of its ∼ 10 33 g of mass must eventually pass through this transition. If the energy equation is being solved in the form of Equation (61), ∼ 10 48 erg of spurious energy would be introduced into the model by EOS blending. Since much of this blending happens after the WD model has cooled to a luminosity of L glyph[lessorsimilar] 0 . 1 L glyph[circledot] , this extra energy corresponds to ∆ t glyph[greaterorsimilar] 100 Myr of extra WD cooling time.  \nThe default form of glyph[epsilon1] grav given in Equation (65) does not suffer from this spurious heating, since it is expressed in terms of thermodynamic derivatives from the EOS rather than e and s . For this form of glyph[epsilon1] grav , the differences between e or s do not directly enter the equations. Instead, changes in e with Dln ρ/ D t and Dln T/ D t are tracked with quantities such as c P and χ T , and Figure 59 shows that these agree well for the EOS blend region. Since the implementation of Equation (65) does not involve any derivatives constructed as finite differences, the fact that quantities such as c P agree to within a few percent guarantees that glyph[epsilon1] grav will be consistent across the blend, with no significant spurious energy injected due to blending. Crucially, the release of latent heat described in Section 8.5 requires switching to glyph[epsilon1] grav in the form of Equation (67) only for zones with Γ > 150, so both EOS blending and crystallization simultaneously receive appropriate treatments with different forms of glyph[epsilon1] grav in different stellar regions.', '8.7. Gravitational Settling': "Equation (65) for glyph[epsilon1] grav ignores changes in internal energy e due to composition changes. Garc'ıa-Berro et al. (2008) point out that a self-consistent evolutionary approach to WD cooling including the effects of 22 Ne settling requires accounting for composition changes due to element diffusion in glyph[epsilon1] grav . They adopt pure 12 C or 16 O core compositions with trace 22 Ne and no other isotopes. While this approach is useful for rigorous study of selfconsistent WD evolution with diffusion fully coupled to evolution, it is not well suited for a general treatment of realistic mixed core compositions.  \nMESA splits element diffusion into a separate step before the main structural solve, and hence diffusive effects are not included in glyph[epsilon1] grav . We ensure that the energy associated with 22 Ne settling is not included in glyph[epsilon1] grav by using Equation (65), and we compensate by including an extra heating term glyph[epsilon1] 22 in Equation (54). This term  \nis calculated using velocities saved from the element diffusion step as described in Section 3.5. Our results for the effects of 22 Ne settling on WD cooling agree well with Garc'ıa-Berro et al. (2008) and with Deloye & Bildsten (2002) who adopt a heating term similar to our approach.", '9. SUMMARY': "We explain significant new capabilities and improvements implemented in MESA since the publication of Paper I, Paper II, and Paper III. Progress in the treatment of convective boundaries (Section 2) and element diffusion (Section 3 and Appendix C) will improve studies of their impact on stellar evolution. Advances to MESA in implicit hydrodynamics (Section 4), approximation of 3D RTI effects (Section 5), and coupling with a public version of the STELLA radiative transfer instrument will enhance the modeling of Type IIP SN light curves from post-explosion to post-plateau (Section 6). We integrate these improvements with an exploration of PPISN and black hole formation models (Section 7). We describe energy conservation in MESA and demonstrate improvements relevant to WD cooling (Section 8). Upgrades to estimating the absolute magnitude of a model in a chosen passband (Appendix A), guidance on importing multi-dimensional models into MESA (Appendix B), and new MESA -based software tools (Section D) will strengthen research and education. Input files and related materials for all the figures are available at http://mesastar.org .  \nWe thank Dan Kasen for many discussions, SEDONA calculations and interactions that led to the implementation of the Sobolev optical depth approach to estimating the Fe II 5169 ˚ A line velocities in Type IIP SNe in MESA . Both Ondˇrej Pejcha and Stefano Valenti were generous in sharing and making available photometry, bolometric luminosities and Fe II line velocities of core collapse SNe. We acknowledge Annop Wongwathanarat for both communications regarding his 3D Rayleigh Taylor simulations of core collapse SNe and for providing detailed simulation results that enabled our comparison of the MESA implementation of the Duffell RTI modeling to 3D data. S.B. is grateful to L.B. and B.P. for hospitality at KITP where the public version of STELLA was finalized. Contributions to STELLA -tools by Petr Baklanov, Ulrich Noebauer, and Marat Potashov are gratefully acknowledged. We thank Carl Fields for initiating and developing MESA-Web . It is a pleasure to thank Conny Aerts, Warrick Ball, Matteo Cantiello, Barbara Castanheira, Sean Couch, Luc Dessart, Pavel Denisenkov, Jim Fuller, Falk Herwig, Raphael Hirschi, Christian Il-  \n-  \n<!-- image -->  \n-  \nFigure 58. Magnitude of the energy differences between the HELM and PC EOS for specific internal energy e (left) and entropy s (right) in a 50/50 C/O mixture. Dashed lines show the EOS blending boundaries for 10 ≤ Γ ≤ 20, and the solid black lines show representative profiles for a 1 . 0 M glyph[circledot] WD cooling from T eff = 26 , 000 K to T eff = 17 , 000 K.  \n<!-- image -->  \n<!-- image -->  \nFigure 59. Percent difference between HELM and PC EOS for c P (left) and χ T (right). Dashed lines show the EOS blending boundaries for 10 ≤ Γ ≤ 20, and the hatched region in the right panel shows the previous default blending range 40 ≤ Γ ≤ 80. The new EOS blending region disagree by only a few percent for these quantities.  \n<!-- image -->  \niadis, Thomas Janka, Sam Jones, Alexandra Kozyreva, Takashi Moriya, Marat Potashov, Eliot Quataert, Ken Shen, Matt Turk, Michael Wiescher, and Matt Wood for discussions. We also thank the participants of the 2016 and 2017 MESA Summer Schools for their willingness to experiment with new capabilities, as well as Mitchell Lachat and Jamie Lombardi for providing us with SPH simulations of stellar mergers.  \nThis project was supported by NSF under the Software Infrastructure for Sustained Innovation (SI 2 ) program grants (ACI-1339581, ACI-1339600, ACI1339606, ACI-1663684, ACI-1663688, ACI-1663696), and by NASA under the Theoretical and Computational Astrophysics Networks (TCAN) program grants  \n(NNX14AB53G, NNX14AB55G, NNX12AC72G). This research is funded in part by the Gordon and Betty Moore Foundation through Grant GBMF5076 and was also supported at UCSB by the National Science Foundation under grant PHY 11-25915. The work at Arizona State University was also supported by the NSF under grant PHY-1430152 for the Physics Frontier Center 'Joint Institute for Nuclear Astrophysics - Center for the Evolution of the Elements' (JINA-CEE). Support for this work was provided by NASA through Hubble Fellowship grant # HST-HF2-51382.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-  \n26555. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562; specifically the Comet cluster at the San Diego Supercomputer Center (SDSC) through allocation TG-AST150065. P.M. acknowledges support from NSF grant AST-1517753 to Vassiliki Kalogera at Northwestern University. A.T. is Research Associate at the Belgian Scientific Research Fund (F.R.S-FNRS). This research made extensive use of the SAO/NASA Astrophysics Data System (ADS).  \nSoftware: crlibm (de Dinechin et al. 2007), sedov (Kamm & Timmes 2007, http://cococubed.asu. edu/research\\_pages/sedov.shtml ), STARSMASHER (Gaburov et al. 2010; Lombardi et al. 2011, https: //jalombar.github.io/starsmasher/ ), Python avaliable from python.org, matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), ipython/jupyter (P'erez & Granger 2007; Kluyver et al. 2016), py mesa reader (Wolf & Schwab 2017), MESA -Web available at http: //mesa-web.asu.edu , MESA -Docker (Bauer & Farmer 2017, https://doi.org/10.5281/zenodo.1002851 ), pyMESA (Farmer 2017, https://doi.org/10.5281/ zenodo.846305 )", 'A. COLORS': "We describe MESA 's implementation of bolometric corrections (BCs) for use in estimating the absolute magnitude of a model in a user-chosen filter system. Note this is different than the colors reported by STELLA (Section 6), as the colors module uses pre-computed tables of BCs while STELLA solves the radiative transfer equations on-the-fly (Blinnikov et al. 1998).  \nThe absolute bolometric magnitude ( M bol ) of a star is defined, with reference to the solar absolute bolometric magnitude, as (Torres 2010):  \nM bol = M bol , glyph[circledot] -2 . 5 log 10 ( L/ L glyph[circledot] ) , (A1)  \nwhere M bol , glyph[circledot] is the absolute bolometric magnitude of the Sun, taken as 4.74 (2015 IAU Resolution B2). This can be transformed into the pass band dependent absolute magnitude, M X , for a nominal pass band X , via  \nM X = M bol -BC X , (A2)  \nBC X is the BC for pass band band X and accounts for the flux emitted outside of the wavelength range of the filter system. The derivation of a BC requires an atmospheric model of a star such that a stellar spectrum can be computed over all wavelengths, a computationally costly process. To prevent the requirement of actually having to generate a spectra at each time-step, we make use of pre-computed BC tables. These define the BC as a function of the stellar photosphere; T eff / K, log ( g/ cms -2 ) and the metallicity [M / H], and are derived from pre-computed grids of stellar atmosphere models, (see e.g., Kurucz 1970; Husser et al. 2013). Given the parameters at the stellar photosphere, we interpolate each set of BCs over log ( T eff / K), log ( g/ cms -2 ) and [M / H] using linear interpolation over nearest neighbors and without extrapolation for points outside of the table range.  \nWe provide two sets of pre-processed tables of BCs, though a user may provide their own. From Lejeune et al. (1998) we provide the Johnson-Cousins-Class bands UBVR c I c JHKLL ' M . This table provides the BCs over the parameter range 2 , 000 ≤ T eff / K ≤ 50 , 000, -1 . 02 ≤ log ( g/ cms -2 ) ≤ 5 . 0 and -5 . 0 ≤ [M / H] ≤ 1 . 0, with a variable sampling rate. Figure 60 shows the time evolution of the absolute magnitude of a 1 M glyph[circledot] star with the pass bands defined in Lejeune et al. (1998). We also provide a set of blackbody BCs for the pass bands UBVR c I c , over the range 100 ≤ T eff / K ≤ 50 , 000 in steps of 100 K. As these are blackbody corrections there is no g or [M / H] dependence.  \nThere are many other possibilities for other pass bands or classes of object (Fukugita et al. 1996; Girardi et al. 2002; Bessell 2011; Bessell & Murphy 2012). Thus the tables we provide are not a definitive set, but merely a reasonable starting point for modeling stellar objects. Other astrophysical objects like WDs, exoplanets, or SN light curves require calculating specialized tables. Users may provide BC tables defined in terms of T eff / K, log ( g/ cms -2 ) and [M / H].", 'B. MODEL RELAXATION': "To simplify the process of importing a model into MESA , we have developed simple relaxation routines that allow the construction of a starting model in hydrostatic equilibrium with specified profiles for composition, angular momentum, and entropy. Examples that motivate importing a model into MESA include multi-dimensional simulations of stellar mergers, common envelope evolution, and the effects of SN explosions on nearby companions.  \nFigure 60. Evolution of the absolute magnitude of a 1 M glyph[circledot] star for the bolometric magnitude and magnitude in the filter bands UBVR c I c JHKLL ' M .  \n<!-- image -->  \nThe relaxation process inputs include 1D profiles of composition and angular momentum. The process also requires either an entropy profile or the profiles of pairs of values ( ρ, T ), ( P gas , T ), or ( ρ, e ), from which MESA extracts the entropy using the eos module. Note that in the case where the entropy is not provided directly, the relaxed model will match the entropy computed by the eos module, but not neccesarily the input ( ρ, T ), ( P gas , T ), or ( ρ, e ) profiles. A good match for the input profiles depends on the input data corresponding to a model in hydrostatic equilibrium computed with an EOS that is consistent with MESA 's.  \nRelaxation is done via pseudo-evolution of a stellar model for which mixing, angular momentum transport, and changes in composition from nuclear burning are suppressed, while a quantity of interest is incrementally altered until it reaches the desired value up to a pre-defined tolerance. Throughout this relaxation process, hydrostatic equilibrium is enforced. The starting stellar model can be any MESA model with the required mass, and for most cases a ZAMS star at Z = 0 . 02 works well. The first two steps in the relaxation of a model fix the composition and angular momentum profiles. This is done by directly adjusting the variables for composition and angular momentum of each cell until the desired values are reached. Since the entropy is a derived quantity in MESA , the third step relaxes the entropy indirectly via the energy equation. This is achieved by adding a heating term that injects energy in regions where the entropy is below the target value, and removes energy in regions where the entropy is above the target value. This specific heating rate is  \nglyph[epsilon1] relax ( m ) = ( 1 -s ( m ) s target ( m ) ) e ( m ) τ , (B3)  \nwhere e ( m ), s ( m ) and s target ( m ) are the specific internal energy, current entropy, and target entropy respectively at the mass coordinate m . The timescale for the relaxation process is specified by τ . The value τ should be chosen to be small enough that energy transport is negligible during the pseudo-evolution. In practice, τ can be chosen to be orders of magnitude smaller than the dynamical timescale of the system.  \nWe verified that using the entropy, composition, and angular momentum profiles of a model computed with MESA as input, the relaxation procedure can reproduce the original model to within 0 . 1%. An example is provided in the test suite under the name relax\\_composition\\_j\\_entropy .  \nWe tested these relaxation routines using the outcome of a stellar merger computed with the STARSMASHER 10 SPH code (Gaburov et al. 2010; Lombardi et al. 2011), configured to use the MESA EOS. Two coeval non-rotating MESA models with ZAMS masses of 20 M glyph[circledot] and 15 M glyph[circledot] are evolved until the 20 M glyph[circledot] star reaches X c = 0 . 34. These models  \nFigure 61. Mass-weighted spherical averages of radial velocity and density from the STARSMASHER simulation of a head-on collision between non-rotating 20 M glyph[circledot] and 15 M glyph[circledot] stars. The dashed line shows the resulting density profile of a MESA model relaxed to the entropy and composition profile of the simulation.  \n<!-- image -->  \n/circledot  \nare then imported into STARSMASHER to simulate a head-on collision, such that the relative velocity of the two stars at infinity is zero. We find that 2 . 18 M glyph[circledot] of material is lost from the system due to the collision.  \nWe compute spherical mass-weighted averages of the composition, ρ , and e . These profiles are input into the MESA relaxation process, along with a zero angular momentum profile since the model is a head-on collision of non-rotating stars. Figure 61 shows that the relaxed model closely follows the input smoothed particle hydrodynamic (SPH) merger model in the central regions, though densities are ≈ 10% larger throughout the inner 25 M glyph[circledot] . Density differences of more than an order of magnitude are present in the outer layers. This is a consequence of these layers not being in hydrostatic equilibrium in the input SPH simulation. The MESA relaxation process matches the entropy rather than density profile of the SPH model assuming hydrostatic equilibrium as discussed above. The relaxed model corresponds to the final configuration if it contracts adiabatically, which is a good approximation as velocities in the SPH model are well below the local sound speed (Pan et al. 2013).", 'C. ELEMENT DIFFUSION IMPLEMENTATION DETAILS': "This appendix provides implementation details not contained in Section 3. Equations (2)-(4) and (9) give the full set of diffusion equations that must be solved to obtain diffusion velocities. For S total species in the plasma (including electrons), Equation (9) provides S -1 equations (one for each ion species), Equation (2) provides S equations (one for each species including electrons), and Equations (3) and (4) each provide one additional equation, for a total of 2 S + 1 independent equations. The 2 S + 1 unknowns are S diffusion velocities w s , S heat flow vectors r s , and the electric field E .  \nThe inputs provided from the MESA model are the number densities n s , temperature T , gradients of each of these quantities d ln n s / d r and d ln T/ d r , species mass in atomic units A s , species mean charge as an average ionization state ¯ Z s , and resistance coefficients K st , z st , z ' st , z '' st (defined in Equation 86 of Paper III). The coefficients are calculated as described in Section 3.3. Together with the mean ionization states, these are the key pieces of input physics that determine the diffusion of all ions. Extra acceleration terms g rad ,s for radiative levitation are either set to zero by default, or calculated as in Hu et al. (2011) when the option to include radiative levitation is enabled.  \nIn the spirit of Thoul et al. (1994), Equations (2)-(4) and (9) are grouped into a single matrix equation:  \nβ rad ,i + α i m p g + ν i k B T dln T d r + k B T ∑ j γ ij dln n j d r = ∑ j ∆ ij W j . (C4)  \nThe vectors capturing the driving terms are  \nα i =    n i A i i = 1 , . . . , S -1 , 0 i = S, . . . , 2 S +1 , (C5)  \nThe vector containing the unknowns is  \nν i =        n i i = 1 , . . . , S -1 , 5 2 n i i = S, . . . , 2 S -1 , 0 i = 2 S, 2 S +1 , (C6)  \nβ rad ,i =    -n i A i m p g rad ,i i = 1 , . . . , S -1 , 0 i = S, . . . , 2 S +1 . (C8)  \nFor i = 1 , . . . , S -1, the right hand side matrix of Equation (C4) is  \nglyph[negationslash]  \nγ ij =    n i δ ij i = 1 , . . . , S -1 , 0 i = S, . . . , 2 S +1 , (C7)  \nW j =        w j j = 1 , . . . , S, r j j = S +1 , . . . , 2 S, q e E j = 2 S +1 . (C9)  \nglyph[negationslash]  \n∆ ij =                    -∑ l = j K il j = i, K ij j = 1 , . . . , S and j = i, ∑ l = j K il z il A l / ( A i + A l ) j = i + S, -K i,j -S z i,j -S A i / ( A i + A j -S ) j = S +1 , . . . , 2 S and j = i + S, n i ¯ Z i j = 2 S +1 . (C10)  \nglyph[negationslash]  \nglyph[negationslash]  \n∆ ij =                                        5 2 ∑ l = j K i -S,l z i -S,l A l A i -S + A l j = i -S, -5 2 K i -S,j z i -S,j A j A i -S + A j j = 1 , . . . , S and j = i -S, -∑ l = j -S K i -S,l [ 3 A 2 i -S + A 2 l z ' i -S,l ( A i -S + A l ) 2 + 4 5 A i -S A l ( A i -S + A l ) 2 z '' i -S,l ] -2 5 K i -S,i -S z '' i -S,i -S          j = i, K i -S,j -S A i -S A j -S ( A i -S + A j -S ) 2 ( 3 + z ' i -S,j -S -4 5 z '' i -S,j -S ) j = S +1 , . . . , 2 S and j = i, 0 j = 2 S +1 . (C11)  \nFor i = 2 S ,  \nFor i = 2 S +1,  \nglyph[negationslash]  \n∆ ij =    n j A j j = 1 , . . . , S, 0 j = S +1 , . . . , 2 S +1 . (C12)  \nIndices i = 1 . . . S -1 capture the S -1 equations (9) for the ions. Indices i = S . . . 2 S -1 capture the S equations (2). Indices i = 2 S, 2 S +1 capture the two constraints in Equations (3) and (4).  \n∆ ij =    n j ¯ Z j j = 1 , . . . , S, 0 j = S +1 , . . . , 2 S +1 . (C13)  \nglyph[negationslash]  \nglyph[negationslash]  \nFor i = S, . . . , 2 S -1, the matrix terms are  \nglyph[negationslash]  \nFor a generic driving term that takes the form of an extra force f s on ions of species s , a term -n s f s appears on the left hand side of Equation (9). This can be accounted for in the matrix setup by adding another vector β f,i to the left hand side of Equation (C4) with the form  \nβ f,i =    -n i f i i = 1 , . . . , S -1 , 0 S, . . . , 2 S +1 . (C14)  \nOne such extra driving force that may be explored with MESA in the future is Coulomb separation in dense matter arising from non-ideal corrections for the ions (Chang et al. 2010; Beznogov & Yakovlev 2013; Diaw & Murillo 2016).  \nThe diffusion velocities are separated into two terms capturing the distinct effects of gravitational settling and ordinary diffusion in the tradition of Equation (11) of Iben & MacDonald (1985):  \nw i = w g i -∑ j σ ij dln C j d r , (C15)  \nwhere C j ≡ n j /n e following the notation of Thoul et al. (1994). These separate terms are constructed by inverting the matrix ∆ ij and then solving Equation (C4) for just one of α , β , ν , and γ ∗ ,j at a time on the left hand side. These results can then be linearly combined to construct w g i and σ ij such that the the full sum in Equation (C15) gives a solution that satisfies the complete set represented by Equation (C4).  \nWhen electrons become degenerate, we drop all S Equations (2) and set the S heat flow vectors to r s = 0. Equation (C4) then represents a system of just S +1 equations and the vectors and matrices simplify considerably, dropping all entries for indices i = S . . . 2 S -1 or j = S +1 . . . 2 S in the definitions given in Equations (C5)-(C13). To avoid discontinuities, we employ a blend that smoothly transitions between the diffusion velocity solutions over a range in η ≡ µ e /k B T , where µ e is the electron chemical potential. By default, the blend is centered around η ≈ 1, with user controls available to adjust the range of this blending region.", 'D. SOFTWARE INFRASTRUCTURE': 'Software is an integral enabler of observation, theory, and computation and a primary modality for realizing the discoveries and innovations expressed, for example, in the astronomy and astrophysics decadal surveys (e.g., National Research Council 1991, 2001, 2011). In this appendix we describe new software stacks at a variety of scales that enhance the research and education infrastructure.', 'D.1. Not A Number': 'Not a Number (NaN) is a numeric data type representing an undefined or unrepresentable value (e.g., Goldberg 1991; Hauser 1996). Examples include 0 / 0 and √ -1 in real arithmetic. In the IEEE 754 floating-point standard (IEEE 2008) there are two types of NaNs: quiet (qNaN) and signaling (sNaN). A qNaN propagates errors resulting from invalid operations or values without triggering a floating point exception. An sNaN precipitates an invalid operation exception whenever an attempt is made to use one as an arithmetic operand. The IEEE 754 standard requires qNaN as the default, while an sNaN can be used to support features such as filling uninitialized memory or other extensions to floating-point arithmetic.  \nNaN and infinity (INF) setting and testing routines are provided within the utils\\_nan.f90 file. A consistent set of interfaces allows for initializing scalars/arrays to NaN values, and testing for qNaN, sNaN, or INF values. Interface overloading allows handling of single, double or quad precision scalars or arrays of rank between 1 and 4. This module provides four generic interfaces. Logical function is\\_nan(x,signal) returns true if x contains NaNs and false otherwise. The optional logical argument signal determines whether qNaN, sNaN or both are tested for. Logical function is\\_inf(x) returns true if x contains INFs and false otherwise. Logical function is\\_bad(x) returns true if x contains NaN or INF values and false otherwise. Routine set\\_nan(x,signal) sets a scalar or array x to NaN values. The optional logical argument signal determines whether a qNaN or sNaN is set.  \nThe library framework of MESA is designed to be interoperable within other software ecosystems. For example, these NaN and INF interfaces are of potential interest to users of MESA or developers of similar software instruments.', 'D.2. MESA-Web': "Stellar evolution software instruments can be complicated to install and use, especially when the aim is primarily pedagogical (e.g., high-school or undergraduate courses). Motivated by the community's expressed need for a lower  \nbarrier to entry for education, a web-based interface to MESA was developed, MESA-Web at http://mesa-web.asu.edu . MESA-Web currently allows choices for the initial mass, metallicity, rotation, mass loss, nuclear reaction network, custom nuclear reactions rates, spatial and temporal resolution, and model output rate.  \nMESA-Web sends the user an email message when their job has completed that contains a URL of a zip file to download. The unzipped output directory contains a MESA history data file holding the time evolution of 57 quantities, as well as a series of MESA profile data files containing information on 56 quantities in each zone of the stellar model at discrete model numbers. Also included in the output is an MP4 formatted video containing a plot dashboard of the abundance profiles, Kippenhahn diagram, Hertzsprung-Russell Diagram, rotational profile, and temperature, density, and pressure profiles.  \nMESA-Web is currently hosted on a 4-core server at Arizona State University and allows jobs to run on a single core for 4 hours of walltime or until the model reaches iron core collapse. Launched in June 2015, MESA-Web has presently served more than 3000 models to over 600 different users at over 40 academic institutions. Efforts to expand MESA-Web 's capabilities include porting the service to a host with enhanced compute resources, simulating core-collapse supernova explosions (see Paper III) and light curves (see Section 6), and binary star evolution (see Paper III).", 'D.3. MESA -Docker': "Docker is a software technology designed to deploy and run applications by using 'containers'. Containers provide much of the virtualization power of traditional virtual machines while requiring far less resource overhead. This allows efficient packaging of an entire operating environment, with all of the necessary libraries and other dependencies for a large software tool such as MESA . The MESA -Docker package (Bauer & Farmer 2017) provides a solution that simplifies the requirements for locally running a full MESA installation with all capabilities available, with only minor overhead associated with running in a container. MESA -Docker will be useful for new users, students with educational projects, and Windows operating system users.", 'D.4. pyMESA': 'pyMESA (Farmer 2017) allows embedding of MESA modules into Python projects. pyMESA currently supports using the equation of state ( eos ), nuclear reaction ( rates ), neutrino ( neu ), atmosphere ( atm ), and opacity ( kap ) packages. This software infrastructure will be useful for users who want to use parts of MESA in their own Python software projects. As an example of these capabilities, Figures 58 and 59 were produced using the pyMESA eos interface to make direct calls to the MESA EOS routines.', 'D.5. MESAstar Model Optimization': 'The MESAstar test suite contains a sample case that shows how to use the simplex optimization algorithm (Nelder & Mead 1965) to find stellar models that minimize a specified χ 2 by automatically adjusting a variety of control parameters. 11 The χ 2 to be minimized can contain both pre-supplied and user-defined terms. Pre-supplied terms include T eff , L , R , g , surface Z , surface Y , and age. An easy-to-use framework allows the user to define other terms to include in the χ 2 . Control parameters include M , Z/X , Y , α MLT , and f ov . Other stellar evolution parameters can be easily added from the extensive set of controls in MESA . We provide a MESA test suite case using this new capability to calibrate a solar model. This can serve as a template for users wishing to use this method to search for models that match the observed properties of specific stars.', 'D.6. http://mesastar.org': 'Reproducibility is bedrock to scientific research. Provenance, as the term relates to software instruments (Van den Bussche & Vianu 2001; Carata et al. 2014), is the ability to record the full history of a result. Scientific research is generally held to be of good provenance when it is documented in detail sufficient to allow reproducibility. The MESA project facilitates provenance by the research community in four ways. One, by curating public releases of the source code, makefiles, test suite, and how the source code was compiled -GNU compilers are redistributed in the MESA Software Development Kit (see Paper II) -at http://mesa.sourceforge.net . Two, by providing bit-for-bit consistency for all results across all the supported platforms (see Paper III). Three, by supporting a user mailing-list to openly share knowledge (see the Manifesto in Paper I). Currently, over 12,000 messages are archived and searchable.  \nFour, by hosting a web-portal at http://mesastar.org to share MESA -oriented software contributions and reposit the MESA files ( inlist , run\\_star\\_extra.f , etc) that specify all the ingredients needed to reproduce a scientific result. Currently, http://mesastar.org offers over 120 MESA -oriented software contributions and inlist repositories.'}