Microphysics
Microphysics copied to clipboard
AMReX-Astro
Move electron chemical potential contributions to energy into network RHS
This issue came up in the context of ECSN and Urca simulations with Maestro/MaestroEx.
These reaction networks do not include the contribution of electron chemical potential to the energy generation RHS, which appears any time Ye is changing. For the ECSN and Urca networks in particular, under some thermodynamic conditions weak reactions are the most important energy generation sources.
Maestro accounts for the electron chemical potential energy generation after the burn by calling the EOS with EXTRA_THERMO defined to get dhdX and dpdX and uses rho_omegadot from the burn to calculate an energy generation contribution to S in addition to rho_Hnuc.
This approach omits the electron chemical potential contribution to energy generation and temperature right hand sides in the network. If energy generation is dominated by reactions that change Ye, the temperature evolution in the network will be incorrect and the rates will be incorrect.
In Castro, this will also give an incorrect temperature after the burn for the same reason, plus the fact that Castro does not add the electron chemical potential energy generation term after the burn so it is omitted entirely.
We can fix this for all networks by writing a subroutine analogous to the temperature RHS or ener_gener_rate to calculate a contribution to ydot(ienuc) as:
#ifdef EXTRA_THERMO
ydot(ienuc) += dedX(1:nspec_evolve) * ydot(1:nspec_evolve)
#endif
And then have temperature_rhs call this subroutine before calculating ydot(itemp).
Thoughts anyone?
@maxpkatz
I worry your over-defining the system here. Changes in electron fraction should also be reflected in the temperature/energy RHS via the mass excesses.
If you're really worried about chemical potentials being the dominant source of energy generation, perhaps consider using an entropy RHS instead of energy. Mueller (1986) provides a good overview of how this can be implemented.
We decided at some point that the sum of (de/dy) * (dy/dt) = 0 identically. So I thought we were in agreement that we don't need the chemical potential term in either code.
Things seem more complicated when there are weak reactions and EOS Coulomb corrections ...
In that equality, our "y" only tracks nuclear species. We also have electrons as another "species" though, and the energy associated with absorbing or emitting an electron is not entirely reflected in the atomic masses. (I know we call this ion masses in our networks, but we include m_e in mion).
When we absorb or emit an electron, we are exchanging energy with the Fermi gas of degenerate electrons, and thus the local internal energy of the fluid. That energy exchange is not just m_e because the electron we capture could be from anywhere in the Fermi distribution with energy greater than the electron capture threshold energy.
This suggests we should account for this energy using the electron chemical potential.
Now that I'm saying this, if we use the electron chemical potential we should possibly not include m_e in mion because that may double count it.
Independent of weak reactions though, there are also de/dabar and de/dzbar terms that come from changes in the Coulomb screening potential for the EOS.
Without either weak reactions or Coulomb screening, treating the ions as an ideal gas then I agree that equality should be true and the contribution to e would come only from ion masses, or equivalently, reaction Q values.
This issue of double counting m_e if we use electron chemical potential may be what @jaharris87 was getting at by saying:
"Changes in electron fraction should also be reflected in the temperature/energy RHS via the mass excesses."
I've looked at the Mueller (1986) paper and it seems like that formulation does make this clearer and there is a path forward there by including electrons as a species and making sure not to double-count m_e.
It is probably also possible to start with that formulation and show what we should do if we integrate a specific energy equation instead as we are currently doing.
Linking the Mueller (1986) paper for reference: http://adsabs.harvard.edu/abs/1986A%26A...162..103M
