MESAHub
Bug in the nuclear reaction rates of reactions with more than two reactants and/or products
Hi, while participating in the Kavli summer program for astrophysics in MPA this year and after helpful discussions with @mathren , @rjfarmer, @2sn and Stephen Justham I believe I found a bug in the reaction rates in mesa-r23.05.1.
It seems that MESA miscalculates the reverse reaction rates of all reactions that involve more than two reactants and/or products. To my understanding, MESA takes the rates of reactions that release energy (Q>0) from the REACLIB database and computes the rate of their reverse reactions from detailed balance. However, it seems that the detailed balance performed by MESA omits the phase space factor that needs to be included for reactions with more than two reactants and two products. From my understanding, MESA handles well the computation of detailed balance of 1+2 ==>3+4 and 1 ==> 2+3 reactions, but gets wrong all the other types of reactions.
In the attached pdf file is a summary of all subclasses of reactions that behave weirdly but in accordance with the detailed balance argument. I plot the reaction rate vs temperature, both on log scales. The REACLIB rates* are plotted in blue, the MESA rates in red and the reverse reaction rates computed considering only the Boltzmann factor exp(-Q/(k_b T), appropriate for 1+2 ==> 3+4 reactions, in green.
This bug causes a drastic change in the equilibrium composition in very high temperatures and densities (e.g., T=7x10^9 K, \rho = 3x10^8 g/cm^3) in nets with many isotopes, as it overestimates some reaction rates by 24 orders of magnitude. I discovered it while running bbq with the nets mesa_330 and mesa_495, where the existence of tritium in the net changed completely the equilibrium composition, particularly by the overestimation of r_neut_nuet_he4_he4_to_h3_li7. So even though there is a relatively small number of reactions that this bug affects, it can have a huge impact on the nucleosynthesis results in high temperatures and densities, and in nets with a large number of isotopes.
*Obtained using the python library pynucastro.
Aldana Grichener
Where the bbq referenced above is @rjfarmer wrapper to the MESA one zone burner https://github.com/rjfarmer/bbq
Btw: I don’t expect this bug to affect most stellar evolution calculations, especially for models that include the standard (commonly used) isotopes nets.
No #526 was related to incorrectly tagging c12->3he4 as a forward rate not a reverse rate.
Is anyone working on a fix for this? Is there any more info about the bug? Thanks so much for reporting this bug, Aldana.
@SDcodenum and me are currently working on fixing the bug, but it might take a little while. What are you trying to simulate with MESA? This bug is probably nothing to worry about in most cases
@AldanaGrichener Glad to hear you’re working on a fix. I’m simulating RSGs up to core collapse using the mesa_204 net. I intend to blow them up with a supernova code, so the nucleosynthesis yields are important to me. Where in the code base should I look to investigate the issue?
Hello @AldanaGrichener, thanks for documenting this so well and digging into the MESA code! If you'd like you can try fix_rates_detailed_balance a potentially working fix. Perhaps you could help verify if it works?
Hi @Debraheem,
Thank you for the fix! I’ve re-created the plots from my original bug report using your fix. For some of the reactions the rate of the reverse reaction from MESA now almost coincides with the rates in REACLIB, but for others, even though the rate from MESA is much closer to REACLIB than before, there is still a couple orders of magnitude difference at high temperatures (e.g., h1_be9 to h2_he4_he4). My understanding is that the REACLIB rates don’t take nuclear partition functions into account while MESA does, but I’m not sure that’s enough to explain the remaining differences.
I’ve attached the plots that compare MESA’s and REACLIB’s reaction rates for the different reactions. I’ve also attached the python scripts I used to generate them, that might be of help. reactions.py extracts all the relevant reaction rates from rates_cache in MESA, which reactionsAnalysis.py then compares to the REACLIB rates.
Thanks @AldanaGrichener for being so helpful! It does indeed seem like a tremendous improvement, but I will double check the partition functions. Thanks for sharing your scripts!
Is there a specific nuclear network you tested your python scripts with? Most default nuclear networks only generate a portion of the rates called by the python scripts. With your network we can implement some minor testing infrastructure that checks these rates in the future.
To prevent delays because of timezones, I believe I know the answer to this one and can answer on @AldanaGrichener 's behalf.
We came across the presumed bug looking at the equilibrium composition reached when using mesa_330.net and mesa_495.net which differs significantly from the equilibrium composition reached with 201 or >800 isotopes (which agree to a reasonable degree with each other). This was done at T=10^9K and rho=3e8 g/cm^3, the attached figure (made by Aldana) shows an example run:
So using mesa_330.net and/or mesa_495.net should reveal issues.
I've made/attached three plots comparing the top 15 isotopes at logT = 9.84, logRho = 8.47 using the bbq burner inlists shared by @AldanaGrichener. With the rate fix, there seems to be good consistency across the three networks I tested, mesa_201.net, mesa_330.net, and mesa_495.net. I will add some documentation for this and merge it in soon. Expect this fix to ship with the next mesa candidate release coming soon.
Note, the equilibrium composition is entirely different than the ones you've arrived at, so this is a significant bug. the fix has some minor spillover into approx21 as well, see https://github.com/MESAHub/mesa/pull/632/files#diff-b0b4e553e85a61a8cec217302c25363051b15e1841cd657054df365f509813d7, but it doesn't appear nearly as significant. However for large networks, I would be more concerned.
Thanks again for raising this issue. I will close the thread when the the branch is merged, unless there are any additional concerns.