PyREC scaffolding (from YREC development, 7be4133)/ sneut.fGlossary

sneut.f

1355 lines · 1 definitions · on GitHub (7be4133)

Python version: sneut.py, 1152 lines, on the right (switch to Python above the code).

Purpose

The thermal neutrino energy-loss rate of stellar matter, from the analytic fits of Itoh et al. (1996, ApJS 102, 411), with its derivatives with respect to temperature, density, mean atomic weight and mean charge. Five processes are summed: pair annihilation, plasmon decay, the photoneutrino process, bremsstrahlung and recombination. The routine is written in a different style from the rest of YREC (lower case, implicit none, c.. comments); it is probably an import of a widely used public routine, our inference.

How it fits

  • Who calls it: NULOSSES, which first computes $\bar A$ and $\bar Z$ from the mass fractions (AZBAR); NULOSSES is called by NEUTRINO, which ENGEB calls only when LNULOS1 is set (engeb.f, lines 1240–1244↗). LNULOS1 is off by default (nmlpar_mod.f, line 246↗; a legacy option in the documentation, "Use Itoh et al. (1996) neutrino loss rates"), but 710 of the 965 namelists in the repository that set it turn it on, with the comment "recommended to be set to TRUE" (settled from the code, Nov 9, 2026).
  • What it calls: IFERMI12 (the inverse Fermi–Dirac integral of order 1/2, to get the electron chemical potential) and ZFERMIM12 (the integral of order −1/2, for its derivative), in the recombination part.
  • Inputs: temp ($T$, K), den ($\rho$, g cm$^{-3}$), abar ($\bar A$), zbar ($\bar Z$).
  • Outputs: snu, the total loss in erg g$^{-1}$ s$^{-1}$, and its derivatives dsnudt, dsnudd, dsnuda, dsnudz. No COMMON blocks.

The method

Below $T = 10^7$ K everything is zero. Otherwise, with $Y_e = \bar Z/\bar A$, $T_9 = T/10^9$ K, $\lambda = T_9/5.9302$ (that is, $kT/m_ec^2$, xl) and $\xi = (\rho Y_e/10^9)^{1/3}/\lambda$ (zeta), each process follows the equations of Itoh et al. that the comments cite:

  • Pair (eqs. 2.5–2.8): $Q_\mathrm{pair} = \tfrac12\bigl[C_+ + C_-\,q_\mathrm{pair}\bigr]\,g(\lambda)\,e^{-2/\lambda}\,f_\mathrm{pair}$, where, with $n$ = xnufam neutrino flavors, $$C_\pm = C_V^2 \pm C_A^2 + (n-1)\bigl(C_V'^2 \pm C_A'^2\bigr),\qquad C_V = \tfrac12 + 2\sin^2\theta_W,\ C_A = \tfrac12,\ C_V' = 1 - C_V,\ C_A' = 1 - C_A,$$ ($C_+$ is tfac1, $C_-$ is tfac2, tfac3 $= C_-/C_+$, tfac4 $= C_+/2$, tfac5 $= C_-/2$), $g(\lambda)$ a polynomial in $\lambda^2$ and $f_\mathrm{pair}$ a rational function of $\xi$ and $\lambda$.

  • Plasma (eqs. 4.1–4.11): transverse and longitudinal plasmon terms with a correction factor $f_{xy}$.

  • Photoneutrino (eqs. 3.2–3.8): a Fourier-type fit in $\tau = \log_{10}(T/10^7)$, with coefficient sets for $10^7$–$10^8$, $10^8$–$10^9$ and above $10^9$ K; negative results are set to zero.
  • Bremsstrahlung (eqs. 5.1–5.26): the weakly degenerate form when $T > 0.01\,T_F$, with $T_F = 5.9302\times10^9\bigl(\sqrt{1 + 1.018(\rho Y_e/10^6)^{2/3}} - 1\bigr)$ K; otherwise the liquid-metal form "with C12 parameters (not too different for other elements)."

  • Recombination (eqs. 6.5–6.14): with the chemical potential $\nu$ from $\mathcal F_{1/2}^{-1}$, only for $-20 \le \nu \le 10$.

The sum is converted from erg cm$^{-3}$ s$^{-1}$ to erg g$^{-1}$ s$^{-1}$ by dividing by $\rho$ ("comment these out to duplicate the Itoh et al. plots"), and every derivative is carried along by the chain rule.

Suspected bugs, checked by running the code

Each bug is marked by what running the code showed: (changes results, large) or (small), the fix moving a run clearly beyond a null change; (stops the run); (wrong output), the fix changing only what is written; (harmless in the runs), within the noise; (never occurred), the code ran but the fault did not; (not reached); (cannot be reached), a guard prevents it; (withdrawn), the suspicion was wrong.

1. (changes results, small) Checked by KC (Nov 10, 2026) A typo in tfac2.

(lines 124–125↗) The parameter is tfac2 = cv*cv - ca*ca + (xnufam-1.0d0) * (cvp*cvp - cap-cap). By the pattern of tfac1 and the definition of $C_-$, the last term should be cvp*cvp - cap*cap. As written it is $C_V'^2 - 2C_A' = C_V'^2 - 1$, so $C_-$ = −1.318 instead of 0.182, and tfac3 $= C_-/C_+$ = −0.784 instead of 0.108. tfac3 enters the pair rate as $1 + $ tfac3$\,q_\mathrm{pair}$ (lines 356–357↗) and the photoneutrino rate as $1 - $ tfac3$\,q_\mathrm{phot}$ (lines 834–835↗), and tfac5 enters both bremsstrahlung forms (line 1015↗, line 1171↗). Taken whenever LNULOS1 is on (off by default, on in most of the repository's namelists) and $\log T$ is above TCUT(5).

A possible fix:

sneut.f

@@ lines 121–128 @@
      5                  cap    = 1.0d0 - ca,
      6                  tfac1  = cv*cv + ca*ca +
      7                           (xnufam-1.0d0) * (cvp*cvp+cap*cap),
+c 11/09/26 KC via Claude
+c code changed to square cap in tfac2, as in tfac1 and C- of Itoh et al.
      8                  tfac2  = cv*cv - ca*ca +
-     9                           (xnufam-1.0d0) * (cvp*cvp - cap-cap),
+     9                           (xnufam-1.0d0) * (cvp*cvp - cap*cap),
      &                  tfac3  = tfac2/tfac1,
      1                  tfac4  = 0.5d0 * tfac1,
      2                  tfac5  = 0.5d0 * tfac2,

How the code itself shows it. The routine cites Itoh et al. (ApJS 102, 411, 1996) for its rates (line 9↗), and its other sums of couplings follow one pattern: tfac1 = cv*cv + ca*ca + (xnufam-1.0d0) * (cvp*cvp+cap*cap) (lines 122–123↗) and tfac6 = cv*cv + 1.5d0*ca*ca + (xnufam - 1.0d0)*(cvp*cvp + 1.5d0*cap*cap) (lines 129–130↗), in which each primed term repeats its unprimed one, squares included. tfac2 is the same sum with the axial term subtracted, cv*cv - ca*ca, and only its primed part, cvp*cvp - cap-cap (lines 124–125↗), is linear in cap. The two are used as a pair: tfac4 and tfac5 are half of tfac1 and tfac2 (lines 127–128↗), and both bremsstrahlung forms take tfac4*fbrem - tfac5*gbrem and tfac4*fliq - tfac5*gliq (line 1015↗, line 1171↗). The fix squares cap as tfac1 and tfac6 do.

Checked by running YREC (Nov 9, 2026, the counts of testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite repeated on Nov 10 with a private start model; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): Of 34 runs (26 of YREC's cases, from a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning, and 8 made-up ones), the four that set LNULOS1 and pass $\log T$ = 7.5 computed the losses, at $\log T$ = 7.50 to 8.07: 33,460,185 times in testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite, 101,222 in examples/class_exercises/exercise_4/Grid_m4p0feh+0p0GS98_std25_TAMS, and 32,030 and 10,622 in the rotating testsuite/Test_m2p0_feh+0p0_vzams150_IES1_TAMS and testsuite/Test_m2p0_feh+0p0_vzams150_IES2_TAMS. In every call the pair, photo, and bremsstrahlung terms, which use tfac3 and tfac5, carried from 1% to all of the loss. With the fix, IES2 keeps its 546 models and end age, and moves $\log L$ by up to $7.95\times10^{-5}$ dex and $\log T_c$ by $3.31\times10^{-5}$, against $1.58\times10^{-8}$ and $1.72\times10^{-8}$ for a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$): a real effect, but a small one. IES1's track is unchanged in the compared columns but for central $X$ ($10^{-11}$), though all three files differ beyond 2% in columns the comparison does not cover. The giant-branch case, where the losses matter most, was not run with the fix.

Definitions

subroutineSNEUT

lines 2–1355 · GitHub

Called by NULOSSES

Flowchart of its calls, in order (2)
IFERMI121196→ZFERMIM121200
Compiler warnings (1)

gfortran -Wall -Wextra, compiling only; nothing was run.

  • line 480: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Three single-precision constants in the recombination part became double precision: the table-12 coefficients c (9.06e-1 and 8.97e-1 to 9.06D-1 and 8.97D-1) and the exponent and factor 4.55 in dd01 = dum**(-4.55D0) and c01; this changes the recombination rate by a few parts in $10^8$. The tfac2 typo (bug 1) is unchanged.