We update the capabilities of the openknowledge software instrument Modules for Experiments in Stellar Astrophysics MESA. The new autodiff module implements automatic differentiation in MESA an enabling capability that alleviates the need for hardcoded analytic expressions or finitedifference approximations. We significantly enhance the treatment of the growth and decay of convection in MESA with a new model for timedependent convection which is particularly important during latestage nuclear burning in massive stars and electrondegenerate ignition events. We strengthen MESAs implementation of the equation of state and we quantify continued improvements to energy accounting and solver accuracy through a discussion of different energy equation features and enhancements. To improve the modeling of stars in MESA we describe key updates to the treatment of stellar atmospheres molecular opacities Compton opacities conductive opacities element diffusion coefficients and nuclear reaction rates. We introduce treatments of starspots an important consideration for lowmass stars and modifications for superadiabatic convection in radiationdominated regions. We describe new approaches for increasing the efficiency of calculating monochromatic opacities and radiative levitation and for increasing the efficiency of evolving the late stages of massive stars with a new operatorsplit nuclear burning mode. We close by discussing major updates to MESAs software infrastructure that enhance source code development and community engagement.
['Stellar physics', 'Stellar evolution', 'Stellar evolutionary models', 'Computational methods', '1621', '1599', '2046', '1965', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Modules for Experiments in Stellar Astrophysics MESA Timedependent Convection Energy Conservation Automatic Differentiation and Infrastructure
{'Modules for Experiments in Stellar Astrophysics (MESA): Time-Dependent Convection, Energy Conservation, Automatic Differentiation, and Infrastructure': "Adam S. Jermyn , 1 Evan B. Bauer , 2 Josiah Schwab , 3 R. Farmer , 4 Warrick H. Ball , 5 Earl P. Bellinger , 4, 6 Aaron Dotter , 7 Meridith Joyce , 8 Pablo Marchant , 9 Joey S. G. Mombarg , 10 William M. Wolf , 11 Tin Long Sunny Wong , 12 Giulia C. Cinquegrana , 13, 14 Eoin Farrell , 15 R. Smolec , 16 Anne Thoul , 17 Matteo Cantiello , 18, 19 Falk Herwig , 20 Odette Toloza , 21, 22 Lars Bildsten , 23, 12 Richard H. D. Townsend , 24 and F.X. Timmes 25  \n1 Center for Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA  \n2 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA  \n3 Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA  \n4 Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany  \n5 School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK  \n6 Stellar Astrophysics Centre, Aarhus University, 8000 Aarhus C, Denmark  \n7 Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755 USA  \n8 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  \n9 Institute of Astrophysics, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium  \n10 Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, Leuven, Belgium  \n11 Department of Physics and Astronomy, University of Wisconsin-Eau Claire, Eau Claire, WI 54701, USA  \n12 Department of Physics, University of California, Santa Barbara, CA 93106, USA  \n13 School of Physics & Astronomy, Monash University, Clayton VIC 3800, Australia  \n14 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)  \n15 School of Physics, Trinity College Dublin, The University of Dublin, Dublin 2, Ireland  \n16  \nNicolaus Copernicus Astronomical Center of the Polish Academy of Sciences, Bartycka 18, PL-00-716 Warszawa, Poland 17 Space sciences, Technologies and Astrophysics Research (STAR) Institute, Universit'e de Li'ege, All'ee du 6 Ao ˆ u t 19C, Bat. B5C, 4000 Li'ege, Belgium  \n18 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA  \n19 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA  \n20 Dept. of Physics and Astronomy, University of Victoria, Victoria, BC V8P5C2, Canada  \n21 Millennium Nucleus for Planet Formation, NPF, Valpara'ıso, Chile  \n22 Departamento de F'ısica, Universidad T'ecnica Federico Santa Mar'ıa, Avenida Espa˜na 1680, Valpara'ıso, Chile  \n23 Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA  \n24 Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA  \n25 School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA", 'ABSTRACT': "We update the capabilities of the open-knowledge software instrument Modules for Experiments in Stellar Astrophysics ( MESA ). The new auto diff module implements automatic differentiation in MESA , an enabling capability that alleviates the need for hard-coded analytic expressions or finite difference approximations. We significantly enhance the treatment of the growth and decay of convection in MESA with a new model for time-dependent convection, which is particularly important during latestage nuclear burning in massive stars and electron degenerate ignition events. We strengthen MESA 's implementation of the equation of state, and we quantify continued improvements to energy accounting and solver accuracy through a discussion of different energy equation features and enhancements. To improve the modeling of stars in MESA we describe key updates to the treatment of stellar atmospheres, molecular opacities, Compton opacities, conductive opacities, element diffusion coefficients, and nuclear reaction rates. We introduce treatments of starspots, an important consideration for low-mass stars, and modifications for superadiabatic convection in radiation-dominated regions. We describe new approaches for increasing the efficiency of calculating monochromatic opacities and radiative levitation, and for increasing the efficiency of evolving the late stages of massive stars with a new operator split  \nnuclear burning mode. We close by discussing major updates to MESA 's software infrastructure that enhance source code development and community engagement.  \nKeywords: Stellar physics (1621); Stellar evolution (1599); Stellar evolutionary models (2046); Computational methods (1965)", '1. INTRODUCTION': "A resurgence of stellar astrophysics research is being fueled by the transformative capabilities in spaceand ground-based hardware instruments providing an unprecedented volume of high-quality measurements of stars, significantly strengthening and extending the observational data upon which stellar astrophysics ultimately rests (National Research Council 2021). Examples include:  \nSeveral individual stars at redshifts of glyph[similarequal] 1 have been discovered by temporary magnification factors of glyph[similarequal] 1000 from microlensing (Kelly et al. 2018; Rodney et al. 2018; Chen et al. 2019). A more persistent and highly magnified star at a redshift of 6.2 has also been discovered with the Hubble Space Telescope (Welch et al. 2022a) by a fortuitous alignment with a foreground galaxy cluster lens caustic (Windhorst et al. 2018). The infrared instruments aboard the James Webb Space Telescope (Gardner et al. 2006; Beichman et al. 2012; Artigau et al.  \n2014; Rieke et al. 2015; Labiano et al. 2021) will search for confirmation and spectral classification of this distant star (Welch et al. 2022b) to define its place on the Hertzsprung-Russell Diagram (HRD), assess how galaxies evolve from their formation (Zackrisson et al. 2011; Robertson 2021), observe the formation of stars (Senarath et al. 2018; Boquien & Salim 2021), and measure the properties of stellar-planetary systems including the Solar System (Sarkar & Madhusudhan 2021; Patapis et al. 2022).  \nIn the late 2020s, kilometer-scale gravitational wave detectors such as Advanced Laser Interferometer Gravitational Observatory (LIGO Scientific Collaboration et al. 2015), Advanced Virgo (Acernese et al. 2015) and Kamioka Gravitational Wave Detector (Akutsu et al. 2021) will routinely detect tens of binary neutron-star mergers with kilonovae annually (Abbott et al. 2018), probe how kilonova r-process nucleosynthetic yields vary with environment (e.g., Barnes et al. 2021), and assess the populations that contribute to the stellar black hole mass distribution, including the presence of any gaps in the distribution (Perna et al. 2019; Zevin et al. 2021; Renzo et al. 2020, 2021; Mandel & Broekgaarden 2022).  \nThe next core-collapse supernova in the Milky Way or its satellites will be a unique opportunity to observe the explosion of a star. The SuperNova Early Warning System is a global network of neutrino experiments sensitive to supernova neutrinos (Al Kharusi et al. 2021) that includes multi-kiloton detectors such as KamLAND (Araki et al. 2005), Borexino (Borexino Collaboration et al. 2018, 2020), SNO+ (Andringa et al. 2016), Daya Bay (Guo et al. 2007), SuperKamiokande (Simpson et al. 2019), and the upcoming HyperKamiokande (Abe et al. 2016), DUNE (Acciarri et al. 2016) and JUNO (JUNO Collaboration 2022). Searching for pre-supernova neutrinos is ongoing, and of interest as they allow tests of stellar and neutrino physics (e.g., Kosmas et al. 2022) and enable an early alert of an impending core-collapse supernova to the electromagnetic and gravitational wave communities (Beacom & Vogel 1999; Vogel & Beacom 1999; Mukhopadhyay et al. 2020; Al Kharusi et al. 2021).  \nSky surveys that probe ever-larger areas of the dynamic sky and ever-fainter transient sources include the Imaging X-ray Polarimetry Explorer (Soffitta et al. 2021), the Compton Spectrometer and Imager (Tomsick & COSI Collaboration 2022), eROSITA (Predehl et al. 2021), Gaia (Gaia Collaboration et al. 2016, 2018, 2021), the Sloan Digital Sky Survey (York et al. 2000; Abdurro'uf et al. 2022), the All-Sky Automated Survey for Supernovae (Chen et al. 2022), Pan-STARRS1 (Flewelling et al. 2020), the Zwicky Transient Factory (Bellm et al. 2019; Dhawan et al. 2022), Gattini-IR  \n(Moore et al. 2016), and the Nancy Grace Roman Space Telescope (Akeson et al. 2019). Roman will measure proper motions of stars several magnitudes fainter than the Gaia mission (Brandt et al. 2021; Dorn-Wallenstein et al. 2021), which is sufficient to probe the main sequence turnoff to distances of glyph[similarequal] 10 kpc and red giants throughout the Galactic halo (Spergel et al. 2015).  \nWide-field spectroscopic surveys in the coming decade will resolve stellar populations and the Milky Way's structure (Bolton et al. 2019) at facilities such as Gaia DR3 (Gaia Collaboration et al. 2021), SDSS-V (Kollmeier et al. 2017), FOBOS (Bundy et al. 2019), Maunakea Spectroscopic Explorer (Marshall et al. 2019), and SpecTel (Ellis & Dawson 2019). For example, FOBOS is a next-generation spectroscopic facility at the W.M. Keck Observatory that will provide multiepoch, high-multiplex, and deep spectroscopic follow-up of panoramic deep-imaging surveys.  \nThe Vera C. Rubin Observatory will conduct a multicolor optical survey of the Southern Hemisphere sky, the Legacy Survey of Space and Time (LSST Science Collaboration et al. 2017; Ivezi'c et al. 2019), to probe dark energy and dark matter (S'anchez et al. 2021; Zhang et al. 2022), explore the transient optical sky (Bianco et al. 2022; Li et al. 2022; Raiteri et al. 2022; Hernitschek & Stassun 2022; Andreoni et al. 2022; Bellm et al. 2022), and build a catalog of solar system objects with an order of magnitude more objects (LSST Solar System Science Collaboration et al. 2020; Schwamb et al. 2021).  \nThe TESS mission (Ricker et al. 2016) is providing systematic measurements of the radii, masses, and ages of 200,000 individual stars sampled at a 2 minute cadence to open a new era of stellar variability exploration (e.g., Huang et al. 2018; Ball et al. 2018; Dragomir et al. 2019; Wang et al. 2019). Within the next decade, the Planetary Transits and Oscillations of Stars mission (PLATO; Rauer et al. 2014) will search for planetary transits across up to one million stars, characterize rocky extrasolar planets around yellow dwarf stars, subgiant stars, and red dwarf stars (Montalto et al. 2021), and investigate the seismology of stars (Miglio et al. 2017; Cunha et al. 2021; Nascimbeni et al. 2022).  \nIn partnership with this ongoing explosion of activity in stellar astrophysics, revolutionary advances in software infrastructure, computer processing power, data storage capability, and open-knowledge software instruments are transforming how stellar theory, modeling, and simulations interact with experiments and observations. Examples include Astropy (Astropy Collaboration et al. 2018, 2022), Athena++ (Stone et al. 2020; Jiang 2021), Castro (Almgren et al. 2020), Dedalus (Burns et al. 2020), emcee (Foreman-Mackey et al.  \n2013), Flash-X (Dubey et al. 2022), GYRE (Townsend & Teitler 2013; Townsend et al. 2018), MAESTROeX (Fan et al. 2019), MESA2Hydro (Joyce et al. 2019), MSG (Townsend & Lopez 2022), Phantom (Price et al. 2018), PHOEBE (Conroy et al. 2020), Starlib (Sallaska et al. 2013), TARDIS (Vogl et al. 2019), TULIPS (Laplace 2022), and yt (Turk et al. 2011).  \nThe previous Modules for Experiments in Stellar Astrophysics software instrument papers (Paxton et al. 2011, MESA I; Paxton et al. 2013, MESA II; Paxton et al. 2015, MESA III; Paxton et al. 2018, MESA IV; Paxton et al. 2019, MESA V), as well as this one, describe new capabilities and limitations of MESA while also comparing to other available numerical or analytic results. We do not fully explore the science implications in this software instrument paper. The scientific potential of these new capabilities will be unlocked by future efforts of the MESA research community.  \nThis MESA VI software instrument paper is organized as follows. Section 2 describes the implementation of automatic differentiation and § 3 introduces a new model for time dependent convection. Section 4 describes improvements to MESA 's implementation of the equation of state (EOS) and § 5 discusses treatments of the energy equation. Section 6 describes treatments of the stellar atmosphere and § 7 introduces new models of starspots and a superadiabatic convection. Section 8 reports improvements to the opacities, § 9 to the element diffusion coefficients, § 10 to the nuclear physics, and § 11 to the physical constants. Section 12 discusses MESA 's infrastructure. Finally, § 13 summarizes MESA VI.  \nImportant symbols are defined in Table 1. Acronyms are defined in Table 2. Components of MESA , such as modules and routines, are in typewriter font e.g., tdc .  \nTable 1 . Important symbols.  \nTable 1 (continued)  \nTable 1 continued  \nTable 1 (continued)  \nNote: Single character symbols are listed first, symbols with modifiers are listed second. Some symbols may be further subscripted, for example, by c (indicating a central quantity), by s (indicating a surface quantity), by a cell index k , or by species index i .  \nTable 2 . Acronyms used in this article.  \nTable 2 continued  \nTable 2 (continued)", '2. AUTOMATIC DIFFERENTIATION': 'MESA solves the equations of stellar evolution implicitly using a Newton-Raphson method, which requires the partial derivatives of each equation with respect to the basic structure variables in each cell (e.g., ρ k , T k ). These derivatives need to be computed accurately, typically to one part in 10 6 , often precluding use of finite differences. These derivatives have historically been computed by hard-coding analytic expressions for each equation. This has accounted for a large fraction of the complexity and sources of error in MESA .  \nWe have largely eliminated this source of error and the associated complexity by using forward-mode operatoroverloaded automatic differentiation (BartholomewBiggs et al. 2000) in the new auto diff module. This functionality provides partial derivatives of expressions automatically with respect to their input variables. The auto diff module provides a number of Fortran derived types for this purpose. For example, we define the type auto diff real star order1 , which contains a floating-point number as well as its first partial derivative with respect to the basic stellar structure variables. The number of partial derivatives is specified at compile-time. If x is a variable of this type, then it contains components x%val representing the value of x and x%d1Array(j) for the value of ∂x/∂η j , where η j is the j -th independent variable.  \nThe auto diff types overload operators to implement the chain rule. This means that a source code line such as f = x * y is equivalent to  \nf%val = x%val * y%val f%d1Array(j) = x%d1Array(j) * y%val + y%d1Array(j) * x%val  \nBasic arithmetic and all special and trigonometric functions used in MESA including functions such as min , max , and abs , are provided. When these functions have discontinuities, we evaluate their derivatives as zero; and where they have discontinuous derivatives, we compute  \ntheir derivatives as the average between the two sides of the discontinuity.  \nUsing auto diff , expressions like  \nF = min ( ρ 1 e T 1 / √ ρ 0 , cosh ( r 2 -r 1 ) ) (1)  \ncan be written as', 'F = min(rho1*exp(T1/sqrt(rho0)),cosh(r2-r1)).': "Together with setup routines that link physical variables (e.g., T1 ) with the independent variables η j , this code automatically provides correct partial derivatives of F .  \nBy contrast, explicitly obtaining the partial derivatives of F requires more complex and error-prone source code:  \nx0 = rho1*exp(T1/sqrt(rho0)) x1 = cosh(r2-r1) F = min(x0,x1) if (x0 < x1) then dF\\_drho1 = x0 / rho1 dF\\_drho0 = -T1 * x0 / (2 * sqrt(rho0)) dF\\_dT1 = x0/sqrt(rho0) dF\\_dr2 = 0 dF\\_dr1 = 0 else dF\\_drho1 = 0 dF\\_drho0 = 0 dF\\_T1 = 0 dF\\_dr2 = sinh(r2-r1) dF\\_dr1 = -sinh(r2-r1) end if  \nThe auto diff module provides overloaded operators that were generated using the SymPy (Meurer et al. 2017) library in Python to compute power series and extract chain-rule expressions. We first optimized these expressions to eliminate common sub-expressions and minimize the number of division operators. We then translated these into Fortran. This functionality is built on top of the CR-LIBM software package (Daramy-Loirat et al. 2006), which enables bit-for-bit identical results across all platforms (see MESA III).  \nThe auto diff module also provides additional auto diff real types for alternative use cases. For convenience, types are provided to support the different hooks in MESA . For operations requiring higher-order derivatives, such as in the EOS (see § 4), additional auto diff real types provide higher-order mixed partial derivatives. The chain-rule expressions rapidly become more complicated for higher-order derivatives, but the basic principle is the same. The auto diff machinery was used to benchmark the Skye EOS (Jermyn et al.  \n2021), with the result that the performance is similar to explicit expressions. Here we provide more detailed benchmarks.  \nTable 3. Ratio of runtimes for evaluating auto diff expressions relative to the same expressions with real(dp) (double precision real ) variables. The real(dp) operations did not calculate any derivatives. Runtimes are averaged over 10 6 trials, performed on the integers 1 -10 6 (cast as real(dp) ), with intermediate results accumulated to prevent the compiler from optimizing away the operations. Label ' * ' refers to the multiplication operator x*x , labels 'Xvar orderY' refer to the number of independent variables X and the maximum partial derivative order Y, label 'f' refers to the function log(cosh(tanh))), and the label 'N' refers to the number of partial derivatives computed. Timing data was obtained on a 2.4 GHz 8-Core Intel Core i9 running on a 2019 Macbook Pro.Table 4. Same as Table 3, but comparing auto diff performance against explicit real(dp) routines that include partial derivatives. Partial derivative expressions were constructed and simplified using Mathematica version 12, then implemented manually in Fortran.  \nTable 3 compares the runtime cost for several operations and several auto diff types to the cost of evaluating the same expressions in real(dp) types calculating no partial derivatives. For operations like multiplication and division, this incurs an overhead of order the number of partial derivatives returned. For more expensive  \noperations, the overhead is much less, as the auto diff expressions are optimized to re-use intermediate results.  \nTable 4 compares three first-order auto diff types and explicit real(dp) routines evaluating the same partial derivatives. There is still overhead for simple operations, but the relative cost no longer scales with the number of derivatives. For sufficiently complex operations, such as f = log(cosh(tanh( x ))), the optimized auto diff functions outperform our explicit routines.  \nFor use in stellar evolution calculations, we find the runtime performance of hand-coded expressions are modestly better than those from auto diff , because most equations do not depend on all of the independent variables. However, this overhead is small compared with the full cost of a timestep in MESA . Moreover, runtime is often significantly reduced by ensuring that all partial derivatives are correct, as inaccurate derivatives result in slow convergence and a larger number of small timesteps. Some parts of the MESA source code do not yet use auto diff , but this is gradually transitioning.  \nFour applications of the profound enabling capability of auto diff are shown in § 3 on time-dependent convection, § 4 on MESA 's implementation of the EOS, § 7.1 on starspots, and § 7.2 on superadiabatic convection. The auto diff module can also be used in run star extras , as well as for software development outside of MESA .", '3. TIME-DEPENDENT LOCAL CONVECTION': 'The mixing length theory (MLT; Biermann 1932; Vitense 1953) has been used to parameterize convection in 1D stellar models for decades. It assumes that convective turbulence is in a steady state in which the energy input by the convective instability balances damping due to turbulent processes and radiative diffusion. This is a good approximation when the composition and structure evolve on time-scales that are long compared to the characteristic time-scales of convection.  \nHowever, during particularly violent episodes of stellar evolution, it is possible for the structure to evolve faster than convection can reach a steady state. This is the case in late-stage nuclear burning in massive stars (preceding core collapse), as well as during electrondegenerate ignition events (e.g., He shell flashes, Ne ignition, etc.). In such cases, the dynamics of convective growth and decay must be incorporated.  \nTo model this, we employ the time-dependent convection (TDC) formalism of Kuhfuß (1986) in the local limit. We build upon the implementation in Smolec & Moskalik (2008), introduced in MESA V to model radial stellar pulsations in the RSP module. More precisely, we use the one-equation version of the Kuhfuß (1986)  \nmodel, both in the RSP module and now for general use in stellar evolution calculations. We caution that combining different mixing models in a stellar evolution calculation might lead to physically inconsistent solutions, because the different models have been developed separately and their underlying assumptions might not be compatible with each other. Examples include combining the newly implemented time-dependent local limit convection model with an overshooting model, or combining TDC with other models for chemical composition gradients, rotation, etc.  \nWe describe the TDC formalism in § 3.1. In § 3.2 and § 3.3 we explain the modifications we have made relative to the implementation in RSP to make TDC numerically stable on long timescales. Section 3.4 then details the TDC solver algorithm. In § 3.5 we identify a change to the implementation which makes TDC agree with MLT in the limit of long timescales. Finally, in § 3.6 we explore the impact of TDC on models of white dwarfs (WDs) accreting He.', '3.1. Formalism': 'Following the Kuhfuß (1986) model, TDC introduces a new variable, the specific kinetic energy in turbulence e t , which evolves according to  \nD e t D t + α P t P t D ρ -1 D t = glyph[epsilon1] q + C -∂L t ∂m . (2)  \nHere P t ≡ (2 / 3) ρe t is the turbulent pressure, α P t is a dimensionless free parameter, glyph[epsilon1] q is the viscous heating of bulk motion, and L t accounts for advection of kinetic energy between mass shells  \nL t = -A αα t ρhe 1 / 2 t ∂e t ∂r . (3)  \nThis expression is the same as the turbulent flux F t in MESAV, but multiplied by A to convert to a luminosity, α t is a convective flux parameter, and h ≡ P/ ( ρg ) is the pressure scale-height. Furthermore,  \nC ≡ αe 1 / 2 t TPQ h √ 6 Y α D ( 8 3 √ 2 3 ) e 3 / 2 t αh -48 σα r α 2 ( T 3 ρ 2 c P κh 2 ) e t (4)  \ngroups together sources and sinks of turbulent kinetic energy, including a source/sink from the superadiabaticity  \nY ≡ ∇ - ∇ ad . (5)  \nThe coefficients α , α D , and α r are free parameters, and Q ≡ ∂ρ -1 /∂T | P is the thermal expansion coefficient.  \nBy default α r = 0, which means that TDC neglects radiative damping of convective motions. We caution that using α r = 0 is an approximation that changes the physical contents and the physical completeness of the model (Kuhfuß 1986, 1987; Wuchterl & Feuchtinger 1998). This choice enables subsequent modifications ( § 3.5) that make TDC reduce to MLT in the limit of long timesteps. Other defaults are α =2, α D =1, and α P t =0. The choice α D =1 is equivalent to the Kuhfuß (1986) choice of C D = (8 / 3) √ 2 / 3 for compatibility with MLT, and α P t =0 implies P t =0; see Table 3 of MESA V and Wuchterl & Feuchtinger (1998).  \nThe turbulent energy is incorporated into the other equations of stellar structure via heat and momentum transport. Specifically, in the momentum equation we include a turbulent pressure term P t . In the luminosity equation we incorporate  \nL = L rad + L conv , (6)  \nwhere L rad is the radiative luminosity and  \nL conv = 4 πr 2 αρc P T w √ 6 Y (7)  \nis the convective luminosity. Here w ≡ √ e t is the turbulent velocity, and the factor of √ 6 arises from a choice of closure constants. Note L t is set to zero in Equation 6, that is, the local limit solution is assumed. Finally, the luminosity enters the total energy equation, which sets the time evolution of the specific internal energy e in each cell.  \nTo implement TDC in MESA , we drop the term glyph[epsilon1] q from the energy equation, simplifying our implementation. We do not expect this term to matter in most cases, because bulk velocities are typically much smaller than convective velocities. When using the Ledoux criterion for convective stability, we further modify TDC in MESA relative to RSP to set Y with the Ledoux gradient ∇ L rather than the adiabatic gradient ∇ ad , as in Kuhfuß (1986).', '3.2. Numerical Stability': 'In MESA , RSP solves equations (2), (6) and (7) implicitly along with other structure equations to evolve e t and e . This approach works well on short (convective/pulsational) timescales, but it is numerically unstable on long (evolutionary) timescales. This poses a challenge, as we want a method that can be used in both limits and smoothly transitions between them.  \nWe conjecture that this numerical instability arises when e t is a solver variable. The superadiabaticity Y sets the time evolution of e t , and hence L conv and L . As L is very sensitive to Y , small errors in Y result in  \nlarge errors in L . These errors are not important over time steps shorter than the thermal timescale of a cell (as is the case in RSP ), as an excess luminosity through one face heats one adjacent cell and cools the other, restoring thermal equilibrium. With much longer time steps, errors in L significantly alter the entropy profile, propagating into Y and producing even larger errors in L with each iteration.  \nAn alternative approach, taken by the MLT implementation in MESA , is to treat the luminosity as a solver variable determined implicitly by the energy equation (e.g., Equation 49). From this, MESA derives the temperature gradient needed to produce that luminosity, and requires that the temperature gradient between cells match that computed by MLT. In effect, this flips the logic around, so that MESA MLT solves for Y given L whereas RSP solves for L given Y . Because L is very sensitive to Y , Y is relatively insensitive to L ; thus, errors in L produce smaller errors in Y , making this approach numerically stable.', '3.3. Flipped Equations': 'To ensure numerical stability over long timescales, we implement TDC in MESA in the same way as MLT, with L as a solver variable. We flip Equation (7) to solve for Y , accounting for the fact that w depends on Y via Equation (2). Doing so requires a few simplifications and a number of new approaches.  \nWe numerically invert Equations (2) and (7) to solve for Y given L . To do this, we note that the time evolution of e t in a single cell is nearly independent of e t in adjacent cells (see Equation 2). The only direct (rather than implicit) coupling between e t,k and e t,k ± 1 arises through L t . For simplicity, we currently set L t = 0. This makes D e t,k / D t independent of e t,k ± 1 except implicitly via the other structure variables. This in turn makes Y k a function only of the local luminosity L k and solver variables in the adjacent cells. We then solve for each Y k using only local information, and preserve the basic structure of the Jacobian in MESA as well as the runtime performance.', '3.4. Numerical Method': 'Our goal is to numerically solve Equations (2) and (7) for Y given L with L t = 0. We first construct machinery to evaluate L given Y , and then perform a numerical root-find to obtain Y given L .  \n3.4.1. L given Y  \nWe use e = w 2 to rewrite Equation (2) as  \n2 D w D t = ξ 0 + ξ 1 w + ξ 2 w 2 ,  \nt (8)  \nwhere ξ 0 , ξ 1 , and ξ 2 are coefficients that we obtain by expanding the definitions of C and P t in Equation (2). We have divided through by w , and so have implicitly excluded one solution branch ( w = 0). We will return to that branch shortly.  \nThe coefficients are given by  \nξ 0 = α h √ 6 c P T ∇ ad Y , (9)  \nξ 2 = -( 8 3 √ 2 3 ) α D αh . (11)  \nξ 1 = -  4 σT 3 ρ 2 c P κ ( 2 α r √ 3 αh ) 2 + 2 3 α P t ρ D ρ -1 D t   , (10)  \nWe fix these to their end-of-step values and solve for w at the end of the time step given the initial value at the start of the time step. This implicit approach is numerically stable, and the required end-of-step values are readily available. The form of the solution to Equation (8) depends on the sign of the discriminant J 2 ≡ ξ 2 1 -4 ξ 0 ξ 2 .  \nWhen J 2 > 0, the system is convectively unstable, with the solution  \nw = -1 2 ξ 2 ( J tanh λ + Jδt 4 + ξ 1 ) , (12)  \nafter a time step δt , where λ is a constant depending on the initial value of w . With long time steps, the solution grows to a plateau w → -( J + ξ 1 ) / (2 ξ 2 ) independent of this initial condition. We show below that this is consistent with MLT.  \nWhen J 2 < 0, the system is convectively stable, with the solution  \nw = 1 2 ξ 2 ( | J | tan λ + | J | δt 4 -ξ 1 ) . (13)  \nThis solution eventually reaches w ( δt ∗ ) = 0 at some time δt ∗ ∼ 1 /J . Beyond that point the system remains fixed at w = 0, which is a valid solution to Equation (2) but which was excluded in the form Equation (8) by dividing through by w . When J 2 < 0 we must additionally check to see if the first root of w occurs before the end of the time step and, if it occurs before, set w = 0 at the end of the step rather than evaluating w with Equation (13).  \nGiven w , we evaluate L conv at the end of the step via Equation (7), and so now have L given Y as desired.', '3.4.2. Numerical Inversion': 'Wenow invert the relation between L and Y by solving  \nL ( Y , ℵ ) = L solver , (14)  \nwhere L solver is the desired luminosity produced by the Newton-Raphson solver, L ( Y , ℵ ) is the relation we constructed using Equations (7) and (8), and ℵ represents additional structure variables.  \nTo solve equation (14), we write it in the form  \nR ( Y ) ≡ L ( Y , ℵ ) -L solver = 0 , (15)  \nand then expand L using Equations (6) and (7) as  \nR ( Y ) ≡ ( L solver -L 0 ∇ L ) -( L 0 + c 0 w ) Y . (16)  \nHere  \nL 0 ≡ 16 πac 3 ( GmT 4 κP ) c 0 ≡ 4 πr 2 α √ 6 ρTc P (17)  \nare positive quantities set by local properties and independent of Y . All quantities in Equation (16) are evaluated at the end of the time step, determining the sign of the solution for Y in advance. The factor L 0 + c 0 w is positive; hence, the sign of Y matches that of the first term, which is independent of Y .  \nWe evaluate L solver -L 0 ∇ L to determine the sign of Y , followed by a change of variables from Y → Z ≡ ln |Y| . This allows more resolution in Y , which can vary by many orders of magnitude across a stellar model. We restrict our search to -100 ≤ Z ≤ 100, covering 10 -43 glyph[lessorsimilar] |Y| glyph[lessorsimilar] 10 43 . We choose such a wide range because we have observed models that enter the extremes of this range, typically involving shocks where both TDC and MLT are suspect. We have not encountered models with |Y| approaching 10 43 , and those with |Y| < 10 -43 are indistinguishable from Y = 0 for the purposes of calculations in MESA , so this window should cover all cases of interest.  \nThe TDC solver identifies and handles a variety of cases. It takes advantage of the fact that d w/ d Y ≥ 0, which follows because the convective velocity always increases as a region becomes more unstable.  \nWe now discuss the different possible solutions. When Y > 0, the root-finding problem is monotonic because d w/ d Y > 0 and d R / d Y < 0. We approach this by performing a bisection search in Z followed by a NewtonRaphson solve. The bisection search ensures that the Newton-Raphson solve starts close to the true solution (we require the range ∆ Z ≤ 1 for termination). The Newton-Raphson solve then rapidly refines the solution to near machine precision and, crucially, imbues the solution with a differentiable dependence on the solver variables, tracked by auto diff ( § 2). Even if the bisection search finds an adequate solution, we still require at least one Newton-Raphson iteration to ensure that the result contains the partial derivatives needed for the MESA Jacobian.  \nWhen Y < 0 and the initial w = 0, the entropy gradient is stable against convection. There is no preexisting turbulence, and so w = 0 for the entire step and L = L rad . This makes Equation (16) linear in Y .  \nFinally, when Y < 0 and the initial w > 0, there can be up to three solutions to Equation (16):  \n- · In one solution, Y becomes large and negative. This forces w → 0 before the end of the time step (e.g., the first root of Equation (13) occurs before time dt passes), and L = L rad .\n- · In the other two solutions, Y becomes small and negative, and w declines but does not reach zero by the end of the step. Here L is carried by a mix of radiation and convection. There are two solutions because there is a tradeoff between the magnitude of Y and the decline of w , which compete in the w Y term in Equation (7).  \nMultiple solutions exist because, for long time steps, both Y and w can evolve significantly in a single step. One could force the time step to be smaller, such that there is just one solution. However, a global time step limit is often undesirable, especially in cases where the precise means by which convection decays (e.g., for a retreating convective boundary on the main-sequence) is not usually of interest. Hence, it is often preferable to select one of the multiple solutions.  \nWe disfavor the solution that decays fastest (e.g., Y is large and negative), as then convection decays on a dynamical timescale, which we suspect is unphysical. Rather, we favor the slower-decaying (e.g., smallermagnitude Y ) solution, which connects smoothly to the Y = 0 limit. These preferences yield this rule: we always select the solution with the smallest |Y| and thus the slowest-decaying convection speed .  \nTo find this solution, it is useful to examine R ( Y ) in a representative case, shown in Figure 1. Each solution is a choice of Y such that R ( Y ) = 0 (Equation 15). Solutions are convective when w > 0 and radiative otherwise.  \nThe first (slowest-decaying) solution is convective, with w ≈ 0 . 2 and Y ≈ 0 . 35. The second solution is also convective, with w ≈ 0 . 05 and Y ≈ 0 . 75. Finally, the third solution is purely radiative, with w = 0 and Y ≈ 0 . 95. The local maximum in R is due to the fact that as Y becomes more negative, w falls but |Y| rises, so the product w Y is not monotonic.  \nHowever, we do not know a priori how many solutions there are. There can be no more than three, but by changing L solver we can make the example shown in Figure 1 have just one (convective) solution. Our approach is to first detect the number of solutions and isolate the one of physical interest.  \nFigure 1. The residual of the luminosity equation R (upper) and w at the end of the time step (lower), plotted as functions of Y in a case where the correct Y < 0 and the initial w > 0. Vertical dashed red lines denote the two special points Y 0 (left) and Y 1 (right). Input parameters were chosen so that there are three solutions to R ( Y ) = 0 (black dots), with L solver = 0 . 1, L 0 = 2, ∇ ad = 1, c 0 = 20, ξ 1 = 0, ξ 2 = -1, w initial = 1, | J | = √ -Y , and δt = t -t 0 = 5, all in arbitrary units.  \n<!-- image -->  \nY  \nThe three solutions must be separated by two special points. The first ( Y 0 ) is the smallest-magnitude Y with w = 0, and the second ( Y 1 ) is the Y such that w > 0 and d R / d Y = 0. These are highlighted in Figure 1. Solutions of the first kind must occur at Y < Y 0 , solutions of the second kind must occur at Y > Y 0 , and at most one solution of the second kind occurs on either side of Y 1 .  \nWe use Y 0 and Y 1 to divide the interval -100 ≤ Z ≤ 100. The discriminant R is monotonic over each subinterval by construction, so in each case we can search for a root using a combination of bisection search and Newton-Raphson refinement. We check the intervals in order, from nearest to Y = 0 to furthest, and terminate the search as soon as a root is found.  \nBecause d w/ d Y > 0, we search for Y 0 using bisection in the interval -100 ≤ Z ≤ 100. We likewise identify Y 1 by a bisection search over -100 ≤ Z ≤ ln |Y 0 | .', '3.4.3. Relation to auto diff': 'TDC returns Y given L and the other solver variables. It additionally returns the partial derivatives of Y with respect to each of those variables. This relies, fundamentally, on the new automatic differentiation feature (see § 2). In particular, we used auto diff to calculate and propagate partial derivatives with respect to 33 variables of stellar structure through a Newton-Raphson solver , producing the partial derivatives of a root-finding procedure with respect to its inputs. The auto diff functionality enables the implementation of TDC.', '3.5. Reduction to Cox MLT': 'We now derive the modifications needed to ensure that TDC in MESA agrees with MLT in the limit of long time steps. While we use the α r =0 approximation in this section for clarity, the need for the correction is not removed by setting α r > 0.  \nIn TDC, the convective luminosity is given by Equation (7). In MLT, the convective luminosity is  \nL conv = 4 πr 2 f 2 ρc P Tv c Λ( ∇-∇ e ) h -1 (18)  \n(Ludwig et al. 1999), where v c is the convective velocity and f 2 is a parameter dependent on the choice of MLT prescription. Finally, ∇ e is the temperature gradient of a convective eddy, which is related to the efficiency parameter  \nΓ ≡ ∇-∇ e ∇ e -∇ L . (19)  \nWe may write  \n∇-∇ e = Γ 1 + Γ ( ∇-∇ L ) = Γ 1 + Γ Y , (20)  \nso  \nL conv = Γ 1 + Γ 4 πr 2 f 2 ρc P Tv c Λ Y h -1 . (21)  \nNext we identify w = √ 3 / 2 v c (in steady state) and Λ = αh , so  \nL conv = 4 πr 2 α Γ 1 + Γ ( 2 3 f 2 2 ) 1 / 2 ρc P Tw Y . (22)  \nThis is nearly the same as Equation (7). In particular, in Cox 1 MLT (Cox & Giuli 1968) f 2 = 1 / 2 and  \nL conv = 4 πr 2 α Γ 1 + Γ ρc P T w √ 6 Y , (23)  \nso the only difference is the term involving Γ.  \nThat term, which controls the convective efficiency, is a genuine difference between TDC-inRSP and MLT. We want TDC in MESA to match the outputs of MLT in the steady state limit, in agreement with Kuhfuß (1986), we modify Equation (7) to include the factor Γ / (1 + Γ), giving  \nL conv = 4 πr 2 α Γ 1 + Γ ρc P T w √ 6 Y . (24)  \nWe evaluate Γ by calling MLT with the same inputs as TDC. We then treat this as fixed during the TDC iterations, which allows us to still use the algorithm described in § 3.4.  \nWith these modifications, the luminosity equations now agree, subject to w = √ 3 / 2 v c in steady state. We now derive the conditions required to make this hold.  \nIn MLT, the convective velocity is given by  \nv 2 c = f 1 Λ 2 gδ ( ∇-∇ e ) h -1 , (25)  \nwhere f 1 is a parameter determined by the choice of MLT and  \nδ ≡ -∂ ln ρ ∂ ln T ∣ ∣ ∣ ∣ P = χ T χ ρ . (26)  \nUsing Λ ≡ αh we can write Equation (25) as  \nv 2 c = α 2 f 1 hg χ T χ ρ ( ∇-∇ e ) . (27)  \nNext, with Equation (20) we find  \nv 2 c = α 2 f 1 hg χ T χ ρ Γ 1 + Γ Y . (28)  \nIn TDC, we have identified the convection speed with v c ≡ √ 2 / 3 w , so we now proceed to prove that this is equivalent to that given by Equation (28). When the TDC discriminant J 2 < 0, then Y < 0 so the system is subadiabatic. Hence, at long times w = 0 and therefore v c = 0, which matches the MLT answer. When the discriminant is positive the system is convectively unstable, so we use Equation (12) and find  \nv c = -√ 2 3 1 2 ξ 2 ( J tanh λ + δtJ 4 + ξ 1 ) . (29)  \nThis solution was constructed assuming not only that J 2 > 0 but also that the relevant root is J > 0. As δt →∞ the tanh term approaches unity. In this limit  \nv c = -√ 1 6 1 ξ 2 ( J + ξ 1 ) , (30)  \nand we also have ξ 1 = 0 because D ρ/ D t = 0 and α r = 0 (Equation 10). Inserting the definition of J and expanding with Equations (9) and (11) we find  \nv 2 c = 1 6 J 2 ξ 2 2 = -2 ξ 0 3 ξ 2 = α 2 c P T ∇ ad Y 8 α D . (31)  \nComparing this with Equation (28), for equality  \nα 2 c P T ∇ ad Y MLT 8 α D = α 2 f 1 hg χ T χ ρ Γ 1 + Γ Y TDC . (32)  \nAs before, to obtain equivalence between MLT and TDC we need to substitute Y Γ / (1 + Γ) for Y in the velocity equation. In addition, we need to have  \nc P T ∇ ad 8 α D = f 1 hg χ T χ ρ . (33)  \nWith some rearranging, and using h = P/ρg , we find  \n1 8 α D f 1 = Pχ T χ ρ ρc P T ∇ ad . (34)  \nWith Γ 3 = 1 + ( P/ρc V T ) χ T ,  \n1 8 α D f 1 = (Γ 3 -1) c V χ ρ c P ∇ ad = Γ 3 -1 Γ 1 ∇ ad = 1 . (35)  \nIn Cox MLT f 1 = 1 / 8, and in TDC by default α D = 1, so the two sides are equal.  \nThe net result is that in the limit of long time steps, TDC and Cox MLT solve the same luminosity equation with the same inputs and so are mathematically identical. We find they agree numerically to around seven decimal places in Y , even when Y glyph[lessmuch] 1. The need for the correction is not removed by setting α r > 0. The asymptotic scaling in the inefficient limit is qualitatively different between (TDC with no correction and α r > 0) and (TDC with correction and α r =0). We further implement the calculation of convective mixing diffusivity and all other derived quantities using v c in the same way in both TDC and Cox MLT.  \nFigure 2 shows the importance of the Y → Y Γ / (1+Γ) correction. In both panels, the solutions for TDC and MLT lie on top of each other. The solution for TDC without the correction of Y in the equations, by contrast, deviates significantly in both panels. This deviation is starkest in the lower panel, which shows a different v c scaling in the inefficient ( ∇ rad -∇ ad glyph[lessmuch] 1) limit.', '3.6. Accreting White Dwarfs': 'WDs accreting He at rates glyph[lessorsimilar] 10 -6 M glyph[circledot] yr -1 undergo He shell flashes (Iben & Tutukov 1989). These flashes can lead to He nova (e.g., V445 Puppis; Ashok & Banerjee 2003), or even double-detonation type Ia supernovae  \n1  \nFigure 2. The values of Y (upper) and v c (lower) plotted as functions of ∇ rad -∇ ad for MLT, TDC, and TDC without the correction Y → Y Γ / (1 + Γ) in the convective velocity and luminosity equations. The time step was chosen to be long enough that TDC reaches equilibrium.  \n<!-- image -->  \n∇  \n-∇  \n(e.g., Shen & Bildsten 2009; Woosley & Kasen 2011; Kupfer et al. 2022). The time-dependent burning is controlled by three timescales: the local nuclear heating time,  \nt heat ≡ c P T glyph[epsilon1] nuc , (36)  \nbeing the characteristic timescale for temperature changes due to nuclear burning; the convective acceleration time,  \nt accel ≡ 4 J = 3 h √ 2 α D c P T ∇ ad Y , (37)  \nbeing the timescale over which convection varies (see Equation 12); and the local dynamical time,  \nt dyn ≡ h c s . (38)  \nIn steady state t accel is proportional to the eddy turnover time  \nt eddy ≡ αh v c , (39)  \nbut in cases of interest t accel and t eddy can be quite different.  \nShen & Bildsten (2009) showed that He shell masses of glyph[greaterorsimilar] 0 . 03 M glyph[circledot] on a ≈ 1 M glyph[circledot] WD can yield a t heat comparable to or shorter than t accel or even t dyn near the base of the convection zone (BCZ). TDC will yield different results than MLT in this limit.  \nWe construct these He flash models by accreting material comprising 99% 4 He and 1% 14 N by mass (as expected for solar metallicity stars that have undergone CNO burning) onto a 1 M glyph[circledot] carbon-oxygen WD at constant log( ˙ M/ M glyph[circledot] yr -1 ) between -7 . 1 and -7 . 4 in steps of 0 . 1 dex. Compressional heating results in a local temperature increase until the He shell ignites. Lower ˙ M results in weaker compressional heating and a more massive He shell at ignition. Because heat is transported from the temperature peak towards both the core and the surface, ignition occurs above the base of the freshly accreted layer. We stop the accretion once a convective zone appears at the ignition site, and continue evolving through the He flash. Both the total accumulated He shell mass and location of ignition are impacted by the included reaction chain 14 N( e -, ν ) 14 C( α, γ ) 18 O (NCO, Hashimoto et al. 1986; Bauer et al. 2017).  \nIn Figure 3, we label our models 1-4 at different ˙ M (with 1 corresponding to -7 . 1, 2 to -7 . 2, etc.) and note the masses enclosed by and exterior to the base of the convection zone (BCZ), which set the pressure P BCZ at ignition. The total accumulated He shell mass ranges from 0 . 03 to 0 . 08 M glyph[circledot] .  \nFigure 3 shows the evolution of T BCZ and P BCZ for models 1-4 with both TDC and Cox MLT. All models initially evolve at nearly constant P BCZ , which increases with He shell mass. As temperature increases in the convection zone, the envelope expands and reduces P BCZ . Concurrently, T BCZ reaches a maximum (Shen & Bildsten 2009). Thicker He shells reach higher peak T BCZ and larger ratios between t heat and t accel . For the contours here, t accel ≈ 10 t dyn (e.g., model 4 reaches t heat /t accel ≈ 0 . 01, and correspondingly t heat /t dyn ≈ 0 . 1.) When models show t heat glyph[lessorsimilar] t accel (models 3 & 4) and start expanding, TDC starts to deviate from Cox MLT, with greater deviations for thicker He shells. TDC shows higher T BCZ than Cox MLT at fixed P BCZ because TDC results in more superadiabatic convection. In contrast, when t heat glyph[greaterorsimilar] t accel (models 1 & 2) TDC and Cox MLT show good agreement in the evolution of P BCZ and T BCZ .  \nt accel glyph[greaterorsimilar] | dln T/ d t | -1 glyph[greaterorsimilar] t heat glyph[greaterorsimilar] t dyn  \nfrom model 1 to model 3, and ultimately to  \nt accel glyph[greaterorsimilar] t dyn glyph[greaterorsimilar] | dln T/ d t | -1 glyph[greaterorsimilar] t heat  \nFigure 3. Evolution of T BCZ and P BCZ during a He flash on a 1 M glyph[circledot] WD with various He shell masses. Solid and dashed lines correspond to models run with TDC and Cox MLT respectively. The values in parentheses indicate the masses in solar units enclosed within and exterior to the BCZ, respectively. Dotted lines are contours of fixed logarithmic ratio between t heat and t accel , from -2 to 1. Dot-dashed line gives the radiation pressure.  \n<!-- image -->  \nThe upper panel of Figure 4 compares several timescales for TDC models 1 & 3. The heating timescale, t heat , trends similarly with | dln T/ d t | -1 , but the latter is larger than t heat by factors of a few, because heat released by nuclear burning is distributed throughout some portion of the convection zone. Due to the sharp dependence of glyph[epsilon1] nuc on T , both timescales decrease sharply with T BCZ until the WD starts to expand. To reduce the noise in | dln T/ d t | -1 displayed in Figure 4, we fit it with a polylogarithmic function. The difference between t heat and | dln T/ d t | -1 decreases with thicker He shells, as heat released by nuclear burning is increasingly trapped locally.  \nAnother relevant timescale is t accel , evaluated at maximum v c . At log( T BCZ / K) glyph[lessorsimilar] 8 . 4, t accel ≈ (3 / 4) t eddy . At log( T BCZ / K) glyph[greaterorsimilar] 8 . 4, t accel evolves more quickly than t eddy , becoming up to 3 (6) times smaller than t eddy in model 3 (4). This is because convection is no longer in steady-state, as t heat glyph[lessorsimilar] t accel for log( T BCZ / K) glyph[greaterorsimilar] 8 . 4 -8 . 5 (see also Glasner et al. 2018). At minimum t heat , the hierarchy of timescales changes from  \n| dln T/ d t | -1 glyph[greaterorsimilar] t heat glyph[greaterorsimilar] t accel glyph[greatermuch] t dyn  \nto  \nFigure 4. Evolution of various time scales in the convection zone (upper), the time derivative of convective velocity at a fixed mass coordinate (middle), and the superadiabaticity at the same location (lower). The upper panel compares the heating timescale (solid line), actual timescale for temperature change (dashed line), convective acceleration timescale (dot-dashed line), and dynamical timescale (dotted line), for models 1 and 3 (TDC only). The first three timescales are evaluated at the BCZ, and the last at maximum convective velocity. Line styles in the middle and lower panels follow that of Figure 3, and in the middle panel we show the gravitational acceleration at the BCZ ( g BCZ , dotted line) for comparison.  \n<!-- image -->  \nin model 4. The fact that convection is not able to reach a steady state on the evolutionary timescale of the He flash explains the difference between TDC and Cox MLT in models 3 & 4 (see Figure 3).  \nWe illustrate the difference in the evolution of v c and Y between TDC and Cox MLT in the middle and lower panels of Figure 4. For each TDC and Cox MLT pair, we locate the mass coordinate at which v c peaks when  \nFigure 5. Profiles of T and P within model 3 at different moments, when both TDC (solid) and Cox MLT (dashed) reach the same T BCZ . Curves of same color correspond to identical T BCZ . Dotted and dot-dashed lines take the same meaning as in Figure 3.  \n<!-- image -->  \nlog( L nuc / L glyph[circledot] ) = 9 (arbitrarily chosen), and evaluate D v c / D t and Y during the initial acceleration phase. Initially, TDC and Cox MLT show good agreement when t heat glyph[greatermuch] t accel (log( T bcz / K) glyph[lessorsimilar] 8 . 4). When t heat glyph[lessorsimilar] t accel (log( T bcz / K) glyph[greaterorsimilar] 8 . 5), TDC shows slower evolution in v c and larger Y than Cox MLT. As v c is lower in TDC, heat is less efficiently transported out of the BCZ, resulting in higher T BCZ and Y near maximum. With a thicker He shell, D v c / D t may become comparable to g (especially for Cox MLT model 4).  \nWe now study the evolution of model 3 in detail. In Figure 5, we compare seven snapshots of the TDC and Cox MLT models in T -P space, when both reach the same T BCZ . The three coolest pairs of curves show good agreement and little superadiabaticity (third panel of Figure 4). For the subsequent three hotter pairs, Y grows up to order unity near peak T BCZ . Once t heat glyph[lessorsimilar] t accel , heat is trapped more locally in TDC compared to Cox MLT. Therefore, TDC reaches the same T BCZ earlier in the evolution, and has a higher P BCZ due to comparably colder outer layers. Likewise, TDC shows less evolution in T -P near the top of the convection zone and more superadiabaticity near the BCZ, again because of stronger heat-trapping near the BCZ.  \nFigure 6 compares TDC and Cox MLT in model 3 as a function of mass coordinate. The two show reasonable agreement in v c when log( T BCZ / K) glyph[lessorsimilar] 8 . 5. At this point, t heat drops below t accel (see Figure 4), which leads to TDC yielding lower v c than Cox MLT. For the same reason, v c near the top of the convection zone appears frozen in TDC for 8 . 6 glyph[lessorsimilar] log( T BCZ / K) glyph[lessorsimilar] 8 . 8.  \nAt fixed T BCZ , TDC shows lower T throughout the convection zone, reflecting a local buildup of heat at the BCZ. Since TDC carries heat out of the BCZ less effi-  \nFigure 7 shows the default MESA EOS boundaries for solar and pure-metal (50% 12 C, 50% 16 O by mass) compositions. Broadly, we prioritize HELM (Timmes & Swesty 2000) at high T and low ρ for handling of the electron-positron plasma. Elsewhere we prioritize Skye (Jermyn et al. 2021), limited by partial ionization at lower T and ρ . We then prioritize FreeEOS (Irwin 2004), then OPAL (Rogers & Nayfonov 2002) and SCVH (Saumon et al. 1995), and finally, when there are no other options, we use an ideal gas with radiation. Blending boundaries between the different EOS prescriptions are set to defaults that have been motivated by specific use cases. For example, we have chosen the boundaries between FreeEOS and Skye such that solar models at the age of the Sun stay fully on FreeEOS and do not encounter the FreeEOS-Skye blend.  \n<!-- image -->  \n/circledot  \nFigure 6. Same as Figure 5, but instead showing T , s , v c , and X ( 12 C) as a function of mass coordinate from top to bottom. In the third panel, c s is shown for comparison. The colors match the colorbar in Figure 5.  \nciently, it also shows a stronger entropy gradient, and for 8 . 6 glyph[lessorsimilar] log( T BCZ / K) glyph[lessorsimilar] 8 . 8, T and s show little evolution near the top of the convection zone.  \nBoth TDC and Cox MLT show appreciable abundance gradients, as 12 C is produced near the BCZ but there is insufficient time for it to be transported outwards. Cox MLT shows higher 12 C abundance overall, as it has more time to reach the same T BCZ and larger v c . As TDC modifies both the T and X i profiles, it may impact the potential for the ignition to develop into a detonation that would result in a thermonuclear transient.  \nIn summary, we see that convection in TDC adjusts more slowly to changes in heating than in Cox MLT. This results in slower, more superadiabatic convection during rapidly-burning phases of evolution. The incorporation of the dynamics of how convection grows and decays is now possible and enabled by default in MESA via TDC.', '4. EQUATION OF STATE': 'MESA models require thermodynamic quantities over a large span of T , ρ , and X i . This involves calling the MESA EOS ∼ 10 4 -10 10 times, depending on the chosen local physics and the number of iterations, cells, and time steps. It would be ideal to have a single EOS that accurately represents the relevant physics in all regimes, obeys all thermodynamic consistency relations to the limits of the arithmetic, and is as efficient in storage and execution as possible. Below we report progress towards this ideal.', '4.1. Skye': 'Skye is an EOS for fully ionized matter (Jermyn et al. 2021). A motivation for developing Skye was eliminating the blend between HELM and PC (Potekhin & Chabrier 2010). There is a blend between HELM  \n<!-- image -->  \nFigure 7. Default EOS boundaries for a solar-like composition (Z = 0.02, left) and for a pure-metal composition (Z = 1, 50% 12 C and 50% 16 O by mass, right). The colors and labels indicate an ideal gas approximation, Skye, FreeEOS, OPAL, SCVH, HELM, and blends between two EOSs. The black curve in the left panel shows a model for the present day Sun.  \n<!-- image -->  \nand Skye that occurs at much higher T and lower ρ (see Figure 7), where the two agree. Skye includes the effects of positrons, relativity, and electron degeneracy (Timmes & Swesty 2000; Baturin et al. 2019), Coulomb interactions (Ichimaru et al. 1987; Potekhin & Chabrier 2010; Potekhin & Chabrier 2000; DeWitt & Slattery 1999; Baiko et al. 2001; Farouki & Hamaguchi 1993), non-linear mixing effects (Caillol 1999; Potekhin et al. 2009; Potekhin & Chabrier 2013; Ogata et al. 1993; Medin & Cumming 2010), and quantum corrections (Hansen & Vieillefosse 1975; Nagara et al. 1987; Carr et al. 1961; Potekhin & Chabrier 2000; Potekhin & Chabrier 2010; Baiko 2019; Baiko & Chugunov 2022). Skye determines the point of Coulomb crystallization in a self-consistent manner, accounting for mixing and composition effects. A defining feature of Skye is the use of analytic Helmholtz free energy terms and automatic differentiation (see § 2) to provide thermodynamic quantities. Skye is thus readily extendable to new physics by including additional terms in the free energy (Jermyn & Timmes 2022).  \nSkye is both a standalone software instrument and integrated into MESA . The two implement the same input physics and options. At times this has required modifications of other parts of MESA . Here we describe the most important of these modifications.', '4.1.1. Crystallization': 'Skye determines the crystallization phase transition by minimizing the Helmholtz free energy, which permits derivatives to be discontinuous at the transition. For instance, the entropy discontinuity reflects the latent heat of crystallization. This posed a challenge in MESA . Con-  \nsider the expression  \nglyph[epsilon1] grav ≡ -T D s D t . (40)  \nThe entropy s ≡ -∂F/∂T | ρ undergoes a discontinuity at the phase transition. If D s /D t is evaluated by finite differences, then no time step will be small enough to produce a converged result for glyph[epsilon1] grav . We could write  \nD s D t = ∂s ∂T ∣ ∣ ∣ ∣ ρ D T D t + ∂s ∂ρ ∣ ∣ ∣ ∣ T D ρ D t , (41)  \nbut this form misses the latent heat of the phase transition because, except for the infinitesimal vicinity of crystallization, thermodynamic derivatives of s contain no information about the transition. At the phase transition, derivatives of s contain a Dirac delta contribution, which cannot be directly implemented in numerical calculations. The choice is between poor convergence (finite differences of s ) or neglecting the latent heat (Equation 41).  \nTo address crystallization, Skye returns a parameter φ that provides a smoothed representation of the phase. Specifically, φ = 1 in the solid phase, φ = 0 in the liquid phase, and near the phase transition φ smoothly interpolates between these limits. The transition in φ is tuned so that the crystallization boundary is numerically resolved and yet spans a small fraction of a stellar model. Using φ , Skye then constructs a smoothed version of the latent heat of crystallization, which is only significant in the transition region. This allows use of Equation (41) to avoid numerical issues near the phase transition, but requires that we include an extra heat source in the en-  \nquation to capture the latent heat:  \nglyph[epsilon1] latent ≡ L T Dln T D t + L ρ Dln ρ D t , (42)  \nwhere L T and L ρ represent the differences between smoothed and original versions of the entropy derivatives T∂s/∂ ln T and T∂s/∂ ln ρ . The original derivatives lack the latent heat, while the smoothed ones contain it, so L T and L ρ produce additional heating. With this procedure, MESA is able to model phase transitions, remain numerically converged, and accurately capture the latent heat of crystallization. This procedure smears only the latent heat of crystallization and does not smear the thermodynamics of the phase transition, which would produce unphysical results such as negative sound speeds.  \nThe Skye EOS approach represents a significant improvement for the MESA latent heat treatment. Previously, MESA relied on a finite difference of the entropy calculated in the PC EOS for solid and liquid phases so that latent heat could be included in glyph[epsilon1] grav via Equation (40), smoothing this quantity near the phase transition for numerical convergence (MESA IV). Another common approach is to include latent heat release with an explicit heating term using l cr ≈ 0 . 77 k B T/ 〈 A 〉 m p based on the calculation of Salaris et al. (2000). Our new approach based on Skye has the advantage that the phase diagram and latent heat release are both calculated from first principles and are self-consistent with the underlying thermodynamics of the EOS. Jermyn et al. (2021) showed that the net latent heat release is commensurate with the Salaris et al. (2000) value.', '4.2. FreeEOS': "We use FreeEOS version 2.2.1 (Irwin 2004) to expand the chemical composition parameter space covered by partial ionization, as compared to the OPAL tables. This replaces the eosPTEH tables of MESA V. FreeEOS minimizes a Helmholtz free energy to span essentially the same thermodynamic range as OPAL.  \nThe FreeEOS tables generated for MESA use the 'EOS1' mode, which is the highest level of physical accuracy provided by FreeEOS. The tables are parameterized by the metal mass fraction Z = 0, 0.02, 0.04, 0.06, 0.08, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, and 1.00. All tables assumed a scaled-solar chemical composition based on Grevesse & Sauval (1998). For Z ≥ 0.80, there is also a set of tables with X ( 12 C) = X ( 16 O) for use with WD interiors. For each Z a range of H mass fraction values between 0 and 1 -Z are provided, allowing for a complementary range of He mass fractions. The tools to generate a new set of MESA EOS tables for an  \narbitrary chemical composition using FreeEOS are provided in MESA DIR/eos/eosFreeEOS builder with the exception of the FreeEOS library, which can be downloaded from the FreeEOS repository.", '4.3. EOS Blends': 'The MESA EOS blends several EOS prescriptions. Each EOS returns fundamental quantities and the partial derivatives of those quantities. The blends of fundamental quantities and derivatives are treated differently because they are used by MESA for different purposes. Fundamental quantities enter into physical equations, and so must be physical (e.g., positive sound speed), while their derivatives are used to construct the solver Jacobian, and so must represent accurate derivatives of the fundamental quantities.  \nThe EOS returns a vector res containing fundamental EOS quantities such as e , s , and c V (see MESA I Table 3), as well as blending fractions for the various EOS components. The EOS also returns corresponding vectors d dlnd and d dlnT of partial derivatives of each of the quantities in res with respect to ρ and T .  \nAt the boundary between a pair of EOS prescriptions (EOS1 and EOS2) we calculate blends of res , d dlnd and d dlnT independently. The EOS at a point in the blending region between EOS1 and EOS2 is evaluated with blending coefficient α ∈ [0 , 1] representing the fraction of EOS1, and 1 -α representing the fraction of EOS2. We construct blending coefficients using the quintic polynomial  \nα = 6 x 5 -15 x 4 +10 x 3 , (43)  \nwhich maps the interval x ∈ [0 , 1] (representing distance across a blend in ρ or T ) onto the interval α ∈ [0 , 1] with zero slope at the blending boundaries. The blending coefficients therefore have non-zero derivatives with respect to ρ and T in blending regions. Quantities in the resulting res blend vector are evaluated as a linear mix using the blending coefficient,  \nres blend = α res 1 +(1 -α ) res 2 . (44)  \nOur choice of a quintic polynomial for the blending coefficient ensures that both α and (1 -α ) are non-negative everywhere in the blending region, and therefore the EOS blending never introduces negative quantities into blends of non-negative values for the EOS res vector. For the derivative vectors, we include additional terms to account for the derivatives of the blending coefficients,  \nd dlnd blend = α d dlnd 1 + ∂α ∂ ln ρ res 1 +(1 -α ) d dlnd 2 -∂α ∂ ln ρ res 2 , (45)  \nand similarly for d dlnT . Including these terms for the blending coefficients in the derivative blends provides correct derivatives for the solver, reducing the number of Newton iterations.  \nSome quantities in the fundamental EOS res vector are themselves derivatives of other EOS quantities, such as c V ≡ ( ∂e/∂T ) ρ . The different blending treatments for EOS quantities and their derivatives mean that thermodynamic identities may be violated in blending regions. Physical equations such as the energy equation must use quantities such as c V from res rather than the theoretically equivalent but numerically different derivative quantities from the d dlnT vector. The latter can lead to unphysical results such as negative heat capacities or negative sound speeds. This inconsistency is unavoidable so long as we must blend between different EOS prescriptions.', '4.4. Thermodynamic Consistency': 'One desirable feature in an EOS is thermodynamic consistency, which ensures that the Maxwell relations hold - e.g., mathematically equivalent forms of the equations of stellar structure are also numerically equivalent within the floating point precision of the arithmetic. Unfortunately, several of the EOS prescriptions in MESA are not fully thermodynamically consistent. This can cause errors in energy conservation, making mathematically equivalent formulations of the structure equations behave differently.  \nHere we report on the current state of thermodynamic consistency in MESA . Figure 8 compares the consistency measure  \ndse ≡ T ∂s/∂T | ρ ∂e/∂T | ρ -1 (46)  \nfor the MESA EOS with PC and OPAL (left, former default) and with Skye and FreeEOS (right, current default). The quantity dse is zero in thermodynamically consistent systems.  \nAs Skye derives all quantities from partial derivatives of a Helmholtz free energy, it is thermodynamically consistent to near machine precision. Without Skye, the corresponding regions of the EOS are covered by PC and HELM. The regions with Skye active show thermodynamic consistency to near machine precision, representing a significant improvement for log( ρ/ g cm -3 ) glyph[greaterorsimilar] 3. The band at log ( T/ K) ≈ 9 -10 in the right panel is due to a blend in the EOS from Skye to HELM, which is required to remedy a floating point loss-of-precision issue in Skye when electron-positron pairs dominate the EOS. In the left panel the PC region shows a stripe of high error due to Coulomb crystallization. FreeEOS is thermodynamically consistent to near machine precision. Our  \ncurrent method of interpolating the MESA FreeEOS tables does not preserve this property. Still, these tables show significant improvement relative to OPAL.', '5. ENERGY EQUATIONS': 'Section 3 of MESA V highlighted the importance of numerical energy conservation in MESA models, and introduced a new form of the energy equation aimed at improving energy conservation. This new form motivated several solver improvements, leading to tighter tolerances for equation residuals and corrections. We now advance that discussion by further explaining the multiple formulations of the energy equations in MESA . We contrast the advantages conferred by each formulation across different applications. We also clarify the meaning of the quantities rel E err and rel run E err reported for MESA runs, and elucidate what these quantities do and do not tell us about the quality of the MESA solutions.  \nAfter reviewing the energy-equation fundamentals and recent implementation improvements in § 5.1, we explore an idealized example problem without any composition changes or EOS complexities in § 5.2. This example motivates a new time-centered approach for evaluating the eps grav form of the energy equation, and demonstrates that a lower value of rel run E err does not always indicate evolution that is more physically accurate. In § 5.3 we describe the additional complexities introduced by thermodynamic inconsistencies that can be present in the EOS, especially in EOS blending regions, and how these manifest in different ways for different energy equation implementations.  \nFinally, in § 5.4 we illustrate the various contributions to energy error terms through the example of a 1 M glyph[circledot] star including both composition changes due to nuclear burning and EOS blends and inconsistencies. This example demonstrates that the quantity rel run E err must be interpreted differently for different forms of the energy equation. When using the dedt form of the energy equation, the energy error reflects the quality of the residuals from the MESA solver, even though larger energy errors associated with the EOS are still present in the model. When using the eps grav form of the energy equation, the energy error reports much larger values that reflect the presence of these EOS errors, even when the quality of solutions may be comparable to or better than the dedt form.  \nConvergence tests and comparisons between multiple forms of the energy equation remain vital for understanding the reliability and accuracy of solutions in different regimes. Significant progress has been made in ensuring that different forms of the energy equation con-  \n<!-- image -->  \nFigure 8. The logarithm of the thermodynamic consistency measure dse for the MESA EOS computed with PC and OPAL (left) and Skye and FreeEOS (right) with X = 0 . 7, Y = 0 . 28, and Z = 0 . 02, with Z comprising equal mass fractions of 12 C and 16 O.  \n<!-- image -->  \nverge to the same result. In degenerate conditions, the eps grav forms generally perform better (i.e., they are closer to the converged answer at a given time resolution). With the dedt form, the numerical energy conservation error often measures the quality of the solution (i.e., the size of the residuals). Focusing on improving that quantity has driven significant solver improvements and motivated the development of an accurate energy accounting infrastructure. This energy accounting work has also motivated improving the eps grav form to account for composition changes, as well as an implicit trapezoidal time-centering scheme. MESA now includes these changes by default when using the eps grav form of the energy equation. Further progress rests on improvements to the quality of the EOS.', '5.1. Fundamentals and Implementations': "MESA has two primary energy equations. One, called the ' eps grav form', 2 is the standard stellar structure energy equation (e.g., Kippenhahn et al. 2012) and is the equation introduced in MESA I. This equation is  \n∂L ∂m = glyph[epsilon1] + glyph[epsilon1] grav , (47)  \nwhere L is the luminosity, glyph[epsilon1] is a specific energy generation source term (e.g., nuclear reactions, neutrinos), and  \nglyph[negationslash]  \nglyph[epsilon1] grav ≡ -D e D t + P D D t ( 1 ρ ) = -T D s D t -∑ i ( ∂e ∂X i ) s,ρ, { X = X i } D X i D t . (48)  \nIn practice, the total Lagrangian time derivative of e is expanded and further manipulated to reach the final form evaluated in MESA (see MESA IV, § 8).  \nThe other, called the ' dedt form', is an energy equation for the time evolution of the total specific energy of a Lagrangian cell,  \nD D t ( e + 1 2 u 2 -Gm r ) = glyph[epsilon1] -∂ ∂m ( L + P A u ) , (49)  \nwhere u is cell velocity and A = 4 πr 2 is the area of the cell face. The relationship between these two forms was derived in MESA IV, Section 8.3 and the dedt form was introduced as a powerful tool in support of improved numerical energy conservation in MESA V, § 3.  \nWhen solutions are numerically converged (i.e., have sufficient space/time resolution to give resolutionindependent results) and the EOS is thermodynamically consistent and provides correct partial derivatives (see § 5.3), these two equations should give identical results. Conversely, the solutions may differ when unconverged.  \nThe error in numerical energy conservation during a step, E step err , is evaluated as the difference between the change in total energy of the model across the time step and the expected change in total energy due to known energy sources and sinks (e.g., nuclear reactions, neutrinos, surface luminosity). Total energy is defined as  \nE tot ≡ ∫ ( e + 1 2 u 2 -Gm r ) d m = ∑ k dm k ( e k + 1 2 u 2 k -Gm k r k ) , (50)  \nwhere dm k is the mass contained within cell k . Additional terms for rotational kinetic energy can also be included in Equation (50) when rotation is enabled, and turbulent energy is included for RSP models.  \nA cumulative sum of the per-step energy errors, E run err , is tracked during a run. When divided by the total energy at the end of the step, E step err and E run err respectively become the quantities rel E err and rel run E err that are reported by MESA . As stated in MESA V, these quantities are primarily meant to represent a measure of the numerical reliability of solutions accepted for MESA evolution steps, rather than a measure of physical validity and completeness of MESA models. In § 5.3 and § 5.4, we focus on further clarifying the meaning of these energy error quantities, which require a different interpretation when using the eps grav form of the energy equation than when using the dedt form.  \nMESA does not solve its discretized, finite-mass form of the stellar structure equations exactly. When a trial solution is accepted, the residual difference between the left- and right-hand sides of the equation becomes an error in numerical energy conservation. Therefore, one necessary step in ensuring good numerical energy conservation is to select tight tolerances for the acceptance of a solution. This requires sufficiently high quality derivatives in the Jacobian that the solver can reach these tolerances in a reasonable number of Newton iterations (see § 3 of MESA V).  \nHowever, even achieving zero residuals is not sufficient to ensure numerical energy conservation. When MESA modifies the stellar model outside of the Newton solve, the resulting changes in total energy must be correctly included in the accounting. When physically appropriate, compensating energy source terms must be included in the equations that are solved during the Newton iterations. For example, mass changes of the stellar model are one such process, and the procedure that ensures numerical energy conservation is described in § 3.3 of MESA V. At that time, this procedure was applied only when using the dedt form of the equation. Now, it is used with all forms of the energy equation, and the less general approach originally used with the eps grav form (MESA III, § 7) has been removed from MESA .  \nThe composition changes associated with element diffusion ( § 3 of MESA IV) and convective premixing ( § 5 of MESA V) are also incorporated in an operator-split  \nmanner (i.e., adjustments to the model made outside of the Newton iterations for the implicit structure solve during an evolutionary step; see also § 10), and so require special accounting. The energy changes due to these composition changes are now tracked and compensating source terms are added to the equations, improving numerical energy conservation.  \nNon-conservation of numerical energy can also occur when the equations being solved are approximated in ways that do not conserve energy. Historically, the default MESA implementation of glyph[epsilon1] grav (MESA I, Equation 12) dropped the term associated with composition changes. While the energy associated with composition changes is dwarfed by the energy released by nuclear reactions (see MESA IV), the integrated energy error introduced by dropping this term is not negligible compared to the value of E tot by the end of the MS.  \nIn MESA V, Figure 25, the ' dLdm -form' calculation (right panel) did not include composition changes in glyph[epsilon1] grav , and so the large values of the relative energy error shown during core He burning effectively quantify the impact of dropping the composition term rather than characterizing the numerical quality of the MESA solution. In this case, the scale of the reported error appears significant because MESA adopts E tot as the reference value for checking cumulative numerical energy conservation. A larger reference value, like the timeintegrated radiated energy of the star, is typically used to justify dropping the composition term from glyph[epsilon1] grav .  \nA continued focus on numerical energy conservation requires equations that are energy conserving, so MESA now includes the composition term in its default implementation of glyph[epsilon1] grav . With ( ρ, T ) as the thermodynamic structure variables, we have  \nglyph[epsilon1] grav = -c V T Dln T D t -[ ρ ( ∂e ∂ρ ) T -P ρ ] Dln ρ D t + glyph[epsilon1] grav ,X i , (51)  \nwhere c V ≡ ( ∂e/∂T ) ρ . As shown in MESA IV, Equation (65), MESA implements the equivalent expression  \nglyph[epsilon1] grav = -c P T [ (1 -∇ ad χ T ) Dln T D t -∇ ad χ ρ Dln ρ D t ] + glyph[epsilon1] grav ,X i , (52)  \nwhere c P ≡ ( ∂e/∂T ) P -( P/ρ 2 )( ∂ρ/∂T ) P and ∇ ad ≡ ( ∂ ln T/∂ ln P ) s . The composition term is  \nglyph[negationslash]  \nglyph[epsilon1] grav ,X i ≡ -∑ i ( ∂e ∂X i ) ρ,T, { X = X i } D X i D t . (53)  \nWhen implemented in MESA , the quantity glyph[epsilon1] grav ,X i is evaluated as a finite-difference approximation to the direc-  \ntional derivative along the change in the composition vector over the time step:  \nglyph[epsilon1] grav ,X i = -1 δt [ e ( ρ, T, { X i } ) -e ( ρ, T, { X i } start )] . (54)  \nThis is analogous to the approach used in evaluating the spatial composition derivatives that enter into the Brunt-Vaisala frequency ( § 3.3 of MESA II). In addition to being simpler to evaluate, this approximation is numerically convenient because it only requires first derivatives of e with respect to composition in order to form the Jacobian. The MESA eos module and its interface with MESAstar have been upgraded either to provide these partial derivatives when available, or to construct approximations to these partial derivatives for the Jacobian based on finite differences using small variations of the composition when analytic derivatives are not available.  \nThe total derivatives of the structure variables in Equation (52) are evaluated as their differences over the time step. In previous implementations of the eps grav form of the energy equation, the thermodynamic quantities that multiply the total derivative quantities were evaluated at the end of the step (in the standard MESA backwards-Euler approach). As a means of further improving numerical energy conservation when using the eps grav form, we have now introduced a higher-order (in time) version of glyph[epsilon1] grav using the implicit trapezoidal rule. This replaces end-of-step quantities with timecentered versions (i.e., averages of the values at the start and end of the step). We refer to this as ' eps grav (centered)' in contrast to the previous implementation, which we indicate as ' eps grav (end of step).' As we shall show in the following sections, including both composition changes and time-centering in the eps grav implementation greatly improves energy conservation, so we now include both of these improvements by default in MESA when using the eps grav form.  \nIn the following sections, we use test cases to demonstrate the performance and physical meaning of numerical energy conservation under the various forms of the energy equation in MESA . We also show that in some circumstances, such as degenerate stars, the eps grav form of the energy equation converges to accurate entropy and temperature evolution substantially faster than the dedt form does, even while reporting larger errors in numerical energy conservation.", '5.2. Results: carbon kh': "As an illustrative test case, we follow an initially lowdensity, 1.3 M glyph[circledot] sphere of pure carbon as it undergoes Kelvin-Helmholtz contraction. The model begins at a  \ncentral density of log( ρ c / g cm -3 ) = 1 and we follow the contraction over a factor glyph[greaterorsimilar] 10 7 increase in ρ c . Nuclear reactions are not considered. For simplicity, we assume that the radiative opacities are given by electron scattering and include standard thermal neutrino losses. This model does not experience convection. We exclusively use the HELM EOS, as the use of a single EOS that is formulated from the Helmholtz free energy avoids most of the EOS inconsistencies that we discuss in § 5.3.  \nThis case is not meant to model a real object, but provides a simple example problem that has neither mass changes nor composition changes. It is nonetheless demanding as the conditions in the star vary tremendously during the evolution as material goes from nondegenerate conditions to conditions of relativistic electron degeneracy, and the dominant energy loss mechanism transitions from radiative diffusion to opticallythin neutrino cooling.  \nWe explore three different versions of the energy equation: the dedt form, the eps grav form (end of step), and the eps grav form (centered). We use a temporal convergence study to illustrate the performance of the different variants of the energy equation. For each equation, we show three time resolutions, and compare against a family of ultra-resolution runs that serve as reference solutions. These ultra-resolution reference runs still show small differences depending on which form of the energy equation is selected, so we also show that smaller level of disagreement to indicate the level of differences that should be interpreted as significant. We interpret the small magnitude of disagreement between ultra-resolution runs as evidence that the different versions of the energy equation are converging to the same result for sufficiently high resolution.  \nFigure 9 shows the trajectory of T c versus ρ c . While this calculation does not consider nuclear reactions, in calculations that do, the T and ρ sensitivity of the nuclear reaction rates motivates obtaining solutions that are converged in these quantities (e.g. Schwab et al. 2015). This example does include thermal neutrinos, which lead to central cooling at high density. The top panel shows that the two eps grav versions agree (to within the line width) at all resolutions, while the dedt form shows visible differences during the evolution after the model has reached its maximum T c . The level of difference from the reference solution is shown in the middle panel. All forms exhibit first-order convergence, where a 1 dex reduction in the time step leads to a 1 dex reduction in the error in T c . However, at a fixed resolution, the eps grav forms show similar performance to each other and superior performance relative to the dedt form.  \nFigure 9. Convergence tests of the evolution of the central quantities for the contracting 1.3 M glyph[circledot] sphere of pure carbon using different forms of the energy equation. The central density is a proxy for time. The top panel shows T c . The middle panel shows the T c difference from a reference solution. The three line styles represent the three levels of resolution as indicated by the legend in the top left corner of the upper panel. The grey region marks the level of agreement between the set of reference solutions. The lower panel shows the time step for each resolution, with the time step of the runs used in the ultra-resolution reference set indicated as the solid grey line.  \n<!-- image -->  \nIn order to understand why the eps grav forms perform better under degenerate conditions, consider an adiabatic change, d e + P d(1 /ρ ) = T d s = 0 at fixed composition. This expression is satisfied exactly for infinitesimal changes and a perfect EOS. When we integrate across a time step, we know the integral of total time derivatives (e.g. D e/ D t or Dln ρ/ D t ) exactly, but approximate the integral over the time step for quantities that are not total time derivatives. The extent to which our scheme will fail to reproduce an adiabatic evolution is the error in approximating these other integrated quantities appearing in the energy equation (e.g., P or c V ). Recall that the usual backwards Euler approach in MESA is effectively like assuming that the non-  \ntal-time-derivative part is constant and equal to the end-of-time step value (e.g., P = P end ).  \nFor nearly adiabatic evolution in electron degenerate conditions, we have a cancellation between large d e and P d(1 /ρ ) terms, but this cancellation ends up incomplete in MESA because the evaluation of the former term is exact while the latter has error. The error is usually small compared to the order of the terms being subtracted, and so imperfect cancellation often will not introduce large errors. But in degenerate material, the scale of the cancelling terms is larger than the thermal energy by roughly the degeneracy parameter η ≡ µ e /k B T , where µ e is the electron chemical potential. Therefore, otherwise small cancellation errors can be amplified by a factor of η for the temperature evolution. 3  \nBy contrast, when we write the glyph[epsilon1] grav form, this cancellation for adiabatic evolution instead occurs in [ ρ ( ∂e/∂ρ ) T -P/ρ ] (Equation 51) which is replaced with c P T ∇ ad χ ρ in the form of Equation (52) that MESA uses for its eps grav implementation. This captures adiabatic density evolution in terms of EOS derivative quantities that are not subject to cancellation errors. Instead, accuracy in this form is limited by the accuracy of our approximations over finite time steps for thermodynamic quantities like ∇ ad appearing in the energy equation. The error associated with time discretization, [ c P T (1 -∇ ad χ T )] end ∆ln T +[ c P T ∇ ad χ ρ ] end ∆ln ρ , is at least a factor of ∼ η smaller than the cancellation error, and in practice can be even better.  \nFigure 10 shows the total energy of the model. The top panel shows that all runs agree in this quantity within the line width, while the middle panel reveals the level of relative difference. We emphasize that even though the dedt runs report by far the best cumulative energy error as a measure of step-to-step internal energy consistency (lower panel of Figure 10), they nevertheless show less accurate evolution of the total energy and temperature relative to the ultra high-resolution reference runs (middle panel of Figure 10). This is because the cumulative energy error reports the degree to which energy is conserved by evolution steps, while the total energy is a function of the global stellar structure which can slowly diverge even with zero reported energy error. This reflects the fact that 'energy error' as reported by MESA is primarily a measure of the internal consistency of the stellar structure solver, and should not be construed as always reflecting globally accurate energy  \nFigure 10. Convergence tests of the evolution of the total energy for the contracting 1.3 M glyph[circledot] sphere of pure carbon using different forms of the energy equation. The top panel shows the negative total energy. The middle panel shows the difference from a reference solution. The grey region marks the level of agreement between the set of reference solutions. The lower panel shows the relative cumulative energy error for each run.  \n<!-- image -->  \nevolution. This reported error is still a useful diagnostic for MESA models, but must be interpreted with caution.  \nAll forms of the energy equation in Figure 10 approximately show first-order convergence in the total energy. The lower panel shows that the different forms exhibit notably different behaviors with increasing time resolution. The dedt form has excellent numerical energy conservation that does not depend on time resolution. The error is roughly the error due to the non-zero residuals in the solution of the energy equation. The eps grav forms display worse performance in this quantity, though the error shrinks as the time step decreases. The 'end of step' form shows first-order convergence, while the 'centered' variant exhibits more rapid, second-order convergence with smaller numerical energy conservation errors at fixed resolution. We would expect these trends to continue until the numerical energy conservation error is no longer dominated by errors due to the temporal  \ndiscretization, at which point it reaches the floor set by non-zero residuals or imperfect EOS thermodynamics.  \nThe pure carbon case shows that the time-centered eps grav form of the energy equation is the best choice for models evolving under degenerate conditions, with the best balance between accurate temperature evolution and step-to-step energy conservation according to Figures 9 and 10. This case was idealized to focus on the effects of finite equation residuals and time discretization. We now move on to discussing the additional complexities introduced by EOS imperfections.", '5.3. Quantifying EOS shortcomings': 'The value of e returned by the EOS is an essential ingredient in evaluating the total energy of the model, and high-quality partial derivatives of EOS quantities are critical for accurate and efficient solver performance. We now discuss three primary EOS issues that influence energy conservation and solver performance. First, an EOS may return low quality partial derivatives that degrade convergence of the implicit solver. We now mitigate this with more careful derivative accounting described in § 5.3.1. Second, an EOS may have internal inconsistencies in its reported thermodynamics. We have mitigated this by upgrading the MESA EOS patchwork with Skye (Jermyn et al. 2021) and FreeEOS (Irwin 2004) as described in § 4. Third, even when individual EOS components yield excellent thermodynamic consistency, the necessity of blending between EOS components to provide continuous coverage across different regimes inevitably introduces additional thermodynamic inconsistency. We have mitigated this last issue by minimizing the number and severity of EOS blends as much as possible, but unavoidable energy inconsistencies remain, and we discuss their implications for energy conservation in § 5.3.2.', '5.3.1. EOS Derivatives': 'In MESA V, we addressed the quality of the EOS derivatives by introducing new options that used bicubic spline interpolation in high-resolution tables of P gas , s , and e . This provided accurate first and second partial derivatives by evaluating analytic derivatives of the interpolating polynomials rather than by interpolating values of tabulated derivatives. While this approach successfully ensured that the partial derivatives corresponded to how the interpolated EOS values actually changed in response to small changes of the parameters, it inevitably led to small, interpolation-related artifacts in partial derivative quantities such as ∇ ad or Γ 1 . In asteroseismic applications that require smooth profiles of the Brunt-Vaisala frequency, this approach proved unsatisfactory.  \nMESA now adopts an approach that separately treats quantities that appear in the equations (and happen to be partial derivatives) and the places where these theoretically equivalent, but numerically different quantities appear in the Jacobian (as partial derivatives of other quantities that appear in the equations). That is, the Jacobian uses the partial derivatives of bicubic spline interpolants, while the equations use the bicubic spline interpolants of partial derivatives. This enables both efficient numerics and smoother solutions at the cost of some additional bookkeeping. A potential pitfall is that negative values for non-negative quantities can be encountered. In practice, we find that we do not encounter negative interpolants for the physical quantities that enter the equations. While we may encounter negative values from the derivatives of the interpolants used for the Jacobian, these only guide the Newton iterations in converging toward a solution. In this scheme, negative derivatives of interpolants cannot introduce physical errors into the equations used for model solutions.', '5.3.2. Thermodynamic Consistency and EOS Blends': 'In order to quantify how models employing the different forms of the energy equation experience inconsistencies in the EOS differently, we establish a measure of the quality of the MESA EOS during the evolution of a model as follows. In the ( ρ, T, { X i } ) basis, the total derivative of the specific internal energy, e , mathematically satisfies  \nglyph[negationslash]  \nD e D t -[ ( ∂e ∂ρ ) T, { X i } D ρ D t + ( ∂e ∂T ) ρ, { X i } D T D t + ∑ i ( ∂e ∂X i ) ρ,T, { X = X i } D X i D t ] = 0 . (55)  \nFor a Lagrangian volume corresponding to a cell k , we evaluate the time integral of the left hand side of Equation (55) across a time step. We replace the sum over individual composition derivatives with a single directional derivative along the direction of composition change over the time step, as in the evaluation of glyph[epsilon1] grav ,X i . We approximate terms that are not the integral of total derivatives using the implicit trapezoidal rule. Weighting by dm k and summing over all cells this gives, for a single step,  \nE step err , eos = n ∑ k =1 dm k { ∆ e -[ ρ ( ∂e ∂ρ ) T ] ∆ln ρ -[ c V T ]∆ln T -∆ e X i } k , (56)  \nwhere overline quantities are the trapezoidal rule estimates corresponding to the average at the start and end of the time step. The value ∆ e X i = e ( { X i } ) -e ( { X i } start ) is the change in specific internal energy due to composition changes alone at a given ρ and T , so  \n∆ e X i = 1 2 [ e ( ρ start , T start , { X i } ) -e ( ρ start , T start , { X i } start ) + e ( ρ, T, { X i } ) -e ( ρ, T, { X i } start ) ] . (57)  \nSumming the per-step errors over a MESA run, then at the n -th time step, we have  \nE run err , eos = n ∑ i =1 E step ,i err , eos . (58)  \nAs the time step is reduced, the error from the temporal discretization shrinks and E run err , eos converges to a measure of the energy error incurred as a result of EOS shortcomings.  \nWe define two other energy errors. By limiting the sum in Equation (56) to those zones that are in an EOS blend during a particular time step, we can isolate the per-step energy error due to the blend and so analogously define E step err , blend and E run err , blend . We define the perstep residual energy error E step err , res as the mass-weighted sum of the energy equation residuals over the model and also track its cumulative value E run err , res .  \nTo understand how these different forms of error might manifest under different treatments of the energy equation, it is helpful to consider the idealized case of an EOS that is a blend of EOS1 and EOS2, identical except for their definition of where the energy zero point lies, so that e EOS1 = e EOS2 + e offset , where e offset is some constant. Physically, either EOS should produce the same evolution, and all EOS derivatives will be the same. However, evolution through the blend between EOS1 and EOS2 will not satisfy Equation (55), and therefore must lead to non-zero values of E step err , eos regardless of which form of the energy equation is used. In particular, the thermodynamic identity c V = ( ∂e/∂T ) ρ is violated in the blending region because the blending coefficient derivatives are not included as part of the blended value of c V (see § 4.3).  \nIn practice, such energy offsets (in addition to other inconsistencies) always occur at the locations of MESA EOS blends because it is impossible to construct a blending region in which two distinct EOS treatments agree exactly. Due to different input physics assumptions in different EOS components, energy offsets have more complexity than simple constant zero point differences across the range of parameters where blending is necessary. As far as possible, we have chosen blending loca-  \ntions to minimize the differences between EOS components and to minimize the residual amount of unavoidable offset (e.g., by adjusting the definitions of internal energy to be as consistent as possible about where the zero point lies). However, no general solution is currently available to completely eliminate inconsistencies for blends between our current EOS components, and EOS blends therefore remain one of the largest potential sources of energy error when present in MESA models.  \nWhen evolving using the dedt form of the energy equation, Equation (49) will lead to e offset being folded into the D e/ D t term of the energy equation for regions of the model evolving through the blend, injecting spurious heating/cooling into those regions. However, the energy error reported by MESA in these regions may still be 0, because MESA evaluates energy error according to the blended e from the EOS. It is therefore possible to have models for which E step err , eos is significantly larger than E err reported by MESA when using the dedt form of the energy equation, as we shall see in the following subsection. On the other hand, in this idealized scenario the eps grav form of the energy equation (Equations 51 and 52) would give the physically correct evolution since it is evaluated in terms of derivatives that are unaffected by the energy offset in the blending region. However, both E step err , eos and rel run E err will report large values under the eps grav form in this scenario, reflecting the thermodynamic inconsistency of e in the blended EOS rather than inaccuracy in the evolution.', '5.4. Results: 1M pre ms to 1TP': 'Having demonstrated the performance of the various forms of the energy equation in idealized circumstances in § 5.2, we now model the evolution of a 1 M glyph[circledot] star from the pre-main sequence to its first He thermal pulse on the AGB. This example includes composition changes due to nuclear reactions, and it uses the current default MESA EOS (blending together FreeEOS and Skye for the ( ρ, T ) regions encountered by this model). A small portion of the envelope of this model encounters the FreeEOS-OPAL/SCVH blend near ZAMS, then the model lies entirely on FreeEOS for most of the first ≈ 5 Gyr of MS evolution (see § 4 and Figure 7), after which the core evolves toward higher density and encounters the FreeEOS-Skye blend. The approach described in § 5.3 allows us to quantify the various sources of energy errors, identifying how much error comes from EOS inconsistencies and blends, and how much is due to residuals of the equation solutions.  \nFigure 11 summarizes the energy errors defined in § 5.3 when using different forms of the energy equation. Because the models are approximately converged and the  \n8  \nFigure 11. Comparison of types of energy errors for a 1 M glyph[circledot] model. Each panel shows a different form of the energy equation. Each line color indicates a different type of energy error. Each line style indicates a different time resolution. The thick black lines show the error in numerical energy conservation that MESA reports as rel run E err . The thinner colored lines show the errors due to the EOS, EOS blends, and equation residuals. The composition term glyph[epsilon1] grav ,X i is included for both of the eps grav panels in this figure.  \n<!-- image -->  \ndifferent forms of the energy equation converge to the same solution, the colored lines at a given time resolution are similar in all panels. In particular, the measure of the energy error associated with the EOS blend has approximately converged (agreement among orange lines in the middle of each panel). However, the measure of the EOS inconsistency in places where this quantity is not dominated by the blend (i.e., where the blue lines are above the orange lines) continues to shrink with increasing time resolution, showing that this measure is not converged and that the true inconsistency is less than that indicated by the highest resolution line.  \nIn the top panel, the dedt form shows a numerical energy error (black lines) that is roughly the energy error associated with the equation residuals (green lines). The EOS energy errors are present in the model, but  \ndo not show up in the reported rel run E err for this form of the energy equation, as explained in § 5.3. In the middle panel, the end-of-step eps grav form shows that the numerical energy conservation error remains above the energy errors associated with the EOS, while these in turn are generally well above the errors associated with the equation residuals. In the bottom panel, the centered eps grav form shows a numerical energy error (black lines) that is roughly the energy error associated with the EOS inconsistencies. In particular, by the end of the main sequence, the numerical energy error is dominated by the blend error. Further resolution increases or improvements to the energy equation will not improve the numerical energy conservation. Progress can only come through improvements to the EOS.', '6. ATMOSPHERE': "The MESA atm module uses an atmosphere model to evaluate the pressure P s and temperature T s at the outermost ( k = 1) cell boundary. These 'model surface' values are in turn incorporated in the outer boundary conditions applied to the stellar model, as specified in Equation (10) of MESA I. Here we describe various improvements and fixes in the atm module since MESA I and MESA II.  \nIn older MESA releases, the choice of atmosphere model was controlled by the which atm option inlist parameter. In recent releases this parameter is renamed atm option , with three possible choices:  \n- · atm option = ' T tau ' : atmosphere based on T ( τ ) relations, as discussed in § 6.1.\n- · atm option = ' irradiated ' : irradiated atmosphere, as discussed in § 6.2.\n- · atm option = ' table ' : tabulated atmosphere, as discussed in § 6.3.  \nIn addition, MESA now offers two complementary approaches to including the atmosphere structure in model data passed as input into pulsation codes. These are described in § 6.4 and § 6.5.", '6.1. Atmospheric T ( τ ) Relations': "Setting atm option = ' T tau ' builds an atmosphere in which temperature at each optical depth τ is specified by a function T ( τ ); P s and T s are obtained by evaluating the atmospheric thermodynamic state at an optical depth τ s corresponding to the nominal model surface. This optical depth can be much smaller or much greater than the optical depth τ ≈ 2 / 3 typically associated with stellar photospheres; that is, the model surface need not correspond to the photosphere.  \nThe T ( τ ) functions have the form  \nT 4 ( τ ) = 3 4 T 4 eff [ τ + q ( τ )] . (59)  \nSelection of T ( τ ) is set by the atm T tau relation inlist parameter, with four possible choices:  \n- · atm T tau relation = ' Eddington ' : the gray, Eddington-approximation relation, where q ( τ ) = 2 / 3.\n- · atm T tau relation = ' solar Hopf ' : the q ( τ ) function described in Equation (A9) of MESA II, which is a fit to Model C of the solar atmosphere by Vernazza et al. (1981), often referred to as VAL C.\n- · atm T tau relation = ' Krishna Swamy ' : the relation from Equation (33) of Krishna Swamy (1966).\n- · atm T tau relation = ' Trampedach solar ' : the relation from Ball (2021), itself a fit to the solar atmosphere simulation by Trampedach et al. (2014).  \nFor a given T ( τ ) relation, the corresponding P ( τ ) throughout the atmosphere is obtained by integrating the hydrostatic balance equation  \n∂P ∂τ = g κ , (60)  \nfrom τ glyph[lessmuch] 1 inward to τ = τ s . In this integration the gravity is set to the constant value g = GM/R 2 , in accordance with the assumption that the atmosphere is geometrically thin and contains negligible mass. The opacity evaluation is controlled by the atm T tau opacity inlist parameter, with three possible choices:  \n- · atm T tau opacity = ' fixed ' -uniform κ throughout the atmosphere, with a value set by the current opacity κ 1 in the outermost cell.\n- · atm T tau opacity = ' iterated ' - uniform κ throughout the atmosphere, with a value obtained from the kap module for ( P s , T s ). As indicated by its name, this choice requires iteration because P s is not known a priori.\n- · atm T tau opacity = ' varying ' - varying κ throughout the atmosphere, with a value obtained from the kap module for ( P, T ) at the local τ .  \nWith the last choice, the dopri5 (5 th -order DormandPrince) differential equation integrator from the num module is employed with specifiable error tolerance and maximum number of steps. With the first and second choices, however, the fact that κ does not depend on τ means that Equation (60) can be integrated analytically to yield  \nP ( τ ) = g κ τ ( 1 + Pextra factor × κL 6 πGMcτ ) . (61)  \nThe second term in the parentheses arises as a constant of integration, and accounts for non-zero radiation pressure in the limit of small τ . The term is obtained from Equation (20.16) of Cox & Giuli (1968). The parameter Pextra factor depends on the assumed angular dependence of the radiation specific intensity, with the default  \nvalue of unity corresponding to isotropic-outward, and a value of 1 . 5 corresponding to radial-outward. While unphysical, setting Pextra factor > 1 . 5 can sometimes be a useful numerical strategy to improve convergence in models that are close to the Eddington limit. However, caution is warranted as this strategy can produce incorrect stellar radii.  \nFor T ( τ ) atmospheres, Table 5 summarizes the mapping between the which atm option parameter choices supported in older MESA releases, and the combinations of atm T tau relation and atm T tau opacity parameter choices that provide the replacement functionality.  \nIn implementing the changes described here, we uncovered two issues that impacted the accuracy and performance of the atm module in older releases.  \nFirst, calls to the eos and kap modules to evaluate ρ ( P, T ) and κ ( ρ, T ) did not use the same tables and/or configuration options as the interior model, leading to possible inconsistencies at the surface where the atmosphere and the interior join. To fix this problem, we implemented a callback system so that the star module can pass appropriately configured EOS and opacity wrapper routines to the atm module.  \nSecond, in cases where atm T tau opacity = ' iterated ' , the partial derivatives of P s and T s with respect to dependent variables ( L 1 , r 1 , ρ 1 , T 1 ) in the outermost cell were incorrectly evaluated; in some cases, this caused the global Newton solver to converge slowly or not at all. To fix this problem, we implemented the correct expressions, which follow from application of the chain rule to the T ( τ ) relation and the hydrostatic solution (Equation 61).", '6.2. Irradiated Atmospheres': "Setting atm option = ' irradiated ' provides functionality similar to ' T tau ' , but adopting the T ( τ ) relation for an externally irradiated atmosphere given in Equation (6) of Guillot & Havel (2011). Equation (60) is integrated analytically, with opacity evaluation controlled by the atm irradiated opacity inlist parameter; the possible choices ' fixed ' and ' iterated ' behave the same as described in § 6.1. The ' iterated ' case replaces the which atm option = ' gray irradiated ' choice described in MESA II, and fixes a bug related to incorrect T eff evaluation.", '6.3. Tabulated White Dwarf Atmospheres': "When atm option = ' table ' , P s and T s are obtained by interpolating in pre-computed atmosphere tables. In addition to the options described in MESA I and MESA II, MESA now provides a new set of atmosphere tables for WDs with He-dominated surfaces (DB WDs).  \nFigure 12. Depth of the surface convection zone (measured as fractional convection zone mass M cvz /M ) for a 0 . 57 M glyph[circledot] , log( g/ cms -2 ) ≈ 8 . 0 WD with a He atmosphere. The curves indicate MESA results with the DB tables and the Eddington T ( τ ) relation, while the points are based on Koester et al. (2020).  \n<!-- image -->  \nThese tables provide P s and T s data over the ranges 5 , 000 K < T eff < 40 , 000 K and 6 . 0 < log( g/ cms -2 ) < 9 . 4. They assume a He dominated composition of n H = 10 -5 n He , and are based on model atmospheres calculated using the Koester (2010) code and evaluated at τ = 25. The limits of this grid are now explained.  \nAt T eff glyph[greaterorsimilar] 40 , 000 K He undergoes its second ionization, and non-LTE effects that are not included in the atmosphere code become important. The lower limit of T eff = 5 , 000 K is imposed by the Koester (2010) code's EOS. Studies of spectroscopic samples of DBs (e.g., Eisenstein et al. 2006; Genest-Beaulieu & Bergeron 2019) have shown that DBs are well described for surface gravities lower than log( g/ cms -2 ) = 9 . 5. The lower limit of log( g/ cms -2 ) = 6 . 0 is imposed by the convergence of the models at T eff = 40 , 000 K, since radiation pressure competes with gravity in this regime.  \nAs a diagnostic for these new tables, we use the depth of the surface convection zone in cooling DB WD models, which is primarily sensitive to the atmospheric boundary condition and the stellar model EOS in the convective region. Figure 12 demonstrates the improvements brought by these new tables, plotting the mass M cvz of the surface convection zone versus T eff . Also shown are results from the atmosphere and envelope models of Koester et al. (2020), as well as the outcome from using atm option = ' T tau ' with atm T tau relation = ' Eddington ' and atm T tau opacity = ' fixed ' .  \nThe new tables reproduce the Koester et al. (2020) calculations much better at low T eff than the Edding-  \nTable 5. Mapping between which atm option parameter choices in releases of MESA up to r11701 (see § 12.4), and the corresponding parameter choices in releases since r12115 (to be used in tandem with atm option = ' T tau ' ).  \nFigure 13. The hydrostatic term ρ -1 d P/ d r (in units of the surface gravity GM/R 2 ) as a function of r/R for GYRE -format models of the present-day Sun that include reconstructed atmospheres with atm option = ' T tau ' , atm T tau relation = ' Eddington ' and atm T tau opacity = ' varying ' . The two curves show releases r11701 (prior to the reconstruction fix) and r22.05.1; only the latter reproduces the correct behavior ρ -1 d P/ d r = -GM/R 2 throughout the reconstructed atmosphere r/R ≥ 1.  \n<!-- image -->  \nton T ( τ ) atmospheres. The latter diverge for T eff glyph[lessorsimilar] 15 , 000 K because the conditions at the surface cross the boundaries of the He opacity table coverage in MESA . The tabulated atmospheres agree with Koester et al. (2020) to T eff ≈ 7 , 000 K. For cooler temperatures, the uncertainty in M cvz is due to the uncertain EOS at the base of the convection zone (Saumon et al. 1995). 4", '6.4. Atmosphere Reconstruction for Pulsation Codes': 'As discussed in MESA II and MESA III, MESA can pass models to the ADIPLS (Christensen-Dalsgaard 2008) or GYRE (Townsend & Teitler 2013) linear pulsation codes, either in-memory during an astero -module optimization or via files written to disk in a variety of for-  \nmats. Often, it is desirable to reconstruct the atmosphere structure from τ = τ s out to τ glyph[lessmuch] 1 before passing them to the pulsation codes. This has no impact on the interior model, but can improve asteroseismic modeling. For T ( τ ) atmospheres ( § 6.1), setting add atmosphere to pulse data = .true. enables this reconstruction. The radial coordinate r throughout the atmosphere is then determined by integrating the τ equation  \n∂r ∂τ = -1 κρ (62)  \noutward from τ = τ s to τ glyph[lessmuch] 1. For this integration, κ is evaluated in accordance with the atm T tau opacity parameter discussed previously, while ρ ( P, T ) is obtained from the eos module for the local τ .  \nIn releases of MESA prior to r12115, the outward integration used an explicit Euler scheme with a default step-size too large to accurately follow r ( τ ). Together with the eos / kap table inconsistency highlighted in § 6.1, this led to departures from hydrostatic balance in the reconstructed atmosphere, as highlighted in Figure D.1 of Christensen-Dalsgaard et al. (2020) for r/R > 1. To fix this issue, release r12115 and later use the dopri5 integrator for the outward integration, with a specifiable error tolerance, step-size and outermost optical depth. Figure 13 demonstrates these improvements by showing hydrostatic balance for GYRE -format models of the present-day Sun calculated using releases r11701 (prefix) and r22.05.1.', '6.5. Atmospheres as Part of the Interior': "The optically-thin outer layers of a star are usually treated separately from the interiors of stars because they do not satisfy the assumptions under which the interior structure equations are derived. Given a T ( τ ) relation, however, it is possible to correct the equation for radiative heat transport so that the temperature stratification produced by solving the stellar structure equations matches the desired T ( τ ) (Trampedach et al. 2014; Mosumgaard et al. 2018a). The radiative temperature gradient ∇ atm of an atmosphere that follows a given  \n4  \nFigure 14. Ratio of the left- and right-hand sides of Equation (59), by which we compare the temperature structure of two models in FGONG format, one with the atmosphere reconstructed as in § 6.4 (integrated) and the other with the atmosphere modeled with the interior as in § 6.5 (interior).  \n<!-- image -->  \nT ( τ ) relation is  \n∇ atm = ∇ rad ( 1 + d q d τ ) . (63)  \nWe can therefore recover any T ( τ ) relation by scaling the canonical ∇ rad by 1 + d q/ d τ . This procedure is enabled with the new option use T tau gradr factor . For the gray, Eddington-approximation relation, q ( τ ) is constant and no correction is necessary.  \nTo include the optically thin layers in the interior model, the surface boundary should be set or relaxed to a smaller τ using the flags set tau factor or relax tau factor and their associated controls. The surface boundary conditions, now at smaller τ , still follow § 6.1. They are evaluated at τ s = tau factor × τ eff rather than τ eff , where τ eff is the optical depth at which the T ( τ ) relation is equal to T eff .  \nThis approach has several advantages. First, like atmospheres reconstructed with atm T tau opacity = ' varying ' , the atmospheric structure is kept consistent with MESA 's EOS and opacity routines. Second, the computation can leverage the parallelization of MESA . Finally, this approach accounts for the fact that g is not exactly constant throughout the atmosphere ( § 6.4), which is assumed by T ( τ ) atmospheres.  \nFigure 14 shows the ratio of the left- and right-hand sides of Equation (59) in a 1M glyph[circledot] ZAMS model with the ' solar Hopf ' T ( τ ) relation when the atmosphere is either reconstructed as in § 6.4 or included as part of the interior. Both models deviate at τ glyph[greaterorsimilar] 2 because convection starts to transport heat.  \nThe reconstructed atmosphere deviates in τ eff < τ glyph[lessorsimilar] 2 because d q/ d τ = 0 just below the photosphere, which is  \nglyph[negationslash]  \ntreated as part of the interior without correcting ∇ rad . When using analytic T ( τ ) relations, this can be corrected by using the use T tau gradr factor option but not when using tabulated atmospheres, for which the correction factors 1 + d q/ d τ have not been provided. Though inconsequential for tables at τ s glyph[greatermuch] 1, where d q/ d τ → 0 anyway, it introduces inaccuracy in the temperature stratification when using tables at small τ . This inaccuracy may be acceptable, depending on the scientific problem. This inaccuracy is not present if d q/ d τ = 0 in the affected regions, which includes the Eddington T ( τ ) relation.", '6.6. Choosing atmosphere options': "There is limited consensus on when to use which atmosphere options, but we nevertheless offer a few guiding remarks.  \nThe most commonly used and current default in MESA is a gray Eddington atmosphere with the surface boundary at τ s = 2 / 3, with κ fixed throughout the atmosphere. If the precise behavior of the stellar atmosphere is not important, this should suffice.  \nThere is a hierarchy of accuracy at the expense of greater computational cost among the choices for atm T tau opacity . On the basis of the self-consistency of κ , ' varying ' is more accurate than ' iterated ' , which is in turn more accurate than ' fixed ' .  \nFor the calculation of pulsations that have significant amplitude near the surface (e.g., solar p modes), it is important to choose an option that allows the atmosphere to be reconstructed for the equilibrium stellar model ( § 6.4).  \nTabulated atmospheres provide boundary conditions that are typically computed using more complete physics (e.g., non-LTE) than can be described by the stellar structure equations. As shown in Figure 14, tabulated atmospheres at τ glyph[lessorsimilar] 2 lead to a stellar model in which the near-surface temperature stratification is equivalent to an Eddington atmosphere. The models in Figure 14 differ by ∼ 50 K at τ eff . This inaccuracy might be outweighed by the benefits of a complete atmosphere model. The correct structure could in principle be recovered by extracting the appropriate T ( τ ) relations (Trampedach et al. 2014) from the detailed atmosphere models, but these are not generally available.  \nIdeally, we would have access to grids of T ( τ ) relations and corresponding bolometric corrections extracted from advanced simulations of stellar atmospheres, with parameters that cover the HR diagram. Until this ideal is realized, stellar modellers must decide which aspects of the atmospheric boundary condition  \nare most important for their calculations and choose appropriate options.", '7.1. Starspots': 'Starspots are common for stars with M glyph[lessorsimilar] M glyph[circledot] . Models of M dwarfs that include starspots and surface magnetism have inflated radii close to those inferred by observation (e.g., Feiden & Chaboyer 2013; Mann et al. 2015).  \nWe thus implement in MESA the treatment of starspots introduced in the SPOTS models of Somers et al. (2020) (also see Cao & Pinsonneault 2022), which are based on the Yale Rotating Stellar Evolution Code ( YREC ; Demarque et al. 2008; Spada et al. 2013) and described in detail in Somers & Pinsonneault (2015) and Somers et al. (2020).', '7.1.1. Starspots formalism': "Somers et al. (2020) parameterize the variance of the surface flux due to magnetic pressure from starspots by modifying the atmospheric boundary condition. Somers & Pinsonneault (2015) characterized the degree of 'spottiness' on the stellar surface by two parameters:  \n- · SPOTF (hereafter f spot ), a coverage fraction, or 'spot filling factor'\n- · SPOTX (hereafter x spot ), representing the temperature contrast between the spotted and unspotted regions at r = R : x spot = T spot /T photosphere .  \nNumerically, values from 0 . 0 to 1 . 0 are permitted for both parameters. Observationally constrained coverage fractions f spot are described in Cao et al. (2022), who find that a value f spot = 0 . 34 is a reasonable fit to observations of sub-solar-mass stars in the λ Ori cluster.  \nThe spot-induced temperature contrast, x spot , is restricted to physically meaningful values of 0 . 5 -1 . 0. A value of x spot = 1 . 0 indicates that the effective temperature in the spotted region, T spot , does not differ from the surrounding, unspotted effective temperature, T amb ( r ) (the 'ambient temperature,' or ATEFF in YREC ). At the surface, T amb ( r = R ) = T photosphere . A value of x spot = 0 . 5, on the other hand, corresponds to the statement that T spot differs from T amb by the maximum degree permitted by magnetic equipartition: namely, when the magnetic pressure contribution constitutes half of the total pressure.  \nThe temperature contrast x spot perturbs the radiative gradient, ∇ rad, spot , in the surface convection zone. This effect can be made depth-dependent via  \nx spot ( r ) = 1 -(1 -x spot ) T amb ( r ) T ( r ) , (64)  \nwhere T ( r ) is the temperature at r and the quantity T photosphere -T spot = T amb ( R ) -T spot is held constant as a function of r . In Equation (64), x spot is a scalar parameter.  \nPer equations (1) through (4) in Somers & Pinsonneault (2015), the 'spotted' luminosity is set to  \nL spotted ≡ L/ [ f spot x 4 spot +(1 -f spot ) ] . (65)  \nYREC models the suppression of convective flux at the stellar surface via a two-part lookup in its atmospheric boundary tables, invoking 'ambient' versus 'spotted' effective temperatures that differ by a factor of  \nα spot = 1 + f spot ( x 4 spot -1) , (66)  \nwhere x spot can optionally be a function of the local temperature at a given depth within the convective envelope.", '7.1.2. MESA implementation': "MESA 's implementation is equivalent, but modifies the surface boundary conditions in terms of pressure rather than T eff , as T eff is strictly an output in MESA . Through the specification of the magnetic pressure term, P B , the temperature contrast is given by  \nx spot ( r ) = P ( r ) -P B P ( r ) , (67)  \nwhere P B is assigned at the start of the step via  \nP B = ( k B N A ρ/µ )(1 -x spot ( R )) T eff . (68)  \nChoosing x spot = 1 corresponds to P B = 0, hence there is no magnetic pressure and no perturbation to P . Choosing x spot = 0 yields the other extreme: P B = P . Once again, a practical limit on P B is set by the assignment x spot = 0 . 5, which corresponds to magnetic equipartition.  \nUsing auto diff ( § 2) we account for the modified pressure term and its partial derivatives at every mass shell. In Equations (67) and (68), P and P B are auto diff variables. Likewise, the radiative temperature gradient due to the presence of spots (i.e., magnetic inhibition of the convective flux),  \n∇ rad, spot = ∇ rad / ( f spot x spot ( r ) 4 +1 -f spot ) , (69)  \nis an auto diff variable, since it is a function of the auto diff quantities ∇ rad , P , x spot ( r ), and scalars.  \nThe quantity P B is assigned once at the beginning of every evolutionary step and held constant over all Newton iterations within that step. The use of auto diff data types ensures that the Newton solver automatically  \nFigure 15. Evolutionary tracks showing the effects of star spots for x spot = 0 . 85. Color indicates mass. Within clusters of identical mass, f spot values are ordered lowest (left-most) to highest (right-most).  \n<!-- image -->  \nreceives correct partial derivatives of the modified radiative temperature gradient with respect to, e.g., depth and other stellar structure variables. The modification to ∇ rad is evaluated once per Newton iteration. To obtain the spotted luminosity, L is adjusted by a factor of α spot (Equation 66).", '7.1.3. MESA Models': 'We demonstrate the starspots functionality by computing a grid of spotted evolutionary tracks for M = 0 . 2 -1 . 3 M glyph[circledot] and Z = 0 . 014. We use the photosphere table option for atmospheric boundary conditions across all tracks, though for the lowest-mass stars (e.g., M glyph[lessorsimilar] 0 . 5 M glyph[circledot] ), a choice of tau 1m1 would be more appropriate. We use the Henyey MLT prescription with α MLT = 1 . 95.  \nFigure 15 shows tracks using f spot = { 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8 } and x spot = 0 . 85. When comparing MESA models to the YREC -based SPOTS models of Somers et al. (2020), we observe the following features:  \n- · Tracks computed with MESA starspots shift smoothly by the same magnitude and in the same direction as a function of f spot .\n- · We also find that the lower the initial value of x spot (i.e., more extreme in terms of magnetism), the greater the impact of an increased coverage fraction.\n- · Below ≈ 0 . 6 M glyph[circledot] , the SPOTS models are cooler than the MESA starspots models. This is due to differences in our choice of the atmospheric boundary condition.  \n- · The degree of radial inflation for an 0 . 3 M glyph[circledot] M dwarf predicted by MESA starspots , using f spot = 0 . 34 and x spot = 0 . 85, is of the order 3%-5%, in agreement with the SPOTS models.  \nThe MESA starspots test case can be found in MESA\\_DIR/star/test\\_suite/starspots .', '7.2. Superadiabatic Convection in Radiation-dominated Regions': "Modeling stars near the Eddington limit is a complex numerical problem. Under such conditions, convective regions with density inversions are expected (Joss et al. 1973; Paxton et al. 2013) and 1D models using standard MLT can develop extended low-density envelopes, becoming red supergiants before finishing their MSevolution (e.g., Sanyal et al. 2015; Sz'ecsi et al. 2015). Three-dimensional radiative hydrodynamic simulations are just starting to explore the physics of energy transport near the Eddington limit (Jiang et al. 2015; Tsang &Milosavljevi'c 2015; Schultz et al. 2020; Goldberg et al. 2022; Moens et al. 2022), and will hopefully provide a way to accurately model these regions in 1D stellar evolution instruments. In the meantime, 1D simulations using MLT near the Eddington limit are expensive, requiring small timesteps down to the point that some calculations become impractical. Enhanced convective energy transport in these regions can inhibit the formation of density inversions and facilitate calculations. One such approach is the use of a density scale height rather than a pressure scale height in MLT (Nishida & Schindler 1967; Stothers & Chin 1973; Maeder 1987).  \nThe MLT++ formalism is a stellar-engineering approach that has been commonly used in MESA to reduce superadiabaticity in regions nearing the Eddington limit (MESA III). Although convection is expected to operate in regions of the star approaching the Eddington limit, the efficiency of convective energy transport is uncertain. By reducing the expected superadiabaticity, MLT++ provides an ad-hoc enhancement. Such an enhancement is supported by results from 3D simulations (Jiang et al. 2015; Schultz et al. 2020), but the method is not calibrated to detailed simulations or observations. The main motivation for MLT++ is to enable computations of massive star evolution up to late stages, and users need to assess if the deviations from a more physical model such as MLT are relevant to their results. One important limitation of MLT++ is that it is a non-local explicit method, which can lead to large step-to-step variations that produce unphysical results and prevent the solver from finding a valid solution. By making use of auto diff ( § 2), we have implemented a fully implicit and local alternative to MLT++, which allows the mod-  \nlling of a larger range of masses and metallicities. Setting the option use superad reduction activates this method.  \nIn hydrostatic equilibrium, the ratio between the radiative luminosity and the local Eddington luminosity is  \nΓ Edd ≡ L rad L edd = 4 aT 4 3 P ∇ . (70)  \nFor a particular model of energy transfer (e.g., TDC as discussed in § 3, or MLT), the expected Eddington factor is  \nΓ Edd , exp ≡ 4 aT 4 3 P ∇ exp , (71)  \nwhere ∇ exp is the temperature gradient predicted by the energy transfer model. As in MLT++, we artificially enhance energy transport in convective regions where the expected Eddington factor is high by adjusting ∇ rad . The difference between the radiative and the Ledoux gradient is reduced to  \n∇ rad , new -∇ L = ∇ rad -∇ L f Γ , (72)  \nwhere ∇ rad , new is the adjusted radiative temperature gradient and f Γ ≥ 1 determines the reduction of ∇ rad -∇ L in the convective region. Such a scaling of ∇ rad can be interpreted as an effective lowering of κ in regions near the Eddington limit. Results from Schultz et al. (2020) suggest that the impact on radiative transfer of a vigorously convecting region supports this choice. The adjusted ∇ rad , new is then used instead of ∇ rad to recompute ∇ according to the convection model.  \nThe functional form of f Γ is arbitrary, and was determined empirically so that Equation (72) can be applied in a large number of cases while minimizing adjustments in the limit Γ Edd , exp → 0. Just as with MLT++, it serves as a stellar-engineering method to circumvent complex evolutionary stages, rather than a specific physical model that accounts for how convection is modified near the Eddington limit. While exploring different options for f Γ , auto diff played a critical role by removing the need to directly specify partial derivatives.  \nOur choice for f Γ is the following:  \nf Γ = 1 + α 1 g (Γ Edd , exp / Γ c -1) + α 2 g (Γ exp / Γ inv -1) √ β × h (( ∇ exp -∇ L ) /δ c ) , (73)  \nwhere β is the ratio of gas to total pressure and Γ inv ≡ 4(1 -β ) / (4 -3 β ) is the Eddington factor at which an  \nFigure 16. Evolution of stars ranging from 10 to 10 6 M glyph[circledot] at solar metallicity, using the implicit method to increase the efficiency of energy transport in regions approaching the Eddington limit.  \n<!-- image -->  \nideal gas with radiation develops a density inversion (Joss et al. 1973, MESA III). The parameters Γ c , α 1 , α 2 and δ c regulate the enhancement of energy transport. The function  \ng ( x ) ≡        0 x < 0 x 2 / 2 0 < x < 1 x -1 / 2 x > 1 (74)  \nis continuous with a continuous first derivative. No correction will be applied if Γ exp < Γ c and Γ exp < Γ inv . If either of those thresholds is exceeded, α 1 and α 2 set the strength of the enhancement in energy transport for each. The 1 / √ β term further enhances the effect in regions dominated by radiation pressure. The function h ( x ) is chosen as  \nh ( x ) =        0 x ≤ 0 6 x 5 -15 x 4 +10 x 3 0 < x ≤ 1 1 x > 1 , (75)  \nsuch that it is equal to zero for x < 0 and equal to one for x > 1, while monotonically increasing in between with a zero derivative at x = 0 and x = 1. This choice ensures significant corrections are only applied in cases where a superadiabaticity comparable to δ c would be expected.  \nFigure 16 illustrates the evolution of stars up to 10 6 M glyph[circledot] at high metallicities (Z = 0 . 0142) using this new approach. The very high mass models are not necessarily meant to represent real stars, but serve as an extreme test of this new approach. In particular, realistic  \nFigure 17. Evolution of a 30 M glyph[circledot] star at solar metallicity. The different tracks use the implicit method to enhance energy transport near the Eddington limit (with different choices for α 1 and α 2 ), MLT++, and the MESA default which includes no energy-transport enhancement. Both simulations with the new method use Γ c = 0 . 5 and δ c = 10 -2 . The MESA default simulation stalls when evolving towards the blue after a RSG phase. Inset plots are made to highlight variations between the runs at TAMS and at the end of the RSG phase.  \n<!-- image -->  \nmodels of supermassive stars need to also take into account general relativistic effects (Chandrasekhar 1964; Fricke 1973), which are not included in these simulations. The calculations all used Γ c = 0 . 5, α 1 = α 2 = 5 and δ c = 10 -2 , which we found to perform consistently across a broad range of masses.  \nA comparison between the new method, MLT++ and a simulation without any enhancement of energy transport is shown for a 30 M glyph[circledot] model in Figure 17. Overall, the new method provides smoother evolution while remaining closer to the result obtained without enhancing energy transport. MLT++ introduces undesirable numerical variations that are particularly visible when the model moves from the blue to the red at TAMS, and when it evolves from the red to the blue after being stripped of most of its H envelope. Lowering α 1 or α 2 , or increasing Γ c or δ c , will produce results closer to those without energy-transport enhancement. Figure 17 also shows two different simulations with the new method using α 1 = α 2 = 5 and 2. The simulations performed with α 1 = α 2 = 2 take almost triple the number of steps and computation wall time. Selection of these parameters requires balancing performance and similarity to the unenhanced behavior. Users need to carefully as-  \nFigure 18. Evolution of a 40 M glyph[circledot] star with a 30 M glyph[circledot] point mass companion with an initial orbital period of 50 days. Models are computed for solar metallicity and using different methods to enhance energy transport in regions near the Eddington limit.  \n<!-- image -->  \nFigure 18 shows the evolution of a 40 M glyph[circledot] model with a 30 M glyph[circledot] point mass companion at an initial orbital pe-  \n<!-- image -->  \n/circledot  \nFigure 19. Mass transfer rate ˙ M transfer as a function of donor mass for the binary system shown in Figure 18. All simulations have two significant phases of mass transfer corresponding to interaction during the MS (Case A) and after TAMS (Case AB).  \nsess if such variations have a meaningful impact in their conclusions.  \nThe stability of this implicit method is particularly useful in simulations of binary systems, where small step-to-step variations in R can result in large changes to ˙ M during Roche lobe overflow. This makes the previous MLT++ method inappropriate.  \nriod of 50 days. Simulations are performed using no enhancement of energy transfer, MLT++, and the implicit method. In all cases, the models experience an initial mass transfer phase during the MS, and a second mass transfer phase right after TAMS. The simulation without enhanced energy transfer stalls during this second mass transfer phase, exemplifying the computational complexity of modelling stars near the Eddington limit.  \nFigure 19 demonstrates that although the MLT++ model evolves beyond detachment, large step-to-step variations lead to large changes in mass transfer rates. In contrast, the new implicit method provides a smooth solution throughout the mass transfer phase. An early example of this implicit method enabling the succesful computation of massive binary evolution are the models of ζ Ophiuchi computed by Renzo & Gotberg (2021).", '8.1. Molecular Opacities': "Molecules contribute significantly to stellar opacity for T glyph[lessorsimilar] 5000K (Alexander & Ferguson 1994; Ferguson et al. 2005). The MESA kap module has been expanded to include low-temperature molecular opacities from Lederer & Aringer (2009) and ÆSOPUS (Marigo & Aringer 2009). Both opacity sets allow for varying CNO levels. CNO-enhanced molecular opacities find applications in models of red giants, AGB stars, and R Coronae Borealis stars (Schwab 2019).  \nLederer & Aringer (2009) provide Rosseland mean opacity tables for 14 metallicities ranging from Z = 10 -5 to Z = 4 × 10 -2 , and three H mass fractions X = 0.5, 0.7, or 0.8. One can specify seven Z-dependent levels of C-enhancement and three Z-dependent levels of N-enhancement, all with Lodders (2003) solar abundances. The tables span 3 . 2 ≤ log( T/ K) ≤ 4 . 05 and -7 ≤ log( ρT -3 6 / g cm -3 ) ≤ 1. These opacities are enabled by setting kappa low T prefix = ' kapCN ' .  \nÆSOPUS (Marigo & Aringer 2009) allows computation of opacity tables for a variety of solar compositions with the optional inclusion of enhancements (and depletions) to C, N and C/O on top of the basic mixture. We provide a set of ÆSOPUS opacity tables. Additional ÆSOPUS opacity tables can be generated through the ÆSOPUS web-interface 5 , and MESA DIR/kap/preprocessor/AESOPUS contains information on preparing the tables for MESA .  \nTo compare MESA and the Monash stellar evolution code (Lattanzio 1986; Frost & Lattanzio 1996; Karakas  \n& Lattanzio 2007), custom ÆSOPUS tables were generated with 3 . 2 ≤ log( T/ K) ≤ 4 . 5 in steps of 0.01 dex and -7 ≤ log( ρT -3 6 / g cm -3 ) ≤ 1 in steps of 0.05 dex. We use steps of 0.5 dex at higher T , the Lodders (2003) solar composition, reference metallicities of Z = 0.01 to 0.10 in steps of 0.01 dex and 0 . 5 ≤ X ≤ 0 . 8 in steps of 0.1. We set the CNO abundance variation factors fc = -1, 0.2, 0.4, 0.6, 1, 1.5, fco = -1, -0 . 5, 0, 0.5, 1, 1.5 and fn = 0, 0.4, 0.7, 1.  \nThese tables were installed in MESA and the Monash code. The physics used in the Monash code is reconstructed as closely as possible in MESA (Cinquegrana et al. 2022; Cinquegrana & Joyce 2022). This includes the basic.net reaction network, treatment of MLT convection (with independently calibrated α MLT parameters of α MLT =1.931 in MESA , α MLT =1.86 in Monash), high temperature and molecular opacities (Iglesias & Rogers 1996; Marigo & Aringer 2009), an Eddington gray atmosphere, and mass-loss approximations (Reimers 1975; Blocker 1995). We use an η Reimers = 0 . 477 efficiency factor, and η Blocker =0.01 in MESA and η Blocker = 0 . 02 in the Monash code (Cinquegrana et al. 2022). Functionality that is not available in both includes the process of defining the border between convective and radiative regions. Here, a relaxation method is used in the Monash code (Lattanzio 1986), and the predictive mixing algorithm in MESA (MESA IV).  \nWe evolve a 3 M glyph[circledot] , Z = 0.014 model with MESA and the Monash code. Both models use the same initial conditions and are terminated at the 11th thermal pulse. Figure 20 compares the L , T eff , and R histories of the two models during the TP-AGB phase. The MESA model reaches higher peak L , smaller peak T eff , and larger peak R than the Monash model. The differences decrease as the evolution proceeds. Both have very similar interpulse periods, glyph[similarequal] 7.2 × 10 4 yr, after the first few pulses.", '8.2. Compton Opacities': 'At temperatures log( T/ K) glyph[greaterorsimilar] 8 . 7, beyond the boundary of the standard opacity tables, the radiative opacity is set by Compton scattering (MESA I, § 4.3). Such conditions are often realized in simulations of massive stars and thermonuclear bursts on neutron stars. MESA now calculates the Compton scattering opacity using the prescription of Poutanen (2017), which improves upon the previous approach from Buchler & Yueh (1976). See § 4 in Poutanen (2017) for a detailed comparison of the two approaches.', '8.3. Conductive Opacities': 'Energy transport via electron conduction plays an important role, especially in degenerate stellar interiors.  \nFigure 20. Evolution of L , T eff and R , respectively, during the TP-AGB phase of 3M glyph[circledot] , Z=0.014 models from MESA (solid) and the Monash code (dashed). The age of the MESA (Monash) stellar model is shown on the bottom (top) x-axis.  \n<!-- image -->  \nThe conductive opacities in MESA are expanded versions of the tables from Cassisi et al. (2007), see MESA II, Appendix A.3. Blouin et al. (2020) improved the conductive opacities for H and He in the regime of moderate Coulomb coupling and moderate degeneracy, primarily relevant for the envelopes of DA and DB WDs. We implement their analytical fits for conductive opacities for H and He in MESA .  \nFigure 21 shows the effect of these new opacities on MESA WD cooling calculations for a 0 . 9 M glyph[circledot] WD with a H atmosphere, confirming the result of Blouin et al.  \nFigure 21. Cooling tracks for 0 . 9 M glyph[circledot] DA WD models using the conductive opacities of Cassisi et al. (2007) and Blouin et al. (2020).  \n<!-- image -->  \n(2020) that the cooling time to reach T eff = 4000K can be reduced by ≈ 2 Gyr. The difference in cooling times is somewhat smaller for WDs with He atmospheres, as also confirmed in MESA models for AM CVn accretors (Wong & Bildsten 2021). Salaris et al. (2022) provide a more detailed study of the impact of these updated opacities on WD cooling timescales. The conductive opacities including the Blouin et al. (2020) corrections are now the default in MESA .', '8.4. Radiative Accelerations from the OP Monochromatic Opacities': "The radiative acceleration g rad ,i of a given species i represents the acceleration caused by the radiation field. MESA III ( § 9) describes the inclusion of radiative levitation from Hu et al. (2011), which includes evaluation of the Rosseland mean opacity κ R and g rad ,i using the Opacity Project (OP) monochromatic opacity tables (Seaton 2005). The computation of both of these quantities requires wavelength integrations of monochromatic opacities according to the local mixture, introducing a significant amount of extra time. MESA V ( § 6.2) briefly describes the steps taken to reduce the time required to evaluate κ R and g rad ,i . We now describe a new implementation from Mombarg et al. (2022).  \nThe OP tables provide monochromatic cross sections σ i , equally spaced in  \nv ( u ) = 15 4 π 4 ∫ u 0 x 4 e -x (1 -e -x ) 3 d x, (76)  \nwith a step size ∆ v , where 10 -3 ≤ u ≡ hν/k B T ≤ 20. The OP tables contain data for H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni. The Rosseland mean opacity for cell k is then given by  \nκ R ,k = 1 µ k   ∑ n 1 ∑ i f i,k σ i ( v n ) ∆ v   -1 , (77)  \nwhere µ k is the mean molecular weight and f i the fractional element abundance of species i . The monochromatic opacities are sampled at 10,000 points in the frequency parameter v n , denoted by index n .  \nIn the new implementation, the OP monochromatic data are converted from ( T, n e ) space to ( T, ρ ) space with  \nlog( ρ/ g cm -3 ) = log( n e / cm -3 ) + log µ -log Ξ -log( N A / g) , (78)  \nwhere Ξ = ∑ i f i Ξ i is the average number of electrons per atom (given by the OP tables) and N A is Avogadro's constant. The tables contain data for 3.5 ≤ log( T/ K) ≤ 8.0 and log( ρ/ g cm -3 ) bounds that depend on log( T/ K). An error is returned for T and ρ beyond these limits.  \nWe select the point in the OP tables with the smallest value of  \n√ (log T OP i -log T k ) 2 0 . 0025 + (log ρ OP i -log ρ k ) 2 0 . 25 , (79)  \nwhere the different denominators reflect the different OP table spacings.  \nFrom this OP data point we select nearest neighbor points and construct a spline interpolant with a maximum degree of three, i.e., a bicubic interpolant. The interpolants are stored such that the opacity can be re-interpolated at the next time step as long as log( T/ K) has changed by less than twice the grid spacing (0.01 dex), log( ρ/ g cm -3 ) has changed less than twice its grid spacing (0.1 dex), and the fractional abundances all satisfy  \n| f i,t -f i,t -1 | f i,t -1 < glyph[epsilon1], (80)  \nwhere t and t -1 indicate the current and previous time steps, respectively. Empirically, we find for MS models that glyph[epsilon1] =10 -4 optimizes the computational efficiency  \nFigure 22. Rosseland mean opacity ( κ R , upper) and the radiative accelerations ( g rad , lower) of 16 O and 40 Ar in a 2 M glyph[circledot] model with an initial Z = 0.02 at a core H mass fraction of X c =0.2. Green dotted curves show the MESA III (Hu) treatment, and black solid curves show the current treatment (Mombarg). Shaded areas indicate convective regions.  \n<!-- image -->  \nwithout sacrificing accuracy in κ R and g rad ,i (Hui-BonHoa 2021).  \nThe radiative acceleration of a species is given by  \ng rad ,i = µκ R µ i c F γ i , (81)  \nwhere µ i is the molecular weight of the species and F is the radiative flux. The factor γ i is given by  \nγ i = ∑ n σ i ( v n )[1 -e -u ( v n ) ] -a i ( v n ) ∑ j f j σ j ( v n ) ∆ v, (82)  \nwhere the index j runs over all species. The correction terms a i ( v n ) are provided by the OP data. The numerator in Equation (82) is precomputed and stored.  \nIn the new implementation, optimized for models with a convective core and radiative envelope, the envelope is divided into two equal regions based on the number of cells (regardless of the outer convection cells). For each region, the γ i factors are pre-computed using an average  \nmixture of the region  \n〈 f i 〉 = ( k 2 ∑ k = k 1 f i,k )( ∑ i k 2 ∑ k = k 1 f i,k ) -1 , (83)  \nwhere k 1 and k 2 are the first and last cell index of a region. Near the boundary of these two regions, 15 cells on each side, the value of log g rad ,i is blended  \nlog g rad ,i = α i log g rad ,i, 1 +(1 -α i ) log g rad ,i, 2 , (84)  \nwhere log g rad ,i, 1 is computed using the average mixture 〈 f i, 1 〉 of the first region, and log g rad ,i, 2 from the average mixture of the second region. The parameter α i varies linearly between 0 and 1 as a function of cell number in the region where blending takes place.  \nFigure 22 compares κ R and g rad ,i from MESA III and the present method. The cell-averaged difference and maximum difference ( | log g Hu rad -log g Mombarg rad | / log g Hu rad ) for 16 O are 0.007 and 0.042, respectively, and for 40 Ar are 0.02 and 0.03, respectively. The differences in g rad for 16 O are typical for most isotopes. The new implementation is about 5 times faster for MS stars.  \nEnabling this capability requires a data file containing the manipulated OP monochromatic data, doi: https://doi.org/10.5281/zenodo.6761478.", '9. DIFFUSION COEFFICIENTS': "MESA III and MESA IV describe the implementation of element diffusion in MESA using the Burgers (1969) equations with diffusion coefficients based on Paquette et al. (1986), and updated with the coefficients of Stanton & Murillo (2016). For strong plasma coupling in the liquid regime (where the multi-component plasma coupling parameter Γ MCP ≡ q 2 e 〈 Z 5 / 3 〉 /a e k B T glyph[greaterorsimilar] 1), these coefficients disagree by a factor of a few with more accurate molecular dynamics (MD) methods (MESA III, Bauer et al. 2020, Caplan et al. 2022).  \nUntil recently, it was unclear how to generalize these results to the arbitrary plasma mixtures needed for stellar models. However, Caplan et al. (2022) show that in the liquid regime (10 glyph[lessorsimilar] Γ MCP glyph[lessorsimilar] 200), diffusion coefficients are approximated to within ∼ 10% accuracy by using a fit to an equivalent one-component plasma (OCP) coefficient calculated with MD, and then scaling with charge as  \nD j = ( Z j 〈 Z 〉 ) -0 . 6 D OCP , (85)  \nwhere D j is the diffusion coefficient of species j , Z j is the charge of species j , 〈 Z 〉 ≡ ∑ j n j Z j /n ions is the average charge of ions in the liquid, and D OCP is the diffusion coefficient for an equivalent OCP with coupling  \nand screening set by the mixture averages. This scaling is justified by the fact that in the liquid regime, ion diffusion can be described in terms of spheres experiencing Stokes-Einstein drift through a viscous liquid (Bildsten & Hall 2001; Daligault 2006), with the effective radii of different ions set by their charge relative to the background plasma (Bauer et al. 2020; Caplan et al. 2022). The D OCP term therefore captures the physics of the global viscosity experienced by all ions, while the scaling with charge captures the different effective radii of ions experiencing Coulomb interactions with the background plasma.  \nSince the diffusion solver in MESA is cast in terms of the Burgers (1969) equations, coefficients must be recast in terms of the binary resistance coefficients K ij between species i and j rather than the net diffusion coefficients D j . At strong plasma coupling (Γ MCP glyph[greaterorsimilar] 10), these coefficients must satisfy the relation 6  \nD j = n j k B T ∑ i K ij . (86)  \nThe resistance coefficients K ij must also be symmetric ( K ij = K ji ), motivating  \nK ij = n i n j ( Z i Z j ) 0 . 6 k B T n ions 〈 Z 0 . 6 〉〈 Z 〉 0 . 6 D OCP , (87)  \nwhere 〈 Z 0 . 6 〉 ≡ ∑ j n j Z 0 . 6 j /n ions . The resistance coefficients of Equation (87) reduce to the desired net diffusion coefficients of Equation (85) under the summation required in Equation (86).  \nWe calculate the OCP diffusion coefficient D OCP using the Caplan et al. (2022) fit to high-resolution MD, which is given in terms of Γ MCP and electron screening length λ e . To calculate the value of λ e as input for the diffusion coefficients, we follow the method described in MESA IV based on Stanton & Murillo (2016) for non-relativistic electrons at ρ < 10 6 g cm -3 . Electrons become relativistic at higher densities, and we therefore switch to the relativistic screening length λ e = ( 2 k F √ α/π ) -1 for densities ρ > 10 6 g cm -3 , where α is the fine-structure constant and k F = (3 π 2 n e ) 1 / 3 .  \nFigure 23 shows the 22 Ne diffusion coefficient in a MESA model of a 0 . 6 M glyph[circledot] C/O WD, which has historically been a key source of uncertainty in the rate at which 22 Ne settles toward the center of the WD (Bildsten & Hall 2001; Deloye & Bildsten 2002; Garc'ıa-Berro et al. 2008). Such  \nFigure 23 also shows the MD coefficients from Caplan et al. (2022) (these closely match the MD coefficients of Hughto et al. 2010) along with diffusion coefficients obtained from Equations (86) and (87) for MESA profiles of our C/O WD model with a similar interior composition. The coefficients based on Stanton & Murillo (2016) vary from the MD results by a factor of two or more, while the residuals for our implementation of the Caplan et al. (2022) coefficients are 10% or less. This represents an order of magnitude improvement on the uncertainties for diffusion coefficients in the strongly coupled regime relative to the previous MESA implementation, while also providing a procedure that is generalizable to compositions other than just C/O/Ne mixtures (Equation 87).  \n<!-- image -->  \nFigure 23. Diffusion coefficient of 22 Ne in the liquid interior profile of a cooling C/O WD model. The residuals in the lower panel are relative to the Caplan et al. (2022) MD for a C/O/Ne plasma mixture shown in the top panel.  \n<!-- image -->  \nC/O mixtures with trace 22 Ne are in the liquid phase for Γ MCP glyph[lessorsimilar] 200 (Caplan et al. 2020; Blouin & Daligault 2021), though distillation of 22 Ne may also occur for sufficient 22 Ne concentration (Blouin et al. 2021).  \nWe have implemented these Caplan et al. (2022) diffusion coefficients as the default coefficients in MESA for Γ MCP > 10, with a smooth transition from the Stanton &Murillo (2016) coefficients over the range 3 < Γ MCP < 10. The diffusion coefficients in the liquid WD regime are now accurate to ≈ 10% for 10 glyph[lessorsimilar] Γ MCP glyph[lessorsimilar] 200, a sub-  \nial improvement compared to our previous methods.  \nWhen crystallization occurs for Γ MCP glyph[greaterorsimilar] 200, we assume that freezing into the solid phase causes the diffusion coefficient to go to zero. We implement this with a smooth turnoff near the crystallization boundary using the smoothed phase parameter φ from Skye ( § 4.1), so that the diffusion velocities are zero for φ ≥ 0 . 5.", '10. NUCLEAR PHYSICS': 'MESA models calculate the energy generation rates and composition changes due to nuclear burning over a large range of T , ρ , and X i . The nuclear evolution of the chemical composition dominates the total cost of a model (memory + CPU) when the number of isotopes is glyph[greaterorsimilar] 30. Here we report progress on nuclear reaction rates, a new operator split burning option, and enhanced reaction rate outputs.', '10.1. Reaction rates updates': "MESA 's default reaction rate for 7 Be( e -, ν e ) 7 Li came from REACLIB (Cyburt et al. 2010). However, REACLIB is only defined for T > 10 7 K, and assumes that all atoms are ionized. While this is a reasonable assumption for reactions that occur deep in the stellar interior, significant Li production can occur in stellar envelopes with T < 10 7 K where the reaction rate then depends on the ionization balance (Schwab 2020). The new default (as of r22.05.1) for this rate is Simonucci et al. (2013) which incorporates ionization contributions.  \nREACLIB defines a reverse reaction as the endothermic direction. However, this direction depends on the nuclear masses assumed during the evaluation. For consistency with the nuclear masses used in MESA , we define the reverse rate based on the nuclear masses from masslib library 5.data . This affects reactions with uncertain nuclear masses; for instance, the reaction 55 Cu( γ, p ) 54 Ni has a Q = 0 . 293 MeV (exothermic) from 'rpsm' (Rauscher 1999, priv. comm) and a Q = -0 . 07256 MeV (endothermic) from 'ths8' 7 (Cyburt et al. 2010). Incorrectly determining which reaction is endothermic leads to large errors when the reverse rate is computed from detailed balance.", '10.2. Operator split nuclear burning': "By default MESA uses a nuclear reaction network that is fully coupled to the stellar hydrodynamics (MESA I), solving for the changes in the composition simultaneously with the changes in stellar structure quantities.  \nThis approach provides consistency between the nuclear physics and stellar hydrodynamics, but can place limits on the maximum time step. When T glyph[greaterorsimilar] 3 × 10 9 K, the time step is limited by the nuclear burning timescale, and the composition enters a dynamic equilibrium state where large forward and reverse reaction rates nearly cancel each other, potentially leading to numerical errors when subtracting large values from one another (MESA III). During the final stages of evolution to corecollapse, time steps of δt < 10 -10 s are common (Farmer et al. 2016). This can make evolution to core collapse a resource-consuming endeavor.  \nMESA has the capability to perform an operator-split procedure ( op split burn ) to compute the composition change and energy generation. Cells with high T use operator splitting. For each such cell, MESA computes the change in composition over the timestep δt with a semi-implicit midpoint rule (Bader & Deuflhard 1983). An operator-split cell can then take multiple sub-steps, allowing the composition to evolve at fixed T and ρ with an adaptive time step. The nuclear energy generated glyph[epsilon1] nuc is calculated by taking the difference between the starting and final compositions, and energy loss due to neutrinos from nuclear reactions glyph[epsilon1] ν is accounted for.  \nThis scheme allows accurate tracking of the net nuclear energy generation rate and composition changes. However, operator-split burning cannot calculate the partial derivatives of these terms with respect to T or ρ for the matrix solver. These partial derivative matrix terms are thus set to zero. This removes the difficulty of the partial derivatives varying in sign and magnitude over short timescales, and enables the solver to more robustly find a solution within the requested tolerances.  \nFigure 24 compares the fully-coupled (unsplit) and operator split solutions in the central T -ρ plane for a 30M glyph[circledot] , solar-metallicity model evolved from the ZAMS to the formation of the iron core. We enable op split burn for any cell with a T > 10 9 K. The evolutions are nearly identical up to T c ≈ 3 × 10 9 K, after which the fully-coupled solution evolves to a slightly larger ρ . These differences are comparable to variations arising from other physics choices, such as the nuclear reaction rate screening prescription (see Appendix A.2 of MESAV). Evolution up to the formation of the iron core took approximately equal number of timesteps. Evolving the models further to the onset of core-collapse with 'gold' tolerances (MESA V), the operator split model only required an additional ≈ 200 time steps, while the fully-coupled model failed to reach the onset of corecollapse after an additional ≈ 140,000 time steps.  \nOperator splitting can provide a significant speedup during Si and Fe burning by reducing the number of  \nFigure 24. Central T -ρ tracks of a 30 M glyph[circledot] star evolved from ZAMS to formation of the iron core.  \n<!-- image -->  \n(  \n)  \ntime steps needed, and may be the only way to make certain problems tractable. It is however less efficient than a fully-coupled model during earlier stages of evolution due to the additional sub-steps taken by each operatorsplit cell. Thus, we suggest considering operator splitting only for models evolving beyond core C-depletion. It is difficult to say which treatment is more accurate; fully-coupled calculations include more physics but can be subject to numerical errors, while our operator-split calculations ignore the physics in the partial derivatives, but provide a more numerically stable solution.", '10.3. Reaction rate output': 'MESA has the new capability to output individual nuclear reaction rates and related quantities to the history or profile data files. Using add raw rates will output the rates of all reactions. Alternatively, raw rate with the name of a specific reaction will report only that reaction. Similarly, add eps neu rates , add eps nuc rates , and add screened rates provide neutrino energy rates, nuclear energy rates, and screened reaction rates.  \nAs an example, Figure 25 shows the total reaction rates of a 1 M glyph[circledot] model evolved to the solar age, and the profiles of the reaction rates at the solar age. Figure 26 further shows the evolution of the solar electron neutrino fluxes and a comparison with the observed solar electron neutrino fluxes (see also Farag et al. 2020) after correction for flavor mixing (Bergstrom et al. 2016). The model agreement with measurements is similar to that obtained using the ASTEC stellar evolution code (Bellinger & Christensen-Dalsgaard 2022).  \n<!-- image -->  \n/circledot  \nFigure 25. Nuclear reactions rates ongoing in the Sun. Upper: Evolution of reaction rates from the pre-main sequence to core H exhaustion. Lower: Profiles at the current solar age as a function of the distance from the solar center. The reactions in the colorbar are ordered from top to bottom according to their frequency.Figure 26. The evolution of solar electron neutrino flux at Earth. The predicted values for the present Sun are in satisfactory agreement with the observed values after correction for flavor effects: Φ(pp) / 10 9 = 5.934 ( MESA ), 5.971 ± 0.035 (Borexino); Φ(pep) / 10 7 = 1.408 ( MESA ), 1.448 ± 0.013 (Borexino); Φ( 7 Be) / 10 9 = 0.476 ( MESA ), 0.480 ± 0.023 (Borexino); and Φ( 8 B) / 10 5 = 4.780 ( MESA ), 5.160 ± 0.110 (Borexino) cm -2 s -1 .  \n<!-- image -->', '11. CONSTANTS': 'The MESA const module provides mathematical, physical, and astronomical constants relevant to stellar astrophysics (MESA I, Section 4.1). The values in this module have been updated to reflect new definitions and conventions. Some microphysics inputs (e.g., opacity, EOS, reaction rates) use slightly different constants in constructing their tables or expressions. Such cases are usually beyond our control, and the updates in this section do not apply to those instances.  \nAs part of the revision of the SI, the values of the physical constants e , h , k B , and N A are now exact (Newell et al. 2018). We have adopted these values and ensured that other constants composed of these and other exact values (e.g., the Stefan-Boltzmann constant) are defined in a consistent manner. For other physical constants, (e.g., G ) we use CODATA 2018 (Tiesinga et al. 2021).  \nMESA follows IAU recommendations for astronomical constants, currently adopting nominal solar and planetary quantities from IAU 2015 Resolution B3 (Mamajek et al. 2015). We follow the recommended procedure of deriving nominal solar and planetary masses from the mass parameters ( GM ) and the adopted value of G . We also adopt the convention that the nominal radii of planetary bodies are the equatorial radii.', '12.1. Migration to GitHub': "MESA development began in 2007 with a Subversion (SVN) repository hosted on SourceForge, and from 2017 -2020 this SVN repository was hosted by Assembla. Beginning in December 2020, development shifted to GitHub, with the root of the new Git repository corresponding to SVN r15140. An archival copy of the SVN development history is available at doi: 10.5281/zenodo.4745225. The shift to GitHub paves the way for a new era of collaborative MESA development and expanded interaction with the astrophysics community.  \nAs Git repositories can have a non-linear commit history, we will no longer denote public releases with a revision number. In SVN, this counted the number of commits from the beginning of MESA 's development. Instead, public releases will now be identified by the date, in the format YY.MM.I, where YY is the final two-digits of the year, MM is the two digit month number, and I is a version (usually 1) to distinguish multiple releases in a given month. While we do not recommend publishing an article based on non-released versions of MESA , if necessary we suggest using the first 7 characters of the Git commit hash as the version number.  \nWith the transition from SVN to Git, the versioning and branching scheme used in the development of MESA changed substantially. A major overhaul of the collection and distribution of continuous integration test results on TestHub (https://testhub.mesastar.org) was necessary. We now highlight these and other changes made since MESA V.  \nWith Git, understanding the relationship of a commit to other commits requires knowledge of the commit graph. To dynamically obtain and update this information, TestHub uses GitHub webhooks and APIs. Upon pushing to the GitHub remote, GitHub now sends a request to TestHub to update its internal tree of commits. This is done by querying the GitHub API for a current list of branches and the head commits of each. TestHub removes any references to branches that are no longer in GitHub and updates the branches table in the TestHub database, adding commits as needed even if they have not had any tests conducted. This allows easy identification of commits that have not yet been tested. We also now store data about commits, such as the authors, commit messages, and commit times directly in the TestHub database for more convenient access.  \nTo view the detailed changes made in a given commit, each commit-specific page on TestHub links to the appropriate page on GitHub. Additionally, most views in TestHub now have a branch selector dropdown so developers can look at commits on a specific branch, and we rely on the GitHub API to order these commits within a branch.  \nFigure 27 shows the total number of individual test instances (i.e., results for a single test case from a single computer) and daily commits submitted since the launch of TestHub in late 2017. Unsurprisingly, the rate of test submissions is strongly correlated with the rate of commits. Substantial deviations occur due to different testing computers coming on- and offline, and because some commits focus on documentation and do not trigger testing runs on automated testing setups. Figure 27 also shows the dates of public releases and key milestones. The switch from SVN to Git was shortly after release r15140, and there is a smooth transition between the eras.  \nThe MESA test suite takes ∼ 100 core hours to complete when run on a modern workstation. To reduce this runtime while still providing diagnostic information, many tests have optional steps that can be skipped, with the next inlist loading a saved model. The total time for a test suite run is roughly halved when these parts are skipped. Skipping optional inlists has drawbacks in requiring saved models, which can generate spurious fail-  \nFigure 27. Daily tests submitted since the instantiation of TestHub, spanning the SVN and Git eras. To smooth over the day-to-day variability, we also show 30-day rolling averages of submitted tests and pushed commits. Below the x-axis are dates of public releases and other milestones such as instrument papers and the SVN-to-Git transition.  \n<!-- image -->  \nFigure 28. Distribution of delay times from when a commit is pushed to when a failure is first reported.  \n<!-- image -->  \nes when skipping inlists. At least one computer runs all inlists on all tests on the main branch.  \nWe now assess how quickly our distributed and automated continuous integration workflow can identify a failure. Figure 28 shows that most failures are detected within an hour of the commit, and more than a third are detected within thirty minutes. This prompt turnaround comes from computing clusters configured to check for new commits every few minutes, and then launch a full test suite run upon detecting a new commit. Some clusters are configured to detect and test commits on any branch of the MESA GitHub repository, allowing  \nfull testing coverage during feature development and ensuring that new features pass the test suite before being merged into the main branch.  \nAwareness of test failure and the commit responsible is useful, but having the detailed output of a failing test is often essential in quickly identifying and rectifying the regression. TestHub now collects information on the runtime, computer architecture, broad failure type (such as a compilation error or a runtime error), and the stdout and stderr logs of the compilation and/or test case in the event of a failure. When present, these logs are accessible by links next to the failure indicator. This allows developers to quickly identify what led to the unexpected behavior without having to request more data from the owner of the computer that submitted it. These updates to TestHub improve the pace, efficiency, and quality of MESA source code development.", '12.3. User Contributions': "MESA has a Zenodo community 8 to encourage users to publicly archive their input and output files. To make user-contributed routines easier to share while minimizing maintenance, we have additionally created the mesa-contrib repository 9 .  \nRoutines in mesa-contrib can be implemented via MESA 's hooks with a few Fortran include statements. If the build environment defines the location of mesa-contrib , then it is included by default as a search path when MESA is compiled. Currently, mesa-contrib contains routines for the atmospheric T ( τ ) relations and corresponding MLT parameters implemented by Mosumgaard et al. (2018b) and for angular momentum transport by the enhanced Tayler-Spruit dynamo described by Fuller et al. (2019). To ease the burden of maintenance, mesa-contrib is only intended to work with the latest public MESA release.", '12.4. NuDocker': 'Reproducibility is a goal of open science and a tenet of scientific research. Provenance, as the term relates to software instruments (Van den Bussche & Vianu 2001; Carata et al. 2014; Stodden et al. 2018), is the ability to record the full history of a result. Scientific research is generally held to be of good provenance when previous results, perhaps decades old, can be reproduced. The aim is to preserve the final knowledge object and the capability to perform the scientific actions that are the foundation of the knowledge object. Thus, we seek  \nto preserve not only the numerical calculation, but accelerate future research by archiving the computational environment.  \nProvenance enables reproducing past simulations and performing new calculations with different physics or numerical options. A challenge is that compilers change, linked libraries evolve, and operating systems progress. Many science results have been obtained with older versions of MESA , which in most cases are not obsolete and therefore remain valuable. NuDocker 10 , 11 provides a solution to the provenance challenge by being able to run older versions of MESA with age-appropriate compilers, libraries, and operating system using light-weight, OSlevel virtualization (e.g., Docker 12 ). NuDocker provides four Docker images that can be launched with one terminal command, and has been tested and used in 14 out-of-the-box MESA versions from r4942 to r22.05.1. A hallmark of Docker virtualization is the minimal performance penalty compared to running natively (Felter et al. 2014, 2015).  \nFigure 29. Reproduction of Figure 13 from Farmer et al. (2015) using release r6794 in NuDocker : HRD of two 8 M glyph[circledot] models; one non-rotating (blue) and one rotating at Ω / Ω crit ∼ 0.2 (red). The evolutions span from H depletion to He depletion.  \n<!-- image -->  \nAs an example of the ability to provide almost decadeold results with NuDocker , Figure 29 reproduces Figure 13 from Farmer et al. (2015) using the same MESA version 6794 from July 2014 and inlists 13 as in Farmer et al. (2015). The virtual containers allow older versions of MESA to be run with bit-for-bit consistency for all versions after 7503 (see MESA III), thereby enabling older versions of MESA to be run on modern hardware, and preserving the required system environment to enable new research.', '13. SUMMARY': "We explain significant new capabilities and improvements implemented in MESA since the publication of MESA I through MESA V. Advances in automatic differentiation ( § 2) and time-dependent convection ( § 3) will open opportunities for future investigations in stellar evolution. Discussion of the current treatment of the energy equation ( § 5), stellar atmospheres ( § 6), and new formalisms for treating starspots and superadiabatic convection ( § 7) will enhance the study of stellar physics. Upgrades to the equation of state ( § 4), opacity ( § 8), element diffusion coefficients ( § 9), nuclear reaction rates ( § 10), and physical constants ( § 11) will increase the robustness of stellar evolution models. The transition to GitHub, upgrades to the MESA TestHub continuous integration framework ( § 12), the opening of a mesa-contrib repository, and NuDocker 's ability to run older versions of MESA will lead to an efficient and distributed model of source code development. Inlists and related materials for all the figures are available at doi: 10.5281/zenodo.6968760.  \nWe thank Dr. Bill Paxton for two decades of his amazing talents, extraordinary energy, and his generous gift of MESA to the astronomy and astrophysics community.  \nWe also thank Amy Mainzer for Figure 4 of MESA III making an appearance in the movie Don't Look Up (first noted by Jared Goldberg), and the participants of the MESA Summer Schools for their willingness to experiment with new capabilities and modalities of delivery. We thank the anonymous referee for a constructive, detailed, and thoughtful report that improved our manuscript. Finally, we thank Susie Groves for her heroic assistance in enabling the completion of this MESA VI instrument paper under challenging circumstances.  \nThe MESA Project is supported by the National Science Foundation (NSF) under the Software Infras-  \ntructure for Sustained Innovation program grants ACI1663684, ACI-1663688, and ACI-1663696. This research was supported in part by the NSF under Grant No. NSF PHY-1748958 for the Kavli Institute for Theoretical Physics. W.H.B acknowledges support from the UK Science and Technology Facilities Council (STFC) through grant ST/R0023297/1. G.C.C acknowledges support by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013 and the Astronomical Society of Australia. R.F acknowledges support of the University of Amsterdam's Helios cluster which was supported by a European Research Council grant 715063, (PI S.E. de Mink) F.H acknowledges funding through an NSERC Discovery Grant, through NSERC project award SAPPJ-202100032 and through the NSF under Grant PHY-1430152 for the JINA Center for the Evolution of the Elements. The Flatiron Institute is supported by the Simons Foundation. A.S.J thanks the Gordon and Betty Moore Foundation (Grant GBMF7392) and the National Science Foundation (Grant No. NSF PHY-1748958) for supporting this work. M.J acknowledges the Lasker Data Science Fellowship awarded by the Space Telescope Science Institute, and thanks Marc Pinnsoneault, Jen van Saders, and Jamie Tayar for many hours of consultation on the Yale Rotating Stellar Evolution Code and its documentation. J.S acknowledges support by NASA through Hubble Fellowship grant # HSTHF2-51382.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555, by the A.F. Morrison Fellowship in Lick Observatory, and by the National Science Foundation through grant ACI-1663688. R.S acknowledges support by the National Science Center, Poland, Sonata BIS project 2018/30/E/ST9/00598. A.T is a Research Associate at the Belgian Scientific Research Fund (F.R.S.-F.N.R.S.). F.X.T acknowledges support by NASA under the Astrophysics Theory Program grant NNH21ZDA001N-ATP, and by the NSF under Grant PHY-1430152 for the JINA Center for the Evolution of the Elements. T.L.S.W thanks support by the Gordon and Betty Moore Foundation through Grant GBMF5076. J.S.G.M acknowledges support by the KU Leuven Research Counsil (grant C16/18/005: PARADISE). O.T was supported by a FONDECYT project 321038. P.M. acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N. This research made extensive use of the SAO/NASA Astrophysics Data System (ADS).  \nFacilities: This work used the Extreme Science and Engineering Discovery Environment (XSEDE; Towns et al. 2014), which is supported by the NSF grant ACI1548562, specifically comet at the San Diego Supercomputer Center through allocation TG-AST180050. We thank Charlie Conroy and the Harvard ITC for providing computational resources for continuous testing of MESA through the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University. J.S acknowledges use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315, and thanks Josh Sonstroem and Brant Robertson for supporting this resource. A.S.J acknowledges use of the rusty supercomputer at the Flatiron Institute, supported by the Simons Foundation, and thanks the Scientific Computing Core for supporting this resource. W.H.B thanks the University of Birming-  \nham's Advanced Research Computing team for support of the BlueBEAR High-Performance Computing service. J.S.G.M thanks the VSC (Vlaams Supercomputer Centrum - Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government - department EWI. T.L.S.W acknowledges use of computational facilities at UC Santa Barbara funded by NSF grant CNS 1725797, and thanks the Center for Scientific Computing for supporting this resource.  \nSoftware: ÆSOPUS (Marigo & Aringer 2009) http://stev.oapd.inaf.it/cgi-bin/aesopus, FreeEOS (Irwin 2004), ipython/jupyter (P'erez & Granger 2007; Kluyver et al. 2016), matplotlib (Hunter 2007), mesaPlot (Farmer 2021) NuDocker https://github.com/ NuGrid/NuDocker, https://doi.org/10.5281/zenodo. 3678601, NumPy (van der Walt et al. 2011).", 'REFERENCES': 'Araki, T., et al. 2005, Phys. Rev. Lett., 94, 081801, doi: 10.1103/PhysRevLett.94.081801  \nAbbott, B. P., Abbott, R., Abbott, T. D., et al. 2018,  \nLiving Reviews in Relativity, 21, 3,  \ndoi: 10.1007/s41114-018-0012-9  \n- Abdurro\'uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35, doi: 10.3847/1538-4365/ac4414  \nAbe, K., Haga, Y., Hayato, Y., et al. 2016, Astroparticle Physics, 81, 39, doi: 10.1016/j.astropartphys.2016.04.003 Acciarri, R., Acero, M. A., Adamowski, M., et al. 2016, arXiv e-prints, arXiv:1601.02984.  \nhttps://arxiv.org/abs/1601.02984  \nAcernese, F., Agathos, M., Agatsuma, K., et al. 2015,  \nClassical and Quantum Gravity, 32, 024001,  \ndoi: 10.1088/0264-9381/32/2/024001  \n- Akeson, R., Armus, L., Bachelet, E., et al. 2019, arXiv e-prints, arXiv:1902.05569.  \nhttps://arxiv.org/abs/1902.05569  \n- Akutsu, T., Ando, M., Arai, K., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A102, doi: 10.1093/ptep/ptab018\n- Al Kharusi, S., BenZvi, S. Y., Bobowski, J. S., et al. 2021, New Journal of Physics, 23, 031201,  \ndoi: 10.1088/1367-2630/abde33  \n- Alexander, D., & Ferguson, J. 1994, The Astrophysical Journal, 437, 879  \nAlmgren, A., Sazo, M., Bell, J., et al. 2020, The Journal of  \nOpen Source Software, 5, 2513, doi: 10.21105/joss.02513  \nAndreoni, I., Coughlin, M. W., Almualla, M., et al. 2022,  \nApJS, 258, 5, doi: 10.3847/1538-4365/ac3bae  \nAndringa, S., et al. 2016, Adv. High Energy Phys., 2016, 6194250, doi: 10.1155/2016/6194250  \nArtigau, \' E., Sivaramakrishnan, A., Greenbaum, A. Z., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, 914340, doi: 10.1117/12.2055191  \nAshok, N. M., & Banerjee, D. P. K. 2003, A&A, 409, 1007, doi: 10.1051/0004-6361:20031160  \nAstropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f Astropy Collaboration, Price-Whelan, A. M., Lian Lim, P., et al. 2022, arXiv e-prints, arXiv:2206.14220.  \nhttps://arxiv.org/abs/2206.14220  \nBaalrud, S. D., & Daligault, J. 2014, Physics of Plasmas,  \n21, 055707, doi: 10.1063/1.4875282  \nBader, G., & Deuflhard, P. 1983, Numerische Mathematik, 41, 373, doi: 10.1007/BF01418331  \nBaiko, D. A. 2019, MNRAS, 488, 5042,  \ndoi: 10.1093/mnras/stz2041  \nBaiko, D. A., & Chugunov, A. I. 2022, MNRAS, 510, 2628, doi: 10.1093/mnras/stab3613  \nBaiko, D. A., Potekhin, A. Y., & Yakovlev, D. G. 2001,  \nPhys. Rev. E, 64, 057402,  \ndoi: 10.1103/PhysRevE.64.057402  \nBall, W. H. 2021, Research Notes of the American  \nAstronomical Society, 5, 7  \nBall, W. H., Chaplin, W. J., Schofield, M., et al. 2018, The Astrophysical Journal Supplement Series, 239, 34, doi: 10.3847/1538-4365/aaedbc  \n- Vogel, P., & Beacom, J. F. 1999, PhRvD, 60, 053003, doi: 10.1103/PhysRevD.60.053003\n- Vogl, C., Sim, S. A., Noebauer, U. M., Kerzendorf, W. E., & Hillebrandt, W. 2019, A&A, 621, A29,  \ndoi: 10.1051/0004-6361/201833701  \n- Wang, S., Jones, M., Shporer, A., et al. 2019, AJ, 157, 51, doi: 10.3847/1538-3881/aaf1b7\n- Welch, B., Coe, D., Diego, J. M., et al. 2022a, Nature, 603, 815, doi: 10.1038/s41586-022-04449-y\n- Welch, B., Coe, D., Zackrisson, E., et al. 2022b, ApJL, 940, L1, doi: 10.3847/2041-8213/ac9d39\n- Windhorst, R. A., Timmes, F. X., Wyithe, J. S. B., et al. 2018, ApJS, 234, 41, doi: 10.3847/1538-4365/aaa760'}