We update the capabilities of the openknowledge software instrument Modules for Experiments in Stellar Astrophysics MESA. RSP is a new functionality in MESAstar that models the nonlinear radial stellar pulsations that characterize RR Lyrae Cepheids and other classes of variable stars. We significantly enhance numerical energy conservation capabilities including during mass changes. For example this enables calculations through the He flash that conserve energy to better than 0.001. To improve the modeling of rotating stars in MESA we introduce a new approach to modifying the pressure and temperature equations of stellar structure as well as a formulation of the projection effects of gravity darkening. A new scheme for tracking convective boundaries yields reliable values of the convective core mass and allows the natural emergence of adiabatic semiconvection regions during both core hydrogen and heliumburning phases. We quantify the parallel performance of MESA on currentgeneration multicore architectures and demonstrate improvements in the computational efficiency of radiative levitation. We report updates to the equation of state and nuclear reaction physics modules. We briefly discuss the current treatment of fallback in corecollapse supernova models and the thermodynamic evolution of supernova explosions. We close by discussing the new MESA Testhub software infrastructure to enhance source code development.
['stars: evolution', 'stars: general', 'stars: interiors', 'stars: oscillations:  including pulsations', 'stars: rotation', 'stars: variables: general', 'Astrophysics - Solar and Stellar Astrophysics']
Modules for Experiments in Stellar Astrophysics MESA Pulsating Variable Stars Rotation Convective Boundaries and Energy Conservation
{'No Header': "Typeset using L A T E X twocolumn style in AASTeX62  \nModules for Experiments in Stellar Astrophysics (MESA): Pulsating Variable Stars, Rotation, Convective Boundaries, and Energy Conservation  \nBill Paxton, 1 R. Smolec, 2 Josiah Schwab, 3, ∗ A. Gautschy, 4 Lars Bildsten, 1, 5 Matteo Cantiello, 6, 7 Aaron Dotter, 8 R. Farmer, 9, 10 Jared A. Goldberg, 5 Adam S. Jermyn, 1 S.M. Kanbur, 11 Pablo Marchant, 12 Anne Thoul, 13 Richard H. D. Townsend, 14 William M. Wolf, 15, 16 Michael Zhang, 17 and F.X. Timmes 16  \n1 Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA  \n2  \nNicolaus Copernicus Astronomical Center of the Polish Academy of Sciences, Bartycka 18, PL-00-716 Warszawa, Poland  \n3 Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA  \n4 CBmA, 4410 Liestal, Switzerland  \n5 Department of Physics, University of California, Santa Barbara, CA 93106, USA  \n6 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA  \n7 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA  \n8 Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA  \n9 Anton Pannenkoek Institute for Astronomy, University of Amsterdam, NL-1090 GE Amsterdam, The Netherlands 10 GRAPPA, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands  \n11 Department of Physics, SUNY Oswego, NY 13126, USA  \n12  \nDepartment of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA 13 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  \n14  \nDepartment of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA 15 Department of Physics and Astronomy, University of Wisconsin-Eau Claire, Eau Claire, WI 54701, USA 16 School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA 17 Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA", 'ABSTRACT': 'We update the capabilities of the open-knowledge software instrument Modules for Experiments in Stellar Astrophysics ( MESA ). RSP is a new functionality in MESAstar that models the non-linear radial stellar pulsations that characterize RR Lyrae, Cepheids, and other classes of variable stars. We significantly enhance numerical energy conservation capabilities, including during mass changes. For example, this enables calculations through the He flash that conserve energy to better than 0.001%. To improve the modeling of rotating stars in MESA , we introduce a new approach to modifying the pressure and temperature equations of stellar structure, and a formulation of the projection effects of gravity darkening. A new scheme for tracking convective boundaries yields reliable values of the convective-core mass, and allows the natural emergence of adiabatic semiconvection regions during both core hydrogen- and helium-burning phases. We quantify the parallel performance of MESA on current generation multicore architectures and demonstrate improvements in the computational efficiency of radiative levitation. We report updates to the equation of state and nuclear reaction physics modules. We briefly discuss the current treatment of fallback in core-collapse supernova models and the thermodynamic evolution of supernova explosions. We close by discussing the new MESA Testhub software infrastructure to enhance source-code development.  \nKeywords: stars: evolution - stars: general - stars: interiors - stars: oscillations - stars: rotation stars: variables: general', '1. INTRODUCTION': "Corresponding author: F.X. Timmes  \nfxtimmes@gmail.com  \n- ∗ Hubble Fellow  \nOne of the foundations upon which modern astrophysics rests is the fundamental properties of stars throughout their evolution. The advent of transforma-  \ntive capabilities in space- and ground-based hardware instruments is providing an unprecedented volume of highquality measurements of stars, significantly strengthening and extending the observational data upon which all of stellar astrophysics ultimately depends. For example, the Parker Solar Probe will provide new information on the flow of energy, structure, and dynamics of the closest star (Parker 1958a; Feng et al. 2010; Cranmer & Winebarger 2018; Gombosi et al. 2018) and the Daniel K. Inouye Solar Telescope will provide high temporal and spatial resolution with adaptive optics to reveal the nature of the the outer layers of the Sun (Parker 1958b; Snow et al. 2018; McComas et al. 2018).  \nThe exceptional precision of stellar brightness measurements achieved by the planet-hunting space telescopes Kepler/K2 (Borucki et al. 2010; Howell et al. 2014) and CoRoT (Auvergne et al. 2009) ushered in a new era in stellar photometric variability investigations. The Transiting Exoplanet Survey Satellite is building upon this legacy by surveying most of the sky in roughly month-long sectors covering four 24 · × 24 · areas from the ecliptic poles to near the ecliptic plane (Ricker et al. 2016). The mission will produce light curves for about 200,000 nearby late-type stars sampled at a 2 minute cadence to open a new era of stellar variability exploration (e.g., Dragomir et al. 2019; Huang et al. 2018; Ball et al. 2018; Shen et al. 2018; Wang et al. 2019). The Characterizing Exoplanets Satellite will complement these surveys by providing a unique, large sample of high precision photometric monitoring of selected bright target stars (Broeg et al. 2013; Gaidos et al. 2017).  \nThe Gaia Data Release 2, containing about 1.7 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. 2018a; Bailer-Jones et al. 2018; Lindegren et al. 2018; Luri et al. 2018). This stellar census will revolutionize a range of questions related to the origin, structure, and evolutionary history of stars in the Milky Way (e.g., Gaia Collaboration et al. 2018b,c; Riess et al. 2018). The infrared instruments aboard the James Webb Space Telescope (Gardner et al. 2006; Beichman et al. 2012; Artigau et al. 2014; Rieke et al. 2015) will search for the first and second generation stars (Rydberg et al. 2013; Kelly et al. 2018; Windhorst et al. 2018), assess how galaxies evolved from their formation (Zackrisson et al. 2011), observe the formation of stars from the initial stages of collapse onwards (Senarath et al. 2018), and measure the physical and chemical properties of stellar-planetary systems (Deming et al. 2009). The Laser Interferometer Gravitational-Wave Observatory and Virgo interferom-  \neters have demonstrated the existence of binary stellarmass black hole systems (Abbott et al. 2017a,b,c) and neutron star mergers (Abbott et al. 2017d,e,f), and will continue to monitor the sky with improved broadband detectors for gravitational waves from compact binary inspirals and asymmetrical exploding massive stars.  \nIn partnership with this ongoing explosion of activity in stellar astrophysics, community-driven software instruments are transforming how stellar theory, modeling, and simulations interact with observations (e.g., Turk et al. 2011; Foreman-Mackey et al. 2013; Ness et al. 2015; Choi et al. 2016; Astropy Collaboration et al. 2018). Modules for Experiments in Stellar Astrophysics ( MESA ) was introduced in Paxton et al. 2011 (Paper I) and significantly expanded its range of capabilities in Paxton et al. 2013 (Paper II), Paxton et al. 2015 (Paper III), and Paxton et al. 2018 (Paper IV). These prior papers, as well as this one, are software instrument papers that describe the capabilities and limitations of MESA while also comparing to other available numerical or analytic results.  \nThis instrument paper describes the major new advances to MESA for variable stars, numerical energy conservation, rotation, and convective boundaries. We do not fully explore the science results and their implications in this paper. The scientific potential of these new capabilities will be unlocked in future work via the efforts of the growing, 1,000-strong MESA research community.  \nMillions of variable stars have been discovered in the Milky Way and Magellanic Clouds, the Local Group (e.g., Optical Gravitational Lensing Experiment , OGLE, Udalski et al. 2015; MACHO Project , Alcock et al. 2003; Palomar Transient Factory , Soraisam et al. 2018) and beyond (e.g., Conroy et al. 2018). Figure 1 shows the broad classifications of these pulsating stars. Pulsating stars such as RR Lyrae and the brighter δ Cephei (the classical Cepheids) are common, and a strong direct relationship between their luminosities and pulsation periods established Cepheids (Leavitt 1908; Freedman et al. 2001; Majaess et al. 2009; Riess et al. 2016, 2018) and RR Lyrae in infrared bands (Clementini et al. 2001; Benedict et al. 2002; Klein et al. 2014; Muraveva et al. 2018a,b) as key distance indicators. New classes of variable stars are still being discovered: Blue LargeAmplitude Pulsators (BLAPs) are a new family of pulsating variable stars (Pietrukowicz et al. 2017). BLAPs are rare; only 14 variable stars are attributed by OGLE to this class after examining glyph[similarequal] 10 9 stars. They vary in brightness by glyph[similarequal] 20% on glyph[similarequal] 30 min timescales (Pietrukowicz et al. 2013). An important new addition to MESA is the capability to model radially-pulsating variable stars.  \nFigure 1. Classes of pulsating variable stars in the Hertzsprung-Russell (HR) diagram, including regions driven by the He II bump ( δ Ceph, δ Sct, RR Lyrae) and Fe bump ( β Ceph, SPB) in the opacity. Backslash ( \\ ) fills represent pressure modes and slash (/) fills represent gravity modes. The zero age main-sequence (ZAMS, black dashed curve) is labeled with the locations of selected masses. The classical instability strip for radial pulsations is shown by the gray dashed curves. Evolution of a 2.1 M glyph[circledot] MESA model (at Z = 0 . 02) from ZAMS to a white dwarf (WD) is shown by the purple curve. Figure design from Papics (2013).  \n<!-- image -->  \nNumerical energy conservation is rarely discussed by stellar evolution software instrument papers, or shown in science papers as part of establishing robustness of the solutions obtained with the software instrument. Yet stellar evolution calculations generally use low-order, implicit time integration with potentially poorly conditioned matrices whose matrix elements contain limitedprecision partial derivatives that can severely limit the quality of solutions. The cumulative effect of such errors can be substantial (Reiter et al. 1995). We implement a set of changes in MESA which, when applicable, can significantly improve the energy conservation properties of stellar evolution models at both global and local levels. This can reduce cumulative errors in energy conservation to 1% or less for applications such as the evolution of a 1 M glyph[circledot] model from the pre-main sequence to the end  \nof core He-burning or a core-collapse supernova from soon after explosion to shock breakout.  \nRotation modifies a star's structure (von Zeipel 1924a; Tassoul 2000; Maeder & Meynet 2000). We present a new approach in MESA for calculating the factors that modify the pressure and temperature equations of stellar structure within the shellular approximation. A rotating star is also oblate, with a larger radius at its equator than at its poles. As a result, the equator has a lower surface gravity and thus a lower effective temperature T eff (von Zeipel 1924b; Chandrasekhar 1933). Hence, the equator is 'gravity darkened', the poles 'gravity brightened', and this effect can play an important role in the classification of stars. The new extensions to MESA open a pathway for correcting T eff and L for aspectdependent effects.  \nStars transport energy by convection, whether within a core, an envelope, or throughout the interior. These convection regions showcase the interplay between composition mixing, gradients, and diffusion, and the transport of energy through the radial exchange of matter. It is necessary to ensure that convective boundaries are properly positioned because their placement can strongly influence the evolution of the stellar model (Gabriel et al. 2014; Salaris & Cassisi 2017, Paper IV). We implement an improved algorithm for correctly locating the convective boundaries and naturally allowing the emergence of adiabatic semiconvection regions during core H and He burning.  \nThe paper is organized as follows. Section 2 introduces a new capability to model large-amplitude radiallypulsating variable stars. Section 3 highlights energy conservation in MESA . Section 4 describes new rotation and gravity darkening factors, Section 5 explores a new treatment of convective boundaries, and Section 6 examines the parallel performance of MESAstar . Appendix A reports updates to the equation of state (EOS) and nuclear reaction modules. Appendix B details properties of the rotation factors. Appendix C discusses the current treatment of fallback in core-collapse supernovae (SN), and the thermodynamic evolution from massive star explosions. Appendix D introduces the MESA Testhub for source code development.  \nImportant symbols are defined in Table 1. Acronyms are defined in Table 2. Components of MESA , such as modules and routines, are in typewriter font e.g., eos .  \nTable 1 . Important symbols. Single character symbols are listed first, symbols with modifiers are listed second, and symbols for the RSP convection model are listed third. Some symbols may be further subscripted, for example, by c (indicating a central quantity), by a cell index k , or by species index i .  \nTable 1 (continued)Table 2 . Acronyms used in this paper.  \nTable 2 (continued)", '2. RADIAL STELLAR PULSATIONS': "Cepheids, RR Lyrae, and other classes of variable stars are observed to brighten and dim periodically. They can be modeled as radially symmetric, large amplitude, nonlinear oscillations of self-gravitating gas spheres. Software instruments for precision asteroseismology such as GYRE (Townsend & Teitler 2013; Townsend et al. 2018) model the small amplitude, linear oscillations of stars. Software instruments such as RSP , described below, are necessary to model the time evolution of large amplitude, self-excited, nonlinear pulsations over many cycles to produce luminosity and radial velocity histories that can be compared to observations.  \nEarly nonlinear radial pulsation models considered purely radiative envelopes (e.g., Christy 1964; Stellingwerf 1975; Castor et al. 1977; Aikawa & Simon 1983). Later, radiation hydrodynamic treatments followed with implicit adaptive grids (Dorfi & Drury 1987; Dorfi & Feuchtinger 1991). While these purely radiative models qualitatively reproduced light and radial velocity curves, it was clear that convection driven by partial ionization of H and He carries most of the flux in the envelopes of RR Lyrae and Cepheids. Prescriptions for coupling convection with pulsations were developed (e.g., Stellingwerf 1982; Kuhfuß 1986) that reside, with modifications, in modern software instruments (Bono & Stellingwerf 1994; Yecko et al. 1998; Koll'ath et al. 2002; Smolec & Moskalik 2008). Models from these software instruments can reproduce the overall morphology of light and radial velocity curves of classical pulsators (e.g., Feuchtinger et al. 2000; Marconi et al. 2015), features of specific objects (e.g., Keller & Wood 2006; Marconi et al. 2013; Smolec et al. 2013), and dynamical phenomena such as the Hertzsprung progression (e.g., Hertzsprung 1926a; Bono et al. 2000). Unsolved problems include doublemode pulsations (Koll'ath et al. 2002; Smolec & Moskalik 2010) and the cyclic modulations of RR Lyrae light curves (e.g., the Blazhko effect, Blaˇzko 1907; Szab'o et al. 2010). For background material we refer the reader to  \nGautschy & Saio (1995, 1996), Buchler (2009), and Marconi (2017).", '2.1. Radial Stellar Pulsations RSP': "RSP is a new functionality in MESAstar that models large amplitude, self-excited, nonlinear pulsations that stars develop when they cross instability domains in the HR diagram (see Figure 1). RSP is closely integrated with the MESA environment. Instead of calling the standard MESAstar routine to evaluate equations and solve for a new model using Newton-Raphson iterations (see Section 3), a separate routine does the same for RSP using a different set of equations and a different NewtonRaphson solver. The different equations include timedependent convection in a form appropriate for modelling nonlinear pulsations, and the different solver uses a band diagonal matrix approach since the equations as currently implemented do not fit into a three-block stencil needed for the standard block tridiagonal solver. Moreover, instead of calling the usual MESAstar routine to get a starting model, a separate routine creates an RSP model envelope that is consistent with the RSP set of equations. RSP uses the same MESA opacity and equation of state (EOS) modules, inlist structure, profile and history output files, photo files for saving and restarting runs, run star extras extensions, and hooks for using externally supplied routines.  \nRSP follows Smolec & Moskalik (2008), where the momentum and specific internal energy equations are  \nD u D t = -A ∂ ∂m ( p + p t + p av ) + U q -Gm r 2 , (1)  \nD e D t +( p + p av ) D V D t = -∂ ∂m [ A ( F r + F c )] -C , (2)  \nwhere D / D t is the Lagrangian time derivative. The generation of a specific turbulent energy, e t , is described by the one-equation Kuhfuß (1986) model  \nD e t D t + p t D V D t = -∂ ∂m ( A F t ) + glyph[epsilon1] q + C . (3)  \nThe latter two equations are added to give an equation for the specific internal and turbulent energies  \nD D t ( e + e t ) + ( p + p t + p av ) D V D t = -∂ ∂m [ A ( F r + F c + F t )] + glyph[epsilon1] q . (4)  \nDefinitions for all terms entering these equations are given in Table 1. RSP solves Equations (1), (3), and (4). The diffusion approximation is used for the radiative flux F r and its numerical implementation follows Stellingwerf (1975). Numerical implementation of the  \nsuperadiabtic gradient follows Stellingwerf (1982). All equations are discretized on a Lagrangian mesh.  \nSeveral quantities enter the convection model. For the momentum equation these are the turbulent pressure p t (Table 1 lists the relationship with the specific turbulent energy e t ) and viscous momentum transfer rate U q . For the turbulent energy equation these are the work done by turbulent pressure, the divergence of the turbulent flux F t , and the viscous energy transfer rate glyph[epsilon1] q . The convective coupling term C = S -D -D r appears with opposite sign in the internal and turbulent energy equations. Generation of the turbulent energy is driven by the source function S , while turbulent dissipation D and radiative cooling D r contribute to its decay. Radiative cooling of convective eddies follows Wuchterl & Feuchtinger (1998). Details of the turbulent convection model are discussed in Kuhfuß (1986), Wuchterl & Feuchtinger (1998) and Smolec & Moskalik (2008).  \nTable 3. Free parameters of the RSP convection model, their base values, and associated MESA controls that multiply the base values.  \nThese terms in the convection model depend on the free parameters listed in Table 3. If radiative cooling and turbulent pressure are neglected, the time-independent version of the Kuhfuß (1986) convection model reduces to standard mixing length theory provided base values are used for α s , α c and α d (associated controls set to 1). Base values for α p and γ r follow Yecko et al. (1998) and Wuchterl & Feuchtinger (1998), respectively. Experience suggests α t glyph[similarequal] 0 . 01, α m glyph[lessorsimilar] 1, and α glyph[lessorsimilar] 2 are useful starting choices.  \nPeriods of pulsation modes depend weakly on the values of these free parameters. Pulsation growth rates and light and radial velocity curves are, however, sensitive to the free parameters. Calibration with multiple observational constraints is unlikely to yield a unique set of parameters that gives satisfactory results across the HR diagram for all pulsation modes. We stress that  \nparameter surveys are an essential part of any science application of RSP .  \nIn Equations (1) and (2), p av is the artificial viscosity pressure (Richtmyer 1957) for numerically handling shocks that may develop during pulsations. We adopt the Stellingwerf (1975) two-parameter formulation as the default. The Tscharnuter & Winkler (1979) artificial pressure-tensor form, which was implemented in Paper III, can also be used in RSP .  \nThe numerical scheme to solve discrete versions of the equations is based on the intrinsically energy conserving method given by Fraley (1968). Details of the numerical implementation, along with RSP 's lineage (Stellingwerf 1975; Kovacs & Buchler 1988), are discussed in Smolec & Moskalik (2008). During the nonlinear integration, u = 0, L = constant and e t = 0 at the inner boundary (See Figure 2). The latter condition holds also for the outermost boundary, i.e., outermost cell is radiative. External pressure is fixed and zero by default.  \nFigure 2. Grid structure in RSP . The inner boundary at the base of the static envelope is defined by a chosen temperature ( RSP T inner ). The model surface has a fixed temperature T outer , derived from T eff , and pressure ( RSP Psurf ). The anchor temperature ( RSP T anchor ) is usually located where H ionizes, shown by the blue curve. The envelope is divided into N Lagrangian mass cells ( RSP nz ). Between the anchor and the surface are N outer cells ( RSP nz outer ), each with a constant mass. Between the inner boundary and the anchor the mass of each cell increases.  \n<!-- image -->", '2.2. RSP in Action': "RSP performs three operations: building an initial model; conducting a linear non-adiabatic (LNA) stability analysis on that model; and integrating the timedependent nonlinear equations.  \nSince the energy density of radial pulsations drops rapidly going inward from a star's surface, a full stellar model reaching to the center is frequently not necessary. The use of RSP is currently restricted to cases in which pulsations are determined by the structure of the envelope and are independent of the detailed structure of the core. RSP begins by building a chemically homogeneous envelope from given stellar parameters ( M , L , T eff , X , and Z ). These parameters can be freely chosen and need not originate from a MESAstar model. It is not yet possible to directly import an envelope from MESAstar into RSP primarily because of the different treatments of convection (a version of mixing length theory in MESAstar versus detailed time evolution of turbulence in RSP ). Tighter integration of MESAstar and RSP is a future project.  \nSpecifications for the initial model include the number of cells and the temperature at the base (see Figure 2). This inner boundary temperature is defined by a chosen temperature ( RSP T inner glyph[similarequal] 2 × 10 6 K) that should be set hot enough so that the eigenvector amplitudes generated in the following stability analysis go to zero, and cool enough to exclude regions of nuclear burning and justify the assumption of chemical homogenity. The model is divided into inner and outer regions at a specified anchor temperature. In the outer region, cells have the same mass; in the inner region cell masses grow by a constant factor so that the innermost cells are significantly larger than the ones at the surface. The anchor temperature should be in the part of the model driving the pulsations. For example, for pulsations in the classical instability strip a value of T anchor =11,000K is typical. In the case of Z -bump pulsations a higher temperature would be appropriate. Proper choice of the number of outer cells and placement of the anchor are necessary to ensure that the driving region is well resolved.  \nThe initial model builder iteratively constructs an envelope in hydrostatic equilibrium that satisfies the RSP equations. Starting from the outer radius determined by L and T eff this process involves selection of a cell mass to be used in the outer part of the envelope and a scale factor that is used to progressively increase cell masses in the inner region. Those choices must match the desired number of cells, both N and N outer , and also satisfy the surface boundary conditions and the required temperatures at the anchor location and at the inner boundary. The model builder is a complex multistage iterative procedure that works well for the range of cases presented in the following but may fail when applied outside of that range.", '2.2.2. Stability Analysis': 'The LNA analysis is performed on the initial model using a full linearization of the RSP equations (for details see Smolec 2009). These include time-dependent convection, moving beyond the frozen-in convection approximation made in software instruments like GYRE . This yields the eigenmodes, periods, and growth rates. The eigenvectors are used to perturb the initial model for the time evolution.  \nFigure 3 shows amplitudes of radial displacements and differential work for the first three eigenmodes of the classical Cepheid model in the MESA test suite . A common resolution for exploratory model surveys, N = 150 and N outer = 40, is adopted. For T glyph[greaterorsimilar] 5 × 10 5 K, the displacements and differential work of all three radial eigenmodes are negligible, indicating that the extent of the computational domain is sufficient.  \nFigure 3. RSP LNA analysis of a classical Cepheid model. Shown are the displacement amplitudes (black) and the differential work (red) done by the lowest three radial eigenmodes. The thickest curves are for the fundamental (F) mode, medium thickness for the first overtone (1O) mode and the thinnest curves for the second overtone (2O) mode. The dots indicate the cell locations. The grey areas show the extent of the convection zones around H ionization, first He and second He ionization, respectively.  \n<!-- image -->', '2.2.3. Evolution in the Linear Regime': 'The initial static model is perturbed with a linear combination of the velocity eigenvectors of the three lowest order radial modes. More specifically, the velocity eigenvectors are scaled to have a surface value of 1. RSP fraction 1st overtone and RSP fraction 2nd overtone multiply the 1O and 2O eigenvectors, respectively. The F-mode eigenvector is then multiplied by (1 RSP fraction 1st overtone -RSP fraction 2nd overtone ). The linear combination  \nFigure 4. Comparison of the fractional growth rate Γ during the initial cycles of the time integration. Horizontal lines show the LNA predictions. An RR Lyrae model ( M =0.65M glyph[circledot] , L =45L glyph[circledot] , T eff =7,100K, X =0.75, Z =0.0014) was initialized with a 0.1 km s -1 amplitude pure F-mode (circles), 1O-mode (squares), or 2O-mode (triangles) and evolved.  \n<!-- image -->  \nof these three scaled eigenvectors is then multiplied by the surface velocity RSP kick vsurf km per sec .  \nThe time integration commences with a constant time step ( RSP target steps per cycle ) and continues for a specified number of pulsation cycles ( RSP max num periods ). A new cycle begins when the model passes through a maximum radius. Controls allow filtering out secondary maxima in the radial velocity curve.  \nFigure 4 shows Γ the fractional growth of the kinetic energy per pulsation period near the start of a time integration, where  \nΓ = 2( E i +1 k , max -E i k , max ) / ( E i +1 k , max + E i k , max ) , (5)  \nand E i k , max is the maximum kinetic energy of the envelope during pulsation cycle i . Agreement between these three time integrations and the corresponding LNA analyses is satisfactory. Similarly, the pulsation periods match the linear values during the low-amplitude phase of development. Consistency between the time integrations and LNA analyses form the basis for interpreting the nonlinear results.', '2.2.4. Different Perturbations, Different Periods': 'Which of the perturbed modes attains large-amplitude pulsations in the nonlinear regime may depend on the initial conditions (e.g., Smolec 2014). Figure 5 shows the results of longer time integrations for the RR Lyrae model shown in Figure 4. The upper triplet of panels, case (a), is for a 1O-mode initialization with a 4.5 km s -1 amplitude. The middle triplet panel, case (b), is for a F-mode initialization with 4.5 km s -1 amplitude. The lower triplet, case (c), is for an F-mode initialization with a 9.5 km s -1 amplitude.  \nFigure 5. Fractional growth rate Γ, period P , and amplitude of radius variation ∆ R during 4,000-cycle integrations of the same RR Lyrae model as in Figure 4 with three different initial conditions labelled (a), (b), and (c).  \n<!-- image -->  \nFor case (a), the pulsations converge towards a single, 1O-mode pulsation. After a glyph[similarequal] 500 cycle transient phase the pulsation period and radius amplitude barely change and Γ glyph[similarequal] 0. For case (b), the model has not converged to a single-periodic mode after 4 , 000 cycles. Despite the pure F-mode initialization, at glyph[similarequal] 300 cycles the pulsation switches toward the 1O mode. This does not prove the model cannot pulsate in the F-mode, as case (c) demonstrates. After a transient phase with beating F and 1O modes, the 1O mode decays and the single-periodic Fmode pulsation grows to saturation.  \nFigure 5 is an example of two different single-mode solutions whose selection depends on the initial conditions. Two stars can have the same physical parameters but pulsate in different modes depending on their evolutionary history.  \nTable 4. Convective parameter sets referred to in the text as A, B, C, or D. Note that the controls multiply base values (see Table 3).  \nFigure 6. Bolometric light curve of the T eff =6,000K Cepheid model of with convective parameters of set A but with varying eddy-viscous dissipation α m .  \n<!-- image -->', '2.2.5. Convection Parameter Sensitivity': 'The final state in the nonlinear regime is usually a single-periodic oscillation. The shape of the light and radial velocity curves may depend on the values of the eight free parameters listed in Table 3. In Table 4 set A corresponds to the simplest convection model. Set B adds radiative cooling, set C adds turbulent pressure and turbulent flux, and set D includes these effects simultaneously. The parameter α m has little effect on the shape of the light curve but strongly affects its amplitude. Figure 6 shows the effect of varying α m on the T eff =6,000K Cepheid model. The free parameter α m may thus be used to match the observed amplitude. For sets A-D, α m was adjusted so that models with different sets have similar amplitudes.  \nFigure 7 shows the effect of these parameter sets on the shapes of the bolometric light, photosphere radial velocity, and radius variation curves for a saturated Fmode RR Lyrae and two saturated F-mode classical Cepheid models. The pulsation periods, radial velocity curves, and radius variation curves show only small differences. For the RR Lyrae models, there are differences in the fine structure of the light curves. For example, the bump before minimum light is weaker when turbulent pressure and turbulent flux are included (sets C and D). The shape of the light curve near maximum light also differs for both the RR Lyrae and Cepheid models.  \nFigure 8 shows the convective luminosity profiles for the models of Figure 7. Depending on pulsation phase, one convective region (darker hues) extends from the surface cells down to cell glyph[similarequal] 90. This convective region is associated with partial ionization of H and He. Another convective region (lighter hues) lies deeper in the envelope, centered at cell glyph[similarequal] 110, and is associated with the second ionization of He. In most of the models these two convective regions merge at pulsation phase glyph[similarequal] 0.5 during maximum contraction when both convective regions are at their strongest and most extended. In the cooler models, the first convective region carries nearly all of the luminosity throughout a pulsation cycle. In the hotter models, this convective region becomes very weak at pulsation phase glyph[similarequal] 0.8 (before maximum expansion) and is barely resolved as it progresses deeper into the envelope. This behavior is pronounced in the RR Lyrae models with radiative cooling (sets B and D), as cooling contributes to damping the turbulent energy and hence the near disappearance of the convective region.', '2.2.6. Artificial Viscosity Sensitivity': 'There are two parameters, α cut and C q , that control the Stellingwerf (1975) artificial viscosity, entering into the definition of p av (see Table 1). Figure 9 shows  \n∆  \nFigure 7. Bolometric light (upper panel), radial velocity (middle panel) and radius curves (lower panel) for an RR Lyrae F-mode model (left; M =0.65M glyph[circledot] , L =45L glyph[circledot] , T eff =6,700K, X =0.75, Z =0.0014) and two F-mode classical Cepheid models ( M =4.15M glyph[circledot] , L =1,400L glyph[circledot] , X =0.73, Z =0.007) at T eff =6,000K (middle panel) and T eff =5,700K (right panel). The mass and luminosity for the Cepheid models are close to the values derived for OGLE-LMC-CEP-227 (Pilecki et al. 2018). Each curve corresponds to a set of convective parameter values listed in Table 4. The mean magnitude of the bolometric light curves is set to zero. Light curves are vertically offset by 0.3 mag, radial velocity curves by 10 km s -1 and radius curves by 0 . 2 R glyph[circledot] (RR Lyrae) or 1 R glyph[circledot] (Cepheids).  \n<!-- image -->  \nthe effect of these parameters on light curves for the T eff =6,000K Cepheid model. The value of C q plays a minor role. For C q = 4, the upper panel shows a different light curve develops only if α cut =0.0, corresponding to artificial viscosity acting for very small compressions, which leads to excessive dissipation that quenches the pulsation amplitude. For α cut ≥ 0.01, the light curves are similar and roughly have the same pulsation amplitude. When α cut =0.1, artificial viscosity turns on only for strong shocks, seen as the wiggle on the ascending branch of the light curve. This choice ( α cut =0.1) numerically captures shocks without excessive dissipation, barely affecting the light curve shape and amplitude. While an artificial viscosity modifies the velocity structure in the envelope at each epoch, we find that these differences are smaller than the differences for the bolometric light curves. The lower panel shows that the Tscharnuter & Winkler (1979) form of artificial viscosity yields light curves with the same amplitude and qualitatively the same shape. Small differences are apparent at a shock-prone phase shortly before the maximum brightness.', '2.2.7. Spatial and Temporal Sensitivity': 'Figure 10 shows the sensitivity of the bolometric light curve to the total number of cells N , number of cells above the anchor N outer , anchor location T anchor , and inner boundary location T inner for an RR Lyrae model ( M =0.65M glyph[circledot] , L =50L glyph[circledot] , T eff =7,000K, X =0.75, Z =0.0014) with convective parameter set A.  \nFor classical pulsators, T inner is typically placed at 2 × 10 6 K, and common choices for T anchor are 11,000 K or 15,000 K. The upper panel shows that the light curves are weakly sensitive to the choice of T inner and T anchor for this RR Lyrae model. Light and radial velocity curves are usually the most sensitive to N and N outer . The lower panel shows this effect is small for this RR Lyrae model. Section 2.4.4 shows a case with a much larger sensitivity.  \nThe default value of 600 time steps per pulsation cycle works well for most cases, but smaller time steps are recommended for models that include radiative cooling, turbulent pressure, turbulent flux, or develop violent pulsations (e.g., the chaotic models of Section 2.4.3). We  \nFigure 8. Evolution of convective luminosity L c /L for the models shown in Figure 7. Cell number on the y-axis serves as a spatial coordinate, with cell 0 marking the stellar surface. The radiative interior of each model is not shown.  \n<!-- image -->  \nstress that there is no unique choice of grid or time step that will work for all applications or guarantees convergence. All nonlinear modeling of variable stars should be accompanied by sensitivity and convergence tests.', '2.3. Current Limitations and Plans for the Future': 'RSP in its present form covers most of the classical instability strip including δ Cepheids, RR Lyrae, High Amplitude δ Scuti and SX Phoenicis stars (see Figure 1), where a single or just a few dominant radial modes are observed. RSP also has applications outside of the classical instability strip as we show below for BLAPs. For stars close to the main sequence, linear growth rates are very small and thus, as we show below, long time integrations are necessary to approach full-amplitude nonlinear pulsations.  \nRSP is currently of limited use for strongly nonadiabatic pulsations with large L/M ratios, including Luminous Blue Variables, Mira-type variables, and type II Cepheids. For the latter, only the shortestperiod BL Her class variables can be reliably modeled (see Section 2.4.3). For the longer-period classes of W Vir and RV Tau variables, either static envelopes  \ncannot be constructed or nonlinear integrations break down at the onset due to violent relaxations of the outermost layers. In the extended envelopes of these variable stars the radiation-diffusion approximation is inadequate due to the low optical depth. The inclusion of pulsation-driven mass loss may also be necessary to study pulsations of these variable stars (Smolec 2016).  \nInclusion of turbulent pressure and flux may lead to convergence difficulties when constructing the static initial envelope. Cooler stellar envelopes with higher Mach number convection are also numerically more difficult than hotter envelope models. These difficulties are rooted in the static grid structure shown in Figure 2. Future developments of RSP should include a more versatile initial model builder, adaptive remeshing during the time integration, and a radiation-hydrodynamic treatment of the radiative energy and flux.', '2.4. Applications of RSP': 'We now apply RSP to variable stars. Examples include the light curves of classical pulsators, modeling of spe-  \nFigure 9. Sensitivity of the bolometric light curve shape of Cepheid model on artificial viscosity: α cut (upper panel) in the Stellingwerf (1975) formulation, and α cut in the Tscharnuter & Winkler (1979) formulation (lower panel). The light curve with C q =4.0 and α cut =0.1 is shown as a gray curve.  \n<!-- image -->  \ncific objects, and models for the dynamics of modulated or chaotic pulsations.', '2.4.1. RR Lyrae variables': "We consider two sequences of RR Lyrae-type models. The first sequence has M =0.65M glyph[circledot] , L =45L glyph[circledot] and [Fe/H] = -1.0 ( X =0.75, Z =0.0014) with T eff varying in 100 K steps for convective sets A-D. Figure 11 shows a gallery of I -band light curves from this sequence. The upper panel shows F-mode pulsators (commonly known as RRab stars) and the lower panel shows 1O-mode pulsators (known as RRc stars). The latter have smaller amplitudes and are less nonlinear in shape (i.e., more sinusoidal) than the F-mode light curves. Models with different convective settings differ the most near minimum and maximum light. For example, F-mode models with convective set B develop a bump preceding minimum light that is absent in the light curves with convective set D. On the other hand, F-mode models with cooler T eff from convective set D develop broad, doublepeaked light curve maxima that are absent in models from convective set B.  \nTo compare the overall morphology of I -band light curves from these sequences with OGLE observations, we perform a Fourier decomposition of the synthetic  \nFigure 10. Sensitivity of the bolometric light curve of an F-mode RR Lyrae model to the grid, labeled as N/N outer /T anchor /T inner . The light curve for the default grid is shown by the gray curve. The upper panel shows the effects of different T anchor / T inner . The lower panel shows the effects of different N/N outer combinations.  \n<!-- image -->  \nlight curves  \nI ( t ) = A 0 + ∑ k A k sin(2 πkft + φ k ) , (6)  \nwhere f is the pulsation frequency, and A k and φ k are amplitudes and phases, respectively. We then construct the amplitude ratios R k 1 and epoch-independent phase differences ϕ k 1 (Simon & Lee 1981):  \nR k 1 = A k A 1 , ϕ k 1 = φ k -kφ 1 . (7)  \nObservationally derived values of R k 1 and ϕ k 1 are taken from the OGLE catalog (Soszy'nski et al. 2014) and shown in Figure 12 by gray dots for RRab (Fmode) stars and blue dots for RRc (1O-mode) stars. The observations show that the Fourier parameters follow progressions with pulsation period, traced by the highest density of data points, but with significant scatter. Fourier parameters from the model I -band light curves are shown with colored symbols. Left panels are for the first sequence of models. Right panels are for the second sequence, which has M =0.65M glyph[circledot] , convective set B, T eff varying in 100 K steps, and either L =(40, 45, 50)L glyph[circledot] at [Fe/H] = -1.0 or L =45L glyph[circledot] at [Fe/H] = -1.5.  \nThe left panels show that F-mode pulsators with convective sets A and B progress similarly. Models with  \n∆  \n[mag]  \nI  \n∆  \nFigure 14 shows I -band light curves and radial velocity curves for the ∆ T eff =300K models. Due to the shift  \n<!-- image -->  \n-0.5  \n0.0  \n0.5  \n1.0  \n1.5  \n2.0  \n0.0  \n0.5  \n1.0  \n1.5  \n2.0  \nPulsation phase  \nFigure 11. I -band light curves of F-mode (upper panel) and 1O-mode (lower panel) pulsators across the instability strip for M =0.65M glyph[circledot] , L =45L glyph[circledot] , [Fe/H] = -1.0, and convective sets B (blue) and D (red). Light curves are labelled with their T eff and period, and offset vertically to facilitate comparisons (by 0.5 in the upper panel and 0.4 in the lower panel).  \nconvective sets C and D also progress similarly but are qualitatively different than those with convective sets A and B. These differences are more pronounced for cooler, longer period, models. The overall match of the Fmode decompositions with the OGLE RRab stars (gray points) is reasonable but shows some systematic differences. For example, the model values of ϕ 31 are larger  \n1O-mode (RRc)  \n7,300K  \n7,200K  \n7,100K  \n7,000K  \n6,900K  \nsets:  \nB D  \n0.65M  \n[Fe/H]=  \n, 45L  \n1.0  \n-  \n/circledot  \n/circledot  \n0.2783d  \n0.2915d  \n0.2914d  \n0.3056d  \n0.3055d  \n0.3206d  \n0.3206d  \n0.3367d  \nthan the observed values. However, the model physical parameters except T eff are fixed, while this is not the case for the OGLE RRab stars.  \nThe match between the 1O-mode models and the OGLE RRc stars (blue points) is worse. The amplitudes and the amplitude ratios are systematically too large. The Fourier phases are systematically too small. This may indicate different convective parameters are needed to reproduce the observed light curve shapes of F-mode and 1O-mode pulsators. Note the 1O-mode instability strip with convective set C is smaller than the F-mode instability strip at the luminosity considered. Models from set C, where the 1O-mode is linearly unstable and the integration is initialized with a 1O-mode velocity perturbation, all switch to an F-mode pulsation.  \nThe right panels in Figure 12 show that the light curve shapes are sensitive to the physical parameters. By varying the luminosity and metallicity in a narrow range, the model sequences match the OGLE values. However, no sequence considered follows the OGLE progression, because RR Lyrae in the Galactic bulge are characterized by mass, luminosity and metallicity distributions that cannot be reproduced with a single sequence.", '2.4.2. Classical F-mode Cepheids': "Cepheids display a feature known as the Hertzsprung progression (Hertzsprung 1926b). A secondary bump in their light and radial velocity curves appears near minimum light on the descending branch when P glyph[similarequal] 5 d. The bump moves towards earlier phases on the descending branch as the period increases and is coincident with maximum light when P glyph[similarequal] 10 d. The bump then moves onto the ascending branch for longer periods and disappears at P glyph[similarequal] 20 d. The bump is driven by a 2:1 resonance between the F-mode and a damped 2O-mode (e.g., Simon & Schmidt 1976; Buchler et al. 1990). This behavior is reflected in Cepheids' Fourier parameters, which follow more complex progressions with pulsation period than the RR Lyrae stars.  \nWe consider 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, and 10.0 M glyph[circledot] models with X =0.736, Z =0.008 and convective parameter sets A-D. Further details are in the test suite example 5M cepheid blue loop . Figure 13 shows the evolutionary tracks during the core He burning, blueloop phase for a selection of masses. Blue and red edges of the instability strip computed with RSP for convective sets B and D are shown. Edges for convective sets A and C largely overlap those for convective sets B and D, respectively. Nonlinear pulsation models are computed with a ∆ T eff =300K or 500K offset from the blue edge.  \nFigure 12. Comparison of the peak-to-peak amplitude and low-order Fourier decomposition parameters of I -band light curves of Galactic bulge RRab stars (gray dots) and RRc stars (blue dots) with synthetic light curves (symbols). In the left panels physical parameters are fixed, except for T eff and different convective parameter sets. In the right panels, the convective model is fixed and physical parameters are varied. Observational data are from Soszy'nski et al. (2014).  \n<!-- image -->  \nin the location of blue edge models, models from convective sets B and D have different T eff and consequently different pulsation periods. Figure 15 compares Fourier parameters of the Cepheid models with those derived from LMC Cepheid light curves (left panel, Soszy'nski et al. 2015; Ulaczyk et al. 2013) and Galactic F-mode Cepheid radial velocity curves (right panel, Storm et al. 2011). The I -band light curves are more sensitive to the convective parameters than the radial velocity curves for both ∆ T eff offsets. As with the RR Lyrae models, the Fourier parameters for convective sets A and B are sim-  \nr to each other, as are those for sets C and D. The radial velocity Fourier parameters follow tighter progressions than the I -band light curves.  \nThe Fourier phases in the left panels of Figure 15 are systematically larger than the observationally-inferred values for P glyph[lessorsimilar] 10 d, with the difference being larger for the cooler ∆ T eff =500K models. For P glyph[greaterorsimilar] 10 d the model Fourier phases are systematically smaller. Large discrepancies for the radial velocity curves are absent in the right panels of Figure 15, except for the amplitudes and R 21 ratio at the shortest periods. The projection  \nFigure 13. Evolutionary tracks in the HR diagram, shaded to show the core He mass fraction, Y c . Blue and red edges of the instability strip for convective sets B and D are shown, along with the locations where non-linear Cepheid models are computed (symbols).  \n<!-- image -->  \n(or p ) factor is the ratio of the pulsation velocity to the radial velocity deduced from spectral line-profile observations, dependent on at least rotation and gravity darkening (see Section 4), and plays a role in the amplitudes. We use p =1.3, close to the average value of determinations based on eclipsing binary Cepheid systems (Pilecki et al. 2018) and interferometric methods (e.g., Breitfelder et al. 2016).  \nThis brief survey is not exhaustive as only a few masses and two model sequences are explored. The M -L relation is also important for Cepheids as evolutionary tracks depend on overshooting, rotation, and metallicity.", '2.4.3. Type II Cepheids': "Type II Cepheids are more similar to RR Lyrae stars than to classical Cepheids due to their lower masses ( M glyph[similarequal] 0.5 M glyph[circledot] ), Pop II chemical composition, and evolutionary history. Their masses are similar to RR Lyrae but they cross the instability strip at larger luminosities and pulsate with longer periods ( P > 1 d). Type II Cepheids are F-mode pulsators except for a few doublemode stars pulsating simultaneously in the F and 1O modes (Smolec et al. 2018; Udalski et al. 2018) and two recently discovered 1O-mode pulsators (Soszy'nski et al. 2019). BL Her variables are a subclass of Type II Cepheids with P glyph[lessorsimilar] 4 d.  \nNonlinear radiative models of Type II Cepheids revealed a variety of complex dynamics including perioddoubled and deterministic chaos pulsations (e.g., Buchler & Kovacs 1987; Kovacs & Buchler 1988; Buchler  \nFigure 14. I -band light curves (upper panel) and radial velocity curves (lower panel) for the Cepheid models with a ∆ T eff =300K offset from the blue edge and convective sets B and D. Light curves are labelled with their pulsation periods and offset vertically by 0.5 mag or 35 km s -1 to facilitate comparison. Radial velocity curves follow the same order.  \n<!-- image -->  \n& Moskalik 1992). With convective pulsation models, modulated pulsations were also found (Smolec & Moskalik 2012). Buchler & Moskalik (1992) discovered perioddoubled pulsations in their survey of BL Her models and predicted that they should be observed in BL Her variables. The period-doubling is caused by a 3:2 resonance of the F and 1O modes and nonlinear phase synchro-  \nFigure 15. Amplitudes and Fourier parameters for I -band light curves (left panel) and radial velocity curves (right panel). Observations are marked with gray dots. Light curve data are from Soszy'nski et al. (2015) and Ulaczyk et al. (2013) and radial velocity curve data are from Storm et al. (2011). For the light curves the peak-to-peak amplitude is shown, while for radial velocity curves the Fourier amplitude is scaled with a p -factor of 1.3. Models are plotted with colored symbols.  \n<!-- image -->  \nnization (Moskalik & Buchler 1990); pulsations repeat only after two cycles of the F-mode.  \nSoszy'nski et al. (2011) report a period-doubled BL Her star, OGLE-BLG-T2CEP-279. We adopt M =0.6M glyph[circledot] , L =184L glyph[circledot] , T eff =6050K, X =0.76, Z =0.01, and convective parameter set A. These are nearly the same physical parameters Smolec et al. (2012) chose for a T2CEP279 model survey. The observed and model light curves are shown in Figure 16. The model amplitudes of the bump at minimum light of the period doubled cycle are larger and the shape of the light maximum is more pro-  \nounced relative to the observed features. The model period, P model =2.6976d, is longer than the observed period, P obs =2.3993d. Still, the light curves are qualitatively similar, devoid of fine-tuning, and demonstrate that RSP can be used to model specific stars.  \nFigure 17 shows the radius variation for two models with M =0.55M glyph[circledot] , L =136L glyph[circledot] , X =0.76, Z =0.0001 and convective set A but with a reduced eddy-viscosity, α m =0.05 (yielding unrealistic light curves). The two models differ only in T eff =5,932K (upper panel) and T eff =6410K (lower panel). Figure 18 shows the return  \nFigure 16. Comparison of the observed I -band light curve of OGLE-BLG-T2CEP-279 (blue dots) with the pulsation model (red curve). The period doubling effect is recognizable upon comparing consecutive maxima and minima. The model is vertically offset by 0.6.  \n<!-- image -->  \nFigure 17. Radius variation in two models differing only in their T eff that show periodic modulation (upper panel) and deterministic chaos (lower panel).  \n<!-- image -->  \nmap of maximum radii for these two models, plotting the maximum radii for each pulsation cycle R n max versus the preceding R n -1 max .  \nFor the cooler T eff =5,932K model, the upper panel of Figure 17 shows a cyclic modulation in the envelope of the period-doubled pulsation. The modulation period of glyph[similarequal] 57 d is longer than the pulsation period of glyph[similarequal] 2.3 d. The return map in the upper panel of Figure 18 is constructed from glyph[similarequal] 8,000 pulsation cycles and shows two loops, corresponding to alternating smaller and larger maximum radii. Since the modulation period is not commensurate with the pulsation period, the return maps develop a locus of points that form the closed lobes. Light curve modulation is common in RR Lyrae stars (Blazhko effect) and periodic pulsation modulation was  \nFigure 18. Maximum radius return maps for the models in Figure 17.  \n<!-- image -->  \nrecently discovered in BL Her variables (Smolec et al. 2018).  \nFor the hotter T eff =6,410K model, the radius variation appears irregular in the lower panel of Figure 17. The return map in the lower panel of Figure 18 reveals a strange attractor, an example of deterministic chaos in nonlinear models. Tracing time series from models is simple, but tracing chaotic dynamics in observations is difficult (e.g., Plachy et al. 2018). Chaotic dynamics is reported in a few type-II Cepheids with longer periods, in the RV Tau variable star range, and in semi-regular variable stars (e.g., Buchler et al. 1996; Kollath et al.  \n1998; Buchler et al. 2004; Plachy et al. 2018). While the T eff =6,410K model has a shorter period, in the BL Her range, such models may provide insight into chaotic dynamics in pulsating stars (see Plachy et al. 2013; Smolec & Moskalik 2014).", '2.4.4. Binary Evolution Pulsators': "Very-low-mass stars do not enter the classical instability strip within a Hubble time. However, mass loss from the more massive component in a close interacting binary can lead to a low-mass star that then evolves through the instability strip. The Binary Evolution Pulsator (BEP, OGLE-BLG-RRLYR-02792) is the prototype of this new class of pulsators (Pietrzynski et al. 2012; Smolec et al. 2013). The BEP's variability is similar to an F-mode RR Lyrae pulsator but with a dynamical mass of glyph[similarequal] 0.26 M glyph[circledot] .  \nFigure 19. I -band light (upper panel) and radial velocity (lower panel) curves of the BEP OGLE-BLG-RRLYR-02792 as a function of pulsation phase. Observations are shown as circles (Pietrzynski et al. 2012) with the radial velocity multiplied by a projection factor of p =1.2. Model curves are shown at three resolutions labeled by ( N/N outer ).  \n<!-- image -->  \nFollowing Smolec et al. (2013), we adopt M =0.26M glyph[circledot] , L =33L glyph[circledot] , T eff =6,910K, X =0.7, Z =0.01, and convective set A. Figure 19 shows I -band light curves and  \nradial velocity curves for the BEP models. Results for a coarse grid ( N/N outer =150/40, gold curves) qualitatively match these observations and have a period of 0.6373 d close to the observed 0.6275 d period. We also show results for a medium grid (300/120, orange curves) and a fine grid (300/120, red curves). The differences are most pronounced around maximum and minimum light. The shape of the light and radial velocity curves approach convergence only on grids with glyph[greaterorsimilar] 600 cells. The amplitude of the model curves are sensitive to convective parameters and can be fine-tuned for a better match with observations.", '2.4.5. Blue Large Amplitude Pulsators': 'The origin of BLAPs, introduced in Section 1, is unknown. They have been modeled as glyph[similarequal] 0.3 M glyph[circledot] shell Hburning stars that are progenitors of low mass WDs and glyph[similarequal] 1.0 M glyph[circledot] stars undergoing core He-burning (Pietrukowicz et al. 2017; Romero et al. 2018; Wu & Li 2018; Byrne & Jeffery 2018). Though mass loss in a close interacting binary must be invoked for both hypotheses, none of the BLAPs are known to be in a binary. Figure 1 shows that BLAPs are located near sdBVs in the HR diagram. The latter have non-canonical abundance profiles that are strongly affected by radiative levitation (see Section 6.2). Following the linear study of Romero et al. (2018), we adopt Z =0.05, to account for the increased envelope metallicity caused by radiative levitation.  \nWe explored nonlinear models with M =1.0M glyph[circledot] , T eff =30,000 K, L =430L glyph[circledot] , and three convection parameter sets. These envelope models used N =280, N outer =140, T anchor =2 × 10 5 K, and T inner =6 × 10 6 K. The LNA analyses show the F and 1O modes are unstable. Figure 20 compares the OGLE-BLAP-011 and 1O-mode ( P glyph[similarequal] 35 min) model I -band light curves. The qualitative agreement is reasonable. Figure 20 shows the model radial velocity curves have amplitudes of glyph[similarequal] 200 km s -1 , a pronounced temporal asymmetry, and light maxima that are narrower than light minima.  \nThese BLAP models lie far off the classical instability strip so that the initial model builder (see Section 2.2.1), which is optimized for classical pulsators, failed to relax the initial models to complete hydrostatic equilibrium. An option is to switch off the relaxation process and commence the time integration with a near-hydrostaticequilibrium initial model. The price to pay is the LNA growth rates are only indicative, not accurate.  \nFigure 21 shows the evolution of Γ for the ( α , α m )=(1.5,0.1) model in Figure 20. Before glyph[similarequal] 10 d, Γ fluctuates around a mean glyph[similarequal] 7.5 × 10 -3 . The fluctuations diminish with time but Γ remains above the LNA value of 6.25 × 10 -3 up to glyph[similarequal] 30 d. Between 30-50 d Γ diminishes and approaches  \nFigure 20. I -band light (upper panel) and radial velocity (lower panel) curves for OGLE-BLAP-011. Observations are shown as circles (Pietrukowicz et al. 2017). RSP models for three convective sets are shown under the glyph[similarequal] 1.0 M glyph[circledot] core He-burning hypothesis. Pulsation periods of the models range between 34.27 and 34.34 min and the observed period is 34.87 min.  \n<!-- image -->  \nzero, the sign that the nonlinear pulsation saturates at its terminal amplitude, yielding the results of Figure 20. In contrast to Γ, the period of the nonlinear pulsations remain close to the LNA period. In cases where the initial model cannot be relaxed, the initial Γ should not be expected to match the LNA analysis.', '2.4.6. High Amplitude Delta Scuti': "High-amplitude δ Scuti (HADS) pulsators are defined to have V -band light curve amplitudes greater than 0.1 mag. HADS lie close to the main-sequence (MS; see Figure 1), where growth rates are usually much smaller than those for RR Lyrae stars or classical Cepheids. This implies long time integrations are needed to drive nonlinear pulsations to saturation.  \nWeconsider a stellar model with M =2M glyph[circledot] , L =30L glyph[circledot] , T eff =6,900K, X =0.7, Z =0.01, and convective set A. This represents a star evolving towards the Hertzsprung gap. The LNA analysis of the initial model reveals that  \nFigure 21. Evolution of the RSP growth rate Γ starting from the unrelaxed initial model for convection set b in Figure 20. The horizontal line shows the LNA growth rate computed from the unrelaxed initial model.  \n<!-- image -->  \nthe F and 1O modes are linearly unstable, with growth rates of 1 × 10 -6 and 6 × 10 -5 , respectively.  \nFigure 22. Evolution of the relative cumulative energy error from model 10,000 ( glyph[similarequal] 15 cycles) to glyph[similarequal] 150 million (at glyph[similarequal] 250,000 cycles). The inset shows the per-step relative energy error over a cycle (600 steps).  \n<!-- image -->  \nSection 3 emphasizes the importance of numerical energy conservation. Figure 22 shows the evolution of the relative cumulative error in the energy for a 250,000 cycle integration ( glyph[similarequal] 150 million time steps, at 600 steps per cycle). The relative cumulative error grows from glyph[similarequal] 3 × 10 -11 after about 15 cycles to glyph[similarequal] 3 × 10 -7 after 250,000 cycles. The inset figure shows that the per-step  \nrelative error in the energy scatters around -2.5 × 10 -15 but is systematically different than zero.  \nFigure 23 shows that after 70,000 cycles the asymmetric V -band light curves have an amplitude of glyph[similarequal] 0.2 mag. The dominant pulsation period is 0.127 d, close to the LNA period of 0.1269 d for the 1O-mode. The amplitude, though, varies cyclically over about four pulsation cycles, reflecting the presence of the F-mode. Continuing the integration to 130,000 cycles leads to a more saturated light curve with an amplitude of glyph[similarequal] 0.3 mag and an unchanged dominant pulsation period. Extending the integration to 250,000 cycles does not lead to any significant changes.  \nFigure 23. V -band light curves of the δ Sct model at three stages during a very-long ( glyph[similarequal] 150 million time step) integration. For each stage, we plot the light curve for 5 pulsation cycles. The legend gives the cycle numbers when the snapshots were taken.  \n<!-- image -->  \nThe robustness in the light curves between 130,000 and 250,000 cycles might suggest that the long-term behavior is that of a double-mode HADS model. However, additional integrations with different initial perturbations, chosen to more adequately sample the neighboring amplitude phase space, are necessary to assess the possible mode selection. Figure 24 shows the evolution of these integrations in the amplitude-amplitude diagram using the analytical signal method (e.g., Koll'ath et al. 2002). The amplitude behavior of the model integration shown in Figure 23 is traced out by the red line. After initially rapid evolution, all the amplitude trajectories develop an arc along which evolution slows markedly. The model in Figure 23 and neighboring trajectories bend to the left toward smaller F-mode amplitudes. Their likely final state is that of a single-periodic 1O pulsation. In contrast, the two right-most trajectories bend to the right towards larger, more dominant,  \nFigure 24. Evolution of fractional radius amplitude, δR/R , for the 1O mode, A 1 and F mode, A 0 . Eight trajectories begin at locations marked with a cross. Circles mark 1,000 d intervals along trajectories. The red curve corresponds to the model in Figure 23 with the three open circles marking the amplitudes after 70,000, 130,000, and 250,000 cycles. The inset on the upper right zooms into the slowing amplitude evolution of the red curve late in the simulation.  \n<!-- image -->  \n0  \nF-mode amplitudes. Despite the limited sampling of phase space in Figure 24, we cautiously conclude that the most likely long-term outcome of this δ Sct model is a single-periodic 1O-mode or F-mode pulsator depending on the initial conditions (see Section 2.2.4).", '3. ENERGY CONSERVATION': "For the following discussion we define the total energy E of a model to be the sum of the internal, potential, and kinetic energies, ignoring rotational and turbulent energy which are currently not included in the energy accounting. To support improved numerical energy conservation 1 , MESAstar provides an option to use what we call the dedt -form of the energy equation:  \nD e D t = glyph[epsilon1] -∂ ∂m ( L + p A u ) -D D t ( 1 2 u 2 -Gm r ) . (8)  \nThis form was introduced in Paper IV and provides an alternative to the dLdm -form of the energy equation (Equation (11) in Paper I). When the time derivative terms are combined, the result is more easily recognizable as an equation for the time evolution of local specific total energy (left hand side) due to local source terms 2  \n1 Note that we are discussing numerical issues in the code rather than questions of the physical completeness and validity of the equations. We will often use the term numerical energy conservation to make this distinction explicit.  \n( glyph[epsilon1] ) and local fluxes between cells (the ∂/∂m term):  \nD D t ( e + 1 2 u 2 -Gm r ) = glyph[epsilon1] -∂ ∂m ( L + p A u ) . (9)  \nThe error in numerical energy conservation E error is the extent to which the time- and mass-integrated Equation (9) is not satisfied when solved in a discretized, finite-mass form. This section discusses recent efforts to improve numerical energy conservation in MESAstar .  \nRecall (from Paper I, Section 6.3) the generalized Newton-Raphson scheme used by MESAstar to solve the stellar equations  \n0 = glyph[vector] F ( glyph[vector]y ) = glyph[vector] F ( glyph[vector]y i + δglyph[vector]y i ) = glyph[vector] F ( glyph[vector]y i ) + [ d glyph[vector] F dglyph[vector]y ] i δglyph[vector]y i + O ( δglyph[vector]y 2 i ) , (10)  \nwhere glyph[vector]y i is the trial solution for the i -th iteration, glyph[vector] F ( glyph[vector]y i ) is the residual, δglyph[vector]y i is the correction, and [ d glyph[vector] F/dglyph[vector]y ] i is the Jacobian matrix. The residual is the left over difference between the left- and the right-hand sides of the equation we are trying to solve, while the correction is the change in the primary variable that is calculated by Newton's rule. The solver generates a series of trial solutions until it produces one that is acceptable according to given convergence criteria. In MESAstar the trial solution is not accepted until the magnitudes of all corrections and residuals become smaller than specified tolerances 3 . If no acceptable trial solution has been found in the allowed maximum number of iterations, the solver rejects the attempt and forces a retry with a smaller time step.  \nIf the numerical accuracy of the partial derivatives forming the Jacobian matrix is not excellent, the reduction in magnitude of the residuals can stall after a few iterations. For this reason, MESAstar has provided a means to use a tight tolerance on residuals for an initial sequence of iterations and then switch to a much relaxed tolerance if no acceptable solution has been found. The benefit of this is that residuals will be driven down when possible, but if the residuals stall at a level above the initial tolerance, the system will still be able to take a step as long as the corrections can be adequately reduced.  \nprocesses such as accretion (see Section 3.3). Importantly, glyph[epsilon1] does not include glyph[epsilon1] grav , the specific rate of change of gravothermal energy, as that source term is not present when using a total form of the energy equation (see Paper IV, Section 8).  \n3 The GARching STellar Evolution Code (GARSTEC) is the only other stellar evolution code we are aware of that considers residuals as well as corrections in deciding when to accept a trial solution (Weiss & Schlattl 2008). Several other codes consider corrections but not, as far as we can tell, residuals (Roxburgh 2008; Christensen-Dalsgaard 2008; Demarque et al. 2008; Scuflaire et al. 2008; Faulkner 1968).  \nThe cost of relaxing the tolerance for residuals of the total energy equation is the creation of numerical energy conservation errors. To obtain good numerical energy conservation we must be able to drive down residuals to low levels, and to do that we must have numerically accurate partial derivatives. This has motivated a major effort to improve partials, and we can now require the solver to keep iterating until it reduces the residuals to a low level that gives good numerical energy conservation. The most significant changes to improve the numerical accuracy of partials were in the eos module and are discussed in Appendix A.1.", '3.1. Gold Tolerances': "To improve energy conservation, a new standard 'gold tolerances' is now the default in MESAstar . This uses tight tolerances that apply even after an arbitrarily large number of iterations. As a result, steps with poor residuals will be rejected, thereby ensuring that if a run succeeds with gold tolerances enabled while using the total energy equation given above, it will have good energy conservation. To show example improvements from this new strategy, in Table 5 we report the results of calculations of a 1 M glyph[circledot] model during the main sequence and the He flash with gold tolerances and compare to the old approach. The 1M glyph[circledot] models are evolved from the ZAMS until the core H mass fraction reaches a value of 10 -6 . The He flash models start at off-center He ignition and terminate when the core He mass fraction drops to 10 -3 . The cumulative energy error is the sum of the energy conservation errors at each of the steps. The relative cumulative energy error is the cumulative energy error divided by the final total energy. This is now much less than 1% during the He flash, a notoriously difficult evolutionary phase from the numerical perspective. The evolution of the errors in energy conservation, as well as the number of iterations required by the solver and the adopted time step, are shown in Figure 25 for the He flash runs. This figure clearly show the superiority of the model adopting gold tolerances and the dedt -form of the equation in terms of energy conservation.  \nNot all cases in the MESAstar test suite are currently able to use gold tolerances. This is primarily because of remaining problems with the numerical accuracy of certain partials, especially in the PC EOS (Potekhin & Chabrier 2010) that is used by WD models (see Appendix A.1) and on the boundary between the OPAL EOS (Rogers & Nayfonov 2002) and the SCVH EOS (Saumon et al. 1995). There are also problems with the numerical accuracy of some partials associated with nuclear reactions at the high temperatures encountered during late stages of evolution, such as Si burning. To  \nFigure 25. Evolution of a 1 M glyph[circledot] model during the He flash and core He burning using gold tolerances and the dedt -form of the energy equation (left panel), and without gold tolerances and using the dLdm -form of the energy equation (right panel). The first panel shows a Kippenhahn plot, where green hatched regions are convective and blue shading shows region of nuclear energy generation. The second panel shows the absolute error in energy conservation divided by the total stellar energy, both per step and cumulative. The third panel shows the time step, while in the fourth panel we report the number of iterations required by the Newton-Raphson solver. Gold tolerances keep the cumulative errors small during the He flash at the price of a larger number of iterations during the first ≈ 200 time steps.  \n<!-- image -->  \naddress these situations, there are controls to allow gold tolerances to be turned off automatically for steps that either require use of the PC EOS or that have extremely high temperatures. To provide feedback, the value of the relative cumulative energy error is monitored and a warning message is written if it exceeds a specified value (2% is the default setting).", '3.2. Definition of glyph[epsilon1] nuc Source Term': 'Table 5. Relative cumulative energy error for 1 M glyph[circledot] runs.  \nPrevious MESA papers have not given a precise definition of glyph[epsilon1] nuc . Motivated by numerical convenience, MESA  \nformerly exploited the fact that the glyph[epsilon1] source term contained the sum of glyph[epsilon1] grav and glyph[epsilon1] nuc and included the response of the internal energy to composition changes due to nuclear reactions in glyph[epsilon1] nuc instead of in glyph[epsilon1] grav . However, since the dedt -form does not include glyph[epsilon1] grav , this is no longer an appropriate choice.  \nIn the current approach, glyph[epsilon1] nuc is evaluated in the net module as a sum over reactions. Schematically,  \nglyph[epsilon1] nuc = N A ∑ reactions ( Q i -Q ν,i ) R i , (11)  \nwhere R i is the molar rate of reaction i , Q i is the change in rest mass energy between the products and reactants, and Q ν,i is the per-reaction average energy of the neutrino (if present). We note the equivalence  \n∑ reactions Q i R i ≡ -∑ isotopes M i c 2 ˙ Y i , (12)  \nwhere M i is the rest mass of isotope i and ˙ Y i is the rate of change of the molar fraction. The approx family of MESA nuclear networks exploit this equivalence and do not strictly follow Equation (11). The right hand side of Equation (12) is a common nuclear physics definition of glyph[epsilon1] nuc (e.g., Equation 11 in Hix & Meyer 2006). The MESA definition of glyph[epsilon1] nuc differs in subtracting off the nuclear neutrino losses, thus enforcing the assumption that they free-stream out of the star.', '3.3. Mass Changes': 'The methods described above perform well for numerical energy conservation when the stellar mass is constant, but extensions are necessary for cases where the mass changes. For energy accounting, we must specify the amount of energy we expect the new mass to introduce and the amount we expect departing mass to remove. To that end we assume that mass being added has the same specific energy as the surface of the model at the start of the step. For mass being removed, we assume that it leaves with a specific energy between what it had at the start of the step and the value at the surface at the start of the step. The exact amount depends on the amount of energy that leaks out of the material as it approaches the surface during the time step. For low rates of mass loss there will be adequate time for the material to adjust so that it leaves with the initial surface value, but for high rates there may not be enough time for adjustment, so that it leaves with a specific energy closer to its starting value. The details of this are presented below.  \nIn addition to providing accurate accounting for the total energy of the model, it is important to ensure  \nthat the energy changes from mass loss or gain are distributed properly within the model. As a guide for this we use the analytic calculations of Townsley & Bildsten (2004). Our new procedure improves upon these by also calculating the distribution of energy in systems with long thermal times, allowing MESA to handle the limit of rapid accretion. We confirm the numerical energy conservation of this method using the > 20 cases in the MESAstar test suite that have mass changes and fully support gold tolerances and the dedt -form of the energy equation. Using the new scheme, each of these completes the test run with a cumulative error in total energy < 2 %. In addition, test cases that depend on the internal distribution of accretion heating continue to yield the expected results.', '3.3.1. Methodology': 'Because MESA works on a Lagrangian mesh, it handles accretion and mass loss in a two-stage process (Paper III, Section 7). In stage I the masses of certain cells are increased or decreased as needed to give the desired end-of-step total mass, but no attempt is made to ensure energy conservation at this stage. In stage II the model thus produced is evolved in time by an amount δt . This separation of stages means that the time step is only taken for a model of fixed mass. However, because the energy of the model changes in stage I, a correction must be added to stage II to make the overall step consistent with energy conservation. Thus, we introduce a new source term glyph[epsilon1] ˙ M that accounts for the heating associated with mass changes.  \nThe change in the mass of cell k during stage I results from the difference between the outward 4 mass flux F m through each cell face. This flux obeys  \ndm mid ,k -dm start ,k = δt ( F m,k +1 -F m,k ) (13)  \nand  \nM mid -M start = -δtF m, 1 = δt ˙ M. (14)  \nDuring stage I the temperature, density, and velocity of each cell are held fixed but the composition is updated to track the flow of material between cells. The subscript start is used for quantities at the start of stage I. The subscript mid is used for quantities at end of stage I, which is the start of stage II. No subscript is used for quantities evaluated at the end of the time step (after stage II). There is no mass change during stage II, so dm k = dm mid ,k . In the following we write dm k rather than dm mid ,k .  \nThe change in energy for cell k during stage I is  \nE mid ,k -E start ,k = dm k E mid ,k -dm start ,k E start ,k (15) = dm k ( E mid ,k -E start ,k ) +( dm k -dm start ,k ) E start ,k ,  \nwhere E k is the total specific energy of cell k , given by the sum of specific potential, kinetic, and internal energies. Neglecting changes in specific energy owing to changes in composition, the difference in E k across stage I is  \nE mid ,k -E start ,k = ( Gm k, C r k, C ) start -( Gm k, C r k, C ) mid . (16)  \nwhere m k, C and r k, C are the mass coordinate and radius at the center of mass of the cell (Paper IV).  \nWe now introduce glyph[epsilon1] k, I , the effective source term in cell k during stage I. This is defined by writing the change in energy in flux-conservative form as  \nE mid ,k -E start ,k = δt ( dm k glyph[epsilon1] k, I + F e,k +1 , I -F e,k, I ) , (17)  \nwhere F e,k, I is the outward flux of energy across face k owing to work and material passing through face k . Inserting Equation (15) and rearranging we obtain  \nglyph[epsilon1] k, I ≡ E mid ,k -E start ,k δt + E start ,k δt ( 1 -dm start ,k dm k ) -F e,k +1 , I -F e,k, I dm k . (18)  \nThe energy flux is  \nF e,k, I = p face ,k A k ˆ v k, I + E face ,k F m,k , (19)  \nwhere  \nˆ v k, I ≡ r mid ,k -r start ,k δt , (20)  \nand the face values E face are interpolated 5 from the cell values E start . The radial coordinate after state I is calculated using the updated cell masses, holding cell densities fixed. The change in total energy of the model during stage I is  \n∑ k ( E mid ,k -E start ,k ) = -δt ( p face , 1 A 1 ˆ v 1 , I + E face , 1 F m, 1 )  \n+ δt ∑ k dm k glyph[epsilon1] k, I . (21)  \nThe term in parentheses on the right-hand side accounts for the energy of new material entering or leaving the model and the work done in the process. The additional glyph[epsilon1] k, I term implies the need for a corrective source term glyph[epsilon1] ˙ M which must be added during stage II so that energy is properly accounted for.  \nTo determine this new source we consider the ratio of the thermal time scale  \nto the mass-change time scale  \nτ th ,k ≡ ∣ ∣ ∣ ∣ x m c p T ( ∇ rad -∇ ad ) L ∣ ∣ ∣ ∣ start ,k , (22)  \nτ m ,k ≡ ∣ ∣ ∣ ∣ x m F m ∣ ∣ ∣ ∣ start ,k , (23)  \nwhere x m,k is the mass above face k . When this ratio is small the thermodynamic state of material is a function primarily of depth (Sugimoto 1970; Sugimoto & Nomoto 1975). This means that as material moves from cell to cell it is a good approximation to the time evolution to suppose that it adopts the state of whichever cell it is in (Paper III). In the opposing limit ( τ th glyph[greatermuch] τ m ) the entropy of material adjusts minimally as it moves from cell to cell.  \nWe account for the effects of the thermal and masschange time scales by tracking the heating of material as it moves from cell to cell in stage I and estimating what fraction of that heat is released as part of L versus carried by the material. Consider the path taken by an infinitesimal fluid element as it moves from face k to face k +1 due to accretion. Over the course of this adjustment the material changes state and releases some heat. We take this heat to be given by  \ndq k = -δtglyph[epsilon1] k, I dm k dm pass , k , (24)  \nwhere  \ndm pass , k = dm start ,k + δt (max(0 , F m,k +1 ) -min(0 , F m,k )) (25)  \nis the total mass of material that at any point during stage I was inside cell k . This evenly distributes the heating which occurs in a cell over all material that starts in, ends in, or passes through the cell. When the thermal time is long, however, the fluid does not have a chance to release all of this heat before it has finished crossing the cell. We estimate the fraction of this heat it releases as the leak fraction  \nf leak ,k ≡ min ( 1 , τ m ,k τ th ,k ) . (26)  \nThis parameterises the extent to which material flowing through a cell follows the implied energy gradient ( f leak ,k ≈ 1) versus evolving adiabatically ( f leak ,k glyph[lessmuch] 1).  \nWe define δE k,j to be the amount of energy that the material which ends up in cell j had not leaked by the time it reached face k . This is zero for k = 1 and for k > j , and for all other faces is given by  \nδE k +1 ,j = (1 -f leak ,k ) ( δE k,j + dq k M k,j ) , (27)  \nwhere M k,j is the amount of material which ends in cell j which passes through cell k during stage I. The heat which is actually released in cell k is then given by  \nglyph[negationslash]  \nglyph[epsilon1] ˙ M,k = δE k,k + ∑ j = k f leak ,k ( δE k,j + dq k M k,j ) δtdm k . (28)  \nThis is our new corrective source term. The same procedure may be used in the case of mass loss, but with δE N +1 ,j = 0 instead of δE 1 ,j and -1 in the subscript on the left-hand side of Equation (27) rather than +1. In the limit of long thermal times most heat is retained and the resulting evolution is adiabatic. In the limit of short thermal times we recover the results of Paper III.  \nAlong the lines of Equation (17), the energy change of cell k during stage II may now be written in fluxconservative form as  \nE k -E mid ,k = δt [ dm k ( glyph[epsilon1] ˙ M,k + glyph[epsilon1] ) + F e,k +1 , II -F e,k, II ] , (29)  \nwhere  \nF e,k, II = L k + p face ,k A k ˆ v k, II (30)  \nˆ v k, II ≡ r k -r mid ,k δt . (31)  \nand  \nMESA solves Equation (29) when using the dedt -form with mass changes. Because the time evolution is implicit all sources are evaluated at the end of stage II. Nuclear burning is evaluated as if material spends the whole step in the cell in which it ends stage I. While this is usually a good approximation it may break down when δt ˙ M is large relative to the masses of cells in burning regions.  \nWhen ˙ M ≥ 0 the above procedure for redistributing energy is conservative, so that  \n∑ k glyph[epsilon1] k, I δtdm k + ∑ k glyph[epsilon1] ˙ M,k δtdm k = 0 . (32)  \nThe first sum represents the effect of stage I; the second sum represents the effect of stage II. When ˙ M < 0 the equality is instead  \n∑ k glyph[epsilon1] k, I δtdm k + ∑ k glyph[epsilon1] ˙ M,k δtdm k = -δE 1 , 0 , (33)  \nwith the additional term δE 1 , 0 being the energy carried out of the star. This term is explicitly accounted for in the MESA energy budget, so in both cases energy is properly accounted for.  \n3.3.2. ResultsFigure 26. Comparison of glyph[epsilon1] ˙ M and the analytic accreting heating expression from Townsley & Bildsten (2004) (TB 04) for a 0 . 33 M glyph[circledot] He WD accreting He and N at a rate of 10 -5 M glyph[circledot] yr -1 . The first panel shows both as functions of depth x m . The second panel shows the same but integrated inward with respect to mass. The third panel shows the ratio of the TB 04 expression to glyph[epsilon1] ˙ M . The final panel shows the thermal time scale τ th and the mass-change time scale τ m . All material to the left of the red vertical line is material which is new in this time step. To the right of the vertical black line MESA transitions to a Lagrangian mesh which continues over the remainder of the core.  \n<!-- image -->  \nTo examine the behavior of glyph[epsilon1] ˙ M we modeled a 0 . 33 M glyph[circledot] WDaccreting He and N at a rate of 10 -5 M glyph[circledot] yr -1 . The first panel of Figure 26 shows a profile of glyph[epsilon1] ˙ M along with  \nthe analytic accretion heating calculations of Townsley & Bildsten (2004). For the most part the two agree closely. The second panel shows the mass integral of the same inward from the surface while the third shows their ratio. Around x m ≈ 10 26 g, τ th becomes long relative to τ m and the two prescriptions differ because that of Townsley & Bildsten (2004) is only applicable where τ th glyph[lessmuch] τ m . The glyph[epsilon1] ˙ M term handles both limits.  \nFigure 27. Comparison of energy conservation for a 0 . 3 M glyph[circledot] He WD model accreting an H/He mixture at a rate of 10 -10 M glyph[circledot] yr -1 . The upper panel shows the relative cumulative error and the lower shows the error in each step.  \n<!-- image -->  \nTo demonstrate improved energy conservation during mass changes with the dedt -form we model accretion onto a 0 . 3 M glyph[circledot] He WD with an initial log( T c / K) = 6 . 7. The accretion rate is fixed at 10 -10 M glyph[circledot] yr -1 . Nuclear reactions are disabled throughout the run. This is repeated with the dLdm -form. The relative cumulative error in energy conservation is shown in the upper panel of Figure 27, and the relative error in each step is shown in the lower panel. Near the beginning of the run there is a period where the dLdm -form performs better, however at those early times both forms do a good job, conserving energy to one part in ∼ 10 5 . At later times the dedt -form produces less error, staying below one part in 10 4 while the dLdm -form yields cumulative errors greater than one part in 10 3 .  \nWe prefer using the dedt -form because of its improved energy conservation and handling of long thermal times.  \nFigure 28. Comparison of the surface luminosity for a 0 . 3 M glyph[circledot] He WD model accreting a H/He mixture at 10 -10 M glyph[circledot] yr -1 (same as Figure 27). The upper panel shows the surface luminosity and the lower shows the relative difference between the dLdm -form and the dedt -form, all as functions of time since the start of accretion.  \n<!-- image -->  \nWe now explore the consequences for stellar evolution of these different prescriptions. The upper panel of Figure 28 shows the luminosity for the same case as Figure 27, again for both the dedt -form and the dLdm -form. The difference is shown in the lower panel. The largest differences are at early times as the models adjust to the accretion. After that, both yield results similar to a few percent, and the relative difference only improves as the luminosity increases. Figure 29 shows T c for the same two runs as a function of time. The differences are small. This is because the core lies beneath the region where cell masses are adjusted significantly, so the precise handling of mass changes only matters for the core insofar as the core temperature is sensitive to the luminosity.  \nMass changes are ubiquitous in binary stellar evolution. To demonstrate the effect of the dedt -form we model the evolution of a binary system with an 8 M glyph[circledot] primary and a 6 . 5 M glyph[circledot] secondary with an initial orbital period of 3 days. This is the same example as in Figure 4 of Paper III. The relative cumulative error in energy conservation is shown in the upper panel of Figure 30. The model is also run with the dLdm -form. The errors are shown independently for the primary and the secondary. Figure 31 compares the mass transfer rate for  \nFigure 29. Comparison of the central temperature for a 0 . 3 M glyph[circledot] He WD model accreting a H/He mixture at 10 -10 M glyph[circledot] yr -1 (same as Figure 27). The upper panel shows the central temperature and the lower shows the relative difference between the dLdm -form and the dedt -form, all as functions of time since the start of accretion.  \n<!-- image -->  \nFigure 30. Comparison of energy conservation for the case of Figure 4 of Paper III: a binary system with mass transfer from the 8 M glyph[circledot] primary to the 6 . 5 M glyph[circledot] secondary and an initial orbital period of 3 days. The relative cumulative energy error computed using both the dedt - and dLdm -forms is shown as a function of time after the start of accretion.  \n<!-- image -->  \nthis system computed using each energy equation. The differences are typically of order 1 %.  \nFigure 31. Comparison of mass transfer rates for the case of Figure 4 of Paper III: a binary system with mass transfer from the 8 M glyph[circledot] primary to the 6 . 5 M glyph[circledot] secondary and an initial orbital period of 3 days. The upper panel shows the mass transfer rate computed using both the dedt - and dLdm -forms as a function of time after the start of accretion. Both are similar to the corresponding curves in Figure 4 of Paper III. The lower panel shows the relative difference between ˙ M computed with the dedt - and dLdm -forms. The spike at the end is due to mass transfer terminating at slightly different times.  \n<!-- image -->', '4. ROTATION': 'For a rigidly-rotating star in hydrostatic equilibrium, surfaces of constant pressure (isobars) coincide with the equipotential surfaces defined by the Roche potential Ψ,  \nΨ( r, θ ) = Φ( r, θ ) -r 2 Ω 2 sin 2 θ 2 , (34)  \nwhere Φ is the standard Newtonian potential, θ is the polar angle, and Ω is the angular frequency of rotation. In one-dimensional stellar evolution calculations the effects of rotation on the stellar structure are usually captured by a simple modification of the stellar equations. These retain their regular form, but include two correction factors,  \nf P = 4 πr 4 Ψ Gm Ψ S Ψ 1 〈 g -1 〉 , f T = ( 4 πr 2 Ψ S Ψ ) 2 1 〈 g 〉〈 g -1 〉 , (35)  \nwhere r Ψ is the volume-equivalent radius of an isobar, and m Ψ and S Ψ are the mass inside that isobar and its  \nsurface area respectively (Kippenhahn & Thomas 1970; Endal & Sofia 1976). The effective gravity is g = |∇ Ψ | , while 〈 g 〉 and 〈 g -1 〉 are surface averages over the equipotential. 6  \nThis approach is still applicable to the case of a differentially-rotating star under the assumption of shellular rotation, in which shells are rigidly-rotating, isobaric surfaces with rotation frequency Ω Ψ (Endal & Sofia 1976; Zahn 1992). The averages in Equation (35) are performed over each isobar (Meynet & Maeder 1997).  \nUntil now, MESA used the method of Endal & Sofia (1976), which considers deviations of the Roche potential from spherical symmetry (Kopal 1959), to compute the f P and f T factors. One issue is that to ensure numerical stability, this approach requires a floor on the correction factors ( f P = 0 . 75 and f T = 0 . 95, Paper II), corresponding to a maximum rotation rate of 60% of critical rotation (the angular frequency at which the centrifugal force would match gravity at the stellar equator).  \nAs stars are centrally condensed and rotational frequencies close to critical are typically reached only in the outermost layers, Ψ is well approximated by the potential of a point mass in rapidly rotating layers (e.g., Maeder 2009). This justifies using the Newtonian potential as Φ = -Gm Ψ /r for the calculation of f P and f T , such that they are only functions of the fraction of critical rotation ω ≡ Ω Ψ / Ω crit , Ψ = Ω Ψ / √ Gm Ψ /r 3 e . Here r e is the equatorial radius of the isobar and Ω crit , Ψ is the rotational frequency at which the centrifugal force is equal to gravity at the equator of the isobar.  \nWe describe a new implementation of centrifugal effects in MESA , which makes use of analytical fits to the Roche potential of a point mass, improving the calculation of rotating stars to ω ≈ 0 . 9.', '4.1. The Roche Potential of a Rigidly-Rotating, Single Star': "For a point mass m Ψ , the dimensionless Roche potential Ψ ' = Ψ / ( Gm Ψ Ω) 2 / 3 can be written in terms of the dimensionless radius r ' = r/ ( Gm Ψ / Ω 2 ) 1 / 3 as  \nΨ ' = -1 r ' -r ' 2 sin 2 θ 2 . (36)  \nNote that  \nr ' e = r e ( Gm Ψ / Ω 2 ) 1 / 3 = ω 2 / 3 , (37)  \n6 In Paper II we defined the volume equivalent radius as r P , while here we adopt the symbol r Ψ as used by Endal & Sofia (1976). This change is to prevent confusion between r Ψ and the polar radius of an isobar, which we denote as r p .  \nFigure 32. Equipotential lines of the dimensionless Roche potential Ψ ' given by Equation (36). Solid lines show the equipotentials for which ω is equal to 0 . 25 , 0 . 5 , 0 . 75 and 1. Dashed lines show ellipses with the same polar and equatorial radii as the actual equipotentials.  \n<!-- image -->  \nsuch that evaluating Equation (36) for θ = 0 ( r ' = r ' p ) and θ = π/ 2 ( r ' = r ' e ) provides the ratio of the polar radius r p to r e as a function of ω ,  \nr e r p = 1 + 1 2 ω 2 . (38)  \nFigure 32 shows how the equipotential surfaces change with increasing ω . For ω glyph[lessorsimilar] 0 . 5 the equipotentials are approximately given by oblate spheroids, while for ω close to unity a cusp develops at the equator. For this critically rotating surface, the polar radius is 2 / 3 of the equatorial one.  \nBy determining the asymptotic behavior of the Roche potential in the limit ω → 0, and, when possible, also in the limit ω → 1, we have constructed analytic fits to properties of interest. These are described in Appendix B, and include fits for the equatorial radius r e ( r Ψ , ω ), the centrifugal corrections f P ( ω ) and f T ( ω ), and the volumes and surface areas of Roche equipotentials, V Ψ ( r e , ω ) and S Ψ ( r e , ω ). Previous versions of MESA approximated the specific moment of inertia of isobaric surfaces as that of a thin spherical shell with radius r Ψ , i rot = 2 r 2 Ψ / 3, but now default to a fit of the form i rot ( r e , ω ). Figure 33 shows the resulting fits for i rot , f P and f T . The new implementation results in values of the specific moment of inertia that are larger by a factor of two as ω → 1.  \n4.2. Implementation in Stellar Evolution Instruments  \nFigure 33. The upper panel shows the specific moment of inertia for a shell of material at different values of ω , normalized by 2 / 3 r 2 Ψ . The model implemented previously in MESA is shown with a dot-dashed green line, and the new ω -dependent model for i rot is shown with a solid orange line. As ω → 1, the moment of inertia becomes that of a ring with radius r e . The lower panel shows the f P and f T factors as a function of ω = Ω / Ω crit .  \n<!-- image -->  \nTo include these fits into a stellar evolution calculation that uses the shellular approximation, a value of ω must be determined for a given specific angular momentum j rot , m Ψ , and r Ψ . From Ω Ψ = ω Ω crit , Ψ = ω √ Gm Ψ /r 3 e and j rot = i rot Ω Ψ , we find  \nj rot √ Gm Ψ r Ψ = ω i rot r 2 e √ r e r Ψ . (39)  \nFor a given j rot , m Ψ and r Ψ , the left-hand side can be directly evaluated, while the right-hand side is a monotonic function of ω for 0 ≤ ω ≤ 1. We compute the solution to this equation for each cell in the stellar model using a bisection method. Given ω , we then use the computed fits to determine the values required by the structure equations: r e , r p , i rot , Ω Ψ , f P and f T . As in the previous implementation of rotation in MESA , we evaluate these quantities explicitly at the beginning and at the end of a step. The analytical nature of the fits allows the possibility of a fully-coupled and implicit implementation in the future.  \nFigure 34 shows a 3 M glyph[circledot] solar metallicity model with different initial rotational velocities, evolved using the current and previous implementations of rotation in MESA . Both methods agree for rotation rates ω < 0 . 5.  \nFigure 34. Stellar evolution calculations with initial mass 3 M glyph[circledot] and solar metallicity, from ZAMS up to terminal age main sequence (TAMS). Different colors indicate different initial rotation rates, defined in terms of the ratio of the rotational frequency to the critical value at the surface at ZAMS. Solid and dotted lines indicate calculations done with the current and the previous implementations of rotation in MESA , respectively.  \n<!-- image -->  \nDifferences at higher rotation rates are due to the aforementioned floor on f P and f T . In our current approach we also define a floor on f P and f T in terms of a maximum value ω max , beyond which the effects are truncated to f P ( ω max ) and f T ( ω max ). We find that calculations using the new strategy are numerically stable near critical rotation, with the simulations shown in Figure 34 being performed with ω max = 0 . 9. This is in comparison to the previous method, that for rapidly-rotating models set a floor on f P and f T corresponding to their values at ω ≈ 0 . 6. Therefore, MESA can now consistently calculate shellular rotation models closer to critical rotation.", '4.3. Gravity Darkening Corrections': 'Rotating stars are subject to gravity darkening (von Zeipel 1924a,b). The variation of flux over the surface and the distorted stellar shape imply that the observed properties of the star vary with the angle between the rotation axis of the star and the line of sight (LOS). Here we describe our approach to calculating geometric factors that allow the intrinsic surface quantities L and T eff to be corrected for projection effects along a given LOS. By intrinsic we mean the total L emitted by the star and the T eff associated with this L given the total surface area of the star and the Stefan-Boltzmann law. This is done in two steps. First, we solve the gravity darkening problem for an arbitrary surface element of  \na rotating star, and second, calculate the projection of the gravity-darkened surface along the LOS.', '4.3.1. The Gravity Darkening Model': 'Figure 35. Solution of the ELR model for 0 ≤ ω ≤ 1; each curve plots the variation of ˜ T eff as a function of cos( θ ) for a different value of ω . Recall that cos( θ ) = 1 corresponds to the pole and cos( θ ) = 0 to the equator.  \n<!-- image -->  \nWe use the gravity darkening model of Espinosa Lara & Rieutord (2011, hereafter ELR), where it is assumed that the radiative flux is directed antiparallel to the effective surface gravity. At a point on the stellar surface with polar angle θ , we find the value of the scaled photosphere radius ˜ r = R/R e by solving  \n1 ˜ r + ω 2 2 ˜ r 2 sin 2 θ = 1 + ω 2 2 . (40)  \nWe then solve  \ncos ϑ +lntan( ϑ/ 2) = 1 3 ω 2 ˜ r 3 cos 3 θ +cos θ +lntan( θ/ 2) , (41)  \nfor the modified angular variable ϑ . Using ELR Equation (31), we use this value to obtain the local T eff .  \nFigure 35 shows the variation of ˜ T eff = T eff [ L/ (4 πσR 2 e )] -over a range of θ for a series of curves with different values of ω . When ω = 0, ˜ T eff = 1 for all θ . When ω = 1, ˜ T eff varies by nearly a factor of 2 between the pole and the equator.  \n1  \n/  \n4', '4.3.2. Projection Effects and Correction Factors': 'We are interested in the projected -the directional average over the surface along the LOST eff and L . The two parameters governing the problem are ω and the inclination angle, i , of the LOS with respect to the rotation axis of the star: i = 90 · when the LOS is in the plane of the equator. We denote the LOS unit vector ˆ l ( i )  \nand the projected surface area Σ proj . Figure 36 shows a grid of Roche equipotential surfaces for different values of ω and i . The color describes the variation of ˜ T eff over the surface.  \nTo calculate the luminosity projected along the LOS, L proj , requires the surface integral  \nL proj = 4 × ∫∫ d glyph[vector] Σ · ˆ l> 0 Fd glyph[vector] Σ · ˆ l, (42)  \nwhere only the flux projected toward the observer, i.e., d glyph[vector] Σ · ˆ l > 0, is kept. The emergent specific intensity from each surface element is assumed to be isotropic. Once L proj and Σ proj are known, the projected T eff can be obtained from the Stefan-Boltzmann law  \nT eff , proj = ( L proj σ Σ proj ) 1 / 4 . (43)  \nAs noted by Georgy et al. (2014), the ratio of L proj /L , and likewise the ratio of T eff , proj /T eff , are geometric factors that depend only on ω and i . We define gravity darkening coefficients  \nC L ( ω, i ) ≡ L proj ( ω, i ) /L (44)  \nand  \nC T ( ω, i ) ≡ T eff , proj ( ω, i ) /T eff . (45)  \nC T and C L are tabulated in MESA for the valid domain of ω and i ; values are readily obtained via bicubic interpolation. The projected L and T eff as viewed from the pole and the equator can be output in the MESA history file; other inclination angles can be accessed via run star extras .  \nFigure 37 shows how C T and C L vary over the ( ω , i )-plane. These can be compared with the upper panels of Figure 2 in Georgy et al. (2014). The variation of C L is greater than that of C T due to the 1 4 -power relationship between L and T eff . At ω = 1, where the geometric factors are the largest, C L varies from -20% to +50% while C T varies by only about ± 2 . 5%.  \nFigure 37 shows a slight, but noticeable, decrease in C T for ω glyph[similarequal] 1 whereas the comparable figures from Georgy et al. (2014) do not. When we calculate the coefficients using oblate spheroids instead of Roche equipotential surfaces, we see no such decrease in C T for nearcritical rotation.  \nFigure 38 demonstrates the effect of gravity darkening on 3 of the 3 M glyph[circledot] tracks from Figure 34 in the HR diagram. In each panel of Figure 38 the non-rotating track is shown for reference as the dotted line. The magnitude of the gravity darkening effect increases with ω and is more substantial for L than for T eff . Relative to the intrinsic track, the polar projection is brighter and hotter, and the equatorial projection is cooler and fainter.  \nFigure 36. A grid of Roche equipotential surfaces for a range of rotation rates ( ω ) and inclination angles ( i ). The color of the surface corresponds to the variation in ˜ T eff with red corresponding to 0.87 and pale yellow to 1.25; the stars with ω = 0 have ˜ T eff = 1.  \n<!-- image -->  \n<!-- image -->  \nFigure 37. The variation of C T (left panel) and C L (right panel) in the ( ω, i )-plane. Note the color scale is different in the two panels.  \n<!-- image -->', '5. CONVECTIVE BOUNDARIES AND SEMICONVECTION REGIONS': "The correct treatment of convective boundaries continues to be a challenging problem. In this section, we discuss three approaches: the 'sign-change' algorithm (Paper I, Paper II), an improved approach called 'predictive mixing' (Paper IV), and a new 'convective premixing' scheme that addresses several remaining issues.  \nEarly versions of MESA located convective boundaries by searching for sign changes in the discriminant y , defined by y = y S ≡ ∇ rad -∇ ad when the Schwarzschild criterion is used to assess convective stability, or by y = y L ≡ ∇ rad -∇ L when the Ledoux criterion is used; here ∇ rad , ∇ ad and ∇ L are the radiative, adiabatic and Ledoux temperature gradients, respectively. As demonstrated in Paper IV, this sign-change algorithm can fail at convective boundaries that exhibit composition dis-  \n<!-- image -->  \n<!-- image -->  \nFigure 38. Three rotating tracks from Figure 34 showing the effect of gravity darkening in the HR diagram of a 3 M glyph[circledot] model from ZAMS to TAMS. The orange line plots the intrinsic values; the yellow line plots the polar projection; the red line plots the equatorial projection. The dotted line shows the evolution of the non-rotating track.  \n<!-- image -->  \ncally, the failing cases exhibit ∇ rad appreciably larger than ∇ ad on the convective side of the boundary. Instead, as argued by Gabriel et al. (2014), physical consistency within local mixing length theory (MLT) dictates that ∇ rad = ∇ ad should hold on the convective side.  \nIn Paper IV we introduced the predictive mixing scheme for treatment of convective boundaries. This scheme improves on the sign-change algorithm by allowing each convection region to expand during a time step until its boundaries satisfy ∇ rad = ∇ ad on their convective side. The expansion is achieved by modifying convective diffusivities in the cells on the radiative side of a boundary. In Paper IV, we applied predictive mixing in five scenarios: a growing convective core in a 1 . 5 M glyph[circledot] star on the MS; a retreating convective core in a 16 M glyph[circledot] star on the MS; growing convective cores in 1 M glyph[circledot] and 3 M glyph[circledot] stars during the core Helium burning (CHeB) phase; and an evolving convective envelope in a 1 M glyph[circledot] star on the MS. Predictive mixing is able to achieve the desired ∇ rad = ∇ ad outcome in most of these cases.  \nHowever, two cases shown in Paper IV continue to exhibit ∇ rad > ∇ ad on the convective side of a boundary of a convection region, highlighting the need for further work and motivating us to develop a new scheme for treating convective boundaries. This scheme, which we dub 'convective premixing', draws inspiration from earlier work by Castellani et al. (1985) and Mowlavi & Forestini (1994). In the following, we first discuss why predictive mixing sometimes fails. Then, we describe the new convective premixing scheme in detail (Section 5.2), and demonstrate its application in various evolutionary scenarios (Sections 5.3-5.5).", '5.1. The Failure of Predictive Mixing': "The 16M glyph[circledot] MS scenario presented in Figure 4 of Paper IV exhibits a convective shell above the abundancegradient region with ∇ rad > ∇ ad at its lower boundary. If predictive mixing is applied to the shell, the lower boundary advances downward to merge with the core while the upper boundary remains fixed in position. The end result is that the entire abundance-gradient region mixes into the core, delivering significant quantities of fresh H. Such behavior is unphysical, and in Paper IV we made the pragmatic choice to avoid this outcome by disabling predictive mixing for the convective shell.  \nAdditionally, the 1 M glyph[circledot] CHeB scenario presented in Figure 6 of Paper IV exhibits ∇ rad > ∇ ad at the upper boundary of its convective core. This issue cannot be resolved by predictive mixing, as discussed in Paper IV.  \nIn both of these scenarios, the problem encountered is an unavoidable consequence of the design of the predictive mixing scheme, which manipulates the diffusion coefficients at cell faces and then relies on MESA 's abundance solver (Paper I, Section 6.2) to update the model. In contrast, the new convective premixing scheme directly updates the abundances, as we now describe.", '5.2. The Convective Premixing Scheme': 'The convective premixing (CPM) scheme is applied at the start of each time step, before any structural or compositional changes that arise due to the evolution of the star. It proceeds by finding the cells where y > 0 on one face (convective) and y < 0 on the other face (radiative). For each of these initial boundary cells, the algorithm considers whether y on the radiative face would change if the adjacent cell outside the convection region is mixed completely with the rest of the convec-  \nFigure 39. Profiles of ∇ rad , ∇ ad , ∇ L , and X as a function of mass coordinate, in the inner part of the 16 M glyph[circledot] MS star at X c = 0 . 15. The panels show the separate runs described in the text. Gray (gold) shading indicates convection (semiconvection) regions.  \n<!-- image -->  \n/circledot  \n/circledot  \ntion region. This putative mixing is performed at constant cell pressure and temperature, and involves recalculating abundances, densities, opacities and the various temperature gradients ( ∇ rad , ∇ ad , ∇ L ) throughout the convection region plus the adjacent cell.  \nIf the radiative face of the boundary cell becomes convective during this putative mixing, then the mixing is committed to the model, overwriting the composition profile throughout the entire (newly extended) convection region. Then, the next adjacent cell outside the convection region is considered for incorporation. This process continues iteratively until the radiative face of the current convective boundary remains radiative during the putative mixing.  \nAt a given point in its evolution, a star typically exhibits multiple convective boundaries (for instance, the left panel of Figure 39 shows three). The order in which CPM processes these boundaries is determined by evaluating a characteristic mixing timescale dr/v conv for the initial boundary cells; here, v conv is the convective velocity on the y > 0 face of the boundary cell, and dr is its radial extent. The boundary with the smallest time scale is processed first, then the boundary with the next-smallest time scale, and so on.  \nIn CPM the mixing is treated as instantaneous. Given that the convective diffusivity D conv is well in excess of 10 10 cm 2 s -1 even near convective boundaries (see, e.g., Figures 11 and 29 of Paper II), the characteristic mixing timescale is τ ∼ ∆ r 2 /D conv glyph[lessorsimilar] 10 4 yr for a typical region with extent ∆ r glyph[lessorsimilar] R glyph[circledot] . This is small compared to typical nuclear timescales, and so the assumption of instantaneous mixing seems warranted for all but the most rapid evolutionary phases.  \nDuring its iterations, CPM naturally handles the transition of cell faces inside the convection region from convective to radiative. This typically happens either at the boundary opposite to the advancing one (causing that boundary to retreat), or at a point inside the convection region (causing the region to split). Because mixing through a face ceases when it transitions from convective to radiative, newly-transitioned faces should be very close to convective neutrality. To improve how closely neutrality is achieved, our CPM implementation divides the cell outside the advancing boundary into a number of virtual sub-cells. Each of these sub-cells is mixed into the convective region in turn, until all have been incorporated. The number of sub-cells is automatically adjusted to ensure that, during each sub-cell mix, at most a single face within the current convective region transitions to radiative.', '5.3. Evolution of a Retreating Convective Core on the Main Sequence': "We evolve a 16 M glyph[circledot] star from the ZAMS to the TAMS. This is the scenario considered in Section 2.3 of Paper IV, and illustrates the behavior of MS stars with retreating convective cores. Here, and for the scenarios presented in the following sections, we assume an initial He mass fraction Y = 0 . 28 and an initial metal mass fraction Z = 0 . 02, and we ignore rotation and mass loss. Figure 39 plots the profiles of ∇ rad , ∇ ad , ∇ L , and X in the inner parts of the star, at a point nearing the TAMS ( X c = 0 . 15). The left panel illustrates a run using the predictive mixing scheme applied at the boundary of the convective core, while the right panel shows a run using  \nFigure 40. Profiles of ∇ rad , ∇ ad , ∇ L , and X as a function of mass coordinate, in the inner part of the 16 M glyph[circledot] MS star at X c = 0 . 15. The panels show the outcome of an artificial switch from the predictive mixing scheme to the CPM scheme. The initial state (a) is shown in the upper left panel, and the subsequent panels, running in alphabetical order (b),. . . ,(f), show how this state evolves during the CPM iterations. Gray (gold) shading indicates convection (semiconvection) regions.  \n<!-- image -->  \nFigure 41. Mass coordinates of the convective core boundary ( m c ) and the top of the abundance-gradient region ( m a ) as a function of MS age, for the 16 M glyph[circledot] star. Different line colors show the separate runs discussed in the text.  \n<!-- image -->  \nthe CPM scheme; in both calculations, the Ledoux criterion is used to assess convective stability, so that y = y L .  \nIn the left panel, the convective shell discussed in Section 5.1 can clearly be seen at the top of the abundance-gradient region, spanning mass coordinates 6 . 2 glyph[lessorsimilar] m/ M glyph[circledot] glyph[lessorsimilar] 7 . 1. In the right panel, the shell is absent and the abundance-gradient region is wider and shallower, extending all the way out to m/ M glyph[circledot] ≈ 7 . 2. Moreover, the region is very close to adiabatic ( ∇ rad = ∇ ad ) , in contrast to the super-adiabatic stratification ( ∇ rad > ∇ ad ) seen in the left panel.  \nSchwarzschild & Harm (1958) were the first to predict the appearance of adiabatically stratified abundancegradient regions outside the convective cores of massive MS stars, labeling them 'semiconvection regions'. The key feature of a semiconvection region is that the adiabatic stratification is continuously maintained by a gradual adjustment of opacity due to the changing abundance profile.  \nTo illustrate how this adjustment naturally arises within the CPM scheme, we simulate a single evolutionary time step of the 16 M glyph[circledot] star where we artificially switch from the predictive scheme to CPM. The starting configuration is the profile shown in the left panel of Figure 39. This panel is reproduced as the upper-left panel of Figure 40; the subsequent panels then show the evolving profiles at selective intermediate stages during the CPM iterations. A clear narrative emerges from these panels: as the lower boundary of the convective shell advances inward, the cell faces near the upper boundary of the shell transition from convective to radiative, causing that boundary to retreat inward. Because there are no  \nabundance gradients within the shell itself, the Ledoux and Schwarzschild criteria give the same condition for this transition to occur: ∇ rad = ∇ ad . The overall effect is that the shell propagates inward as a whole, leaving behind it a 'wake' with an adiabatic stratification. Eventually, the propagating shell merges with the core, leading to a final state (seen in the lower-right panel of Figure 40) that closely resembles the right panel of Figure 39 (the small differences are because the former has an evolutionary history determined by the predictive mixing scheme, while the latter has a history determined by CPM).  \nThe seminal paper by Schwarzschild & Harm (1958) triggered significant interest in semiconvection regions, with a particular focus on their final stratification. Sakashita & Hayashi (1961) argued that y L = 0 should apply in semiconvection regions, rather than y S = 0 as originally proposed. However, Kato (1966) reasoned that because the former stratification is super-adiabatic ( ∇ rad > ∇ ad ), slow mixing by overstable g-mode oscillations will drive it toward the same y S = 0 outcome. Subsequently, Gabriel (1970) suggested that Kato's mechanism is superfluous, due to the appearance of propagating convective shells that continually adjust the abundance profile to achieve y S = 0. Gabriel's narrative closely mirrors the one we give above, and the correspondence between his Figure 1 and our Figure 40 is striking.  \nA possible source of confusion in this discussion is the relationship between the semiconvection region shown in the right panel of Figure 39, and the semiconvective mixing discussed in Section 4.1 of Paper II. The latter implements the mixing envisaged by Kato (1966); while it will ultimately yield a y S = 0 stratification, it is a fundamentally different mechanism than the propagating convective shells shown in Figure 40. A critical distinction lies in the role played by the convective stability criterion. While our calculations adopt the Ledoux criterion ( y = y L ), repeating them with the Schwarzschild criterion ( y = y S ) leads to results very similar to the ones already shown (although the abundance profiles are rather more jagged). In contrast, Kato's mechanism requires the Ledoux criterion to establish an initially super-adiabatic stratification, which the mechanism then drives toward adiabaticity. In hindsight, it seems prudent to reserve the label 'semiconvection' for the adiabatically stratified regions envisaged by Schwarzschild & Harm (1958), and avoid using it as in Paper II to describe a mechanism that can generate these regions (for additional examples of this conflation, see e.g., Silva Aguirre et al. 2010; Noels et al. 2010; Ding & Li 2014; Moore & Garaud 2016).  \nFigure 42. Period spacings ∆ P plotted as a function of period P , for glyph[lscript] = 2 g-modes of the 16 M glyph[circledot] MS star at X c = 0 . 15. Different symbols show the separate runs discussed in the text.  \n<!-- image -->  \nTo summarize the core and near-core evolution of the 16 M glyph[circledot] star, Figure 41 plots the mass coordinate m c of the convective core boundary as a function of MS age for the separate runs using the predictive mixing and CPM schemes. We also show the mass coordinate m a of the top of the abundance-gradient region outside the core; this coincides with the maximal extent of the core in the predictive mixing case, and to the top of the semiconvection region in the CPM case. During the star's evolution, the opacity change due to progressive H depletion at the center lowers ∇ rad throughout the core, causing its boundary to retreat inward. In the CPM case, this retreat is mirrored by an outward growth of the semiconvection region, a behavior first noted by Schwarzschild & Harm (1958). The growth slowly feeds fresh H from the envelope down into the core, which causes the core to shrink slightly less rapidly than in the predictive mixing case, thereby marginally prolonging the star's MS lifetime.  \nAs the star nears the TAMS, the mass m a -m c of the abundance-gradient region is ≈ 30% larger in the CPM case than with predictive mixing. This region plays a key role in determining the shape of the Brunt-Vaisala frequency near the core, and the two cases should therefore exhibit differences in their g-mode oscillation spectra. We explore this by using release 5.2 of GYRE (Townsend & Teitler 2013; Townsend et al. 2018) to evaluate the star's normal-mode frequencies at X c = 0 . 15. Figure 42 plots the resulting period echelle diagram (period P versus period-spacing ∆ P ) for glyph[lscript] = 2 g modes in the period range 0.5-5 d. Both cases show significant departures from the uniform period spacing ∆ P ≈ 10 4 s  \npredicted by asymptotic theory; these departures are caused by the abundance-gradient region, and therefore are an asteroseismic diagnostic of the star's age (e.g., Miglio et al. 2008). As the figure shows, the CPM scheme yields significantly smaller non-uniformity in ∆ P for P glyph[greaterorsimilar] 2 d than the predictive mixing scheme. These different outcomes are potentially testable by the TESS mission (Ricker et al. 2016), which will observe many massive stars, some with the long time baselines necessary to detect g-modes with multi-day periods. It will first be necessary to explore whether the differences persist when the additional effects of core overshoot and rotation are included.", '5.4. Evolution of the Convective Core During Core He Burning': "We now evolve a 1 M glyph[circledot] star through the CHeB phase; this is the scenario considered in Section 2.4 of Paper IV, and illustrates the behavior of He-burning stars with growing convective cores. Figure 43 plots the profiles of ∇ rad , ∇ ad , and Y in the inner parts of the star, at three points during its evolution corresponding to Y c = 0 . 9, 0.6 and 0.3. The left panels show a run using the predictive mixing scheme, while the right panels show a run using the CPM scheme; in both calculations, the Ledoux criterion is used to assess convective stability.  \nThe left panels, which reprise Figure 6 of Paper IV, show how a local minimum in ∇ rad has developed by Y c = 0 . 6. This minimum is held just above the ∇ ad threshold using the predictive superad thresh control, with the result that the core continues to grow slowly without splitting. With the CPM scheme, a growing semiconvection region develops above a smaller convective core.  \nThe emergence of semiconvection regions during the CHeB phase may have first been proposed by Schwarzschild & Harm (1969), but it was Castellani et al. (1971) and Eggleton (1972) who considered the possibility in detail. Later, Castellani et al. (1985) described a 'concatenated convective mixings' scheme for simulating the formation of the semiconvection region; their Figure 3, which can be regarded as a CHeB analog to our Figure 40, reveals how the core splits to form a convective shell, which then propagates outward leaving a wake with an adiabatic stratification. Mowlavi & Forestini (1994) demonstrated how the Castellani et al. (1985) scheme can be generalized to work in other evolutionary phases, and together these two papers provided the original inspiration for the CPM scheme.  \nTo summarize the core and near-core evolution of the 1 M glyph[circledot] star, Figure 44 plots the mass coordinate m c of the convective core boundary as a function of CHeB age, for  \nFigure 43. Profiles of ∇ rad , ∇ ad , and Y as a function of mass coordinate, in the inner part of the 1 M glyph[circledot] star. The panels correspond to different stages during the CHeB phase: Y c = 0 . 9 (upper panel), Y c = 0 . 6 (middle panel), and Y c = 0 . 3 (lower panel). The left panels show the run using the predictive mixing scheme, and the right panels the run using CPM. Gray (gold) shading indicates convection (semiconvection) regions.  \n<!-- image -->  \n/circledot  \n/circledot  \nFigure 44. Mass coordinates of the convective core boundary ( m c ) and the top of the abundance-gradient region ( m a ) as a function of CHeB age, for the 1 M glyph[circledot] star. Different line colors show the separate runs discussed in the text.  \n<!-- image -->  \nthe separate runs using the predictive mixing and CPM schemes. For the CPM run, we also show the mass coordinate m a of the top of the abundance-gradient region outside the core; there is no abundance-gradient region in the predictive mixing run. The figure shows that the semiconvection region forms in the CPM run at an age ≈ 25 Myr. The top of the semiconvection region then closely tracks the outward growth of the core boundary from the predictive mixing run, through to an age ≈ 95 Myr. At this juncture, oscillations start to occur in m c for the CPM run. These oscillations are due to breathing pulses, which disrupt the semiconvection region and introduce discontinuities in the abundance profile. Because m a becomes ill-defined when the discontinuities appear, we do not plot it beyond this point. After two final, large-amplitude pulses, the star reaches the end of the CHeB at an age ≈ 125 Myr, slightly later than the final age for the predictive mixing run.  \nDebate continues as to whether core breathing pulses during the CHeB phase are physical or numerical (see, e.g., Salaris & Cassisi 2017, and references therein). In the present context, we note that the instantaneous mixing assumed in the CPM scheme will tend to exacerbate the pulses, because it does not account for the finite time required for He ingested at the top of the semiconvection region to be transported down to the core. We are currently considering improvements that address this shortcoming, but in the meantime we recommend that the CPM scheme be used with caution in the late stages ( Y c glyph[lessorsimilar] 0 . 15) of the CHeB phase.", '5.5. Evolution of a Growing Convective Core on the Main Sequence': 'We now evolve a 1 . 5 M glyph[circledot] star from the ZAMS to the TAMS; this is the scenario considered in Section 2.4 of Paper IV, and illustrates the behavior of MS stars with initially growing convective cores. Figure 45 plots the profiles of ∇ rad , ∇ ad , ∇ L and X in the inner parts of the star when X c = 0 . 30. The left panel illustrates a run using the predictive mixing scheme, while the right panel shows a run using the CPM scheme; in both calculations, the Ledoux criterion is used to assess convective stability.  \nThe left panel shows a small super-adiabatic region just outside the core that is stabilized against convection by the abundance gradient ( ∇ ad < ∇ rad < ∇ L ). If slow mixing via the Kato (1966) mechanism were allowed to proceed, this region would eventually transform into a semiconvection region. The right panel shows that the CPM scheme naturally reproduces this semiconvection region. Ledoux (1947) suggests the existence of this semiconvection region (though not labeled as such). The mean molecular weight profile µ ∝ mp 7 / 5 he obtains for the abundance-gradient region outside the core yields an adiabatic stratification ( ∇ rad = ∇ ad ).  \nTo bring the present analysis to a close, Figure 46 plots the mass coordinate m c of the convective core boundary as a function of MS age, for the separate 1 . 5 M glyph[circledot] runs using the predictive mixing and CPM schemes. We also show the mass coordinate m a of the top of the abundance-gradient region outside the core. In the predictive mixing case, this region first appears at an age ≈ 1 . 5 Gyr, when the core boundary reverses direction and begins to retreat inward. In the CPM case, the abundance-gradient region is present from the start, and coincides with the semiconvection region up until an age ≈ 1 . 75 Gyr. At this point, ∇ rad outside the core drops below ∇ ad , and the semiconvection region becomes radiative; then, without any further ingestion of H from the envelope, m a remains fixed.', '6. PARALLEL PERFORMANCE': 'Here we provide updates on the parallel performance of MESAstar . For each test, simulations with different numbers of parallel threads were performed on the same computer with no other CPU-intensive tasks taking place. The tests were performed on one Intel Xeon E5-2699V4 processor with 22 physical cores. Although this processor allows hyperthreading (i.e., two threads running on one physical core), we restrict these tests to one thread per physical core because tests found that enabling hyperthreading results in a performance penalty rather than a benefit.  \n<!-- image -->  \n/circledot  \n/circledot  \nFigure 45. Profiles of ∇ rad , ∇ ad , ∇ L , and X as a function of mass coordinate, in the inner part of the 1 . 5 M glyph[circledot] MS star at X c = 0 . 30. The panels show the separate runs described in the text. Gray (gold) shading indicates convection (semiconvection) regions.Figure 46. Mass coordinates of the convective core boundary ( m c ) and the top of the abundance-gradient region ( m a ) as a function of MS age, for the 1.5 M glyph[circledot] star. Different line colors show the separate runs discussed in the text.  \n<!-- image -->', '6.1. Parallel Scaling of MESAstar': "MESA uses the OpenMP application programming interface to parallelize certain operations. Among these we distinguish three categories. The first category is operations that are parallel per cell in the stellar model, including the EOS, opacity, and nuclear network. The next category concerns parts of the problem that are a mixture of parallel and serial execution. These include matrix manipulations, the Newton-Raphson solver, and atomic diffusion calculations. The final category is those operations that are serial, such as the main evolve loop and the adjustment of the total mass of the star via mass loss or accretion.  \nThree MESA test suite cases are considered: black hole , 7M prems to agb , and 1M pre ms to wd . 7 The black hole test uses 6 structure variables and a nuclear network with 22 isotopes. The 7M prems to AGB test uses 4 structure variables and a nuclear network with 10 isotopes. The 1M pre ms to wd test uses 4 structure variables and a nuclear network with 8 isotopes and includes rotation.  \nFigure 47. Speed up for the three test cases described in the text when run on different numbers of parallel threads from 1 to 22. Lines are illustrations of Equation (46) for different values of p .  \n<!-- image -->  \nAs a metric, we define 'speed up' to be the ratio of the measured run-time on 1 thread to that on N threads. The theoretical speed up,  \n1 (1 -p ) + ( p/s ) , (46)  \nis predicted by Amdahl (1967), where, for our purposes, s is the number of parallel threads and p is the fraction of the code that will benefit from parallel execution for a given test case.  \nFigure 47 shows the speed up for each of the three cases. In general, we expect cases with more variables (both structure and network) to benefit more from parallel execution. Indeed, the black hole case shows the best scalability with p ≈ 0 . 96, followed by 7M prems to AGB with p ≈ 0 . 90, and 1M pre ms to wd with p ≈ 0 . 81. These results are consistent with the numbers of variables included in the respective tests.  \nFor each case, Figure 48 shows the relative fraction for the three categories described above and the speed up. The parts of MESAstar that are parallel-per-cell scale nearly linearly with the number of threads while the mixed category scales weakly and the serial component is essentially flat. The serial portion becomes an increasing fraction of the total runtime as the other categories decrease with an increasing number of parallel threads. This may become a greater issue with the move towards many-core processors in the future.", '6.2. OP Mono Opacities and Radiative Levitation': 'Paper III (Section 9) describes the inclusion of radiative levitation in MESA via the work of Hu et al. (2011). These capabilities were originally developed as part of STARS (Eggleton 1971; Pols et al. 1995) and evaluate the opacity and radiative acceleration using the OP monochromatic opacity tables (Seaton 2005). Due to the differing approaches of STARS and MESA , a number of derivatives were being calculated but not used in the MESA implementation of the opacity routines. By eliminating the evaluation of these unused quantities, and by pre-computing some frequently re-used stimulated emission factors, we achieved at least a factor of 5 reduction in the time required to evaluate opacities and radiative accelerations. These optimizations translate into an improvement in total runtime relative to previous versions of MESA when making use of the OP monochromatic opacities or radiative levitation capabilities without any compromise to the numerical results.  \nWe demonstrate in Figure 49 the difference in execution time for runs of the radiative levitation test suite case using MESA versions before and after the changes described herein. The computational expense  \nof calculating the radiative accelerations continues to dominate the runtime of such models, accounting for more than 65% of run-time even when using 20 parallel threads, meaning these capabilities can benefit from progress toward many-core architectures.', '7. SUMMARY': "We explain significant new capabilities and improvements implemented in MESA since the publication of Paper I, Paper II, Paper III and Paper IV. The addition of the RSP radial pulsation functionality in MESAstar (Section 2) provides a new capability to model radiallypulsating variable stars. Advances to MESA in numerical energy conservation (Section 3), rotation factors and gravity darkening (Section 4 and Appendix B), and convective boundaries (Section 5), will open opportunities for future investigations in stellar evolution. Improvements in the computational efficiency of MESA on current generation multicore x86 instruction set architectures (Section 6) will inform future development directions. Upgrades to the EOS and nuclear reaction physics (Appendix A) will increase the robustness of stellar evolution models. Discussion of the current treatment of fallback and comparisons of the thermodynamic evolution of supernova models from different software instruments (Appendix C) will enhance the study of massive star explosions. Introduction of the MESA Testhub software infrastructure (Appendix D) for webbased, automated, daily examination of the MESAstar and MESAbinary test suites will lead to more efficient source code development. Input files and related materials for all the figures are available at http://mesastar.org and https://doi.org/10.5281/zenodo.2582656.  \nWe thank Warrick Ball and Evan Bauer for their sustained engagment with the MESA project. Both Richa Kundu and Susmita Das graciously shared their variable star calculations. It is a pleasure to thank Conny Aerts, Sean Couch, Franck Delahaye, Luc Dessart, Ebraheem Farag, Carl Fields, Chris Fontes, Chris Fryer, Jim Fuller, Falk Herwig, Thomas Janka, Sam Jones, Sanjib Gupta, Joyce Guzik, Max Katz, John Lattanzio, Abhijit Majumder, Wendell Misch, Joey Mombarg, Viktoriya Morozova, Sterl Phinney, Eliot Quataert, Rene Reifarth, Toshio Suzuki, Katie Mussack Tamashiro, Dean Townsley, Todd Thompson, Suzannah Wood and Mike Zingale for discussions. We also thank the participants of the 2018 MESA Summer School for their willingness to experiment with new capabilities. The improvements discussed in Section 5 were in large parts discussed during RHDT and AT stays at the Munger residence.  \nFigure 48. Breakdown of the fraction of the execution time (upper row) and the speed up (lower row) of three different categories listed in the text. Each column is a single test case. The line in the lower panels is Equation (46) for p = 0 . 9 and p = 1.  \n<!-- image -->  \nThe MESA project is supported by the National Science Foundation (NSF) under the Software Infrastructure for Sustained Innovation program grants (ACI1663684, ACI-1663688, ACI-1663696). This research benefited from interactions that were funded in part by the Gordon and Betty Moore Foundation through Grant GBMF5076 and was also supported at UCSB by the NSF under grant 17-48958. This research was also supported at ASU by the NSF under grant PHY1430152 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-26555. The Center for  \nComputational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.  \nP.M. acknowledges support from NSF grant AST1517753 and the Senior Fellow of the Canadian Institute for Advanced Research (CIFAR) program in Gravity and Extreme Universe, both granted to Vassiliki Kalogera at Northwestern University. S.M.K thanks the Indo-US Science and Technology Forum for financial support. A.T. is a Research Associate at the Belgian Scientific Research Fund (F.R.S-FNRS). R.F is supported by the Netherlands Organisation for Scientific Research (NWO) through a top module 2 grant with project number 614.001.501 (PI de Mink). R.S. acknowledges support from the IdP II 2015 0002 64 grant of the Polish Ministry of Science and Higher Education and from SONATA BIS grant, 2018/30/E/ST9/00598, from the National Science Center, Poland. R.H.D.T. acknowledges support from the NSF under grant AST-1716436.  \nFigure 49. Execution time of the radiative levitation test suite case. The same test was run using MESA versions before and after the changes described herein to demonstrate the difference in performance obtained by these improvements.  \n<!-- image -->  \nM.Z. was supported by the Heising-Simons Foundation through grant #2017-274. J.A.G. is supported by the National Science Foundation Graduate Research Fellowship under grant number 1650114. A.S.J. acknowledges support from the Gordon and Betty Moore Foundation under Grant GBMF7392. This work was in part carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. This paper is based upon work supported by the National Aeronautics and Space Administration (NASA) under Contract No. NNG16PJ26C issued through the WFIRST Science Investigation Teams Program. This research made extensive use of the SAO/NASA Astrophysics Data System (ADS).  \nSoftware: gnuplot (Williams & Kelley 2015), ipython/jupyter (P'erez & Granger 2007; Kluyver et al. 2016), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), and Python frompython.org.", 'A. UPDATES TO PHYSICS MODULES': 'A.1. Equation of State  \nThe EOS is evaluated by the eos module. Figure 50 shows the default coverage in the ρ -T plane. The new PTEH option extends the eos coverage to lower densities (to 10 -18 g cm -3 ) and higher metallicities ( Z up to 1.0) than allowed by OPAL ( ρ glyph[greaterorsimilar] 10 -10 g cm -3 and Z ≤ 0.04 only). Previous versions of eos use HELM to provide approximate results for the low density or Z > 0 . 04 cases now covered by PTEH. The PTEH tables are created using the approach of Pols et al. (1995) as implemented by Paxton (2004) in a program derived from Eggleton (1971). PTEH includes a solution to the Saha equations to obtain the dissociation and ionization stages H + , H, H 2 , He, He + , and He ++ and assumes full ionization of C, N, O, Ne, Mg, Si, and Fe.  \nIn addition to PTEH, there are two other new eos options, DT2 and ELM 8 . These are motivated by the need for more numerically accurate partials as discussed in Section 3. In this context, the desired numerical accuracy of the partials is achieved by evaluating analytic partials of the interpolating polynomials rather than by interpolating values of tabulated partials (as is done with OPAL and SCVH data in MESA ). As a result, the partials correspond to how the interpolated eos values will actually change in response to small changes of the parameters, whereas interpolated values of partials will be less accurate predictors of the response to such variations. All three new options use bicubic spline interpolation in high resolution tables of log( p gas /erg cm -3 ), log( e /erg g -1 ), and log( s /erg g -1 K -1 ) to obtain first and second partial derivatives of these quantities with respect to log( ρ /g cm -3 ) and log( T /K). The options DT2 and ELM use tables holding values derived from OPAL/SCVH and HELM respectively. Figure 51 shows ∇ ad in a region of the ρ -T plane covered by DT2 and PTEH.', 'A.1.1. Blending': 'The control structure for blending the various EOS sources in Figure 50 considers a particular source EOS (PTEH, PC, HELM, etc.) to determine what fraction f of the final result comes from that source. A recursive calling structure is used:  \nlevel 1: If f PTEH =0, then return result of level 2. If f PTEH =1, then evaluate PTEH and return. Otherwise, call level 2, blend result with that of PTEH, and return.  \nFigure 50. The ρ -T coverage of the EOS used by the eos module. PTEH is from Pols et al. (1995), HELM is from Timmes & Swesty (2000), PC is from Potekhin & Chabrier (2010), OPAL is from Rogers & Nayfonov (2002), SCVH is from Saumon et al. (1995), and the low-density cold region in the lower left is treated as an ideal neutral gas. The region between SCVH and PC is currently problematic from input physics and numerical perspectives and treated as an ideal gas (see Chabrier et al. 2019 for a recent treatment that is not yet in MESA ). The blue curve shows the profile of a 25 M glyph[circledot] star that has reached an iron core infall speed of 1,000 km s -1 and the purple curve shows the profile of a 0.8 M glyph[circledot] WD.  \n<!-- image -->  \nFigure 51. Adiabatic gradient in the ρ -T plane for X = 0 . 70 and Z = 0 . 02. Regions undergoing H 2 , H, He, or He + dissociation/ionization are colored blue and labelled. Previous versions of eos truncated the H 2 dissociation region at log( ρ /g cm -3 )= -10, the limit of the SCVH data in eos . Other ionization bands occur at a high enough temperature to be mainly covered by OPAL.  \n<!-- image -->  \nlevel 2: If f PC =0, then return result of level 3. If f PC =1, then evaluate PC and return result to level 1. Otherwise, call level 3, blend result with that of PC, and return the blended result to level 1.  \nlevel 3: If f PTEH , Z = 0 (i.e., Z glyph[lessorsimilar] 0.04), then call level 4 with DT2 and return the result to level 2. If f PTEH , Z = 1 (i.e., Z glyph[greaterorsimilar] 0.04), then call level 4 with PTEH and return the result to level 2. Otherwise (i.e., Z near 0 . 04), call level 4 with each of DT2 and PTEH, blend the result in Z and return the blended result to level 2.  \nFigure 52. Relative difference in χ ρ, gas between the eos derivative and a Richardson limit based numerical derivative across the ρ -T plane. The left panel shows the results for the OPAL/SCVH and HELM options, and the right panel shows the results for the PTEH, DT2, and ELM options. Relative differences for the new eos options are glyph[similarequal] 10 -7 , except for a region between SCVH and PC. Note the large difference in accuracy between the old OPAL/SCVH options using interpolated partials and the new DT2 option using partials of interpolating polynomials.  \n<!-- image -->  \nFigure 53. Thermodynamic consistency metric dpe ; see Equation (A2) for the definition. The left panel shows the results for OPAL/SCVH and HELM, and the right panel shows the results for PTEH, DT2, and ELM. As expected, HELM gives machine precision consistency and so is superior in this to ELM. OPAL and SCVH have an advantage over DT2 because they use interpolated values of thermodynamically consistent partials -which degrades their numerical accuracy as shown in Figure 52.  \n<!-- image -->  \nlevel 4: If f PTEH / DT2 =0, then return result of level 5. If f PTEH / DT2 =1, then evaluate PTEH/DT2 9 and return result to level 3. Otherwise, call level 5, blend the result with that of PTEH/DT2 and return the blended result to level 3.  \nlevel 5: If f ELM =0, then return result of level 6. If f ELM =1, then evaluate ELM and return result to level 4. Otherwise, call level 6, blend the result with that of ELM, and return the blended result to level 4.  \nlevel 6: Evaluate HELM and return result to level 5.  \nThe blends use a quintic polynomial with zero slope at boundaries. The partial derivatives of the blend polynomial are included in the calculation of the final result to maintain high numerical consistency in the blend region and across EOS region boundaries.  \nThere remain challenges to providing a broad coverage EOS given the current need to combine multiple sources, some of which may not provide the necessary thermodynamic information. For example, when the log( ρ /g cm -3 ) blend region extended from 3.0 to 3.1, a negative χ ρ, gas resulted because log( p gas /erg cm -3 ) from DT2 were slightly greater than log( p gas /erg cm -3 ) from ELM. This could give a drop in p gas as log( ρ /g cm -3 ) transitions from DT2 to ELM. Extending the blend region from 2.98 to 3.12 is enough to ensure χ ρ, gas > 0 in the DT2 to ELM transition.', 'A.1.2. Numerical Accuracy of Partials': 'To check the numerical accuracy of partials using the new options, we have compared their results with iteratively acquired high precision numerical derivatives (Ridders 1982; Press et al. 1992). For example, Figure 52 shows the relative difference between the eos derivative and a Richardson iterative numerical derivative across the ρ -T plane for χ ρ, gas . The right panel shows that new options give relative errors of glyph[similarequal] 10 -7 , while the left panel shows that the previous options have larger errors, particularly in the OPAL/SCVH regions. As mentioned in Section 3, those larger errors limit the ability of the MESAstar Newton-Raphson solver to reduce residuals and that in turn can lead to increased errors in numerical energy conservation.', 'A.1.3. Thermodynamic Consistency': "The first law of thermodynamics is an exact differential, which implies  \np = ρ 2 ∂e ∂ρ ∣ ∣ ∣ ∣ ∣ T,Y i + T ∂p ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i , ∂e ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i = T ∂s ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i , -∂s ∂ρ ∣ ∣ ∣ ∣ ∣ T,Y i = 1 ρ 2 ∂p ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i . (A1)  \nAn EOS is thermodynamically consistent if these relations are satisfied. Thermodynamic inconsistency may manifest itself as an artificial buildup or decay of the entropy during what should be an adiabatic flow. Models that are sensitive to the entropy may suffer inaccuracies if thermodynamic consistency is systematically violated over sufficiently long timescales. Equation (A1) may be recast in a form suitable for evaluating numerical inconsistencies  \ndpe = ρ 2 p ∂e ∂ρ ∣ ∣ ∣ ∣ ∣ T,Y i + T p ∂p ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i -1 , dse = T   ∂s ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i / ∂e ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i   -1 , dsp = -ρ 2   ∂s ∂ρ ∣ ∣ ∣ ∣ ∣ T,Y i / ∂p ∂T ∣ ∣ ∣ ∣ ∣ ρ,Y i   -1 . (A2)  \nIdeally, dpe , dse , and dsp are zero. Figure 53 shows the first thermodynamic consistency quantity, dpe , across the ρ -T plane for an older ( MESA r8845) and current version of eos . The other two thermodynamic consistency metrics, dse and dsp , show similar magnitudes. In general, the thermodynamic consistency with DT2 and ELM is reduced relative to the older options directly using OPAL/SCVH and HELM. This is because the thermodynamic consistency relations can only be approximated by bicubic splines. If possible, Hermite interpolation of the Helmholtz free energy would make use of partial derivatives from the EOS and guarantee thermodynamic consistency (e.g., Timmes & Swesty 2000). However, a bicubic Hermite interpolation produces discontinuities in second derivative quantities (e.g., ∂p gas /∂ρ | T,Y i ) that are problematic for MESA 's Newton-Raphson solver and a biquintic Hermite interpolation requires partial derivatives that are unavailable from some constituents of the EOS patchwork. So for now, we compromise by providing the previous options that give better thermodynamic consistency and the new options that provide better numerical accuracy of partials. The impact of the thermodynamic inconsistencies of the new approach must be evaluated on a case-by-case basis.", 'A.2.1. Nuclear Reaction Rates': "The Joint Institute for Nuclear Astrophysics (JINA) REACLIB library, which provides the default nuclear reaction rates for MESA , has been updated from the jina reaclib results v2.2 snapshot to the default snapshot 10 . This update includes changes to the fitting formula for a few neutron capture rates that previously returned erroneous values at T glyph[greaterorsimilar] 8 × 10 9 K. We have modified the default snapshot to include a missing, temperature-independent 26 Al → 26 Mg weak reaction rate. Reaction rates between the 26 Al ground and meta-stable states now use Gupta & Meyer (2001). Nuclear partition functions use JINA winvne v2.0.dat table 11 and JINA masslib library 5.data 12 provides atomic masses.  \nA cell's temperature may exceed log( T/ K)=10.0 in shocks and explosive burning, which is beyond the range of validity for the fits to the reaction rates and the partition functions. Previously, MESA extrapolated for log( T/ K) > 10.0, leading to erroneous reaction rates. Reaction rates are now set equal to their log( T/ K)=10.0 values when log( T/ K) > 10.0.  \nMESA now includes the option to use the electron-capture and β -decay rates from Suzuki et al. (2016), which cover sd -shell nuclei with A = 17 -28. The primary application for these tables is the evolution of high-density oxygen-neon cores (e.g., Miyaji et al. 1980; Miyaji & Nomoto 1987; Jones et al. 2013). Compared to the on-the-fly weak rate approach described in Paper III, these tabulated rates are less computationally expensive and more accurate at high temperatures ( T glyph[greaterorsimilar] 10 9 K), but less accurate at low temperatures ( T glyph[lessorsimilar] 10 8 K) and do not allow explicit updating of nuclear physics.", 'A.2.2. Reaction Rate Screening': 'The plasma coupling parameter of two reactants Γ i,j and the ion sphere radius a i are  \nΓ i,j = Z i Z j e 2 0 . 5 ( a i + a j ) k B T , a i = [ 3 Z i 4 π ( Z i n i + Z j n j ) ] 1 / 3 , (A3)  \nwhere Z i is the atomic charge of isotope i , e is the electron charge, k B is the Boltzmann constant, and n i is the ion number density of species i . MESA applies screening factors to correct nuclear reaction rates for plasma interactions (e.g., Salpeter 1954). Previous versions defaulted to using one set of expressions for the weak screening regime (Γ i,j ≤ 0.3, Dewitt et al. 1973; Graboske et al. 1973) and another set of expressions for the strong screening regime (Γ i,j ≥ 0.8, Alastuey & Jancovici 1978; Itoh et al. 1979). A linear blend of the weak and strong screening factors is used in the 0.3 < Γ i,j < 0.8 intermediate regime. Silicon burning reactions in the cores of massive stars often operate in this intermediate regime. The numerical blending provides a smooth and continuous function that equals the non-blended values and their derivatives at the edges. A new default for screening, based on Chugunov et al. (2007), includes a physical parametrization for the intermediate screening regime and reduces to the familiar weak and strong limits at small and large Γ values. We extend the Chugunov et al. (2007) one-component plasma results to a multi-component plasma following Itoh et al. (1979), where the Z i are replaced with the average charge ¯ Z .  \nFigure 54 compares three MESA nuclear reaction rate screening options on the evolution of a 25 M glyph[circledot] model. The differences are small from H-burning to the onset of core collapse. However, the Chugunov et al. (2007) implementation takes ≈ 20% fewer time steps, retries, and backups which indicates a numerically smoother solution.', 'B. ANALYTICAL APPROXIMATIONS TO THE ROCHE GEOMETRY OF A SINGLE STAR': "In this appendix we compute various properties of the Roche potential of a single star. These are used for the computation of centrifugal effects in stellar structure, as discussed in Section 4. Through these derivations we denote dimensionless properties using an apostrophe, with distances being normalized as r ' = r/ ( Gm Ψ / Ω 2 ) 1 / 3 and the potential as Ψ ' = Ψ / ( Gm Ψ Ω) 2 / 3 .  \nFigure 54. Evolution in the central ( ρ -T ) plane of a 25 M glyph[circledot] star until the onset of Fe core-collapse with three different nuclear reaction rate screening options. 'Previous' denotes the previous MESA default option, 'Chugunov' is the new default option, and 'None' applies no screening correction to any nuclear reaction rate. Locations when a fuel undergoes central ignition are labelled.  \n<!-- image -->  \n(  \n)", 'B.1. Volume of Roche Equipotentials': "Following the diagram in Figure 55, the dimensionless volume equivalent radius r ' Ψ can be computed in terms of ω as  \n4 3 πr ' 3 Ψ ≡ V ' Ψ = 2 ∫ ω 2 / 3 0 dV ' dx ' dx ' , (B4) dV ' dx ' =2 πx ' y ' , y ' = √ ( 1 ω 2 / 3 + ω 4 / 3 2 -x ' 2 2 ) -2 -x ' 2 ,  \nwhere the expression for y ' ( x ' ) can be derived directly from the Roche potential. The integral can be solved analytically for the case of ω = 1, providing the volume of a critically rotating star in terms of its equatorial radius (Kopal 1959),  \nV Ψ ( ω = 1)= 4 π 3 r 3 e [ 3 √ 3 -4 + 3 ln ( 3( √ 3 -1) √ 3 + 1 )] glyph[similarequal] 4 π 3 (0 . 8149 r e ) 3 . (B5)  \nIn the opposite case where ω → 0, the equipotential is well approximated by an ellipsoid of revolution with  \nr p ( ω ) = r e ( 1 -ω 2 2 + O ( ω 4 ) ) , (B6)  \nsuch that its volume is  \nV Ψ ( ω ) = 4 π 3 r 3 e ( 1 -ω 2 2 + O ( ω 4 ) ) . (B7)  \nBy numerically integrating Equation (B4), a simple polynomial approximation to the volume can then be constructed, which is consistent with the value at critical rotation given by Equation (B5) and the asymptotic behaviour at small  \nFigure 55. Quantities used for the integration of the f P and f T factors, which require the computation of volume and surface areas of equipotentials.  \n<!-- image -->  \nω given by Equation (B7):  \nV Ψ ( ω ) glyph[similarequal] 4 π 3 r 3 e ( 1 -ω 2 2 +0 . 1149 ω 4 -0 . 07376 ω 6 ) . (B8)  \nThis expression has an error < 0 . 25% for 0 ≤ ω ≤ 1. In all the asymptotic expressions considered in this work ω appears in series of even powers, so we do not include odd terms in any fit. Also, since the value for ω = 1 is fixed, Equation (B8) is a fit with only 1 free parameter.  \nStellar evolution instruments typically use the radial coordinate r Ψ , so it is useful to have polynomial fits for r Ψ ( ω ) as well. In the limit of ω → 0,  \nr Ψ = r e ( 1 -ω 2 6 + O ( ω 4 ) ) . (B9)  \nA polynomial fit that matches this and the value at critical rotation 0 . 8149 r e from Equation (B5) is  \nr Ψ = r e ( 1 -ω 2 6 +0 . 01726 ω 4 -0 . 03569 ω 6 ) , (B10)  \nwhich has an error < 0 . 15% for 0 ≤ ω ≤ 1. Similarly, the expression  \nr e = r Ψ ( 1 + ω 2 6 -0 . 0002507 ω 4 +0 . 06075 ω 6 ) (B11)  \nhas an error < 0 . 2% in the same range.", 'B.2. Surface Area of Roche Equipotentials': "The computation of the dimensionless surface area S ' Ψ is given by  \nS ' Ψ = 2 ∫ ω 2 / 3 0 dS ' dx ' dx ' , dS ' dx ' = 2 πx ' √ ( dy ' dx ' ) 2 +1 . (B12)  \nIn the limit ω → 0, dS ' /dx ' can be approximated as  \nwhich upon integration provides the approximate form of the surface area of a slowly rotating equipotential,  \ndS ' dx ' ( ω ) = 2 πx ' r ' e √ r ' 2 e -x ' 2 [ 1 -1 2 ( x ' r ' e ) 2 ω 2 + O ( ω 4 ) ] , (B13)  \nS Ψ ( ω ) = 4 πr 2 e ( 1 -ω 2 3 + O ( ω 4 ) ) . (B14)  \nThis result is consistent with that for an oblate spheroid. Using Equation (B14) and the numerically computed value S Ψ ( ω = 1) = 8 . 832 r 2 e , a fit to S Ψ that has an error < 0 . 01% in the range 0 < ω < 1 is  \nS Ψ ( ω ) glyph[similarequal] 4 πr 2 e ( 1 -ω 2 3 +0 . 08525 ω 4 -0 . 04908 ω 6 ) . (B15)", 'B.3. Surface Averages of Gravity': "The surface average of the dimensionless gravity 〈 g ' 〉 is computed as  \n〈 g ' 〉 = 2 S ' Ψ ∫ ω 2 / 3 0 g ' dS ' dx ' dx ' , (B16)  \nwith an equivalent expression for 〈 g '-1 〉 . For reference, the dimensionless gravity is given by g ' = g/ ( Gm Ψ Ω 4 ) 1 / 3 . In the limit of ω → 0, it can be shown that  \ng ' ( ω ) = 1 r ' 2 e { 1 + [ 1 -2 ( x ' r ' e ) 2 ] ω 2 } , (B17)  \nwhich combined with Equations (B13) and (B16) results in  \nS Ψ 〈 g 〉 = 4 πGm Ψ ( 1 -2 3 ω 2 + O ( ω 4 ) ) . (B18)  \nBy numerically computing this integral, a simple fit that matches the computed value at ω = 1 and has an error < 0 . 35% for 0 ≤ ω ≤ 1 is  \nS Ψ 〈 g 〉 = 4 πGm Ψ ( 1 -2 3 ω 2 +0 . 3045 ω 4 +0 . 001382 ω 6 ) . (B19)  \nSimilarly the average of the inverse gravity, in the limit of ω → 0, can be shown to be  \nS Ψ 〈 g -1 〉 = 4 πr 4 e Gm Ψ ( 1 + O ( ω 4 ) ) . (B20)  \nHowever, this integral diverges for ω = 1, as for a critically rotating star the effective gravity at its equator becomes zero. Although the expression for S Ψ 〈 g '-1 〉 cannot be integrated in the limit ω → 1, by comparing it to the numerical results we have verified that it is approximately given by S Ψ 〈 g -1 〉 ∝ -ln(1 -ω 4 ). Combining this information with Equation (B20) we have found the following fit with an error < 0 . 85% in the range 0 ≤ ω ≤ 0 . 9999:  \nS Ψ 〈 g -1 〉 = 4 πr 4 e Gm Ψ A ( ω ) , A ( ω ) = 1 -0 . 1076 ω 4 -0 . 2336 ω 6 -0 . 5583 ln(1 -ω 4 ) . (B21)", 'B.4. Moments of Inertia': "The specific moment of inertia i rot is needed to determine Ω Ψ from the specific angular momentum j rot and volume equivalent radius r Ψ . To compute i rot , consider a shell of material extending from Ψ to Ψ + d Ψ. At each point in its surface, its thickness is given by  \nds = |∇ Ψ | -1 d Ψ = g -1 d Ψ . (B22)  \nAssuming a constant density ρ in the shell (as in the shellular approximation),  \ndm ' =2 ρ ' d Ψ ' ∫ ω 2 / 3 0 g '-1 dS ' dx ' dx ' = ρ ' d Ψ S ' Ψ 〈 g '-1 〉 , (B23)  \ni ' rot = 2 ∫ ω 2 / 3 0 x 2 g '-1 dS ' dx ' dx ' S ' Ψ 〈 g '-1 〉 , (B24)  \nwhere dm ' = dm/m Ψ and ρ ' = ρ/ [ m Ψ ( Gm Ψ / Ω 2 ) -3 / 2 ] . From Equation (B23),  \nwith the dimensionless specific moment of inertia defined as i ' rot = i rot / ( Gm Ψ / Ω 2 ). Preserving the fit for S Ψ 〈 g -1 〉 given by Equation (B21) and using Equations (B13) and (B17), in the limit of ω → 0  \ni rot ( ω ) = 2 3 r 2 e 1 + 1 5 ω 2 + O ( ω 4 ) A ( ω ) . (B25)  \nEquation (B22) implies that i rot ( ω = 1) = r 2 e , as in the extreme of critical rotation g = 0 at the equator and almost all mass between two close equipotentials lies in a ring of radius r 2 e . Using this information, we construct the following fit which has an error < 0 . 9% for 0 ≤ ω ≤ 0 . 9999:  \ni rot = 2 3 r 2 e × B ( ω ) A ( ω ) , B ( ω ) = 1 + 1 5 ω 2 -0 . 2735 ω 4 -0 . 4327 ω 6 -3 2 × 0 . 5583 ln(1 -ω 4 ) , (B26)  \nwhere the 3 / 2 factor in the last term ensures the desired result as ω → 1.", 'B.5. Computation of f P And f T': 'Using all the fits constructed so far, f P and f T can be evaluated directly. However, to provide a more compact expression, we keep only the fit for S Ψ 〈 g -1 〉 and use Equations (35), (B10) and (B18) to determine the behavior of the remaining terms when ω → 0,  \nf P ( ω ) = 1 -2 3 ω 2 + O ( ω 4 ) A ( ω ) , f T ( ω ) = 1 + O ( ω 4 ) A ( ω ) . (B27)  \nThe following fits for f P and f T are derived, which have errors < 0 . 8% and < 1 . 6% in the range 0 ≤ ω ≤ 0 . 9999 respectively:  \nf P ( ω ) = 1 -2 3 ω 2 -0 . 06837 ω 4 -0 . 2495 ω 6 A ( ω ) , f T ( ω ) = 1 + 0 . 2185 ω 4 -0 . 1109 ω 6 A ( ω ) . (B28)', 'B.6. Determination of ω': 'For given values of r Ψ , m Ψ , and j rot , ω can be determined from the implicit equation  \nj rot √ Gm Ψ r Ψ = ω i rot r 2 e √ r e r Ψ . (B29)  \nThe left hand side can be directly evaluated, while the right hand side is a monotonic function of ω for 0 ≤ ω ≤ 1.  \nWe compute a fit to the right hand side of Equation (B29). Equations (B9) and (B25) can be used to determine the form of this term in the limit ω → 0. In the limit of ω → 1 all material is concentrated in a equatorial ring, such that j rot ( ω = 1) = √ Gm Ψ r e . Using this information, we find the following fit which has an error < 0 . 8% in the range 0 ≤ ω ≤ 0 . 9999:  \nj rot √ Gm Ψ r Ψ = 2 3 ωC ( ω ) A ( ω ) , C ( ω ) = 1 + 17 60 ω 2 -0 . 3436 ω 4 -0 . 4055 ω 6 -0 . 9277 ln(1 -ω 4 ) . (B30)  \nThis allows the computation of partial derivatives of ω with respect to r Ψ and j rot ,  \n∂ω ∂r Ψ = -3 4 j rot r 3 / 2 √ Gm Ψ [ d dω ( ωC ( ω ) A ( ω ) )] -1 , ∂ω ∂j rot = 3 2 1 √ Gm Ψ r Ψ [ d dω ( ωC ( ω ) A ( ω ) )] -1 . (B31)', 'C. CORE COLLAPSE SUPERNOVA EXPLOSIONS': "Paper IV described modeling the evolution of core-collapse supernova (SN) ejecta up to shock breakout with MESA , and with STELLA beyond shock breakout. Modifications since Paper IV have focused on fallback in weak explosions of red supergiant (RSG) stars. In these weak explosions, the total final explosion energy is positive, but insufficient to unbind all material. Thus, some material falls back onto the central object during the subsequent evolution. To quantify and remove this material, we introduce two new user controls. First, we implement a new criterion to select which material is excised from the model during the ejecta evolution. 13 At each time step, MESA calculates the integrated total energy from the innermost cell to cell j above it:  \nE j = j ∑ i =inner ( e i -Gm i r i + 1 2 u 2 i ) dm i . (C32)  \nIf E j < 0, then there is a bound inner region, and MESA continues this sum outward until it reaches a cell k with local positive total energy ( e k -Gm k /r k + u 2 k / 2 > 0). MESA removes material inside this cell, making cell k the new innermost cell. Second, to remove any slow-moving, nearly hydrostatic material left near the inner boundary, a minimum innermost velocity can be specified at handoff to STELLA . All material below the innermost cell that has a velocity above this specified velocity is not included in STELLA input files. 14 A velocity cut between 100 - 500 km s -1 has little effect on light curve properties and the photospheric evolution of Type IIP SNe, and can greatly reduce numerical artifacts that may arise from an inward-propagating shock hitting the inner boundary in STELLA . Such a cut can also lead to a factor of 10 or more reduction in the number of time steps required to produce a light curve. While this scheme is useful in quantifying and excising fallback material, it is not a satisfactory treatment of fallback. For a more complete description of these fallback criteria, see Appendix A of Goldberg et al. (2019).  \nNext, we look at the evolution of material deep within the ejecta of a Type IIP supernova, at such large optical depths that the outcome is not sensitive to any particular treatment of radiation transport in different software instruments. For this restricted regime, we compare results using MESA+STELLA to quantities derived from the open-source 1D gray radiation hydrodynamics code SNEC (Morozova et al. 2015), both using the same MESA Type IIP ejecta model at shock breakout. This should yield meaningful density and velocity comparisons nearly everywhere in the ejecta and meaningful temperature comparisons very deep within the ejecta.  \nWe explode (with E exp =10 51 erg) the 99em 19 RSG progenitor model from Paper IV, which has a ZAMS mass of 19 M glyph[circledot] , and a mass of 17.8 M glyph[circledot] and a radius of 603 R glyph[circledot] at time of explosion. We excise the inner 1.5 M glyph[circledot] and explode the remaining 16.3 M glyph[circledot] of ejecta using a thermal bomb as described in Paper IV. The mass of radioactive 56 Ni is set to M Ni =0.03M glyph[circledot] . The resulting ejecta is influenced by the inclusion of the Duffell (2016) prescription for mixing due to the Raleigh-Taylor instability in MESA . We then pass our MESA models at shock breakout to both STELLA and SNEC . To mimic the additional line opacities from metals, SNEC employs an opacity floor. We use SNEC 's default opacity floor recommended in Bersten et al. (2011), which is set to κ floor , core =0.24cm 2 g -1 for the metal-rich core material, κ floor , envelope =0.01cm 2 g -1 for the envelope, and proportional to metallicity for the intermediate region.  \nFigure 56 shows density and temperature profiles as a function of mass coordinate, with color corresponding to the number of days after shock breakout. Both instruments agree in the density evolution of the expanding ejecta. The temperature evolution is also in agreement in the very optically thick inner ejecta. Temperature differences at the surface reflect the differing treatments of opacity and radiation transfer between SNEC and STELLA .  \nFigure 57 focuses on the deep core during the first days of the evolution post shock-breakout at the location of the reverse shock generated at the H/He boundary of the initial model. Density, velocity, and temperature profiles over the first 20 days are plotted with the corresponding MESA profile at shock breakout, which serves as the common input in both STELLA and SNEC . By day 20 the reverse shock has reached the inner boundary in both models. Although there are slight differences within the first day, both instruments agree on global properties of the reverse shock and its effects on the temperature and density profiles out to day 20.", 'D. MESA TESTHUB': 'Development of the MESA source code is a collaborative process with multiple commits each day from developers working on separate parts of the codebase. In addition to serving as starting templates for science projects, the test cases in MESAstar and MESAbinary exist to detect when changes to the codebase cause unintended deviations from expected behavior (i.e., bugs). The number of test cases grows with time and is currently more than 100, with a total runtime on the order of 10 hours on multicore workstations. Since MESA is committed to supporting reproducibility by giving bit-for-bit identical results on a variety of different hardware and software platforms (Paper III, Section 10), the test suite must be checked on a representative sample of host systems; just as it takes a team to create MESA , it takes a team to test it.  \nTo prevent slowdowns in development that would be caused by running the test suite on multiple hosts before every commit, we have developed the MESA Testhub (testhub.mesastar.org). The Testhub is a web application that collects and organizes the results of test suite submissions via a companion Ruby gem called mesa test . Every day, submissions from multiple computers and clusters with diverse hardware, operating systems, and compilers check out the most recent revision of MESA , run the test suite, and upload their results to the Testhub. Each day, a summary e-mail is sent to developers detailing which, if any, revisions had failing test cases submitted in the previous 24 hours. With  \nFigure 57. Density (upper panel), velocity (middle panel), and temperature (lower panel) profiles of the 99em 19 progenitor, exploded with 10 51 erg, up to 20 days after shock breakout. The reverse shock originating at the H/He boundary makes its way back through the expanding ejecta. The MESA profile at shock breakout (thick gray line) is used as the input for subsequent evolution in STELLA (solid curves) and SNEC (dashed curves). Numbers in the legend correspond to the day post shock-breakout for each profile.  \n<!-- image -->  \n/circledot  \nFigure 56. Density (upper panel) and temperature (lower panel) profiles in STELLA (solid colored lines) and SNEC (dashed lines) on the plateau of a Type IIP SN model with M Ni = 0 . 03 M glyph[circledot] at 5 to 70 days after shock breakout in MESA . Lighter colors indicate later days.  \n<!-- image -->  \n/circledot  \na quick check of the daily email from the Testhub, developers are able to detect cases that fail to give the expected output with bit-for-bit identical results on different computers, or take more than the specified number of time steps, retries, or backups. With daily coverage, we can promptly diagnose issues as soon as they arise by looking at changes from a only handful of commits while maintaining a brisk development pace. The addition of the Testhub has yielded a significant improvement in the pace and quality of MESA development.', 'REFERENCES': "Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, PhRvL, 118, 221101,  \ndoi: 10.1103/PhysRevLett.118.221101  \n- -. 2017b, ApJ, 851, L35, doi: 10.3847/2041-8213/aa9f0c\n- -. 2017c, PhRvL, 119, 141101,  \ndoi: 10.1103/PhysRevLett.119.141101  \n- -. 2017d, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9\n- -. 2017e, ApJ, 851, L16, doi: 10.3847/2041-8213/aa9a35\n- -. 2017f, ApJ, 850, L40, doi: 10.3847/2041-8213/aa93fc\n- Aikawa, T., & Simon, N. R. 1983, ApJ, 273, 346, doi: 10.1086/161373\n- Alastuey, A., & Jancovici, B. 1978, ApJ, 226, 1034, doi: 10.1086/156681\n- Alcock, C., Alves, D. R., Becker, A., et al. 2003, ApJ, 598, 597, doi: 10.1086/378689\n- Amdahl, G. M. 1967, in Proceedings of the April 18-20, 1967, Spring Joint Computer Conference, AFIPS '67 (Spring) (New York, NY, USA: ACM), 483-485. http://doi.acm.org/10.1145/1465482.1465560\n- Artigau, ' E., Sivaramakrishnan, A., Greenbaum, A. Z., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, 914340\n- Astropy Collaboration, Price-Whelan, A. M., Sip'ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411, doi: 10.1051/0004-6361/200810860\n- Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58, doi: 10.3847/1538-3881/aacb21\n- Ball, W. H., Chaplin, W. J., Schofield, M., et al. 2018, The Astrophysical Journal Supplement Series, 239, 34, doi: 10.3847/1538-4365/aaedbc  \nBeichman, C. A., Rieke, M., Eisenstein, D., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, 84422N  \nBenedict, G. F., McArthur, B. E., Fredrick, L. W., et al. 2002, AJ, 123, 473, doi: 10.1086/338087  \nBersten, M. C., Benvenuto, O., & Hamuy, M. 2011, ApJ, 729, 61, doi: 10.1088/0004-637X/729/1/61  \nBlaˇzko, S. 1907, Astronomische Nachrichten, 175, 325,  \n- doi: 10.1002/asna.19071752002 Bono, G., Marconi, M., & Stellingwerf, R. F. 2000, A&A, 360, 245. https://arxiv.org/abs/astro-ph/0006229 Bono, G., & Stellingwerf, R. F. 1994, ApJS, 93, 233, doi: 10.1086/192054 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402 Breitfelder, J., M'erand, A., Kervella, P., et al. 2016, A&A, 587, A117, doi: 10.1051/0004-6361/201527030 Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, in European Physical Journal Web of Conferences, Vol. 47, European Physical Journal Web of Conferences, 03005 Buchler, J. R. 2009, in American Institute of Physics Conference Series, Vol. 1170, American Institute of Physics Conference Series, ed. J. A. Guzik & P. A. Bradley, 51-58 Buchler, J. R., Koll'ath, Z., & Cadmus, Jr., R. R. 2004, ApJ, 613, 532, doi: 10.1086/422903 Buchler, J. R., Kollath, Z., Serre, T., & Mattei, J. 1996, ApJ, 462, 489, doi: 10.1086/177167 Buchler, J. R., & Kovacs, G. 1987, ApJL, 320, L57, doi: 10.1086/184976 Buchler, J. R., & Moskalik, P. 1992, ApJ, 391, 736, doi: 10.1086/171384 Buchler, J. R., Moskalik, P., & Kovacs, G. 1990, ApJ, 351, 617, doi: 10.1086/168500 Byrne, C. M., & Jeffery, C. S. 2018, MNRAS, 481, 3810, doi: 10.1093/mnras/sty2545 Castellani, V., Chieffi, A., Tornambe, A., & Pulone, L. 1985, ApJ, 296, 204, doi: 10.1086/163437 Castellani, V., Giannone, P., & Renzini, A. 1971, Ap&SS, 10, 355 Castor, J., Davis, Jr., C. G., & Davison, D. K. 1977, , Tech. Rep. LA-6664, TRN: 77-008571, Los Alamos Scientific Lab., NM (USA) Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, 872, 51, doi: 10.3847/1538-4357/aaf99f Chandrasekhar, S. 1933, MNRAS, 93, 462,  \ndoi: 10.1093/mnras/93.6.462"}