Modules for Experiments in Stellar Astrophysics (MESA): Pulsating Variable Stars, Rotation, Convective Boundaries, and Energy Conservation

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.

One of the foundations upon which modern astrophysics rests is the fundamental properties of stars throughout their evolution. The advent of transformative 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).
The 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. 2018;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).
The 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 interferometers 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.
In 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.
This 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 implica-tions 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. 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 MESA model (at Z = 0.02) from ZAMS to a white dwarf (WD) is shown by the purple curve.
Millions 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;Freed-man et al. 2001;Majaess et al. 2009;Riess et al. 2016Riess et al. , 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 Large-Amplitude 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 10 9 stars. They vary in brightness by 20% on 30 min timescales (Pietrukowicz et al. 2013). An important new addition to MESA is the capability to model radially-pulsating variable stars.
Numerical 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 model from the pre-main sequence to the end of core He-burning or a core-collapse supernova from soon after explosion to shock breakout.
Rotation 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.
Stars 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.
The 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.
Important 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. Table 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.

Name
Description Appears   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. Early 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áth 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áth et al. 2002;Smolec & Moskalik 2010) and the cyclic modulations of RR Lyrae light curves (e.g., the Blazhko effect, Blažko 1907;Szabó et al. 2010). For background material we refer the reader to Gautschy & Saio (1995, Buchler (2009), and Marconi (2017).

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 Newton-Raphson 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.
RSP follows Smolec & Moskalik (2008), where the momentum and specific internal energy equations are where D/Dt is the Lagrangian time derivative. The generation of a specific turbulent energy, e t , is described by the one-equation Kuhfuß (1986) model The latter two equations are added to give an equation for the specific internal and turbulent energies Definitions 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). Several 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 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). These 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 0.01, α m 1, and α 2 are useful starting choices.
Periods 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 parameter surveys are an essential part of any science application of RSP.
In 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.
The 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).  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 Touter, 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 Nouter cells (RSP nz outer), each with a constant mass. Between the inner boundary and the anchor the mass of each cell increases.

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.

Building an Initial Model
Since 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.
Specifications 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 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,000 K 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.
The 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.

Stability Analysis
The LNA analysis is performed on the initial model using a full linearization of the RSP equations. 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. Figure 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 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. Figure 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.

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 of these three scaled eigenvectors is then multiplied by the surface velocity RSP kick vsurf km per sec.
The 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. Figure 4 shows Γ the fractional growth of the kinetic energy per pulsation period near the start of a time integration, where and 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.

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 Figure 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).
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.
For case (a), the pulsations converge towards a single, 1O-mode pulsation. After a 500 cycle transient phase the pulsation period and radius amplitude barely change and Γ 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 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. Figure 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. Table 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).

Control
Set A Set B Set C Set D

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,000 K 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. Figure 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. Figure 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 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 110, and is associated with the second ionization of He. In most of the models these two convective regions merge at pulsation phase 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 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.

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 the effect of these parameters on light curves for the T eff = 6,000 K 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. 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. 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.65 M , L = 50 L , T eff = 7,000 K, X = 0.75, Z = 0.0014) with convective parameter set A.

Spatial and Temporal Sensitivity
For 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.
The 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 stress 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.

Current Limitations and Plans for the Future
RSP in its present form covers most of the classical instability strip including δ Cepheids, RR Lyrae, High    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 (RR Lyrae) or 1 R (Cepheids).
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.
RSP 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 cannot 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).
Inclusion 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.

Applications of RSP
We now apply RSP to variable stars. Examples include the light curves of classical pulsators, modeling of specific objects, and models for the dynamics of modulated or chaotic pulsations.

RR Lyrae variables
We consider two sequences of RR Lyrae-type models. The first sequence has M = 0.65 M , L = 45 L 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.
To compare the overall morphology of I-band light curves from these sequences with OGLE observations, we perform a Fourier decomposition of the synthetic light curves where f is the pulsation frequency, and A k and φ k are amplitudes and phases, respectively. We then construct the amplitude ratios R k1 and epoch-independent phase differences ϕ k1 (Simon & Lee 1981): Observationally derived values of R k1 and ϕ k1 are taken from the OGLE catalog (Soszyński 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. The left panels show that F-mode pulsators with convective sets A and B progress similarly. Models with convective 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, Pulsation phase 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 than the observed values. However, the model physical parameters except T eff are fixed, while this is not the case for the OGLE RRab stars. The 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.
The 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.

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 5 d. The bump moves towards earlier phases on the descending branch as the period increases and is coincident with maximum light when P 10 d. The bump then moves onto the ascending branch for longer periods and disappears at P 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.
We consider 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, and 10.0 M 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 = 300 K or 500 K offset from the blue edge. Figure 14 shows I-band light curves and radial velocity curves for the ∆T eff = 300 K models. Due to the shift in the location of blue edge models, models from convec-tive 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ński 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 similar 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.
The Fourier phases in the left panels of Figure 15 are systematically larger than the observationally-inferred values for P 10 d, with the difference being larger for the cooler ∆T eff = 500 K models. For P 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 (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).
This 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.

Type II Cepheids
Type II Cepheids are more similar to RR Lyrae stars than to classical Cepheids due to their lower masses (M 0.5 M ), 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ński et al. 2019  tatively similar, devoid of fine-tuning, and demonstrate that RSP can be used to model specific stars. Figure 17 shows the radius variation for two models with M = 0.55 M , L = 136 L , 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,932 K (upper panel) and T eff = 6410 K (lower panel). Figure 18 shows the return map 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 . For the cooler T eff = 5,932 K model, the upper panel of Figure 17 shows a cyclic modulation in the envelope of the period-doubled pulsation. The modulation period of 57 d is longer than the pulsation period of 2.3 d. The return map in the upper panel of Figure 18 is constructed from 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 recently discovered in BL Her variables (Smolec et al. 2018).
For the hotter T eff = 6,410 K 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. 1998;Buchler et al. 2004;Plachy et al. 2018). While the T eff = 6,410 K model has a shorter period, in the BL Her range, such models may provide insight into chaotic dy- Observations are marked with gray dots. Light curve data are from Soszyński 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.

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 0.26 M .
Following Smolec et al. (2013), we adopt M = 0.26 M , L = 33 L , T eff = 6,910 K, X = 0.7, Z=0.01, and convective set A. Figure 19 shows I-band light curves and radial 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) Pulsation phase 16750 16800 10 11 12 30250 30300 30350 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 600 cells. The amplitude of the model curves are sensitive to convective parameters and can be fine-tuned for a better match with observations.

Blue Large Amplitude Pulsators
The origin of BLAPs, introduced in Section 1, is unknown. They have been modeled as 0.  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.
We explored nonlinear models with M = 1.0 M , T eff =30,000 K, L = 430 L , and three convection parameter sets. The LNA analyses show the F and 1O modes are unstable. Figure  light curves. The qualitative agreement is reasonable. Figure 20 shows the model radial velocity curves have amplitudes of 200 km s −1 , a pronounced temporal asymmetry, and light maxima that are narrower than light minima.
These 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. Figure 21 shows the evolution of Γ for the (α,α m )=(1.5,0.1) model in Figure 20. Before 10 d, Γ fluctuates around a mean 7.5×10 −3 . The fluctuations diminish with time but Γ remains above the LNA value of 6.25×10 −3 up to 30 d. Between 30-50 d Γ diminishes and approaches zero, 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.

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.
We consider a stellar model with M = 2 M , L = 30 L , T eff = 6,900 K, 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 the F and 1O modes are linearly unstable, with growth rates of 1×10 −6 and 6×10 −5 , respectively.
Section 3 emphasizes the importance of numerical energy conservation. Figure  relative cumulative error in the energy for a 250,000 cycle integration ( 150 million time steps, at 600 steps per cycle). The relative cumulative error grows from 3×10 −11 after about 15 cycles to 3×10 −7 after 250,000 cycles. The inset figure shows that the per-step relative error in the energy scatters around −2.5×10 −15 but is systematically different than zero. Figure 23 shows that after 70,000 cycles the asymmetric V -band light curves have an amplitude of 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 0.3 mag and an unchanged dominant pulsation period. Extending the integration to 250,000 cycles does not lead to any significant changes. The 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áth 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, F-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).
Evolution of fractional radius amplitude, δR/R, for the 1O mode, A1 and F mode, A0. 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.

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: This 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 ( ) and local fluxes between cells (the ∂/∂m term): The 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.
Recall (from Paper I, Section 6.3) the generalized Newton-Raphson scheme used by MESAstar to solve the stellar equations where y i is the trial solution for the i-th iteration, F ( y i ) is the residual, δ y i is the correction, and [d F /d 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.
If 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. The 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.

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 model during the main sequence and the He flash with gold tolerances and compare to the old approach. The 1 M 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.
Not 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 address 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).

Definition of nuc Source Term
Previous MESA papers have not given a precise definition of nuc . Motivated by numerical convenience, MESA formerly exploited the fact that the source term contained the sum of grav and nuc and included the response of the internal energy to composition changes due to nuclear reactions in nuc instead of in grav . However, since the dedt-form does not include grav , this is no longer an appropriate choice.
In the current approach, nuc is evaluated in the net module as a sum over reactions. Schematically, where 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 where M i is the rest mass of isotope i andẎ 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 nuc (e.g., Equation 11 in Hix & Meyer 2006). The MESA definition of nuc differs in subtracting off the nuclear neutrino losses, thus enforcing the assumption that they free-stream out of the star.

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 sur-face 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.
In addition to providing accurate accounting for the total energy of the model, it is important to ensure that 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.

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 Ṁ that accounts for the heating associated with mass changes.
The 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 and During 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 .
The change in energy for cell k during stage I is where 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 where m k,C and r k,C are the mass coordinate and radius at the center of mass of the cell (Paper IV). We now introduce 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 where 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 The energy flux is wherev k,I ≡ r mid,k − r start,k δt , and the face values E face are interpolated 5 from the cell values E start . The change in total energy of the model during stage I is The 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 k,I term implies the need for a corrective source term Ṁ which must be added during stage II so that energy is properly accounted for.
To determine this new source we consider the ratio of the thermal time scale to the mass-change time scale where 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 τ m ) the entropy of material adjusts minimally as it moves from cell to cell.
We 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 where is the total mass of material that at any point during stage I was inside cell k. 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 This 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 1).
We 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 where 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 This 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. Along the lines of Equation (17), the energy change of cell k during stage II may now be written in fluxconservative form as where F e,k,II = L k + p face,k A kvk,II andv k,II ≡ r k − r mid,k δt .
MESA 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Ṁ is large relative to the masses of cells in burning regions. WhenṀ ≥ 0 the above procedure for redistributing energy is conservative, so that k k,I δtdm k + k Ṁ ,k δtdm k = 0.
The first sum represents the effect of stage I; the second sum represents the effect of stage II. WhenṀ < 0 the equality is instead with 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.

Results
x m [g] To examine the behavior of Ṁ we modeled a 0.33 M WD accreting He and N at a rate of 10 −5 M yr −1 . The first panel of Figure 26 shows a profile of Ṁ along with the 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 τ m . The Ṁ term handles both limits. To demonstrate improved energy conservation during mass changes with the dedt-form we model accretion onto a 0.3 M He WD with an initial log(T c /K) = 6.7. The accretion rate is fixed at 10 −10 M 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 .
We prefer using the dedt-form because of its improved energy conservation and handling of long thermal times. We now explore the consequences for stellar evolution of these different prescriptions.   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.
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.
Mass 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 primary and a 6.5 M 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 this system computed using each energy equation. The differences are typically of order 1 %.   Figure 4 of Paper III: a binary system with mass transfer from the 8 M primary to the 6.5 M 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.  Figure 4 of Paper III: a binary system with mass transfer from the 8 M primary to the 6.5 M 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Ṁ computed with the dedt-and dLdm-forms. The spike at the end is due to mass transfer terminating at slightly different times.
For a rigidly-rotating star in hydrostatic equilibrium, surfaces of constant pressure (isobars) coincide with the equipotential surfaces defined by the Roche potential Ψ, where Φ 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, where r Ψ is the volume-equivalent radius of an isobar, and m Ψ and S Ψ are the mass inside that isobar and its surface 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 This 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).
Until 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).
As 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.
We 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.
(38) Figure 32 shows how the equipotential surfaces change with increasing ω. For ω 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.
By 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 = 2r 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.

Implementation in Stellar Evolution Instruments
To 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 As ω → 1, the moment of inertia becomes that of a ring with radius re. The lower panel shows the fP and fT factors as a function of ω = Ω/Ωcrit. j rot , m Ψ , and r Ψ . From Ω Ψ = ωΩ crit,Ψ = ω Gm Ψ /r 3 e and j rot = i rot Ω Ψ , we find For 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. Figure 34 shows a 3 M 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. Differences 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 maxi- mum 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.

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 a rotating star, and second, calculate the projection of the gravity-darkened surface along the LOS. We 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 radiusr = R/R e by solving We then solve cos ϑ + ln tan(ϑ/2) = 1 3 ω 2r3 cos 3 θ + cos θ + ln tan(θ/2), (41) for the modified angular variable ϑ. Using ELR Equation (31), we use this value to obtain the local T eff . Figure 35 shows the variation ofT eff = T eff [L/(4πσR 2 e )] −1/4 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.

Projection Effects and Correction Factors
We are interested in the projected -the directional average over the surface along the LOS-T 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 vectorl(i) and 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 ofT eff over the surface. To calculate the luminosity projected along the LOS, L proj , requires the surface integral where only the flux projected toward the observer, i.e., d Σ ·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 As 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 and C 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. Figure 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%. Figure 37 shows a slight, but noticeable, decrease in C T for ω 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. Figure 38 demonstrates the effect of gravity darkening on 3 of the 3 M 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.

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 "pre- dictive mixing" (Paper IV), and a new "convective premixing" scheme that addresses several remaining issues.
Early versions of MESA located convective boundaries by searching for sign changes in the discriminant y = ∇ rad − ∇ ad (or y = ∇ rad − ∇ L , when the Ledoux stability 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 discontinuities. Typically, 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.
In 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 star on the MS; a retreating convective core in a 16 M star on the MS; growing convective cores in 1 M and 3 M stars during the core Helium burning (CHeB) phase; and an evolving convective envelope in a 1 M star on the MS. Predictive mixing is able to achieve the desired ∇ rad = ∇ ad outcome in most of these cases.
However, 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).

The Failure of Predictive Mixing
The 16 M 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.
Additionally, the 1 M 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.
In 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.

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 convection 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.
If 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.
At 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.
In 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 10 4 yr for a typical region with extent ∆r R . 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.
During 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.

Evolution of a Retreating Convective Core on the Main Sequence
We evolve a 16 M 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 the CPM scheme; in both cases, the Ledoux stability criterion is adopted.
In 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 m/M 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 ≈ 7.2. Moreover, the region is very close to adiabatic (∇ rad = ∇ ad ) and thus convectively neutral according to the Schwarzschild stability criterion, in contrast to the super-adiabatic stratification (∇ rad > ∇ ad ) seen in the left panel. Schwarzschild & Härm (1958) were the first to predict the appearance of convectively neutral (∇ rad = ∇ ad ) abundance-gradient regions outside the convective cores of massive MS stars, labeling them 'semiconvection regions'. The key feature of a semiconvection region is that convective neutrality is continuously maintained by a gradual adjustment of opacity due to the changing abundance profile.
To illustrate how this adjustment naturally arises within the CPM scheme, we simulate a single evolutionary time step of the 16 M 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 abundance gradients within the shell itself, the Ledoux and Schwarzschild stability 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).
The seminal paper by Schwarzschild & Härm (1958) triggered significant interest in semiconvection regions, with a particular focus on the criterion adopted for convective neutrality. Sakashita & Hayashi (1961) argued that the Ledoux criterion ∇ rad = ∇ L should apply in semiconvection regions, rather than ∇ ad = ∇ rad as originally proposed. However, Kato (1966) reasoned that because the resulting regions are super-adiabatic (∇ rad > ∇ ad ), slow mixing by overstable g-mode oscillations will eventually drive them toward ∇ rad = ∇ ad . Subsequently, Gabriel (1970) suggested that Kato's mechanism is unnecessary, due to the appearance of propagat- ing convective shells that continually adjust the abundance profile to achieve ∇ rad = ∇ ad . Gabriel's narrative closely mirrors the one we give above, and the correspondence between his Figure 1 and our Figure 40 is striking.
A 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 refers to the slow mixing in super-adiabatic regions that was proposed by Kato (1966); while this mixing will yield a region satisfying ∇ rad = ∇ ad , if given enough time, it is fundamentally a different mechanism. A critical distinction lies in the role played by the convective stability criterion. While our calculations adopt the Ledoux criterion, repeating them with the Schwarzschild criterion 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. In hindsight, it seems prudent to reserve the label 'semiconvection' for the adiabatically neutral regions envisaged by Schwarzschild & Härm (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).
To summarize the core and near-core evolution of the 16 M 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 & Härm (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.
As 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-Väisälä 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 = 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 predicted 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 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.

Evolution of the Convective Core During Core He Burning
We now evolve a 1 M 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 cases, the Ledoux stability criterion is adopted.
The 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.
The emergence of semiconvection regions during the CHeB phase may have first been proposed by Schwarzschild & Härm (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.
To summarize the core and near-core evolution of the 1 M star, Figure 44 plots the mass coordinate m c of the convective core boundary as a function of CHeB age, for the 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.
Debate 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 0.15) of the CHeB phase.

Evolution of a Growing Convective Core on the Main Sequence
We now evolve a 1.5 M 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 cases, the Ledoux stability criterion is adopted.
The 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 µ ∝ m p 7/5 he obtains for the abundance-gradient region outside the core yields an adiabatic stratification (∇ rad = ∇ ad ).
To 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 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.

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.

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.
Three 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.
As 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, is 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. Figure 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.
For 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.

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. We 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 of 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. 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 soft-ware 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. The 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 (ρ 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.  (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 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 WD.
In 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:  Figure 51. Adiabatic gradient in the ρ − T plane for X = 0.70 and Z = 0.02. Regions undergoing H2, H, He, or He + dissociation/ionization are colored blue and labelled. Previous versions of eos truncated the H2 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. ideal gas Figure 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 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.
level 5: If f OPAL/SCVH =0, then return result of level 6. If f OPAL/SCVH =1, then evaluate OPAL/SCVH and return to level 4. Otherwise, call level 6, blend the result with that of OPAL/SCVH, and return the blended result to level 4.
level 6: Evaluate HELM and return result to level 5.
The 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.
There 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 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 p = ρ 2 ∂e ∂ρ Ideally, 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,Yi ) 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 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.
A 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.
MESA 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 10 9 K), but less accurate at low temperatures (T 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 where 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 Itoh et al. (1979), where the Z i are replaced with the average chargeZ. Figure 54 compares three MESA nuclear reaction rate screening options on the evolution of a 25 M 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 and the potential as Ψ = Ψ/(Gm Ψ Ω) 2/3 .

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 dV dx = 2πx y , y = 1 ω 2/3 + where 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), In the opposite case where ω → 0, the equipotential is well approximated by an ellipsoid of revolution with r p (ω) = r e 1 − ω 2 2 + O(ω 4 ) , such that its volume is By 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 ω given by Equation (B7): This 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. Stellar 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, r Ψ = r e 1 − ω 2 6 + O(ω 4 ) .

B.3. Surface Averages of Gravity
The surface average of the dimensionless gravity g is computed as with 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 g (ω) = 1 r 2 e 1 + 1 − 2 x r e 2 ω 2 , which combined with Equations (B13) and (B16) results in By numerically computing this integral, a simple fit that matches the computed value at ω = 1 and has an error < 0.35% for 0 ≤ ω ≤ 1 is S Ψ g = 4πGm Ψ 1 − 2 3 ω 2 + 0.3045ω 4 + 0.001382ω 6 .
Similarly the average of the inverse gravity, in the limit of ω → 0, can be shown to be However, 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: Gm Ψ A(ω), A(ω) = 1 − 0.1076ω 4 − 0.2336ω 6 − 0.5583 ln(1 − ω 4 ).

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 ds = |∇Ψ| −1 dΨ = g −1 dΨ.
Assuming a constant density ρ in the shell (as in the shellular approximation), dm = 2ρ dΨ 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: If 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.
Next, 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.
We explode (with E exp = 10 51 erg) the 99em 19 RSG progenitor model from Paper IV, which has a ZAMS mass of 19 M , and a mass of 17.8 M and a radius of 603 R at time of explosion. We excise the inner 1.5 M and explode the remaining 16.3 M of ejecta using a thermal bomb as described in Paper IV. The mass of radioactive 56 Ni is set to M Ni = 0.03 M . 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.24 cm 2 g −1 for the metal-rich core material, κ floor, envelope = 0.01 cm 2 g −1 for the envelope, and proportional to metallicity for the intermediate region. Figure 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. Figure 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  , 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. 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.
To 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 a 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.