I noticed for later versions of MESA the option use_flux_limiting_with_dPrad_dm_form has ben removed. It would be great to get it back in order to have an approximate way to model the transition from diffusion to free-streaming radiative transport. Neglecting such physical prescription makes the radiative transport in the optically thin zones less reliable.

In my current scientific cases, I use MESA 23.05.01 to inject energy in stellar envelopes to model common envelope events. When the outer envelope of the primary expands, part of its outer regions reach velocities higher than the escape velocity. But, as such regions become optically thin, the current radiative transport leads to an overestimation of the luminosity. The energy of such escaping regions is being radiated at superluminal rates, they loose their kinetic energy and become bound. These runs use the HLLC Riemann solver, do not have any artificial viscosity that could damp the kinetic energy of such regions. I did some diagnostics by calculating the flux limiter for the profiles this behaviour appears and it looks that the surface luminosity is being overestimated by several order of magnitud (depending on the configuration of the common envelope).