rates.f

Computes the nuclear reaction rates at one shell, for the abundance updates. Its header calls it "a stripped-down version of ENGEB": the same reactions, screening and constants as ENGEB, but without the energy generation. The rates are those of Bahcall's Neutrino Astrophysics (1989), with the alpha captures from Caughlan and Fowler (1988), as the comments say.

• Who calls it: MIX (for each shell, before the composition is mixed and burned), ROTMIX (for each shell of its substeps, since 7be4133), and XTIME (at the shell IEND, for the time step).
• Inputs: $\log_{10}\rho$ (DL), $\log_{10}T$ (TL), the mass fractions X, Y, XHE3, XC12, XC13, XN14, XO16, XO18, and the shell index IU.
• Outputs: the arrays HR1…HR13 (the rates of the 13 reactions) and HF1, HF2 (two branching fractions) at element IU. KEMCOM reads them (all but HR9, HR13 and HF2) to update the abundances; no routine computes with HF2, the Be7 branching $F_1$ (unused?). ROTMIX calls RATES itself, shell by shell, for its own substeps.
• Modules: CROSS_MOD (the derivative scalings QS0E, QQS0EE and the N15+p branching factors FO16, FC12, set in REMAP, shared with ENGEB); NMLPAR_MOD (the cross-section factors SSTANDARD(17), set in PARMIN or by REMAP; WEAKSCREENING, 0.03 by default in PARMIN; the temperature cuts TCUT); CONSTANTS_MOD (CLN $=\ln 10$, CC13 $=1/3$, CC23 $=2/3$).

1. Composition sums. Over the 13 isotopes (n, H1, D, H3, He3, He4, C12, C13, N14, O16, O18, Ne20, Mg24; those of n, D, H3, Ne20, Mg24 set to zero) with mass fractions $X_i$, bare nuclear masses $A_i$ (ANUC) and charges $Z_i$: $$\mu_I^{-1} = \sum_i \frac{X_i}{A_i}\ (\texttt{AMU}),\quad \mu_e^{-1} = \sum_i \frac{Z_iX_i}{A_i}\ (\texttt{EMU}),\quad \texttt{ZET} = \sum_i\frac{Z_i^2X_i}{A_i},\quad \texttt{XTR} = \sum_i\frac{Z_i^{1.58}X_i}{A_i}.$$
2. Below the temperature cut $\log T \le$ TCUT(1), all rates are zero (with the default 6.5 this stops pp burning and $^3$He production too, not only deuterium burning; see engeb.f, question 2). On that path (lines 346–351↗) F1 and F3 are not set, so HF1(IU), HF2(IU) would get the values kept (by SAVE) from the previous call; no caller takes it: MIX and ROTMIX call RATES only above TCUT(1) and zero the rates and HF1, HF2 of the colder shells themselves (mix.f, line 181↗, mix.f, lines 219–237↗; rotmix.f, lines 88–103↗), and XTIME's EQBURN does not read HF1 (settled from the code, Nov 9, 2026).
3. Electron degeneracy (for Salpeter's screening): with $\rho/\mu_e$ (RWE), $p_F^2/(m_ec)^2 = 1.017677\times10^{-4}(\rho/\mu_e)^{2/3}$ and $E_F/kT = 5.92986\,T_9^{-1}\bigl(\sqrt{1+p_F^2/m_e^2c^2}-1\bigr)$. The factor $f'/f$ (FPRF) is 1 for $E_F/kT \le 0.01$, 0 for $\log_{10}(E_F/kT) \ge 1.5$, and a sixth-order polynomial in $\log_{10}(E_F/kT)$ between.
4. Screening. $\Lambda_0 = 5.9426\times10^{-6}\,T_9^{-3/2}\sqrt{\rho/\mu_I}$ (XXL), $\tilde\zeta = \sqrt{(\texttt{ZET} + f'/f\cdot\mu_e^{-1})\,\mu_I}$ (ZCURL), $\bar Z = \mu_I/\mu_e$ (ZBAR). For reaction $j$ with charge product $Z_1Z_2$ (ZPRD): $$U_\mathrm{wk} = Z_1Z_2\,\Lambda_0\,\tilde\zeta .$$ If $U_\mathrm{wk} \le$ WEAKSCREENING, weak screening is used, $U = U_\mathrm{wk}$. Otherwise the intermediate form $$U_\mathrm{int} = 0.38\,\Lambda_0^{0.86}\,\frac{\texttt{XTR}\;\zeta_{86,j}}{\mu_I^{-1}\,\tilde\zeta^{0.58}\,\bar Z^{0.28}}$$ (Graboske et al. 1973, with their coefficients Z86), or, for $U_\mathrm{wk} > 2$, the strong form $$U_\mathrm{str} = 0.624\,\bar Z^{1/3}\Lambda_0^{2/3}\Bigl[\zeta_{53,j} + 0.316\,\bar Z^{1/3}\zeta_{43,j} + \frac{0.737\,\zeta_{23,j}}{\bar Z\,\Lambda_0^{2/3}}\Bigr]$$ (USTR, with the coefficients Z53, Z43, Z23) when it is smaller than $U_\mathrm{int}$ or when $U_\mathrm{wk} \ge 5$. The rate is multiplied by $e^{U}$.
5. Reactions 1–7 (pp, He3+He3, He3+He4, p+C12, p+C13, p+N14, p+O16), with $T_9 = T/10^9$ K: $$r_j = \rho\,\Bigl[T_9^{-2/3} + q_{1}T_9^{-1/3} + s'_j\,(q_{2} + q_{3}T_9^{1/3}) + s''_j\,(q_{4}T_9^{2/3} + q_{5}T_9)\Bigr]\exp\bigl(q_{6}T_9^{-1/3} + q_{7} + (q_{8}T_9)^2 + U_j\bigr)\;S_j,$$ with $q_{k} =$ Qk(j) (the terms of eq. 3.14 of Neutrino Astrophysics, as the comments explain), $s'_j$, $s''_j$ = QS0E(j), QQS0EE(j) (user scalings of the S-factor's first and second derivatives) and $S_j$ = SSTANDARD(j). Rates below $10^{-30}$ are set to zero. All are skipped when $X = 0$.

6. Be7. Electron capture: $\lambda_e = 1.752\times10^{-10}\,T_9^{-1/2}\,[1 + 0.004(1000\,T_9 - 16)]\,\mu_e^{-1}\,S_{15}$. Proton capture: $\lambda_p = 3.128813\times10^5\,X\,S_{16}\,e^{-10.2625\,T_9^{-1/3}}$ times the bracket of step 5 with Q(8), and its own weak or intermediate screening. Both are multiplied by $\rho\cdot 8.584981\times10^{22}$, and $F_1 = \lambda_e/(\lambda_e+\lambda_p)$ is the fraction of Be7 that captures an electron.

   Three comments disagree with the code. "F3 IS THE FRACTION OF THE N14 THAT IS BURNED BY P, ALPHA" (line 516↗) should read N15: F3 = C12ALPHA/(C12ALPHA + O16GAMMA) (line 598↗) splits the two branches of N15 + p. The comment above the screening loop says WEAKSCREENING is passed "in the FLUX COMMON block" (lines 433–434↗), but it comes from NMLPAR_MOD (nmlpar_mod.f, line 151↗; /WEAK/ at rp_1). And "DL AND DT" (line 331↗) are DL and TL, the routine's arguments (line 6↗). 8. Alpha captures (only when $\log T \ge$ TCUT(4)): C13$(\alpha,n)$O16 (RATE(8)), C12$(\alpha,\gamma)$O16 (RATE(10)), N14$(\alpha,\gamma)$F18 (RATE(11)) and the triple-alpha (RATE(12), $\propto\rho^2 T_9^{-3} e^{-4.4027/T_9}$), as fits of the Caughlan and Fowler form. RATE(9) and RATE(13) are always zero, and every rate at or below $10^{-5}$ is set to zero.

All constants are DATA statements (double precision since rp_1), with earlier sets kept as comments (BP00, 7/21/03). The routine duplicates much of ENGEB, so a change to one must be copied to the other.

Physics Options: WEAKSCREENING (default 0.03) is the threshold below which weak screening is used for the reaction rates. The light-element cross sections (XSLI6, XSLI7, XSBE91–XSBE93, keV b, with the switches LXLI6–LXBE93) are turned into the factors SLI6…SBE93 by PARMIN.