Question about compute_rev_ratio

[haakoan]

Question about compute_rev_ratio for 2-to-2 reverse rates in reaclib_support.f90

I am running late-stage stellar burning simulations in 3D using the FLASH code, for core-collapse supernova progenitors.
Initial conditions were constructed with MESA and I used approx21 for the MESA simulations and the FLASH
simulations.
When comparing the end state of my MESA and FLASH simulations I found differences
in the final compositions right above the iron core. In general, I see an enhancement in
sulfur in FLASH.

To understand where this difference comes from, I investigated the rates
relevant for this regime (temperatures around 1e9 K)

Concretely at T = 2.0e9 for the rate r_ar36_ap_k39
I find a ratio of MESA / FLASH: 4.53e-3 / 6.37e-4 = 7.11
Based on my understanding of these rates and how they are calculated,
the factor 7.11 is almost exactly the A-ratio 1 / (39/144)^(3/2) = 1/0.141.
The ratio stays essentially constant at ~7 across T = 2-4e9 K.
I see a similar pattern for r_si28_ap_p31 vs r_p31_pa_si28
(approx 6.9 at around T=1e9 K which matches 1 / (31/112)^(3/2) for the Si28/P31 pair).

I am not very familiar with MESA, but I am trying to understand the
code and find where the difference is coming from.

In rates/private/reaclib_support.f90 around line 214 (current main):

tmp = product(mass(1:Ni))/product(mass(Ni+1:Nt))
rates% inverse_coefficients(1,i) = pow(tmp, 1.5d0*(Ni-No)) * (product(g(1:Ni))/product(g(Ni+1:Nt)))

The code cites Reichert et al. 2023, which does indeed have the Ni-No factor in the A-products.
Other papers, for example Smith et al. 2023 (arXiv:2210.09965), does not have this factor (see their appendix B)

For a 2→2 reaction like K39(p,α)Ar36 (Ni = No = 2), pow(tmp, 1.5·(Ni−No))
evaluates to pow(tmp, 0) = 1.

If I'm reading things right, this line was changed in #632 (fix for #575) from
a constant exponent 1.5d0 to 1.5d0*(Ni-No).

The similar code pynucastro does it differently.
Smith et al. 2023 appendix B equation B16 writes the detailed-balance
ratio as two separate factors with different exponents:

• A-ratio: (∏ A_reactants / ∏ A_products)^(3/2) — exponent always 3/2, independent of (r-p)
• Thermal factor: (m_u kB T / 2π ℏ²)^(3(r-p)/2) — exponent scales with (r-p)

In pynucastro/rates/derived_rate.py, this is built up as a combination of
A_nuc^2.5 and 1/A factors that collapse to A^(1.5) regardless of
net_stoich, plus a separate net_stoich-dependent block that supplies the
thermal factor. The two factors are kept separate.

Question:

• Is the current pow(tmp, 1.5d0*(Ni-No)) line intentional, in the sense
  that there's a convention or a separate correction applied elsewhere that
  I haven't found? I did try to trace things through
  compute_some_inverse_lambdas in reaclib_eval.f90 and through
  net_approx21.f90 for the (α,p) cases I care about, and I didn't see an
  additional factors getting applied.

[Debraheem]

Hello Haakoan, thanks for raising this question. I've moved this discussion to #974.