PaperlayerView exact revision on arXiv

Source-linked AI summary

Modules for Experiments in Stellar Astrophysics (MESA): Time-Dependent Convection, Energy Conservation, Automatic Differentiation, and Infrastructure

Adam S. Jermyn, Evan B. Bauer, Josiah Schwab, R. Farmer, Warrick H. Ball, Earl P. Bellinger, Aaron Dotter, Meridith Joyce, Pablo Marchant, Joey S. G. Mombarg, William M. Wolf, Tin Long Sunny Wong, Giulia C. Cinquegrana, Eoin Farrell, R. Smolec, Anne Thoul, Matteo Cantiello, Falk Herwig, Odette Toloza, Lars Bildsten, Richard H. D. Townsend, F. X. Timmes

arXiv:2208.03651v2astro-ph.SRastro-ph.EPastro-ph.HEastro-ph.IM

TL;DR

MESA needs improved treatments for rapidly changing convection, equation-of-state consistency, energy accounting, and stellar-modeling infrastructure. This paper introduces automatic differentiation and broad physics and software updates, resulting in expanded capabilities, greater model robustness, and opportunities for future stellar-evolution investigations.

  • Problem

    MESA requires improved treatment of convection that can evolve faster than steady state and equation-of-state issues affecting energy conservation and solver performance.

    Source: §3 · p. 7, §5.3 · p. 23

  • Method

    The paper updates MESA through automatic differentiation, time-dependent convection, equation-of-state and energy-equation improvements, expanded stellar physics, efficiency methods, and infrastructure upgrades.

    Source: §13 · p. 43

  • Results

    The updates expand MESA capabilities, enhance stellar-physics studies and model robustness, and improve software development and community engagement.

    Source: §13 · p. 43

  • Takeaways & Limitations

    MESA provides a broader and more robust framework for investigating stellar evolution across rapidly changing convection, degenerate conditions, massive-star evolution, and other stellar-physics regimes.

    Source: §3 · p. 7, §5.2 · p. 22, §13 · p. 43

    Operator splitting is less efficient during earlier evolution, and its accuracy relative to fully coupled calculations is difficult to determine.

    Source: §10.2 · p. 39

Abstract

from arXiv · show

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 late-stage 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 nuclear burning mode. We close by discussing major updates to MESA's software infrastructure that enhance source code development and community engagement.

1. INTRODUCTION

MESA VI responds to rapidly expanding stellar observations and advances in open-knowledge software by introducing new capabilities while documenting their scope and limitations. The paper organizes updates spanning automatic differentiation, stellar physics, energy treatment, opacities, nuclear physics, and infrastructure, while deferring broader scientific implications to future community work.

Motivation

  • New observational facilities are generating an unprecedented volume of high-quality stellar measurements, strengthening the empirical foundation of stellar astrophysics.Examples include gravitational-wave detectors, neutrino networks, sky surveys, spectroscopic facilities, and stellar-variability missions.Source: §1 · pp. 2–3
  • Advances in software infrastructure, computing, data storage, and open-knowledge instruments are transforming interactions among stellar theory, modeling, simulations, experiments, and observations.The paper places MESA within an ecosystem that includes Astropy, Athena++, Castro, Dedalus, GYRE, and other community tools.Source: §1 · pp. 3–4

Scope

  • MESA VI presents new capabilities and limitations without fully exploring their scientific implications, which are left to future efforts by the MESA research community.The software-instrument paper compares capabilities with available numerical or analytic results, while future community work is expected to unlock their scientific potential.Source: §1 · p. 4

Organization

  • The paper covers automatic differentiation, time-dependent convection, the equation of state, energy equations, atmospheres, starspots, opacities, diffusion, nuclear physics, physical constants, infrastructure, and a final summary.These topics are developed in Sections 2–12 and summarized in Section 13.Source: §1 · p. 4

2. AUTOMATIC DIFFERENTIATION

MESA’s new auto_diff module replaces much hard-coded derivative computation with forward-mode operator-overloaded automatic differentiation, reducing complexity and error while supporting higher-order derivatives. Benchmarks show performance close to explicit expressions, and the capability enables developments across several MESA physics modules.

Motivation and design

  • Automatic differentiation largely eliminates the complexity and error sources of hard-coded analytic derivatives required by MESA’s implicit Newton-Raphson solver.The solver requires derivatives accurate to about one part in 10^6, often precluding finite differences.Source: §2 · p. 5

Implementation

  • The auto_diff module uses compile-time-sized Fortran types, overloaded operators, chain-rule expressions, and broad function support to compute partial derivatives automatically.Discontinuous functions use zero derivatives, while discontinuous derivatives are evaluated as averages between the two sides; optimized operators are generated with SymPy and support bit-for-bit reproducibility through CR-LIBM.Source: §2 · pp. 5–6

Applications

  • Additional auto-diff types support higher-order mixed partial derivatives and applications to time-dependent convection, the equation of state, starspots, and superadiabatic convection.The module can also be used in run_star_extras and software development outside MESA.Source: §2 · pp. 6–7

Performance

  • For sufficiently complex operations, optimized auto-diff functions outperform manually implemented derivative routines, although simple operations retain overhead.Relative cost against explicit routines no longer scales with the number of derivatives for the tested first-order types.Source: Table 4 · p. 6, §2 · p. 7
  • In stellar-evolution calculations, hand-coded derivatives are modestly faster, but their small auto-diff overhead is outweighed by the full timestep cost and improved convergence from correct derivatives.Incorrect derivatives can cause slow convergence and more small timesteps; some MESA source code is still transitioning to auto_diff.Source: §2 · p. 7

3. TIME-DEPENDENT LOCAL CONVECTION · 3.1. Formalism · 3.2. Numerical Stability

MESA extends time-dependent convection beyond steady-state MLT for rapidly evolving stellar structures, using a local Kuhfuß formalism that evolves turbulent kinetic energy. The implementation modifies the RSP approach to improve numerical stability across short convective and long evolutionary timescales.

3. TIME-DEPENDENT LOCAL CONVECTION

  • Time-dependent convection is needed when stellar structure evolves faster than convection reaches steady state, including late-stage massive-star burning and electron-degenerate ignition events.Steady-state MLT is appropriate when composition and structure evolve slowly relative to convective timescales.Source: §3 · p. 7

3.1. Formalism

  • MESA implements the one-equation local Kuhfuß (1986) formalism, extending its prior RSP use to general stellar-evolution calculations.The model evolves specific turbulent kinetic energy et and incorporates turbulent pressure, convective luminosity, and energy transport into the stellar-structure equations.Source: §3 · p. 7, §3.1 · pp. 7–8
  • The default αr = 0 neglects radiative damping, an approximation that changes the model’s physical completeness but enables modifications making TDC reduce to MLT at long timesteps.Other defaults are α = 2, αD = 1, and αPt = 0, with αPt = 0 implying Pt=0.Source: §3.1 · p. 8
  • The formalism couples turbulent kinetic energy to heat and momentum transport through turbulent pressure and convective luminosity, with the local-limit assumption setting Lt to zero.The resulting luminosity contributes to the total energy equation and determines the evolution of specific internal energy.Source: §3.1 · p. 8
  • MESA drops viscous bulk-motion heating and uses the Ledoux gradient in Y when applying the Ledoux stability criterion.The omitted heating term is expected to be unimportant because bulk velocities are typically much smaller than convective velocities.Source: §3.1 · p. 8

3.2. Numerical Stability

  • Implicitly solving et with the structure equations works on short convective or pulsational timescales but is numerically unstable on long evolutionary timescales.The desired method must operate in both limits and transition smoothly between them.Source: §3.2 · p. 8
  • Long-timestep errors amplify because small errors in superadiabaticity Y produce large luminosity errors that alter entropy, feed back into Y, and increase subsequent luminosity errors.Shorter-than-thermal-timescale steps can restore equilibrium, whereas long steps allow the error propagation to affect the entropy profile.Source: §3.2 · p. 8
  • MESA’s MLT-inspired stable strategy treats luminosity as the solver variable, solving for Y given L because Y is relatively insensitive to luminosity errors.This reverses the RSP logic of solving for L given Y and reduces error amplification.Source: §3.2 · p. 8

3.3. Flipped Equations

MESA implements time-dependent convection with luminosity L as the solver variable for long-timescale numerical stability. The flipped formulation solves for Y while preserving local solvability, Jacobian structure, and runtime performance.

3.3. Flipped Equations

  • MESA uses L as the solver variable and flips Equation (7) to solve for Y, accounting for w's dependence on Y through Equation (2).Source: §3.3 · p. 8
  • The implementation numerically inverts Equations (2) and (7) to determine each Y_k from local luminosity L_k and adjacent-cell solver variables.Source: §3.3 · p. 8
  • Setting L_t = 0 removes direct neighboring-cell coupling, preserving MESA's basic Jacobian structure and runtime performance.Source: §3.3 · p. 8

3.4. Numerical Method

The numerical method constructs L from Y by implicitly evolving convective velocity, then inverts this relation with branch-aware root finding. Automatic differentiation propagates derivatives through the Newton-Raphson solve, enabling TDC derivatives with respect to stellar-structure variables.

3.4. Numerical Method

  • The overall procedure solves Equations (2) and (7) by constructing L given Y and numerically root-finding Y given L.This provides the stated numerical route for solving the coupled relations with L_t = 0.Source: §3.4 · p. 8, §3.4.2 · p. 9

3.4.1. L given Y

  • TDC implicitly evolves w using end-of-step coefficients, handles growth or decay according to the discriminant, sets w = 0 after a stable-branch zero crossing, and evaluates L.This produces L given Y while maintaining numerical stability.Source: §3.4.1 · p. 9

3.4.2. Numerical Inversion

  • The inversion searches in Z ≡ ln |Y| over −100 ≤ Z ≤ 100, then combines bisection with Newton-Raphson refinement to obtain a near-machine-precision, differentiable solution.The search window covers 10^-43 ≲ |Y| ≲ 10^43, while at least one Newton-Raphson iteration preserves partial derivatives for the MESA Jacobian.Source: §3.4.2 · p. 9
  • For Y < 0 with initial w > 0, up to three solutions can occur, so the solver identifies solution intervals and selects the smallest-|Y| root for slowest-decaying convection.The smallest-|Y| choice connects smoothly to Y = 0 and avoids the solution where convection decays on a dynamical timescale.Source: §3.4.2 · p. 10

3.4.3. Relation to auto diff

  • Automatic differentiation propagates partial derivatives through the Newton-Raphson solver for Y with respect to 33 stellar-structure variables, enabling the TDC implementation.TDC returns Y given L and the other solver variables together with their partial derivatives.Source: §3.4.3 · p. 11

3.5. Reduction to Cox MLT

The TDC formulation is corrected so that its long-time-step limit matches Cox MLT mathematically and numerically, including in inefficient convection. The correction also preserves the existing TDC iteration algorithm while aligning derived convective quantities.

3.5. Reduction to Cox MLT

  • The required correction replaces Y with YΓ/(1 + Γ) by adding the factor Γ/(1 + Γ) to the TDC convective-luminosity equation.The Γ term controls convective efficiency and is the genuine difference between TDC-in-RSP and MLT; Γ is evaluated with MLT inputs and held fixed during TDC iterations.Source: §3.5 · pp. 11–12
  • Equivalence requires w = 3/2vc in steady state, with Cox MLT f1 = 1/8 matching the default TDC choice αD = 1.For subadiabatic states, both approaches give zero convective velocity at long times; for unstable states, the corrected velocity equation recovers the MLT solution.Source: §3.5 · pp. 11–12
  • TDC convective mixing diffusivity and other derived quantities are calculated from vc consistently with Cox MLT, while the correction remains necessary when αr > 0.Without correction and with αr > 0, the inefficient-limit asymptotic scaling differs qualitatively from corrected TDC with αr = 0.Source: §3.5 · pp. 11–12
  • TDC and Cox MLT become mathematically identical at long time steps and agree numerically to around seven decimal places in Y, even when Y ≪1.Figure 2 shows corrected TDC and MLT solutions overlapping, whereas uncorrected TDC deviates significantly, especially in the inefficient-convection limit.Source: §3.5 · p. 12

3.6. Accreting White Dwarfs

MESA’s time-dependent convection (TDC) treatment changes helium-flash evolution in accreting white dwarfs when nuclear heating varies faster than convection can adjust. Compared with Cox MLT, TDC produces slower, more superadiabatic convection and modifies thermal, velocity, and abundance profiles during rapid burning.

Accreting White Dwarfs

  • When the heating timescale is shorter than the convective acceleration time, TDC diverges from Cox MLT, producing higher TBCZ at fixed PBCZ and greater deviations for thicker helium shells.For models 3 and 4, TDC and Cox MLT diverge during expansion; when theat ≳ taccel in models 1 and 2, their evolution agrees well.Source: §3.6 · p. 13, Figure 3 · p. 13
  • TDC convection cannot reach steady state during rapidly burning helium flashes because taccel evolves faster than teddy once theat ≲ taccel.In models 3 and 4, taccel becomes up to 3 and 6 times smaller than teddy, respectively, at high TBCZ.Source: §3.6 · pp. 13–14, Figure 4 · p. 14
  • TDC traps heat near the BCZ, yielding colder outer layers, higher PBCZ at fixed TBCZ, stronger superadiabaticity, and reduced or frozen convective velocities near the convection-zone top.These differences become prominent after theat drops below taccel, while TDC and Cox MLT initially show reasonable agreement in convective velocity.Source: §3.6 · pp. 14–15, Figure 5 · p. 14, Figure 6 · p. 15
  • Both treatments develop carbon abundance gradients, but Cox MLT produces more overall 12C because it has more time for outward transport, while TDC alters profiles relevant to possible detonation.The modified thermal and composition profiles may affect whether ignition develops into a thermonuclear transient.Source: §3.6 · p. 15
  • TDC is enabled by default in MESA, incorporating the slower growth and decay of convection during rapidly burning phases.The models accrete helium-rich material onto a 1 M⊙ carbon-oxygen white dwarf, with accumulated helium shells spanning 0.03–0.08 M⊙.Source: §3.6 · pp. 13, 15

4. EQUATION OF STATE

MESA’s EOS combines multiple prescriptions to cover broad thermodynamic regimes, with Skye and FreeEOS extending physical coverage and improving consistency. New treatments address crystallization, EOS blending, and thermodynamic errors affecting energy conservation and solver convergence.

4. EQUATION OF STATE

  • MESA evaluates its EOS extremely frequently across broad T, ρ, and Xi ranges, motivating continued progress toward accurate, thermodynamically consistent, and efficient coverage.The EOS may be called approximately 10^4–10^10 times depending on local physics, iterations, cells, and time steps.Source: §4 · p. 15
  • MESA prioritizes HELM, Skye, FreeEOS, OPAL, and SCVH across thermodynamic regimes, with ideal gas plus radiation as a fallback and composition-dependent blend boundaries.The default boundaries are selected for specific use cases, including keeping solar models at the Sun’s age entirely on FreeEOS.Source: §4 · p. 15, Figure 7 · p. 16

4.1. Skye

  • Skye models fully ionized matter with positrons, relativity, degeneracy, Coulomb interactions, mixing, quantum corrections, and self-consistent crystallization using automatic differentiation of Helmholtz free energy (Jermyn et al. 2021).Its free-energy formulation makes Skye readily extendable by adding further physical terms.Source: §4.1 · p. 16

4.1.1. Crystallization

  • MESA captures crystallization latent heat with a smoothed phase parameter and an additional heat source, preserving numerical convergence without smearing the phase-transition thermodynamics.The Skye-based phase diagram and latent heat are calculated self-consistently from first principles, and the net release is commensurate with the Salaris et al. (2000) value.Source: §4.1.1 · pp. 16–17

4.2. FreeEOS

  • FreeEOS version 2.2.1 replaces MESA’s eosPTEH tables and expands partial-ionization coverage beyond OPAL across a broad chemical-composition parameter space.MESA provides EOS1 tables over Z = 0–1 with hydrogen and helium variations, carbon–oxygen tables for high-Z white-dwarf interiors, and tools for generating arbitrary compositions.Source: §4.2 · p. 17

4.3. EOS Blends

  • EOS blends preserve non-negative fundamental quantities and include blending-coefficient derivatives in Jacobian derivatives, providing correct solver derivatives and reducing Newton iterations.Fundamental quantities and derivatives are blended differently because physical equations require physical values while the solver requires accurate derivatives; this can violate thermodynamic identities in blend regions.Source: §4.3 · pp. 17–18

4.4. Thermodynamic Consistency

  • Skye and FreeEOS are thermodynamically consistent to near machine precision, while Skye improves consistency significantly for log(ρ/g cm−3) ≳3 and FreeEOS improves over OPAL.Skye’s Helmholtz-free-energy formulation enables near-machine-precision consistency; FreeEOS remains inconsistent after table interpolation but improves substantially relative to OPAL.Source: §4.4 · p. 18

5. ENERGY EQUATIONS

This section examines MESA’s energy-equation implementations and accounting, including the effects of EOS derivative quality, EOS blending, and equation residuals. It assesses how energy-equation formulation shapes accuracy and the interpretation of reported energy errors, including during 1 M⊙ evolution.

5. ENERGY EQUATIONS

  • The reported rel E err and rel run E err primarily diagnose numerical reliability, not physical validity, and their meaning differs between dedt and eps grav.For dedt, energy error mainly reflects solver residual quality; for eps grav, larger values can reflect EOS inconsistencies even when the solution is comparably accurate or better.Source: §5 · pp. 18–19, §5.1 · p. 20

5.2. Results: carbon kh

  • The two energy equations should agree for converged solutions with a thermodynamically consistent EOS, whereas degenerate conditions favor eps grav because it avoids cancellation errors amplified by the degeneracy parameter η.In the carbon test, all forms converge toward the same result at sufficiently high resolution, but eps grav reaches accurate entropy and temperature evolution faster under degeneracy.Source: §5.1 · pp. 19, 21, §5.2 · p. 22

5.1. Fundamentals and Implementations

  • MESA now accounts for composition changes in all energy-equation forms and uses implicit trapezoidal time-centering by default for eps grav, improving numerical energy conservation.The accounting includes operator-split composition changes from diffusion and convective premixing, with compensating source terms added to the solved equations.Source: §5.1 · pp. 20–21

5.2. Results: carbon kh

  • The centered eps grav form best balances accurate temperature evolution and step-to-step energy conservation for degenerate conditions, outperforming dedt in the carbon-contraction test.The eps grav variants outperform dedt for central-temperature convergence; the centered version achieves second-order convergence in numerical energy conservation, while dedt can report smaller errors despite less accurate global evolution.Source: §5.2 · pp. 21–23, Figure 9 · p. 22, Figure 10 · p. 23

5.3. Quantifying EOS shortcomings

  • EOS internal inconsistencies and unavoidable inconsistencies from blending remain important energy-accounting limitations, motivating explicit thermodynamic-consistency measures and continued EOS improvements.MESA upgrades its EOS patchwork with Skye and FreeEOS, while minimizing but not eliminating blend-related inconsistencies.Source: §5.3 · p. 23, §5.3.2 · p. 24

5.3.1. EOS Derivatives

  • MESA mitigates EOS derivative problems by using spline derivatives in the Jacobian while using spline-interpolated derivatives in the equations, enabling efficient numerics and smoother solutions.Negative interpolant derivatives may guide Newton iterations but cannot introduce physical errors into the model equations.Source: §5.3 · p. 23, §5.3.1 · p. 24

5.3.2. Thermodynamic Consistency and EOS Blends

  • EOS blends can violate thermodynamic consistency even when blended EOS components differ only by an energy zero point, producing nonzero EOS energy errors regardless of the energy-equation form.Blending-coefficient derivatives are omitted from the blended cV, violating cV = (∂e/∂T)_ρ.Source: §5.3.2 · p. 24
  • MESA minimizes but cannot eliminate EOS-blend offsets, which remain a major potential energy-error source; dedt can inject spurious heating or cooling, whereas eps grav preserves the physical evolution in the idealized case.Energy diagnostics may report zero for dedt despite a larger actual EOS error, while eps grav reports large errors reflecting EOS inconsistency rather than inaccurate evolution.Source: §5.3.2 · pp. 24–25

5.4. Results: 1M pre ms to 1TP

  • In the 1 M⊙ evolution, EOS-blend errors approximately converge, whereas non-blend EOS inconsistency measures shrink with increasing time resolution and remain unconverged.The model evolves from the pre-main sequence to its first He thermal pulse using FreeEOS blends with OPAL/SCVH and Skye.Source: §5.4 · p. 25
  • Figure 11 shows that dedt hides EOS energy errors in reported rel run E err, while centered eps grav becomes dominated by EOS inconsistencies and blends by the end of the main sequence.Further resolution or energy-equation improvements cannot reduce this error; progress requires improving the EOS.Source: §5.4 · pp. 25–26, Figure 11 · p. 25

6. ATMOSPHERE

MESA expands its atmosphere capabilities with selectable T(τ), irradiated, and tabulated models, plus two approaches for incorporating atmospheric structure into pulsation and interior models. These updates improve consistency, accuracy, and flexibility across stellar modeling applications.

6.1. Atmospheric T(τ) Relations

  • T(τ) atmospheres support four relations and fixed, iterated, or varying opacity treatments, while corrected callbacks and derivatives improve atmosphere–interior consistency and Newton convergence.The model surface optical depth τs need not equal the photospheric value τ ≈2/3; the varying-opacity option uses dopri5 integration.Source: §6.1 · pp. 26–27

6.2. Irradiated Atmospheres

  • Atmosphere selection includes an irradiated option using the Guillot & Havel (2011) T(τ) relation, with fixed or iterated opacity and a corrected Teff evaluation.The iterated irradiated option replaces the former gray-irradiated choice.Source: §6.2 · p. 27

6.3. Tabulated White Dwarf Atmospheres

  • For a 0.57 M⊙, log(g/cm s−2) ≈8.0 helium-atmosphere white dwarf, new DB tables reproduce Koester et al. (2020) calculations much better at low Teff than Eddington T(τ) atmospheres, which diverge for Teff ≲15,000 K, while tables agree to Teff ≈7,000 K.Source: Figure 12 · p. 27, §6.3 · pp. 27–28

6.4. Atmosphere Reconstruction for Pulsation Codes

  • Atmosphere reconstruction outward from τs to τ ≪1 now uses dopri5 with configurable controls, correcting hydrostatic-balance errors from older explicit-Euler integrations.The reconstruction affects pulsation-code input rather than the interior model and is enabled with add atmosphere to pulse data = .true.Source: §6.4 · p. 28, Figure 13 · p. 28

6.5. Atmospheres as Part of the Interior

  • Including atmospheres as part of the interior keeps EOS and opacity treatments consistent, enables parallel computation, and accounts for gravity variation across the atmosphere.The method corrects the radiative gradient using 1 + dq/dτ and places the surface boundary at a controlled smaller optical depth.Source: §6.5 · p. 29

6.6. Choosing atmosphere options

  • Opacity treatments form an accuracy–cost hierarchy: varying is more accurate than iterated, which is more accurate than fixed, while Eddington with fixed opacity remains the default.Atmosphere reconstruction is particularly important for pulsations with substantial near-surface amplitudes.Source: §6.6 · p. 29

7. CONVECTION IN THE OUTER LAYERS OF STARS

This section adds MESA treatments for starspots and an implicit, local method to enhance convection near the Eddington limit. The implementations reproduce spotted-star inflation and improve numerical stability and evolution in massive-star and binary calculations.

7.1.1. Starspots formalism

  • The starspots formalism represents surface magnetism through coverage fraction fspot and temperature contrast xspot = Tspot/Tphotosphere, with observationally reasonable fspot = 0.34.Both parameters range from 0.0 to 1.0, while physically meaningful xspot values are restricted to 0.5−1.0; xspot = 0.5 corresponds to magnetic equipartition.Source: §7.1.1 · p. 30

7.1.2. MESA implementation

  • MESA implements starspots by modifying pressure-based surface boundary conditions through magnetic pressure PB, with auto diff supplying derivatives throughout the stellar structure calculation.PB is fixed at each evolutionary-step start and held constant during Newton iterations; the spotted luminosity is adjusted by αspot.Source: §7.1.2 · pp. 30–31

7.1.3. MESA Models

  • MESA starspots predict 3%-5% radius inflation for a 0.3 M⊙ dwarf with fspot = 0.34 and xspot = 0.85, agreeing with SPOTS models.Tracks shift smoothly with fspot, while lower xspot produces a larger coverage-fraction impact; below ≈0.6 M⊙, SPOTS models are cooler because of atmospheric-boundary differences.Source: §7.1.3 · p. 31, Figure 15 · p. 31

7.2. Superadiabatic Convection in Radiation-dominated Regions

  • The new convection treatment targets density inversions and impractical timesteps in radiation-dominated massive stars, where standard MLT can produce extended low-density envelopes.Its stability is especially useful for binaries because MLT++ step variations can cause large mass-transfer-rate changes during Roche-lobe overflow.Source: §7.2 · pp. 31, 33
  • Near the Eddington limit, MESA adds a fully implicit, local alternative to MLT++ that reduces superadiabaticity by adjusting ∇rad and recomputing convection.The enhancement scales ∇rad − ∇L using fΓ, acts when the expected Eddington factor is high, and is an engineering method rather than a calibrated physical model.Source: §7.2 · pp. 31–32
  • The implicit method produces smoother 30 M⊙ evolution than MLT++ while remaining closer to unenhanced behavior; MLT++ shows undesirable numerical variations.For a 40 M⊙ binary, the unenhanced model stalls during second mass transfer, whereas the implicit method provides a smooth solution throughout mass transfer.Source: §7.2 · pp. 33–34, Figure 17 · p. 33, Figure 18 · p. 33

8. OPACITY

MESA expands opacity physics across molecular, Compton, conductive, and monochromatic-opacity treatments. These updates support composition-dependent stellar modeling, improved degenerate-envelope cooling calculations, and faster radiative-acceleration evaluations.

8.1. Molecular Opacities

  • MESA adds low-temperature molecular opacity tables from Lederer & Aringer (2009) and ÆSOPUS, supporting variable CNO abundances in cool-star models.For 3 M⊙, Z = 0.014 TP-AGB models, MESA reaches higher peak L, smaller peak Teff, and larger peak R than Monash, while interpulse periods converge to ≃7.2×10^4 yr after the first few pulses.Source: §8.1 · p. 34, Figure 20 · p. 35

8.2. Compton Opacities

  • At log(T/K) ≳8.7, MESA now computes Compton-scattering opacity with Poutanen (2017), replacing the previous Buchler & Yueh (1976) prescription for extreme-temperature stellar conditions.These conditions commonly occur in massive-star simulations and thermonuclear bursts on neutron stars.Source: §8.2 · p. 34

8.3. Conductive Opacities

  • MESA implements Blouin et al. (2020) analytical conductive-opacity fits for H and He in moderately coupled, moderately degenerate regimes relevant to DA and DB white-dwarf envelopes.A 0.9 M⊙ hydrogen-atmosphere white-dwarf cooling calculation confirms the result of Blouin et al.Source: §8.3 · p. 35, Figure 21 · p. 35

8.4. Radiative Accelerations from the OP Monochromatic Opacities

  • The new OP monochromatic-opacity implementation is about 5 times faster for main-sequence stars while preserving κR and grad,i accuracy.For 16O, cell-averaged and maximum differences are 0.007 and 0.042; for 40Ar, they are 0.02 and 0.03. Stored bicubic interpolants and precomputed factors reduce repeated wavelength-integration costs.Source: §8.4 · pp. 36–37

9. DIFFUSION COEFFICIENTS

MESA now uses Caplan et al. (2022) molecular-dynamics-informed diffusion coefficients for strongly coupled plasmas, generalizing them to arbitrary mixtures through charge-scaled equivalent one-component-plasma coefficients. The implementation improves liquid white-dwarf diffusion accuracy to approximately 10% and smoothly suppresses diffusion upon crystallization.

Motivation

  • Diffusion coefficients based on Stanton & Murillo (2016) differ from molecular-dynamics results by factors of a few in strongly coupled plasmas.The discrepancy applies when ΓMCP ≳1 and motivates the updated treatment.Source: §9 · p. 37

Method

  • MESA approximates arbitrary-mixture diffusion to within approximately 10% for 10 ≲ ΓMCP ≲ 200 by scaling an MD-calibrated equivalent OCP coefficient with ion charge.The OCP term represents global viscosity, while charge scaling represents species-dependent effective radii in the liquid.Source: §9 · p. 37

Validation

  • Caplan et al. (2022) coefficients reduce residuals to 10% or less versus molecular dynamics, while the previous coefficients differ by a factor of two or more.Figure 23 compares 22Ne diffusion in a 0.6 M⊙ C/O white dwarf; the update is an order-of-magnitude improvement and generalizes beyond C/O/Ne mixtures.Source: Figure 23 · p. 38, §9 · p. 38

Implementation

  • The Caplan et al. (2022) coefficients are now MESA defaults for ΓMCP > 10, with a smooth transition from Stanton & Murillo (2016) coefficients over 3 < ΓMCP < 10.The coefficients are recast as symmetric Burgers-equation resistance coefficients Kij for use by MESA’s diffusion solver.Source: §9 · pp. 37–38

Limits

  • For crystallized matter at ΓMCP ≳ 200, MESA smoothly turns diffusion off using the Skye phase parameter, making diffusion velocities zero for φ ≥ 0.5.The strong-coupling resistance-coefficient relation can neglect thermal-diffusion corrections of up to about 20% when ΓMCP ≪ 10.Source: §9 · pp. 37–38

10. NUCLEAR PHYSICS

MESA updates nuclear physics through revised reaction rates, operator-split burning for high-temperature evolution, and expanded reaction-rate diagnostics. These changes improve consistency, numerical robustness, late-stage evolution efficiency, and comparison of solar-neutrino predictions with observations.

10.1. Reaction rates updates

  • MESA replaces the default 7Be electron-capture rate with Simonucci et al. (2013), incorporating ionization contributions relevant to Li production below 10^7 K.The revision addresses REACLIB’s T > 10^7 K and fully ionized-atom assumptions; reverse rates are also defined consistently with masslib nuclear masses to avoid large detailed-balance errors.Source: §10.1 · p. 38

10.2. Operator split nuclear burning

  • MESA’s operator-split burning evolves high-temperature composition with adaptive semi-implicit sub-steps, computes net nuclear and neutrino energy changes, and improves solver robustness by zeroing unavailable T and ρ derivatives.Fully coupled burning can impose δt < 10^-10 s near core collapse because of fast reactions and numerical cancellation, while operator splitting cannot provide the corresponding partial derivatives.Source: §10.2 · p. 39
  • The operator-split model reached core-collapse onset with only ≈200 additional time steps, whereas the fully coupled model failed after ≈140,000 additional steps.For a 30 M⊙ solar-metallicity model, both approaches were nearly identical until Tc ≈3 × 10^9 K; operator splitting is recommended beyond core C-depletion despite reduced early-stage efficiency and omitted partial derivatives.Source: §10.2 · p. 39, Figure 24 · p. 39

10.3. Reaction rate output

  • MESA can now output individual, screened, nuclear-energy, and neutrino-energy reaction rates to history or profile files.The new controls support either all reactions or a named reaction, enabling detailed inspection of nuclear burning and related quantities.Source: §10.3 · p. 39
  • The solar model’s predicted present-day neutrino fluxes agree satisfactorily with Borexino measurements after flavor correction, with agreement similar to ASTEC.Figure 26 reports Φ(pp)/10^9 = 5.934 (MESA), 5.971 ± 0.035 (Borexino) and Φ(8B)/10^5 = 4.780 (MESA), 5.160 ± 0.110 (Borexino) cm−2 s−1.Source: §10.3 · p. 39, Figure 26 · p. 40

11. CONSTANTS

MESA updates its constants module to reflect revised SI definitions, CODATA 2018 values, and IAU 2015 recommendations for astronomical quantities.

Scope of updates

  • The constants module updates mathematical, physical, and astronomical values, but these updates generally do not apply to microphysics tables or expressions built with separate constants.Opacity, equation-of-state, and reaction-rate inputs may therefore retain slightly different constants.Source: §11 · p. 40

Physical constants

  • Exact SI values for e, h, kB, and NA are adopted, with derived constants defined consistently; other physical constants such as G use CODATA 2018.The Stefan-Boltzmann constant is given as an example of a derived constant treated consistently.Source: §11 · p. 40

Astronomical constants

  • MESA adopts IAU 2015 nominal solar and planetary quantities, derives nominal masses from mass parameters GM and G, and uses equatorial planetary radii (Mamajek et al. 2015).These conventions follow the recommended IAU procedure for astronomical constants.Source: §11 · p. 40

12. INFRASTRUCTURE

MESA’s infrastructure updates support collaborative development, continuous integration, user contributions, and reproducible execution of older software versions. The transition to GitHub, expanded TestHub, mesa-contrib, and NuDocker improve development, sharing, testing, and provenance.

12.1. Migration to GitHub

  • MESA moved development from SVN to GitHub in December 2020, enabling more collaborative development and broader interaction with the astrophysics community.The new Git repository begins at SVN r15140, with the prior development history archived at doi: 10.5281/zenodo.4745225.Source: §12.1 · p. 40

12.2. MESA TestHub

  • Most test failures are detected within an hour, while TestHub links commits, branches, failure logs, runtime data, and architectures to accelerate regression diagnosis.More than a third of failures are detected within thirty minutes; clusters test new commits every few minutes, including on development branches.Source: §12.2 · pp. 41–42, Figure 28 · p. 41

12.3. User Contributions

  • MESA’s Zenodo community and mesa-contrib repository support public archival of inputs and outputs and sharing of user-contributed routines through MESA hooks.mesa-contrib currently includes atmospheric T(τ) relations with corresponding MLT parameters and enhanced Tayler–Spruit angular-momentum transport routines.Source: §12.3 · p. 42

12.4. NuDocker

  • NuDocker runs older MESA versions with age-appropriate environments, reproducing past simulations and enabling new calculations despite changing compilers, libraries, and operating systems.It provides four Docker images and has been tested with 14 out-of-the-box versions from r4942 to r22.05.1.Source: §12.4 · p. 42

13. SUMMARY

MESA VI summarizes substantial capabilities and improvements since MESA I–V, spanning stellar-physics modeling, numerical methods, and software infrastructure. These updates expand future research opportunities, strengthen model robustness, and support more efficient, distributed development.

Core capabilities

  • Automatic differentiation and time-dependent convection add capabilities that open opportunities for future investigations in stellar evolution.The summary identifies these advances as central new capabilities in MESA VI.Source: §13 · p. 43

Stellar physics

  • Updates to the energy equation, stellar atmospheres, starspots, and superadiabatic convection enhance the study of stellar physics.These treatments address physical modeling areas highlighted in the summary.Source: §13 · p. 43

Model robustness

  • Upgrades to the equation of state, opacity, element diffusion coefficients, nuclear reaction rates, and physical constants increase the robustness of stellar evolution models.The summary groups these improvements as robustness-enhancing updates.Source: §13 · p. 43

Software infrastructure

  • GitHub, MESA TestHub continuous integration, mesa-contrib, and NuDocker support a more efficient and distributed model of source code development.NuDocker also enables running older MESA versions, while inlists and figure materials are available at doi: 10.5281/zenodo.6968760.Source: §13 · p. 43
Loading 2208.03651v2…