eqstat2.f

EQSTAT2 is the core of YREC's equation of state: it chooses among, and blends, the Yale Saha solver, the relativistic degenerate gas, the SCV tables and the OPAL 2006 table, depending on temperature, density and the run's settings. Its header calls it "the old EOS routine"; since 2006 the code calls it through EQSTAT.

• Who calls it: EQSTAT (once or five times per call, see there) and CALCAD directly.
• Inputs: TL, PL, X, Z, LDERIV, LATMO; DL as a starting guess.
• What it calls: EQSAHA (Saha ionization), EQSCVE (SCV tables), EQRELV (fully ionized, degenerate gas), and for OPAL OEQOS06 with EQBOUND06 (the OPAL 1995 and 2001 tables, with OEQOS, EQBOUND, OEQOS01 and EQBOUND01, removed since rp_1).
• Outputs: as listed in EQSTAT: $T$, $P$, $\rho$, $\beta$, the ionization fractions, $\mu$ terms, $\eta$, $Q_{\rho T}$, $Q_{\rho P}$, $c_P$, $\nabla_\mathrm{ad}$ and their derivatives.
• Modules: NMLPAR_MOD (the namelist settings: the flag LOPALE06, LSCV, TSCUT, the temperature below which Saha is solved, and LDIFZ, metal diffusion); COMP_MOD (the envelope composition XENV, ZENV, AMUENV, FXENV, set in STARIN); CONSTANTS_MOD (constants: $\ln 10$, $R$ as CGAS, $a/3$ as CA3); IOUNITS_MOD (ISHORT).

1. Composition. AMU $=\sum_i X_i/A_i$ (ions per unit mass, in units of $1/m_\mathrm{H}$) and the number fractions FX of 12 species (Na, Al, Mg, Fe, Si, C, H, O, N, Ar, Ne, He; atomic weights in ATOMW2) are rescaled from the envelope values. Without metal diffusion, any change in $Z$ is taken as carbon ("ASSUME EXCESS Z(METALS) IS IN THE FORM OF CARBON(12)"); with it, all metals scale together by $Z/Z_\mathrm{env}$.

2. Radiation. $1-\beta = aT^4/(3P)$; then BETA14 is redefined as $4(1-\beta)/\beta$ (with $1-\beta$ capped at 1). This is the radiation term of the gas's response to temperature: for an ideal gas with radiation, $Q_{\rho T} = -(4-3\beta)/\beta = -1 - 4(1-\beta)/\beta$, and EQSAHA builds its $Q_{\rho T}$ (QDT = -1 - BETA14 - …, the last term from ionization), $c_P$ and their derivatives on it. After an OPAL result BETA14 is set to $1-\beta$ instead (line 334↗, line 348↗), which would be wrong in these formulas, but no caller reads the BETA14 it gets back: EQSAHA receives it within this call, before the OPAL step, and the callers of EQSTAT and EQSTAT2 only declare it and pass it on (settled from the code, Nov 9, 2026: line 165↗, lines 182–184↗, lines 311–358↗).

3. Fully ionized gas ($\log T \ge$ TSCUT): EQRELV alone, with $\mu_e^{-1} = 0.9921X + 0.4995(1-X)$ and $R/\mu = R(\mu_\mathrm{ion}^{-1} + \mu_e^{-1})$.

4. Partial ionization ($\log T <$ TSCUT): EQSAHA, and if LSCV is set, EQSCVE, which overwrites $\rho$ and the first derivatives where its table covers the point. A block after it would compute the SCV second derivatives numerically, but LSCVD is set to false just before it, so it never runs: unused? (lines 189–241↗; see bug 1).

   Where SCV supplies $\rho$, RMU, EMU and the second derivatives still come from EQSAHA (lines 181–187↗; eqsaha.f, lines 194–195↗), so below the transition band $\mathcal{R}/\mu$ is not $\beta P/(\rho T)$ there: EQSCVE's question 1.

5. Transition band (TSCUT − 0.5 < $\log T$ < TSCUT): EQRELV is also called, and the two are blended with a smooth step. With $w = 2(\mathrm{TSCUT} - \log T)$, between 0 and 1, $$W(w) = w^3(6w^2 - 15w + 10),\qquad q = q_\mathrm{rel} + W\,(q_\mathrm{Saha} - q_\mathrm{rel})$$ for $\log\rho$, the ionization fractions and the derivatives; then $R/\mu = \beta P/(\rho T)$ and $\mu_e^{-1}$ follows. Two gaps in the blend of the second derivatives are known to the developers (KNOWN_ISSUES.md, section 40, "EOS second derivatives"): The differences are not formed for QDTT and QDTP, so they come out as $q_\mathrm{rel} + W\,q_\mathrm{Saha}$ (lines 255–256↗, lines 277–285↗, lines 298–299↗), and the term from $dW/d\log T$ is left out (WT1, commented out, line 276↗). The same section records that the converged models are unaffected; what suffers is the rate and robustness of the Newton iteration.

6. OPAL. If LOPALE06 is set and the point is in the table's rough range ($1870 \le T \le 2\times10^8$ K, $\log\rho \le 7$), OPAL's values replace the above, or are blended with them by the weight $f$ from EQBOUND06: $$q = q_\mathrm{Yale} + f\,(q_\mathrm{OPAL} - q_\mathrm{Yale})$$ for $\log\rho$, $\beta$, the $\mu$ terms, $Q_{\rho T}$, $Q_{\rho P}$, $c_P$ and $\nabla_\mathrm{ad}$. OPAL gives no second derivatives, so those stay from the Yale/SCV result. If OPAL fails (alternate return) or the point is out of its table, the routine exits with the Yale/SCV result. At the end the routine restores Z from its copy Z0 (line 91↗, line 362↗), but nothing changes Z: no statement here, and neither OEQOS06 nor EQSCVE (or EQSCVG) assigns it (settled from the code, Nov 9, 2026: oeqos06.f, line 91↗, oeqos06.f, line 122↗): unused?

Microphysics in YREC describes the layers: OPAL (LOPALE06, the 2006 tables, Rogers & Nayfonov 2002, valid for $3.3 < \log T < 8.3$ with a density edge; "strongly recommended," although off by default), then SCV (LSCV, from three tables for pure H, He and metals, combined in a thermodynamically consistent way; the metal table was made with YREC's own EOS), then YREC's own EOS, with a ramp near the edge of the OPAL tables. With both on, SCV is used only outside OPAL; a full SCV calculation needs OPAL off and SCV on. OPAL assumes a fixed heavy-element mixture, so settling does not change $Z$ in it (a minor effect, since $P$ hardly depends on the metals). YREC's own EOS assumes full ionization above $\log T$ = TSCUT (default 6.0); the documentation warns not to change it, since the Saha equation for H and He fails at high density. YREC's semi-analytic derivatives are used for the structure equations even where a table gives the pressure.

What the code uses. The documentation's last sentence holds for the second derivatives only: where OPAL applies, its own first derivatives ($Q_{\rho T}$, $Q_{\rho P}$, $c_P$, $\nabla_\mathrm{ad}$) replace or are blended with the Yale ones (step 6), and SCV's replace them where its tables cover the point (step 4). Beyond OPAL's range ($T > 2\times10^8$ K or $\log\rho > 7$) only EQRELV applies: full ionization, relativistic degenerate electrons, and, as its one Coulomb term, the Debye–Hückel correction when LDH is set (as in every test-suite run; see EQRELV for its defaults). Of the 51 test-suite runs, 45 set LOPALE06 and 47 LSCV.