  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.
