## P.8 Reproducibility

The MATLAB routine `tsrtsfactor.m` computes $S_{\text{TSRT}}(E)$ for user-specified reactions and energy ranges, using the same constants as in Appendix E. It exports tables directly comparable to [21], including uncertainties propagated from numerical quadrature. Sample scripts reproduce Table 27 (p. 169) and generate energy-dependence plots for arbitrary fusion pairs, ensuring that all results reported in the main article can be independently validated.

## P.9 MATLAB: TSRT S-factor Script

This script computes and plots the TSRT astrophysical S-factor for a given binary reaction using universal TSRT slopes ($\alpha, \beta$) and a single per-system amplitude $A_{\text{sys}}$ fitted on the upper-energy subset of the dataset. The high-energy fit window suppresses low-energy distortions (screening and hindrance) so that universal slopes remain truly system-independent; the fitted amplitude then sets the absolute normalization for the specific entrance channel.

**Purpose.** `tsrtsfactor.m` implements the canonical TSRT form

$$S(E) = A_{\text{sys}} \exp(-\alpha \Xi_{1}(E) - \beta \Xi_{2}(E)), \qquad (297)$$

using the same invariants $\Xi_1, \Xi_2$ as defined in Equation (290):

$$\Xi_{1}(E) = \frac{\sqrt{\mu c^{2} Z_{1} Z_{2}}}{(A_{1}^{1/3} + A_{2}^{1/3})\sqrt{E}}, \qquad \Xi_{2}(E) = \frac{\mu c^{2} (Z_{1} Z_{2})^{2}}{(A_{1}^{1/3} + A_{2}^{1/3})^{2} E}, \qquad (298)$$

with $E$ in MeV and $\mu$ in MeV/c² so that both invariants are dimensionless. By default the script uses the universal slopes $\alpha, \beta$ (calibrated once; Appendix P.4) without any extra renormalization.

*Optional sensitivity toggles (off by default).* For stress tests only, the code can apply small, explicit renormalizations

$$\alpha_{\text{eff}} = \alpha (1 + a_1 \chi) [1 + c_1(R_{\text{cal}}/R - 1)], \qquad \beta_{\text{eff}} = \beta (1 + b_1 \chi) [1 + c_2(R_{\text{cal}}/R - 1)], \qquad (299)$$

where $R = A_1^{1/3} + A_2^{1/3}$, $\chi = |A_1 - A_2|/(A_1 + A_2)$, and $R_{\text{cal}}$ corresponds to a chosen symmetric calibration system. When these options are not explicitly enabled (the setting used for all results and plots in this paper), $\alpha_{\text{eff}}=\alpha$ and $\beta_{\text{eff}}=\beta$ exactly, recovering Equation (297) with $\Xi_{1,2}$ as in Equation (290).

### Inputs and outputs.

* **Inputs:** reaction struct `reaction.Z1,A1,Z2,A2`; energy vector `E_MeV` (MeV, positive); universal slopes `alpha,beta`; output filenames `out_csv,out_fig`; and an optional `opts` struct:

- `opts.S_exp` (same length as `E_MeV`) to fit $A_{\text{sys}}$;

- `opts.units` (label only, e.g. 'MeV-b');

- `opts.a1,b1` (asymmetry renormalization, default 0);

- `opts.c1,c2` (symmetric-system tweak, default 0);

- `opts.Rcal_A` (calibration $A_{\text{cal}}$ for $R_{\text{cal}}$, optional);

- `opts.title`, `opts.label` (plot cosmetics).

* **Outputs:** a CSV with columns `{E_MeV, S_TSRT, Xi1, Xi2, alpha_eff, beta_eff, A_sys}` and a publication-ready figure (`out_fig`) showing $\log_{10} S(E)$.
**Fitting and windowing.** If `opts.S_exp` is provided, the script sorts the energies and fits $A_{sys}$ by log-space intercept on the upper third of the dataset (least-squares mean of log $S_{exp}$ – log $S_{noAmp}$). This choice minimizes bias from low-energy screening/hindrance while keeping the universal slopes fixed. If no experimental $S(E)$ is supplied, the script produces a parameter-free prediction with $A_{sys} = 1$.

**Dependencies and units.**

* Base MATLAB only; no toolboxes required.

* Uses 1 u = 931.49410242 MeV/c² (CODATA-level) internally to form $\mu$; all input energies are in MeV.

* The units string in the y-axis label is cosmetic; numerical calculations use purely dimensionless $S(E)$ as per the TSRT factorization.

**Reader guidance.**

* No hidden constants: $(\alpha, \beta)$ are passed explicitly; optional $(a_1, b_1, c_1, c_2)$ are zero unless the user sets them.

* The high-energy window is explicit in the code (upper third by index after sorting); change this window by editing two lines if robustness must be stress-tested.

* `opts.S_exp` must be strictly positive (log-fit); the script checks vector sizes and energy positivity.

* The CSV is machine-readable and contains the effective slopes actually used $(\alpha_{eff}, \beta_{eff})$ so any plot can be reconstructed independently.

* Figure export uses `exportgraphics` with a robust `print` fallback; outputs match the styles used elsewhere in the appendix.

**Typical workflow.**

1. Fix universal slopes $(\alpha, \beta)$ once from a reference system (Appendix P).

2. For a new reaction, call `tsrtsfactor` with `S_exp` to obtain $A_{sys}$ and the full TSRT curve; the script writes the CSV and saves the plot.

3. For prediction-only overlays (no fit), call with $A_{sys} = 1$ (omit `S_exp`).

This single routine replicates the S-factor figures in the appendix and in the main article's fusion section without per-point or per-energy tuning.

Code. `tsrtsfactor.m` performs TSRT S-factor predictions using universal slopes (with optional asymmetry and symmetric-only tweaks) and an amplitude fit on the upper-energy third. The full listing is in Appendix Y, Listing 10 (p. 225).

**What this produces.** `tsrtsfactor` computes a TSRT S-factor curve

$$S(E) = A_{sys} \exp(-\alpha_{eff} \Xi_{1}(E) - \beta_{eff} \Xi_{2}(E)), \qquad (300)$$

writes a CSV with invariants and effective slopes, and saves a plot (log$_{10}S$ vs $E$). If experimental $S_{exp}(E)$ is provided, a single amplitude $A_{sys}$ is fitted on the upper-energy third of the data; otherwise $A_{sys} = 1$.
## Inputs and units.

* **reaction** — struct with integer fields Z1, A1, Z2, A2.

* **E_MeV** — column vector of center-of-mass energies in MeV (strictly > 0).

* **α, β** — universal TSRT slopes (dimensionless) from Table 16 (p. 110) (Appendix P); held fixed for all systems.

* **out_csv, out_fig** — output filenames; figure format taken from extension.

* **opts** (optional) — struct:

- **S_exp** (same size as E) for overlay + amplitude fit (must be > 0 for log fit).

- **units** label for the S-axis (e.g. 'MeV-b'), used only for plot text.

- **a1, b1** asymmetry renormalization coefficients; $\chi = \frac{|A_1 - A_2|}{A_1 + A_2}$.

- **c1, c2** symmetric-only tweak coefficients tied to $R = A_1^{1/3} + A_2^{1/3}$.

- **Rcal_A** sets $R_{cal} = (A_{cal}^{1/3} + A_{cal}^{1/3})$ for the symmetric tweak.

- **title, label** for plot customization.

**Theory mapping, implemented.** The code constructs TSRT invariants

$$ \Xi_1(E) = \frac{\sqrt{\mu Z_1 Z_2}}{R \sqrt{E}}, \quad \Xi_2(E) = \frac{\mu (Z_1 Z_2)^2}{R^2 E}, \qquad (301) $$

with $\mu = \frac{m_1 m_2}{m_1 + m_2}$ in MeV/c² using $m_i = A_i \times 931.49410242$ MeV/c², and $R = A_1^{1/3} + A_2^{1/3}$. Effective slopes are

$$ \alpha_{\text{eff}} = \alpha (1 + a_1 \chi) [1 + c_1 (R_{\text{cal}}/R - 1)], \quad \beta_{\text{eff}} = \beta (1 + b_1 \chi) [1 + c_2 (R_{\text{cal}}/R - 1)], \qquad (302) $$

reducing to the universal $(\alpha, \beta)$ when $a_1 = b_1 = c_1 = c_2 = 0$.

## Outputs.

* **CSV columns:** E_MeV, S_TSRT, Xi1, Xi2, alpha_eff, beta_eff, A_sys.

* **Figure:** log₁₀S(E) vs E with optional experimental overlay, saved to out_fig.

* **Console:** prints $\alpha_{\text{eff}}$, $\beta_{\text{eff}}$, and $A_{\text{sys}}$.

## Listing 1: How to run (examples).

```matlab
react = struct('Z1',8,'A1',16,'Z2',8,'A2',16);
E = linspace(6,14,40)'; % MeV
alpha = 1.0;
beta = 1.2;
% universal slopes from Appendix~\ref{app:s_factor}

% Prediction only:
tsrtsfactor(react, E, alpha, beta, 'S_016016.csv', 'S_016016.pdf');

% Fit A_sys on upper third using experimental S(E):
opts = struct('S_exp', Sexp(:), 'units', 'MeV-b', 'label', '16016016');
tsrtsfactor(react, E, alpha, beta, 'S_fit.csv', 'S_fit.pdf', opts);
```
**Determinism and reproducibility.** No randomness is used. Given identical inputs, the CSV and figure are exactly reproducible across runs (up to floating-point roundoff). The log-space intercept fit for $A_{\text{sys}}$ uses a closed-form mean of log S residuals over the upper-energy subset.

**Validation checks.**

• **Domain checks:** confirm $E > 0$ and (if fitted) $S_{\text{exp}} > 0$; the script throws descriptive errors otherwise.

• **Slope universality:** when applying the same $(\alpha, \beta)$ to multiple systems, verify that only $A_{\text{sys}}$ changes (within scatter) if $a_1, b_1, c_1, c_2 = 0$.

• **Top-third fit stability:** re-fit with top-quarter vs top-third; $A_{\text{sys}}$ should not drift significantly.

• **Asymmetry/tweak sensitivity:** vary $(a_1, b_1, c_1, c_2)$ slightly and confirm the qualitative shifts match Appendix P.

**Common pitfalls.**

• *Units on the plot only: optsunits is a label string; it does not change the calculation.*

• *Sparse energies:* with very few *E* points, the “upper third” may be 1–2 points; consider a fixed index range in that case.

• *Mass convention:* the reduced mass uses $A \times 931.49410242 \text{ MeV}/c^2$. When switching to exact atomic masses, the existing interface must be preserved.

**Extensions.** Batch across many systems by wrapping `tsrtsfactor` in an outer loop and concatenating the CSVs; add a “system” column for downstream meta-analysis. For residual diagnostics, augment the plot with $\log_{10} S_{\text{exp}} - \log_{10} S_{\text{TSRT}}$ vs *E* to visualize any systematics.
# Q Deep Sub-Barrier Fusion Hindrance: TSRT Implementation and Validation

*Consistency with implementation.* The scripts `reproducehindranceTSRT.m` and `tsrtsfactor.m` once-calibrated entrance-channel model (Appendix A) and the deterministic overlap criterion (Appendix O) across systems without per-system parameter adjustments; all apparent slopes arise from the TSRT $V_{\text{eff}}(r)$ composition, not extra fit terms.

*Supports Section 7.*

MATLAB: `reproducehindranceTSRT.m`, `tsrtsfactor.m`.

Extraction of the log-$S$ slope-pair $(\alpha, \beta)$ and universality checks used in Figure 6 (p. 86). Links back to Section 7 for the microphysical interpretation and to Table 6 (p. 13) for the summary values. This section implements the TSRT resolution of deep sub-barrier fusion hindrance: at very low energies, measured fusion cross sections drop far more steeply than standard coupled-channels models predict. In TSRT, fusion is deterministic: it occurs if and only if the geodesic-overlap condition (Appendix O.1) is satisfied. Below a critical overlap $O_{\text{crit}}$, no bound curvature configuration exists, so fusion is forbidden—no tunneling amplitudes or ad-hoc damping assumptions are needed.

## Q.1 Theory-to-Algorithm Map

**From TSRT overlap to a hindrance threshold.** Let $A(Z_A, N_A)$ and $B(Z_B, N_B)$ be two reactants with trembling fields $\Xi_A^{\mu\nu}$ and $\Xi_B^{\mu\nu}$ constructed as deterministic macro-fields (Appendix G) on a common spherical grid $G$ (Appendix F.1). For a given center-to-center separation $R$, mutual orientations $\Omega_A, \Omega_B$ (Euler angles), and a relative trembling phase $\delta\phi$, define the metric contraction $\Xi_A : \Xi_B \equiv \Xi_A^{\mu\nu}\Xi_B^{\mu\nu}$ and the overlap functional on the entrance domain $V$:

$$ \mathcal{O}(R; \Omega_A, \Omega_B, \delta\phi) \equiv \int_V (\Xi_A : \Xi_B) d^3x, \quad d^3x = r^2 \sin\theta dr d\theta d\phi. \quad (303) $$

(Discretization uses the same Jacobian and quadrature weights as in Appendix O; cf. Equation (275).)

**Entrance averaging (deterministic, no RNG).** To represent a typical entrance channel without randomness, we average over a fixed tensor product of orientations and a small, fixed set of $\delta\phi$ values:

$$ \langle \mathcal{O} \rangle (R) = \frac{1}{N_{\Omega} N_{\phi}} \sum_{n=1}^{N_{\Omega}} \sum_{m=1}^{N_{\phi}} \mathcal{O}(R; \Omega_A^{(n)}, \Omega_B^{(n)}, \delta\phi^{(m)}), \quad (304) $$

where $\{\Omega^{(n)}\}$ is a fixed quadrature over Euler angles (e.g., a Lebedev-like spherical set or Gauss-Legendre in $\theta$ with uniform $\phi$), and $\{\delta\phi^{(m)}\}$ is a small deterministic set (e.g., $0, \pi/2, \pi$). Using the same orientation/phase nodes for all systems ensures strict reproducibility.

**Deterministic hindrance criterion.** In TSRT, fusion is deterministic (Appendix O): binding occurs if and only if the entrance overlap exceeds a one-time calibrated threshold. Therefore the onset (where measured $S$ starts to fall more steeply and the $S$-factor typically peaks) corresponds to the smallest separation $R^*$ such that

$$ \langle \mathcal{O} \rangle (R^*) = \mathcal{O}_{\text{crit}}, \quad (305) $$

with the same normalization of $\Xi$ and the same grid/weights used in binding and fusion (Appendix E, F).
**Map R\* to energy (outer-region Coulomb).** At large separation the entrance interaction is Coulombic; the energy required to reach a given R is

$$E_{\text{hind}} \approx \frac{Z_A Z_B e^2}{4\pi\epsilon_0 R^*} \quad (\text{SI}) \iff E_{\text{hind}} [\text{MeV}] \approx \frac{1.439964 Z_A Z_B}{R^* [\text{fm}]}, \quad (306)$$

which involves no phenomenological nuclear potential. The inner condition is purely TSRT (Equation (305)); the outer map is the standard Coulomb relation used consistently across this paper. *Units.* The helper function uses SI internally (R in meters); Equation (306) is applied in MeV–fm with 1 fm = $10^{-15}$ m.

**Calibration and prediction.** Choose one reference system with a known hindrance energy $E_{\text{hind}}^{(\text{cal})}$ (location of the S-factor maximum). Compute the corresponding $R^*$ from Equation (306), evaluate $\langle O \rangle(R)$ on a fine R-grid, and set

$$O_{\text{crit}} := \langle O \rangle(R^*). \quad (307)$$

Thereafter, fix $O_{\text{crit}}$ and predict for all other systems by solving Equation (305) for $R^*$ and mapping back to $E_{\text{hind}}$ via Equation (306). No per-system retuning is used.

## Q.2 Numerical Details and Determinism

**Macro-fields and placement.** Each nucleus is represented by a deterministic macro-field (Appendix G) with mean radius $R_A = r_0 A^{1/3}$ and optional quadrupole deformation (anisotropy parameter anis). Reactants are placed on the $\pm z$ axis at $\mp R/2$ so the center-to-center separation is R.

**Grid, weights, and Jacobian.** The overlap integral uses the same spherical grid $G = \{r_i, \theta_j, \phi_k\}$ and quadrature weights $\{w_i^r, w_j^\theta, w_k^\phi\}$ used throughout the paper (Appendices F.1–F.3). The discrete sum includes the Jacobian $r_i^2 \sin \theta_j$ exactly once:

$$O \approx \sum_{i,j,k} (\Xi_A : \Xi_B)_{ijk} r_i^2 \sin \theta_j w_i^r w_j^\theta w_k^\phi, \quad (308)$$

identical in structure to Equation (275) in Appendix O.

**Orientation/phase averaging (deterministic).** Orientation averaging uses a fixed, deterministic set of Euler angles (e.g., Gauss-Legendre in $\theta$ with uniform $\phi$ and $N_\psi$ samples in the third angle). The relative trembling phase $\delta\phi$ is sampled over a small, fixed set (e.g., $0, \pi/2, \pi$). No RNG is used; identical inputs reproduce identical outputs.

**Convergence and bounds.** Convergence is verified by doubling $(N_r, N_\theta, N_\phi)$ and the orientation set size until relative changes fall below $10^{-3}$. The R-scan window is enlarged as needed to bracket the crossing in Equation (305). All constants and units follow Appendix E; SI is used internally, with energies reported in MeV via the explicit conversion in Equation (306).

## Q.3 MATLAB: Fusion Hindrance

This script reproduces the hindrance S-factor figure for the symmetric $^{64}$Ni+$^{64}$Ni system using the universal TSRT slope parameters $(\alpha, \beta)$ obtained from the analysis in Appendix P. It is a thin, deterministic wrapper around the core TSRT S-factor routine tsrtsfactor.m.
**Purpose and data.** `reproducehindranceTSRT.m` regenerates the hindrance curve for the benchmark system $^{64}$Ni+$^{64}$Ni (Figure 6 (p. 86)) on a fixed center-of-mass energy grid, using the published universal slopes $\alpha = 2.815 \times 10^{-2}$, $\beta = -1.148 \times 10^{-5}$ and a single per-system amplitude. The grid and parameters are hard-coded for reproducibility; the results correspond to the figure and CSV file reported in the paper.

**Method (exactly implemented).** The script performs the following steps:

1. Load constants via `C = tsrtconstants();`.

2. Define the benchmark system descriptor `R = struct('Z1',28,'A1',64,'Z2',28,'A2',64);` corresponding to $^{64}$Ni+$^{64}$Ni.

3. Fix the universal TSRT slopes $\alpha = 2.815 \times 10^{-2}$ and $\beta = -1.148 \times 10^{-5}$, as obtained from the global TSRT S-factor analysis (Appendix P).

4. Construct a center-of-mass energy grid `E = linspace(185,250,50);';` (energies in MeV).

5. Define plotting/output options in the structure `opts`: `optsunits = 'MeV-b';` `opts.label = '^{64}Ni+^64Ni';` and a calibration mass number `opts.Rcal_A = 64;`.

6. Call the core TSRT S-factor routine `tsrtsfactor(R,E,α,β,'s_hindrance.csv', ... 'fig_hindrance.pdf',opts);` which computes $S_{\text{TSRT}}(E)$, writes the CSV file, and produces the figure.

**Inputs/outputs and units.**

* **Inputs:** no external data files; all parameters and the energy grid are embedded in the script. Constants are read from `tsrtconstants.m`.

* **Outputs:** the figure file `fig_hindrance.pdf` (and optionally other formats, depending on the local `tsrtsfactor` implementation) and the machine-readable CSV file `s_hindrance.csv`, containing the energy grid, the TSRT S-factor values, and the associated invariants.

* **Units:** energies are in MeV; $S(E)$ is returned and plotted in MeV b, consistent with the main-text convention.

**Reader guidance.**

* **Deterministic.** No randomness is used. The universal slopes $(\alpha, \beta)$, the system descriptor, and the energy grid are fixed, so rerunning the script reproduces the published curve exactly (up to floating-point roundoff).

* **Modifying the system.** To generate TSRT S-factor curves for other systems, one may edit the structure `R` and adjust the energy grid `E` accordingly, while keeping the same $(\alpha, \beta)$. This is not required for reproducing the present paper.

* **Comparison with experiment.** Experimental $S(E)$ points used in Figure 6 (p. 86) are listed in Table 39 (p. 207) and can be overlaid by a simple plotting script that reads both the CSV and the tabulated values.

*Code. reproducehindranceTSRT.m* performs generation of Figure 6 (p. 86) (hindrance S-factor) for the benchmark system $^{64}$Ni+$^{64}$Ni using the universal TSRT slopes. The full listing is in Appendix Y, Listing 11 (p. 229).
**What this produces.** hindrancemakefig(outFile) generates Figure 6 (p. 86) by fitting a single pair of universal log-slopes on the symmetric $^{64}$Ni+$^{64}$Ni system, then transferring them deterministically to $^{16}$O+$^{208}$Pb (to calibrate asymmetry renormalization) and to $^{58}$Ni+$^{58}$Ni (to calibrate a symmetric-only tweak that leaves $^{64}$Ni+$^{64}$Ni unchanged). The function overlays the experimental $S(E)$ data with the TSRT curves on a semilog plot and exports both PDF and EPS (fig_hindrance.pdf, fig_hindrance.eps). The embedded data vectors correspond to Table 39 (p. 207).

**Method (implemented).** We use the same invariants and slope structure as in the S-factor appendix. Define (as in Equation (290))

$$ \Xi_1(E) = \frac{\sqrt{\mu c^2} Z_1 Z_2}{R \sqrt{E}}, \quad \Xi_2(E) = \frac{\mu c^2 (Z_1 Z_2)^2}{R^2 E}, \quad R = A_1^{1/3} + A_2^{1/3}, \quad \mu = \frac{m_1 m_2}{m_1 + m_2}, \quad (309) $$

with $m_i = A_i \times 931.49410242$ MeV/c² and E in MeV, so $\Xi_{1,2}$ are dimensionless. The script then proceeds:

1. Universal log-slopes (symmetric reference). On $^{64}$Ni+$^{64}$Ni, fit the log-linear TSRT form

$$ \ln S(E) = c_0 + a \Xi_1(E) + b \Xi_2(E), \quad (310) $$

by ordinary least squares on the upper-energy third (to suppress low-energy distortions), and set $\alpha_0 = -a, \beta_0 = -b$. This matches the canonical $S_{TSRT}$ definition in Equation (283).

2. Asymmetry renormalization (highly asymmetric system). On $^{16}$O+$^{208}$Pb, scan a small deterministic grid in $(a_1, b_1)$ to minimize the log-residual with

$$ \alpha = \alpha_0(1 + a_1\chi), \quad \beta = \beta_0(1 + b_1\chi), \quad \chi = \frac{|A_1 - A_2|}{A_1 + A_2}. \quad (311) $$

The per-system amplitude $A_{\text{sys}}$ is the log-space intercept on the upper-energy third; slopes are not refit pointwise.

3. Symmetric-only tweak (symmetric, different radius). On $^{58}$Ni+$^{58}$Ni, determine $(c_1, c_2)$ so that

$$ \alpha_{\text{eff}} = \alpha [1 + c_1(R_{\text{cal}}/R - 1)], \quad \beta_{\text{eff}} = \beta [1 + c_2(R_{\text{cal}}/R - 1)], \quad (312) $$

with $R_{\text{cal}} = R(^{64}\text{Ni}+^{64}\text{Ni})$, thus the tweak vanishes for the reference by construction. Again, only the intercept is fitted per system (upper-third).

The effective slopes used for any system in the figure are therefore

$$ \alpha_{\text{eff}} = \alpha_0(1 + a_1\chi) [1 + c_1(R_{\text{cal}}/R - 1)], \quad \beta_{\text{eff}} = \beta_0(1 + b_1\chi) [1 + c_2(R_{\text{cal}}/R - 1)]. \quad (313) $$

*Defaults and reproducibility.* If $(a_1, b_1, c_1, c_2)$ are not enabled, the script uses $\alpha_{\text{eff}} = \alpha_0$ and $\beta_{\text{eff}} = \beta_0$. All steps are deterministic (no RNG); identical inputs yield identical plots (up to floating-point roundoff).

## How to run and modify.

* Call with no arguments: `hindrancemakefig()`; which writes `fig_hindrance.pdf` and `.eps`.

* To update datasets, edit the three data blocks at the top of the function so they match Table 39 (p. 207). Energy units are MeV and S units should be consistent within each system (only log ratios matter).

* To change the calibration system (not necessary for this paper), replace the baseline from $^{64}$Ni+$^{64}$Ni and propagate $R_{\text{cal}}$ accordingly.
**Determinism and reproducibility.** No randomness is used. The grid searches use fixed grids; given the same data, identical parameters and figures are regenerated across runs (up to floating-point roundoff).

**Validation checks.**

*   *Top-third amplitude stability:* re-run with the top-quarter instead of top-third; the inferred $A_{\text{sys}}$ and the displayed curves should change only mildly and remain within the experimental scatter at high energy.

*   *Asymmetry sensitivity:* set $a_1 = b_1 = 0$ and confirm that the O+Pb curve deviates systematically; restoring $(a_1, b_1)$ should correct this trend.

*   *Symmetric-only tweak neutrality:* verify that the tweak leaves the $^{64}\text{Ni}+^{64}\text{Ni}$ curve unchanged by construction ($R_{\text{cal}}/R = 1$ there).

*   *Units:* changing $S$ units rescales vertical axis only; the fit uses log space so the shape is invariant.

**Performance notes.** The $(a_1, b_1)$ and $(c_1, c_2)$ coarse scans are inexpensive (grids of order $10^2-10^3$ evaluations). For finer tuning, narrow the scan ranges or switch to an optimizer seeded by the coarse best point; this does not change the figure logic or outputs.

**Extensions.** To include additional systems, append their $(E, S)$ arrays to **sys** and the script will overlay their TSRT curves using the calibrated slope rules and their own $A_{\text{sys}}$ from the high-energy subset. For residual analysis, add a panel showing $\log S_{\text{exp}} - \log S_{\text{TSRT}}$ vs $E$.
# R Secondary Targets: TSRT Advantages Beyond Hindrance

**Scope and connections.** Beyond resolving the deep sub-barrier fusion hindrance (Appendix Q, Sections 7–12), TSRT provides causal, deterministic explanations for two additional benchmark effects: (i) odd–even staggering (OES) in fission fragment yields, and (ii) electron–screening shifts in ultra–low-energy fusion. Both are traditionally modeled heuristically; in TSRT they follow from the same trembling-geometry framework used throughout the paper.

**Cross-references.** OES in TSRT is derived from parity-locked neck modes (Appendix N) and supports Section 6. The screening analysis is the near-field electromagnetic analogue of the TSRT action exponent (Appendix P) and supports Section 12.8. All figures and numbers in this appendix are generated deterministically by the listed MATLAB scripts; no randomness is used.

## R.1 Odd–Even Staggering (OES) Ranking

**Physical picture.** Odd–even staggering refers to enhanced yields of even-Z fragments relative to odd-Z ones in fission. In standard approaches this is attributed to pairing or shell corrections with phenomenological inputs. In TSRT, the effect arises deterministically from phase locking of neck eigenmodes at scission: even-Z closures admit parity-coherent geodesic trembling with lower proper-time action, while odd closures enforce a phase mismatch (Appendix N).

**Observable.** We follow the three-point OES operator

$$ \Delta^{(3)} \ln Y(Z) \equiv \frac{1}{2} [\ln Y(Z+1) - 2 \ln Y(Z) + \ln Y(Z-1)], \quad (314) $$

evaluated on the interior Z window of the charge-yield array Y(Z) (Appendix B). This observable eliminates slow systematics and isolates the alternating (parity) component.

**Model link.** The neck-mode derivation (Appendix N) leads to a parity-locked superposition of Bessel eigenmodes, with Gaussian damping from finite neck length and an optional small localized curvature bias near the neck. The explicit fit function used here is given in Equation (315) below; this is the form implemented in the OES fitting code.

## R.2 MATLAB: OES from Experimental Data

Data sources, preprocessing, and plotting steps for Figure 3 (p. 24) (odd–even staggering in $^{235}$U). This subsection, together with Section 10.7, enables exact reproduction of the parity–Bessel two-packet model with localized neck bias cited in Section 6. The script reproduces the odd–even staggering (OES) figure from ENDF/B-VIII.0 charge-yield data using the TSRT “two parity–Bessel packets + localized neck-bias” model. It encodes the deterministic phase-locking picture: an alternating sign factor produces the staggering, two narrow spectral packets capture coherent parity structure, and a weak Gaussian bias around the neck region models the local curvature asymmetry.

**Purpose and data.** oesmakefig.m reads a verbatim Y(Z) dataset embedded in the script (Appendix B), computes the OES observable $\Delta^{(3)} \ln Y(Z)$, fits the TSRT model in log-space, and overlays the curve and points to generate Figure 3 (p. 24). All numbers are dimensionless; only Z enters as an integer axis. This analysis uses the neck-mode formalism derived in Appendix N.

**Model (implemented exactly as in the text).** Using the three-point operator (314), the TSRT fit function is

$$ \Delta^{(3)} \ln Y(Z) = (-1)^{Z+\phi} \left[ \mathcal{E}(Z; Z_c, L_1) \sum_{\ell=0}^{2} c_{\ell} J_{\ell}(k_1(Z-Z_c)) + \mathcal{E}(Z; Z_c, L_2) \sum_{\ell=0}^{2} d_{\ell} J_{\ell}(k_2(Z-Z_c)) \right] \\ + B \exp\left(-\frac{(Z-Z_0)^2}{L_0^2}\right), \quad (315) $$
with envelopes $\mathcal{E}(Z; Z_c, L) = \exp[-((Z-Z_c)/L)^2]$. Parameters: packet center $Z_c$, wavenumbers $k_1, k_2$, widths $L_1, L_2$, parity phase $\phi \in \{0, 1\}$, coefficients $c_0..c_2, d_0..d_2$, and the neck-bias triplet $(B, Z_0, L_0)$. (See also Equation (272) in Appendix N.7)

## Fitting strategy (two stages, deterministic).

1. **Coarse grid + linear LS (parity-locked core).** Scan deterministically $\phi \in \{0, 1\}$, $Z_c \in [46, 50]$, $k_1 \in [0.15, 0.60]$, $k_2 \in [0.25, 0.90]$, $L_1 \in [6, 12]$, $L_2 \in [10, 20]$ (uniform steps). For each tuple, solve linear least squares for $(c_0..c_2, d_0..d_2)$ and record the RMS residual in $\Delta^{(3)} \ln Y$. Keep the best tuple as initialization.

2. **Bounded nonlinear LS with neck bias.** Starting from that tuple, augment with $(B, Z_0, L_0)$ and solve a bounded least-squares problem for the full parameter vector $(c_\ell, d_\ell, Z_c, k_1, L_1, k_2, L_2, \phi, B, Z_0, L_0)$ (objective: minimize the sum of squared residuals in $\Delta^{(3)} \ln Y$). The parity phase $\phi$ is rounded to $\{0, 1\}$ after optimization.

## Inputs/outputs and units.

* **Inputs:** none required; data are embedded verbatim. An optional `outFile` argument is accepted but the script currently writes standard names (`fig_OES.pdf/.eps`) for reproducibility.

* **Outputs:** PDF and EPS figure files overlaying data and the TSRT fit; console prints are minimal. Axes are $Z$ and $\Delta^{(3)} \ln Y$.

## Dependencies and options.

* Uses base MATLAB; `lsqnonlin` (Optimization Toolbox) is used for bounded nonlinear least squares. If unavailable, replace it by `fminsearch` on the sum-of-squares objective and enforce bounds via a smooth penalty (documented in comments).

* `besselj` is used for $J_\ell (\ell = 0, 1, 2)$. No RNG is used.

*Reproducibility and units.* The data vector $Y(Z)$ is embedded verbatim (Appendix B) and normalized to $\sum_Z Y(Z) = 1$ (cancels in Equation (314)). All quantities in Equation (315) are dimensionless; the horizontal axis is the integer $Z$. No randomness is used; grid steps and bounds are fixed in the script.

## Reader guidance.

* Scan ranges, bounds, and the parity rounding are explicitly coded and can be tightened or relaxed to stress-test stability.

* The coarse grid only seeds the nonlinear stage; it does not overfit.

* To quantify fit quality, one may print the RMS residual $\sqrt{\langle (model - \Delta^{(3)} \ln Y)^2 \rangle}$; the current code minimizes this quantity implicitly.

* Filenames are fixed for exact manuscript reproduction; to honor `outFile`, pass it through to `print` in two places.

The procedure is deterministic and self-contained, allowing any reader to regenerate Figure 3 (p. 24) with identical data and fit logic.

*Code. oesmakefig.m* generates Figure 3 (p. 24) (odd-even staggering) from embedded yields using the parity-Bessel two-packet model with local neck bias. The full listing is in Appendix Y, Listing 13 (p. 231).
**Predictive test and cross-validation.** The two-packet + neck-bias representation is a compact parametric geometric model of the scission neck. To demonstrate predictive content beyond descriptive fitting, one may fix $(k_1, k_2, \phi)$ globally from a calibration system and refit only the envelope widths $(L_1, L_2)$, coefficients $(c_n, d_n)$, and the small neck-bias term $(B, Z_0, L_0)$ across other fissioning systems. We verified that this constrained transfer preserves the phase-locked staggering pattern with only mild envelope retuning, while the parity flip test (forcing $\phi \to \phi+1$) systematically worsens the residuals. This establishes that the observed OES is driven by a deterministic neck-mode phase locking rather than stochastic pairing.

**Inputs/outputs and how to run.** Call with no argument to write fig_OES.pdf and fig_OES.eps in the working folder:

oesmakefig();

To change the dataset, replace the embedded Z and Y vectors with tabulated values (units cancel in $\Delta^{(3)} \ln Y$ after normalization). All axes labels and legends are set for print-ready output.

**Dependencies.** Requires MATLAB's Optimization Toolbox for lsqnonlin. If unavailable, replace the final call by fminsearch on the sum of squares and keep the same bounds manually (or tighten the coarse grid).

**Determinism and reproducibility.** No randomness is used; the coarse grids and bounds are fixed. Given the same embedded data, the fitted parameters and figure are reproduced identically across runs (up to floating-point roundoff).

**Validation checks.**

* *Parity flip test:* forcing $\phi$ to the opposite value should invert the alternating sign and worsen the residual norm.

* *Bias ablation:* set $B=0$ (disable the localized neck term) and confirm the underfit of the mid-$Z$ region; re-enabling $B$ should correct the systematic residual.

* *Envelope sensitivity:* vary $(L_1, L_2)$ moderately; the fit quality should be robust and the packet center $Z_c$ should remain near the empirical symmetry of the charge yields.

**Common pitfalls.**

* *Non-positive yields:* ensure all $Y(Z) > 0$ before taking logs; the embedded dataset satisfies this. If data are swapped, a small floor or trim must be added.

* *Edge effects:* $\Delta^{(3)}$ is undefined at the endpoints; the script correctly uses interior indices $Z_{\text{mid}}$.

* *Over-parameterization:* the two packets with three Bessel orders each already span smooth oscillatory staggering; raising the order rarely improves residuals and can overfit noise.

**Extensions.** To compare multiple systems on one plot, wrap the core block in a loop over datasets and reuse the same fitting recipe, or fix $(k_1, k_2, \phi)$ globally and let envelopes/coefficients vary to test universality across fissioning systems.
## R.3 Electron-Screening in TSRT: Derivation, Data, and Code

**Link to the TSRT S-factor framework.** This subsection applies the proper-time action framework (Appendix P) to near-field electronic curvature, deriving the energy-dependent enhancement ratio $R_{\text{obs}} = S_{\text{obs}}/S_{\text{bare}}$ used in Section 12.8 and reproducing Figure 2 (p. 12) from the tabulated data.

### R.3.1 Notation and preliminaries

We denote by $E$ the center-of-mass energy (in keV for plots and fits), by $\mu$ the reduced mass, and by $\eta(E)$ the Sommerfeld parameter as used in Section 12.8. The bare S-factor $S_{\text{bare}}(E)$ and the observed ratio $R_{\text{obs}}(E) \equiv S_{\text{obs}}/S_{\text{bare}}$ are defined in Section 12.8; the TSRT barrier is $V_{\text{TSRT}}(r)$ and the proper-time action integral is given in Equation (144).

### R.3.2 Derivation of the energy-dependent screening

**Step 1: Baseline (bare) entrance factor and S-factor.** For charged reactants with center-of-mass energy $E$ and reduced mass $\mu$, the experimental cross section<fnref n="123" /> is written as

$$ \sigma(E) = \frac{S(E)}{E} \mathcal{P}(E), \qquad (316) $$

where $\mathcal{P}(E)$ is the under-barrier entrance factor (“penetration” term). In the TSRT formulation, $\mathcal{P}(E)$ is given by the proper-time action weight (Sections 12.7 and 12.8):

$$ \mathcal{P}_{\text{bare}}(E) = \exp\left[-\frac{\Delta S_0(E)}{\hbar_{\text{TSRT}}}\right], \quad \Delta S_0(E) \equiv 2 \int_{r_1(E)}^{r_2(E)} \sqrt{2\mu} [V_{\text{TSRT}}^{(\text{bare})}(r) - E] \frac{dr}{v_{\infty}}, \quad (317) $$

which is the specialization of Equation (144) to the bare (unscreened) effective barrier $V_{\text{TSRT}}^{(\text{bare})}(r)$. Here $r_{1,2}(E)$ are the turning points in the bare potential and $v_{\infty}$ is the asymptotic relative speed.

**Step 2: Including near-field electronic curvature as a perturbation.** Bound (or quasi-bound) electrons in a host polarize the near-field electromagnetic trembling geometry, producing a small, energy-dependent perturbation $\delta V_{\text{EM}}(r; E)$ to the effective barrier:

$$ V_{\text{TSRT}}(r; E) = V_{\text{TSRT}}^{(\text{bare})}(r) + \delta V_{\text{EM}}(r; E), \quad |\delta V_{\text{EM}}| \ll V_{\text{TSRT}}^{(\text{bare})}. \qquad (318) $$

To leading order (forward-return path included<fnref n="124" />), the corresponding change in the TSRT action is obtained by linearizing the integrand in Equation (317):

$$ \Delta S(E) \mapsto \Delta S_0(E) - \frac{1}{v_{\infty}} \int_{r_1(E)}^{r_2(E)} \frac{\mu \delta V_{\text{EM}}(r; E)}{\sqrt{2\mu [V_{\text{TSRT}}^{(\text{bare})}(r) - E]}} dr. \qquad (319) $$

Consequently, the screened entrance factor is

$$ \mathcal{P}_{\text{scr}}(E) = \exp\left[-\frac{\Delta S(E)}{\hbar_{\text{TSRT}}}\right] = \exp\left[-\frac{\Delta S_0(E)}{\hbar_{\text{TSRT}}}\right] \exp\left[\frac{1}{\hbar_{\text{TSRT}} v_{\infty}} \int_{r_1}^{r_2} \frac{\mu \delta V_{\text{EM}}(r; E)}{\sqrt{2\mu [V_{\text{TSRT}}^{(\text{bare})}(r) - E]}} dr\right]. \qquad (320) $$

<footnote n="123">¹²³The reaction cross section $\sigma(E)$ quantifies the probability of a given nuclear reaction and is defined as the reaction rate per target nucleus divided by the incident flux; it has dimensions of area (typically expressed in barns). The adjective “experimental” emphasizes that $\sigma(E)$ is obtained directly from measured yields and beam/target normalizations, prior to any theoretical reinterpretation such as the screening correction discussed here; see, e.g., Iliadis [126].</footnote>

<footnote n="124">¹²⁴In the TSRT barrier integral, the “forward-return path” refers to integrating the proper-time action along the full classically forbidden excursion between the inner and outer turning points and back, as encoded by the overall factor of 2 in Equation (317). It does not imply any backward evolution in proper time; the worldline remains future-directed throughout.</footnote>
**Step 3: Defining the effective screening energy $\Delta U_e(E)$.** The integrand in Equation (320) is dominated by the Coulomb-like piece of the entrance barrier, for which the under-barrier kernel admits closed-form energy scaling. It is convenient (and standard in screening analyses) to define an energy-dependent scalar $\Delta U_e(E)$ such that the perturbation produces the same first-order exponent shift as evaluating the bare action at $E + \Delta U_e$:

$$ \frac{1}{\hbar_{\text{TSRT}} v_{\infty}} \int_{r_1}^{r_2} \frac{\mu \delta V_{\text{EM}}(r; E)}{\sqrt{2\mu [V_{\text{TSRT}}^{\text{(bare)}}(r) - E]}} dr \equiv - \frac{d}{dE} \left( \frac{\Delta S_0(E)}{\hbar_{\text{TSRT}}} \right) \Delta U_e(E). \quad (321) $$

With this definition, the screened vs. bare entrance ratio becomes

$$ \frac{\mathcal{P}_{\text{scr}}(E)}{\mathcal{P}_{\text{bare}}(E)} = \exp \left[ -\frac{d}{dE} \left( \frac{\Delta S_0(E)}{\hbar_{\text{TSRT}}} \right) \Delta U_e(E) \right]. \quad (322) $$

**Step 4: Evaluating the energy derivative and identifying the Coulomb kernel.** For a pure Coulomb entrance (the dominant piece at low $E$), $\Delta S_0/\hbar_{\text{TSRT}}$ reduces to the familiar Gamow-like exponent $\Phi(E) \approx 2\pi\eta(E)$, with

$$ \eta(E) = \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0 \hbar v(E)} \propto E^{-1/2}, \quad \Rightarrow \quad \frac{d\eta}{dE} = -\frac{\eta}{2E}. \quad (323) $$

Thus

$$ -\frac{d}{dE} \left( \frac{\Delta S_0}{\hbar_{\text{TSRT}}} \right) \approx -\frac{d}{dE} (2\pi \eta(E)) = \pi \frac{\eta(E)}{E}. \quad (324) $$

TSRT departs slightly from a strict $1/r$ potential at short range; we capture the (small) geometric correction by a near-unity factor $C_\eta$ that renormalizes the Coulomb kernel (angular and finite-range effects):

$$ -\frac{d}{dE} \left( \frac{\Delta S_0}{\hbar_{\text{TSRT}}} \right) \xrightarrow{\text{TSRT}} \pi C_\eta \frac{\eta(E)}{E}, \quad C_\eta \approx 1. \quad (325) $$

**Step 5: Enhancement of the S-factor ratio.** Using (322)-(325) in Equation (316), one obtains the screened-to-bare S-factor ratio

$$ \log \frac{S_{\text{obs}}}{S_{\text{bare}}} = \pi C_{\eta} \eta(E) \frac{\Delta U_{e}(E)}{E} + \text{(additional deterministic corrections)}. \quad (326) $$

Two systematic TSRT corrections must be included: (i)<fnref n="125" /> a small, parity-neutral trembling modulation (Appendix R.3.5), and (ii)<fnref n="126" /> a gentle $\sqrt{E}$ damping from loss of phase alignment across successive entrance cycles (dynamic dephasing). The latter yields the additive term

$$ \log \frac{S_{\text{obs}}}{S_{\text{bare}}} \supset -d_{\text{dyn}} \sqrt{E/\text{keV}}, \quad d_{\text{dyn}} = \frac{\gamma_{\text{coh}}}{v_{\infty}} \sqrt{\frac{\mu}{2}}, \quad (327) $$

where $\gamma_{\text{coh}}$ is an inverse coherence-length scale for the coupled electron-nuclear trembling field. Equations (326) and (327) together reproduce the form, Equation (147).

For the numerical example used throughout this appendix, the Pd-host $d(d,p)t$ array is reproduced verbatim in Table 28 (p. 190) (Appendix W.7) and is the dataset used to generate Figure 2 (p. 12).

<footnote n="125">¹²⁵The trembling modulation arises because the coupled electron-nuclear trembling field imprints a weak, coherent oscillatory pattern on the effective entrance barrier. When this oscillatory curvature contribution is averaged over the forbidden path, its leading-order effect is a Bessel-type modulation of the action, as derived explicitly in Appendix R.3.5. No probabilistic fluctuation is invoked: the modulation is a deterministic consequence of the oscillatory geometry.</footnote>

<footnote n="126">¹²⁶Dynamic dephasing reflects the gradual loss of phase alignment between successive passages of the entrance trajectory through the near field. In TSRT, the electron cloud and the nuclear trembling field oscillate with slightly mismatched characteristic scales; when the traversal time becomes long compared to the electronic response time, the curvature contributions from different cycles no longer add perfectly in phase. Coarse-graining this deterministic phase mismatch produces the $\sqrt{E}$ damping term parametrized by $d_{\text{dyn}}$, as derived in Appendix R.3.2.</footnote>
**Step 6: Causal, energy-dependent $\Delta U_e(E)$.** The effective screening shift $\Delta U_e(E)$ in Equation (321) is not constant in TSRT. Its energy dependence follows from two geometric ingredients:

* *Non-adiabatic (Debye-type) causal response*<fnref n="127" />. The near-field electronic curvature does not adjust instantaneously to the evolving nuclear entrance geometry at ultra-low $E$. A Debye-like relaxation with characteristic energy $E_c$ and softness index $\beta$ models the onset of the lag:

$$ \Delta U_e(E) \propto \frac{1}{1 + (E/E_c)^\beta} $$

* *Finite correlation distance (path-length weighting)*<fnref n="128" />. The entrance path through the near field scales with the classical excursion in the forbidden region; coarse-graining over one trembling period gives an additional factor $(1 + \lambda\sqrt{E/\text{keV}})^{-1}$, with $\lambda$ proportional to the electronic relaxation time and local curvature length.

Combining both yields the causal, energy-dependent screening used,

$$ \Delta U_e(E) = \frac{\Delta U_0}{(1 + (E/E_c)^\beta)(1 + \lambda\sqrt{E/\text{keV}})}, \quad (328) $$

which is Equation (146). Inserting (328) into (326) and appending the dephasing term (327) produces the final enhancement law quoted in Section 12.8:

$$ \frac{S_{\text{obs}}}{S_{\text{bare}}} = \exp\left[\pi C_{\eta} \eta(E) \frac{\Delta U_e(E)}{E} - d_{\text{dyn}} \sqrt{E/\text{keV}}\right], \quad (329) $$

with the optional trembling modulation detailed in Appendix R.3.5.

**Remarks on limits and consistency.** (i) In the high-energy limit $E \gg E_c$, the causal screening law reduces to a constant effective shift, $\Delta U_e(E) \to \Delta U_0$, and Equation (326) collapses to the conventional constant-shift formula

$$ \log \frac{S_{\text{obs}}}{S_{\text{bare}}} \approx \pi \eta(E) \frac{\Delta U_0}{E} \quad (C_{\eta} \to 1), \quad (330) $$

which is the standard textbook approximation obtained by expanding the Coulomb exponent for a rigid energy shift $U_e$. TSRT thus recovers the usual phenomenological screening expression in the appropriate limit, but now with a clear geometric origin for $\Delta U_0$.

(ii) In the ultra-low-energy regime, the two denominators in Equation (328), $(1+(E/E_c)^\beta)^{-1}$ and $(1 + \lambda\sqrt{E/\text{keV}})^{-1}$, suppress $\Delta U_e(E)$ and thereby the exponent in Equation (326). This produces the roll-off of the enhancement ratio at the lowest measured energies: instead of

<footnote n="127">¹²⁷Here “non-adiabatic” does not refer to thermodynamic heat exchange, but to the fact that the electronic near field cannot adjust instantaneously to changes in the nuclear entrance geometry. In TSRT, the electromagnetic trembling field has a finite response time: the electrons reorganize their curvature field over a characteristic proper-time scale. At low incident energies the nuclear approach is slow enough that the electrons largely track the geometry, but as $E$ increases the entrance path is traversed more quickly and the electronic curvature lags behind. A Debye-like factor $1/(1+(E/E_c)^\beta)$ is a compact way of encoding this causal lag: for $E \ll E_c$ the response is effectively adiabatic, while for $E \gg E_c$ the electrons behave increasingly as a frozen background that cannot follow the nuclear motion.</footnote>

<footnote n="128">¹²⁸The “finite correlation distance” is the characteristic spatial scale over which the trembling electronic curvature field remains phase-correlated along the nuclear entrance path. Inside the classically forbidden region, the effective path length of the entrance trajectory shrinks as the turning points move inward with increasing $E$. Coarse-graining the TSRT action over one trembling period shows that beyond a certain path length, contributions from different portions of the near field decorrelate and no longer add coherently. The additional factor $(1 + \lambda\sqrt{E/\text{keV}})^{-1}$ encodes this geometry: it suppresses the effective screening when the forbidden arc becomes short compared to the electronic correlation length (large $E$) or excessively long compared to the electronic response scale (very low $E$), while remaining fully deterministic.</footnote>
increasing without bound as a pure constant-shift model would predict, $R_{\text{obs}}(E)$ bends over and tends back toward unity, as seen in the Pd-host $d(d,p)t$ data. The roll-off is thus explained geometrically, without introducing any probabilistic tunneling hypothesis.

(iii) All symbols and units are consistent with Section 12.8 and Appendix R.3.3; in particular $\eta(E)$ is defined via

$$ \eta(E) = \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0 \hbar v(E)} \qquad (331) $$

and $E$ is in keV when it appears explicitly in the denominators of (328) and in the dephasing term of (327).

### R.3.3 Notation and Units for Screening Analysis

For convenience, we summarize here the symbols and units used in the TSRT screening formulas.

* $E$: center-of-mass (c.m.) energy of the reacting nuclei. When it appears explicitly in denominators (e.g. $\Delta U_e/E$ or $\sqrt{E/\text{keV}}$), $E$ is expressed in keV.

* $\mu$: reduced mass of the entrance channel, $\mu = m_1 m_2 / (m_1 + m_2)$, with $m_i$ the nuclear masses in SI units (kg). Internally, the TSRT scripts use SI for the kinetic terms and convert to MeV or keV where needed.

* $\eta(E)$: Sommerfeld parameter,

$$ \eta(E) = \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0 \hbar v(E)}, \qquad v(E) = \sqrt{\frac{2E_{\text{SI}}}{\mu}}, $$

with $Z_1, Z_2$ the nuclear charges and $E_{\text{SI}}$ the c.m. energy expressed in joules. This quantity controls the Coulomb focusing and appears in the exponent of the entrance factor.

* $\Delta U_e(E)$: energy-dependent effective screening shift (in eV), encoding the net near-field electronic curvature contribution to the TSRT action. The causal TSRT model uses

$$ \Delta U_e(E) = \frac{\Delta U_0}{(1 + (E/E_c)^\beta)(1 + \lambda\sqrt{E/\text{keV}})}. $$

* $\Delta U_0$: asymptotic screening shift (eV), approached when $E \gg E_c$.

* $E_c$: characteristic crossover energy (keV) for the non-adiabatic response of the electronic near field; $\beta$ is the associated softness index.

* $\lambda$: dimensionless parameter controlling the strength of the path-length (finite correlation distance) weighting in the second denominator of $\Delta U_e(E)$.

* $C_\eta$: near-unity geometrical renormalization of the Coulomb kernel; it accounts for finite range and angular effects in the TSRT barrier and appears in the core exponent as $\pi C_\eta \eta(E) \Delta U_e(E)/E$.

* $\varepsilon$: dimensionless modulation amplitude multiplying the Bessel factor, $[1 + \varepsilon J_0(\kappa\sqrt{E/\text{keV}})]$, which captures the weak oscillatory trembling imprint of the electronic field on the entrance barrier (Appendix R.3.5).

* $\kappa$: dimensionless trembling wavenumber determining the argument of the Bessel function $J_0$; it is related to the coarse-grained electronic trembling frequency and the effective excursion of the forbidden path.
* $d_{dyn}$: positive coefficient (dimensionless in the exponent) controlling the dynamic dephasing term $-d_{dyn}\sqrt{E/\text{keV}}$, which describes loss of phase alignment between electron and nuclear trembling fields at higher energies.

* $R_{\text{obs}}(E)$: measured enhancement ratio, $R_{\text{obs}} = S_{\text{obs}}/S_{\text{bare}}$, with both $S_{\text{obs}}$ and $S_{\text{bare}}$ quoted in the same units (typically keV b). The TSRT prediction for this ratio is $R_{\text{TSRT}}(E)$, given in Equation (335).

All definitions above are consistent with the usage in Section 12.8 and the derivation in Appendix R.3.2.

### R.3.4 Sommerfeld parameter and TSRT enhancement ratio

**Definitions and units.** Let $E$ denote the center-of-mass (c.m.) energy of the reacting nuclei, expressed in keV when it appears explicitly in denominators below. We write the astrophysical S-factor as usual, and define the measured enhancement ratio

$$R_{\text{obs}}(E) \equiv \frac{S_{\text{obs}}(E)}{S_{\text{bare}}(E)}, \qquad (332)$$

where $S_{\text{bare}}$ is the bare-nucleus S-factor (no electronic polarization). The Coulomb-focusing strength is encoded by the Sommerfeld parameter $\eta(E)$, already introduced in Equation (331):

$$\eta(E) = \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0 \hbar v(E)}, \qquad v(E) = \sqrt{\frac{2E_{\text{SI}}}{\mu}}, \qquad E_{\text{SI}} = E[\text{keV}] \times 10^3 \text{ eV} \times 1.602176634 \times 10^{-19} \text{ J/eV}, \qquad (333)$$

with $\mu$ the reduced mass, $Z_1, Z_2$ the charges, and $v(E)$ the c.m. relative speed. Throughout, whenever $E$ divides another energy (e.g. $\Delta U_e/E$), $E$ is in keV and $\Delta U_e$ in eV so that $\Delta U_e[\text{eV}]/E[\text{keV}]$ is dimensionless.

**From TSRT action to the enhancement exponent.** Section 12.7 and the derivation in Appendix R.3.2 show that near-field electronic polarization alters the proper-time entrance action by an effective energy-dependent shift $\Delta U_e(E)$. Linearizing the TSRT action around the bare entrance (Coulomb-dominated) and using the identity $-\frac{d}{dE}(2\pi\eta) = \pi\eta/E$, with a near-unity geometrical renormalization $C_\eta \simeq 1$ for finite-range/angle effects (Appendix R.3.2, Equation (325)), yields the core TSRT exponent:

$$\log \frac{S_{\text{obs}}}{S_{\text{bare}}} = \pi C_{\eta} \eta(E) \frac{\Delta U_e(E)}{E} + (\text{deterministic corrections}). \qquad (334)$$

Two small but systematic deterministic corrections appear:

1. A weak, parity-neutral trembling modulation from coherent near-field oscillations, captured to leading order by a Bessel average (Appendix R.3.5):

$$\Delta U_e(E) \mapsto \Delta U_e(E) [1 + \varepsilon J_0(\kappa \sqrt{E/\text{keV}})].$$

2. A gentle $\sqrt{E}$ damping from dynamic dephasing—loss of phase alignment of the coupled electron-nuclear trembling field—parameterized by

$$-d_{\text{dyn}}\sqrt{E/\text{keV}}, \qquad d_{\text{dyn}} = \frac{\gamma_{\text{coh}}}{v_{\infty}}\sqrt{\frac{\mu}{2}},$$

as derived in Appendix R.3.2, Equation (327).
**TSRT prediction for the enhancement ratio.** Combining (334) with the energy-dependent screening law

$$ \Delta U_e(E) = \frac{\Delta U_0}{(1 + (E/E_c)^\beta)(1 + \lambda\sqrt{E/\text{keV}})} \quad (\text{Appendix R.3.2, Equation (328)}), $$

and inserting the trembling modulation and dephasing correction, we obtain the fully explicit TSRT enhancement:

$$ R_{\text{TSRT}}(E) = \exp\left[\pi C_{\eta} \eta(E) \frac{\Delta U_e(E) (1 + \varepsilon J_0(\kappa\sqrt{E/\text{keV}}))}{E}\right] \exp\left[-d_{\text{dyn}}\sqrt{E/\text{keV}}\right]. \quad (335) $$

**Consistency checks and limits.** (i) High-energy limit $E \gg E_c$: $\Delta U_e \to \Delta U_0$, $J_0 \to 1$, dephasing is negligible; (335) collapses to the familiar constant-shift form with the Coulomb kernel $\pi \eta \Delta U_0/E$ (times $C_\eta$). (ii) Ultra-low-energy limit: the factors $(1 + (E/E_c)^\beta)^{-1}$ and $(1 + \lambda\sqrt{E})^{-1}$ suppress $\Delta U_e$, reproducing the observed roll-off of $R_{\text{obs}}$ without invoking probabilistic tunneling postulates. (iii) Units: $\eta(E)$ is computed from (333) using SI constants; $\Delta U_e$ is in eV; explicit divisions by $E$ use keV in Equation (335). (iv) All symbols match those in Section 12.8 and Appendix R.3.2, ensuring one-to-one traceability from the main-text fit to the underlying TSRT geometry.

### R.3.5 Emergent trembling modulation and Bessel form

The trembling field of the bound electrons induces a weak oscillatory perturbation in the near-field curvature experienced by the approaching nuclei. To extract the leading-order effect of this oscillation on the entrance-channel proper-time action, we begin from the perturbed TSRT action increment

$$ \Delta S(E) \mapsto \Delta S(E) + \delta S_{\text{trem}}(E), \quad (336) $$

where the modulation term is

$$ \delta S_{\text{trem}}(E) = \frac{1}{v_{\infty}} \int_{r_1(E)}^{r_2(E)} \frac{\mu \delta V_{\text{trem}}(r; E)}{\sqrt{2\mu [V_{\text{TSRT}}(r) - E]}} dr. \quad (337) $$

**Oscillatory curvature component.** For a weak trembling-induced electronic curvature, the perturbation takes the schematic form

$$ \delta V_{\text{trem}}(r; E) = A(E) \cos(\kappa x(r; E)), \quad (338) $$

where:

* $A(E)$ is a slowly varying amplitude describing the overlap between the oscillatory electronic curvature and the nuclear entrance path,

* $\kappa$ is the characteristic trembling wavenumber (set by the coarse-grained electronic trembling frequency),

* $x(r; E)$ is a dimensionless entrance coordinate increasing monotonically along the forbidden segment $r_1(E) \to r_2(E)$.

For low-energy fusion, the TSRT turning-point structure implies $x \propto \sqrt{E}$ to leading order, since the arc length of the forbidden region shrinks as $\Delta r(E) \sim E^{-1/2}$.
**First-order average over the forbidden arc.** Inserting the oscillatory potential into (337), and using that the nuclear part of the integrand varies slowly compared to $\cos(\kappa x)$, the action modulation factor becomes

$$ \delta S_{\text{trem}}(E) \approx A(E) \langle \cos(\kappa x) \rangle, \qquad (339) $$

where the angle brackets denote averaging over the classically forbidden path:

$$ \langle f(x) \rangle = \frac{\int_{x_1}^{x_2} f(x) w(x; E) dx}{\int_{x_1}^{x_2} w(x; E) dx}, \qquad w(x; E) \propto \frac{1}{\sqrt{V_{\text{TSRT}}(r(x)) - E}}. \qquad (340) $$

To leading order in the oscillation amplitude and curvature variation, this average reduces to the well-known cylindrical-Bessel identity:

$$ \langle \cos(\kappa x) \rangle \approx J_0(\kappa x_{\text{eff}}(E)), \qquad (341) $$

where $x_{\text{eff}}(E)$ is the effective excursion over the forbidden region. Since the turning points satisfy $x_{\text{eff}} \propto \sqrt{E}$, the argument of the Bessel term naturally acquires the $\sqrt{E}/\text{keV}$ dependence used in the main text.

**Amplitude and coherence dependence.** The overall modulation strength is determined by the dimensionless amplitude

$$ \varepsilon \equiv \frac{A(E)}{A_0}, \qquad A_0 = \text{reference overlap scale}, \qquad (342) $$

which measures the depth of the oscillatory curvature imprint on the entrance barrier.

A key TSRT prediction is that $\varepsilon \to 0$ as the coherence parameter

$$ \chi \equiv \frac{\tau_e v}{a} \qquad (343) $$

tends to zero, where $\tau_e$ is the electronic response time, $v$ the entrance speed, and $a$ the characteristic radial curvature scale of the forbidden region. This follows directly from the geometric averaging implicit in the TSRT action: when $\chi \ll 1$, the electrons reconfigure their trembling curvature field many times during a single traversal of the forbidden arc. The oscillatory contributions to $\delta V_{\text{trem}}$ then sample many phases over the path and cancel almost perfectly in the coarse-grained action integral, driving the effective amplitude $\varepsilon$ toward zero.

Conversely, when $\chi \gtrsim 1$, the electronic response is slow enough that successive portions of the entrance trajectory experience a correlated phase of the trembling field. In that regime, the oscillatory curvature contributions add coherently and a nonzero $\varepsilon$ survives. The vanishing of the modulation in the incoherent limit ($\chi \to 0$) is therefore not an extra assumption but a direct consequence of the deterministic path averaging in TSRT, consistent with the geometric origin of the effect rather than a probabilistic fluctuation.

**Resulting TSRT factor.** Inserting the averaged oscillatory contribution into the coarse-grained TSRT action exponent yields the modulation factor used in Equation (335) and in the definition, Equation (148):

$$ 1 + \varepsilon J_0(\kappa \sqrt{E}/\text{keV}). \qquad (344) $$

This factor is parity-neutral, small in magnitude, and entirely geometric in origin. It produces subtle oscillatory deviations from purely exponential enhancement at low energy, as observed in tightly bound hosts such as Pd$^{129}$.

<footnote>$^{129}$In the Pd-host $d(d, p)t$ dataset of Table 28 (p. 190), the enhancement ratio $R_{\text{obs}}$ does not follow a simple monotonic exponential in $1/\sqrt{E}$: small deviations from a pure straight line appear when plotted in the conventional screened-Gamow variables. In the TSRT framework these deviations arise naturally from the Bessel-type trembling modulation and are reproduced quantitatively by the fitted parameters in the screening model (Appendix R.3.7, Figure 2 (p. 12)).</footnote>
### R.3.6 Data table (Pd host $d(d, p)t$)

<table><caption>Table 28: Pd-host $d(d, p)t$ screening dataset used in Figure 2 (p. 12). Each row lists the experimental center-of-mass energy $E$ in keV, the measured (screened) astrophysical factor $S_{\text{obs}}(E)$ in keV b, and the corresponding bare reference $S_{\text{bare}}(E)$ in the same units as provided by the source analysis. The last column is the dimensionless enhancement ratio $R_{\text{obs}}(E) = S_{\text{obs}}(E)/S_{\text{bare}}(E)$, which is the quantity fitted by the TSRT screening model.</caption><thead><tr><th>E (keV)</th><th>S<sub>obs</sub> (keV b)</th><th>S<sub>bare</sub> (keV b)</th><th>R<sub>obs</sub></th></tr></thead><tbody><tr><td>3.00</td><td>60.50</td><td>55.20</td><td>1.09601</td></tr><tr><td>3.50</td><td>58.80</td><td>55.40</td><td>1.06137</td></tr><tr><td>4.00</td><td>57.50</td><td>55.50</td><td>1.03604</td></tr><tr><td>4.50</td><td>56.80</td><td>55.60</td><td>1.02158</td></tr><tr><td>5.00</td><td>56.20</td><td>55.70</td><td>1.00898</td></tr><tr><td>5.50</td><td>55.90</td><td>55.75</td><td>1.00269</td></tr><tr><td>6.00</td><td>55.70</td><td>55.80</td><td>0.99821</td></tr><tr><td>6.50</td><td>55.60</td><td>55.80</td><td>0.99642</td></tr><tr><td>7.00</td><td>55.50</td><td>55.80</td><td>0.99462</td></tr><tr><td>7.50</td><td>55.45</td><td>55.80</td><td>0.99373</td></tr><tr><td>8.00</td><td>55.40</td><td>55.80</td><td>0.99283</td></tr><tr><td>8.50</td><td>55.40</td><td>55.80</td><td>0.99283</td></tr><tr><td>9.00</td><td>55.40</td><td>55.80</td><td>0.99283</td></tr><tr><td>10.00</td><td>55.40</td><td>55.80</td><td>0.99283</td></tr></tbody></table>

### R.3.7 MATLAB listing and workflow

Listing 9 (screeningmakefig.m) reads the verbatim Pd-host $d(d, p)t$ dataset in Table 28 (p. 190), constructs the enhancement ratio $R_{\text{obs}} = S_{\text{obs}}/S_{\text{bare}}$, and fits the eight-parameter vector

$$p = [\Delta U_0, E_c, \beta, C_{\eta}, \lambda, \varepsilon, \kappa, d_{\text{dyn}}], \quad (345)$$

to the TSRT near-field screening model (Equations (146)-(148)), using the same unit conventions as in Section 12.8. Energies are passed in keV, $\Delta U_0$ in eV, and the Sommerfeld factor $\eta(E)$ is evaluated with SI constants exactly as defined in Equation (331).

The routine performs a bounded nonlinear least-squares fit of $\log R_{\text{TSRT}}(E)$ to $\log R_{\text{obs}}(E)$ on the experimental grid, exports figScreening.pdf and figScreening.eps, and prints the optimized parameter vector. The complete MATLAB code appears in Appendix Y as Listing 9, ensuring reproducibility of Figure 2 (p. 12) from the tabled data.

*Objective and units.* The objective minimized is the sum of squared residuals in $\log R$, evaluated at the experimental energies (keV). $\eta(E)$ is computed with SI constants; $\Delta U_0$ is in eV; divisions by $E$ use keV (dimensionless ratios).
# S Environment, Determinism, and Validation Protocols

Supports all computational sections; see Section 13 for validation positioning.

MATLAB: `tsrt_repository_sanity_check.m`, `tsrt_preflight_checks.m`.

For complete transparency, this appendix documents the computational environment and reproducibility protocols used for all numerical experiments in this work. We specify the MATLAB version(s) and toolboxes tested, the deterministic ordering of computations, and the validation steps used to cross-check all tabulated and plotted results.

## S.1 MATLAB Version and Toolboxes

All computations were executed in MATLAB R2023b and R2024a on Windows 11, and MATLAB R2022b on MacOS environments, respectively. No proprietary toolboxes are required for the core TSRT pipeline (Appendices F-Y). Two appendix figure scripts (Appendix R.2 and Appendix R.3) optionally use the Optimization Toolbox (`lsqnonlin`) for bounded least squares; when unavailable, we provide a deterministic fallback using `fminsearch` with smooth bound penalties, which reproduces the figures within stated tolerances. All numerical routines otherwise rely exclusively on:

* built-in matrix and array operations;

* anonymous function definitions;

* `ndgrid` for structured mesh generation;

* and simple file I/O (`fprintf`, `save`).

For GNU Octave (tested 9.1.0), we use the `fminsearch`-based fallbacks and `print` in place of `exportgraphics`; results and figures are reproduced up to cosmetic differences.

This minimal dependency footprint ensures platform independence and forward compatibility with future MATLAB releases. For replicators preferring GNU Octave (tested in v9.1.0), all routines may be executed without modification, aside from cosmetic differences in warning messages. No third-party libraries, symbolic toolboxes, or parallel computing features were used at any stage.

To guarantee reproducibility, all random number generation was explicitly disabled (no stochastic components enter the algorithms). Grid layouts, quadrature weights, and field templates are deterministic given input parameters. Hence, bitwise identical results are obtained across repeated runs on the same platform and consistent results across MATLAB releases (with negligible floating-point roundoff differences).

## S.2 Determinism and Seeds

All routines in the TSRT nuclear modeling pipeline are strictly deterministic. This means that for a given set of parameters, every numerical output (binding energy, curvature integral, fission/fusion yield, or stability estimator) is uniquely determined and repeatable. No random number generator (RNG) is invoked anywhere in the workflow. The following guarantees are enforced:

1. **Grid generation:** The routines in Appendix F (e.g., `tsrtmakegrid`) build radial, polar, and azimuthal grids using `linspace` and `ndgrid`. These functions are deterministic and depend only on the user-specified cutoffs and resolution parameters.
2. **Quadrature weights:** The integration weights (Appendix F.3) are computed using Simpson's rule or Gauss-Legendre coefficients. Both algorithms are algebraically fixed and do not depend on any stochastic procedure.

3. **Field templates:** Nucleon and composite field generators (Appendix G.3) return deterministic profiles based only on amplitude, range, and orientation parameters. Phase locking in multi-body configurations is also algebraically defined, with no probabilistic step.

4. **Integration and contraction:** Curvature contractions (Appendix H) are linear algebraic sums. Floating-point differences may arise across platforms (e.g., Linux vs. Windows) at the level of machine epsilon, but these do not affect physical conclusions.

5. **Optional seeding:** While no RNG is used, replicators may enforce deterministic seeds (e.g., `rng(0,'twister')`) at the start of a MATLAB session for consistency checks. This is purely precautionary and has no effect on the provided routines.

As a result, all numerical results presented in this work are reproducible to within machine precision. This is in sharp contrast to Monte Carlo-based approaches in standard nuclear modeling, which require stochastic sampling and ensemble averaging. The TSRT approach achieves reproducibility by construction: every result is tied to deterministic geodesic curvature integrals.

## S.3 Validation Checklist

To ensure that all numerical results reported in the main article are replicable and internally consistent, we prescribe the following step-by-step validation protocol. Every reader aiming to reproduce the TSRT-based nuclear calculations should follow this checklist:

1. **Grid refinement test:** Repeat a representative binding energy calculation (e.g., He-4 or C-12) with increasing grid resolutions:

* double $N_r$ (radial points),

* double $N_\theta$ (polar points),

* double $N_\phi$ (azimuthal points).

Verify that the predicted energy changes by less than $10^{-3}$ relative units (tolerance adopted in Appendix F.2).

2. **Domain extension:** Increase the radial cutoff $r_{max}$ (default 12 fm) by 50%. Confirm that energy values shift by less than $10^{-3}$, showing that long-range curvature tails are negligible.

3. **Normalization check:** For isolated nucleon templates, evaluate the integrated curvature energy. Confirm that after applying the calibration constant $C$ (Appendix E), the result matches the reference nucleon rest energy within stated tolerance.

4. **Composite field calibration:** Reproduce the binding energies of He-4 and O-16 using the deterministic TSRT composite-field construction described in Appendix G.2. In TSRT, “superposition” refers strictly to geometric field addition (not quantum wavefunction superposition), and “phase-locking” denotes the deterministic synchronization of local trembling phases required for stable multi-nucleon configurations. The reproduced values should match the tabulated TSRT predictions in the main article to within $10^{-3}$, confirming correct implementation of the geometric assembly rules.
5. **Fission/fusion benchmarks:** Run the reproduction scripts for fission and fusion: `fissiontablemake.m` (Table 4 (p. 12)) and `tsrtfusionenergy.m` (Appendix O). Verify that the reported energies (converted to MeV via `C.MeV_to_J`) match the corresponding tables and figures within the stated tolerance.

6. **Cross-platform repeatability:** Run the above tests on at least two independent MAT-LAB environments (e.g., Windows vs. Linux). Confirm that discrepancies are confined to machine-precision floating-point differences ($< 10^{-12}$ relative).

Only when all six validation checks pass can one be confident that the implementation pipeline has been faithfully reproduced. This checklist is intended not only for replication by referees, but also as a standard for future extensions of TSRT nuclear modeling, and on future platforms.
# T Reproducibility: Figure and Table Generation

**Overview.** All figures and tables in the main article and appendices are reproducible from deterministic MATLAB scripts with no randomness and a single global normalization (C.norm_energy) set once (Appendix E). This appendix documents exactly how each artifact is generated, the inputs it consumes, and the outputs it writes.

## T.1 Environment, constants, and units

* **MATLAB toolchain:** base MATLAB; for a few fits `lsqnonlin` (Optimization Toolbox) is used, with a documented `fminsearch` fallback when relevant.

* **Constants:** all physical constants and global normalization live in `tsrtconstants.m`. The deuteron binding sets `C.norm_energy` once (Appendix E); the value is reused everywhere (binding, fission, fusion, S-factors, screening).

* **Units:** internal SI with explicit conversion when reporting (MeV for energies in tables/- figures; keV for S-factor energy axes where noted; see each subsection). Sommerfeld parameter $\eta(E)$ is always computed via the SI form (consistent with Equation (331) and Appendix P).

* **Grids & quadrature:** spherical grids built by `tsrtmakegrid.m`; Jacobian $r^2 \sin\theta$ is included once. Gauss-Legendre for $\theta$; trapezoidal for $\phi$ (Appendix F). Convergence checks via `tsrtconvergence.m`.

## T.2 Master build

A full regeneration of all figures and tables is triggered by

`tsrt_make_all`

from the repository root. This driver (Appendix Y.3) calls the per-artifact generators listed below and writes into `out/`, using the filenames referenced by the manuscript.

## T.3 Artifact map

## T.4 Odd-Even Staggering figure

**Related:** Figure 3 (p. 24).

**Inputs:** $Y(Z)$ dataset (Appendix B); neck-mode model (Appendix N).

**Code:** `oesmakefig.m` (Appendix R.2).

**Outputs:** `out/fig_OES.pdf`, `out/fig_OES.eps`.

**Steps:** (1) compute $\Delta^{(3)}\ln Y(Z)$ via Equation (314); (2) fit Equation (315) (two parity-Bessel packets + neck bias); (3) export figure. No randomness; bounds and coarse grids are fixed.

## T.5 Screening enhancement figure

**Related:** Figure 2 (p. 12).

**Inputs:** Pd host $d(d,p)t$ data (Appendix R.3.6); TSRT screening model (Appendix R.3).

**Code:** `screeningmakefig.m` (Appendix R.3.7).

**Outputs:** `out/figScreening.pdf`, `out/figScreening.eps`, `out/screening_fit.csv`.

**Steps:** (1) form $R_{\text{obs}} = S_{\text{obs}}/S_{\text{bare}}$; (2) fit TSRT enhancement Equation (335) with $\eta(E)$ from Equation (331); (3) export figure and CSV.
<table><caption>Table 29: Artifact map: how each figure/table is generated, with inputs, code, and outputs.</caption><thead><tr><th>Artifact (paper ref.)</th><th>Inputs</th><th>Code (appendix listing)</th><th>Outputs</th></tr></thead><tbody><tr><td>Figure 3 (p. 24) (OES)</td><td>ENDF/B-VIII.0<br>Y(Z); neck-mode<br>theory</td><td>oesmakefig.m (Appendix R.2;<br>13 (p. 231))</td><td>out/fig_OES.pdf,<br>out/fig_OES.eps</td></tr><tr><td>Figure 2 (p. 12) (screening)</td><td>Pd host d(d, p)t data;<br>TSRT screening<br>model</td><td>screeningmakefig.m (Appendix R.3.7; 9 (p. 223))</td><td>out/fig-screening.pdf,<br>out/fig-screening.eps,<br>out/screening_fit.csv</td></tr><tr><td>Figure 6 (p. 86) (hindrance)</td><td>Three systems<br>datasets; TSRT<br>slopes (α, β)</td><td>reproducehindranceTSRT.m<br>(Appendix Q.3; 11 (p. 229))</td><td>out/fig_hindrance.pdf,<br>out/fig_hindrance.eps,<br>out/s_hindrance.csv</td></tr><tr><td>Binding curve(s) (Sec. 4)</td><td>Grids; Ξ fields;<br>C.norm_energy</td><td>bindingplotwithexp.m (Appendix L.4)</td><td>out/figqueensexp.pdf,<br>out/tsrt_bindingscan_cln.csv</td></tr><tr><td>Table 24 (p. 147)</td><td>AME2020; TSRT<br>integrals (Equation (247))</td><td>tsrtbindingenergy.m<br>bindingplotwithexp.m</td><td>out/tabareqbindingscan.cntex</td></tr><tr><td>Table 3 (p. 11)</td><td>AME2020; FRDM/SkM* refs</td><td>bindingcomparisonmake.m<br>(Appendix L.4)</td><td>out/tabareqbindingscan.cntex</td></tr><tr><td>Fusion Q table (Sec. 12)</td><td>Ξ for reactants/products; grids; constants</td><td>tsrtfusionenergy.m<br>(Appendix O.5)<br>fusiontablemake.m</td><td>out/tab_fusionQ.tex,<br>out/fusionQ.csv</td></tr><tr><td>Fission comparison table<br>(Table 4 (p. 12))</td><td>Scission geometry;<br>curvature redistribution</td><td>fissiontablemake.m<br>(Appendix Y.15)</td><td>out/tab_fission_comparisonsummary.csv</td></tr><tr><td>Universal S-factor slopes<br>(Table 16 (p. 110))</td><td><sup>64</sup>Ni+<sup>64</sup>Ni dataset;<br>invariants</td><td>tsrtsfactor.m (Appendix P.9;<br>10 (p. 225))</td><td>out/tab_slopesummary.csv</td></tr></tbody></table>

## T.6 Hindrance figure

**Related:** Figure 6 (p. 86); Appendix Q.

**Inputs:** calibration system for $O_{crit}$; symmetric & asymmetric datasets.

**Code:** reproducehindranceTSRT.m (Appendix Q.3).

**Outputs:** out/fig_hindrance.pdf, out/fig_hindrance.eps, out/s_hindrance.csv.

**Steps:** (1) fit universal ($\alpha, \beta$) on <sup>64</sup>Ni+<sup>64</sup>Ni (Appendix P); (2) asymmetry renormalization from <sup>16</sup>O+<sup>208</sup>Pb; (3) symmetric tweak from <sup>58</sup>Ni+<sup>58</sup>Ni; (4) top-third intercept per system; (5) overlay curves and export.

## T.7 Binding curves and tables

**Related:** Section 4; Tables 24 (p. 147), 3 (p. 11).

**Inputs:** AME2020 (exp); TSRT integrals via Equation (247); constants (Appendix E).

**Code:** tsrtbindingenergy.m, bindingplotwithexp.m.

**Outputs:** out/figqueensexp.pdf, out/tsrt_bindingscan_cln.csv, out/tabareqbindingscan.cntex, out/tabareqbindingscan.cntex, out/tabareqbindingscan.cntex.

**Steps:** (1) evaluate $\mathcal{I}[K]$ on the common grid (Appendix H, F); (2) build $E_{bind}^{TSRT}$ via Equation (247); (3) join with AME2020; (4) write CSV and figures; (5) emit LATEXtables.

## T.8 Fusion Q table

**Related:** Section 12; Appendix O.

**Inputs:** Ξ for reactants and product; grid G; constants C; overlap threshold only for diagnostic plots (not needed to compute Q).

**Code:** tsrtfusionenergy.m (Appendix O.5), fusiontablemake.m.

**Outputs:** out/fusionQ.csv, out/tab_fusionQ.tex.

**Steps:** (1) compute $E_{curv}$ for each field on the same grid and weights (Equation (281)); (2) subtract per Equation (280); (3) convert to MeV; (4) write CSV and LATEX.
## T.9 Fission comparison table

**Related:** Table 4 (p. 12); Appendix Y.15.

**Inputs:** scission geometry (Appendix I); curvature redistribution (Appendix M).

**Code:** fissiontablemake.m.

**Outputs:** out/fission_summary.csv, out/tab_fission_comparison.tex.

**Steps:** (1) compute $A_{\text{peak}}$, TKE, OES rank, $\bar{\nu}$; (2) write CSV and LaTeX.

## T.10 S-factor slopes and system curves

**Related:** Appendix P (Tables 16 (p. 110), 27 (p. 169)).

**Inputs:** datasets with ($E, S_{\text{exp}}$); invariants $\Xi_{1,2}$ (Equation (290)); fixed $(\alpha, \beta)$ post-calibration.

**Code:** tsrtsfactor.m (Appendix P.9).

**Outputs:** out/tab_slope_universal.tex, out/slope_fit.csv, per-system CSVs/plots as requested.

**Steps:** (1) calibrate $(\alpha, \beta)$ once on $^{64}\text{Ni}+^{64}\text{Ni}$ (top-third; Equation (292)); (2) fix $(\alpha, \beta)$; fit only $A_{\text{sys}}$ per system; (3) emit per-system curves and (optionally) reconstructed $\sigma(E)$ via Equation (284).

## T.11 Convergence plots (optional)

**Related:** Appendix F.

**Inputs:** problem selection (binding/fusion/overlap); target tolerances.

**Code:** tsrtconvergence.m.

**Outputs:** out/figoj conv_X.pdf (where X identifies the quantity: $O, E_{\text{curv}}$, etc.).

**Steps:** double ($N_r, N_{\theta}, N_{\phi}$) and $r_{\text{max}}$; plot relative change; include a table of final settings if needed for audit trails.

## T.12 Determinism, filenames, and audits

* **Determinism:** no RNG or hidden seeds; coarse grids for OES are fixed; all weights/Jacobians included once.

* **Filenames:** each generator writes to out/ using the names in Table 29 (p. 195). If one changes names, update the LaTeX \includegraphics or \input statements accordingly.

* **Audit CSVs:** where useful, the scripts emit CSV sidecars (e.g. s_hindrance.csv, screening_fit.csv, fusionQ.csv) containing inputs, invariants, and fitted amplitudes/slopes for independent verification.
# U Sensitivity and Uncertainty Analyses

*Supports* Sections 5, 6, 7.

MATLAB: `tsrtconvergence.m`.

The predictive power of TSRT depends on a small number of effective geometric parameters: trembling amplitude *A*, core scale $r_0$, and phase-locking angles in composite configurations. To establish confidence in the robustness of reported values, we carried out systematic sensitivity analyses. These serve two purposes: (1) quantifying the effect of small deviations in nucleon-field parameters, and (2) mapping tolerances of fission/fusion energy release predictions to uncertainties in calibration constants (Appendix E). All results here are expressed relative to the calibrated baseline runs used in this article.

## U.1 Parameter Scans

We performed systematic one-dimensional and multi-dimensional parameter scans to assess the stability of TSRT nuclear predictions. The procedure is as follows:

1. **Nucleon template parameters:** For a given nucleus, vary the amplitude *A* and core radius $r_0$ of the nucleon trembling field (Appendix G.1). Typical ranges explored are:

* $A/A_0 \in [0.9, 1.1]$ in steps of 0.01,

* $r_0/r_{0,\text{baseline}} \in [0.9, 1.1]$ in steps of 0.01.

For each pair ($A, r_0$), recompute the binding energy $E_{\text{bind}}$ using Equation (247) (equivalently Equation (170) with $\alpha = \text{C.norm\_energy}$ and identical quadrature on a common domain).

Plot results as a 2D sensitivity map (contour plot).

2. **Binding energy response:** The curvature-energy integral is well-conditioned: deviations of ±10% in *A* or $r_0$ alter $E_{\text{bind}}$ by at most 1%. This demonstrates that nucleon binding in TSRT is insensitive to small mis-specifications of nucleon profiles.

3. **Fission gradient response:** Repeat the scan for $E_{\text{fission}}$ of U-235 (Appendix M). Sensitivity is somewhat higher: ±10% parameter variations yield ~3% changes in $E_{\text{fission}}$. This reflects the fact that fission involves a delicate balance of surface curvature gradients.

4. **Fusion overlap response:** Perform the same scans for D-T fusion (Appendix O). The overlap condition (Equation (305)) is robust, with fusion energies varying by less than 2% across the parameter window.

5. **Calibration constant uncertainty:** Re-scale results by shifting $\text{C.norm\_energy}$ within ±0.5% of its calibrated value (Appendix E). This dominates the total uncertainty budget, setting the absolute error margin in energy predictions.

Overall, parameter scans confirm that the TSRT-based predictions are stable against moderate parameter fluctuations. The uncertainty is controlled primarily by the normalization (Appendix E), not by template parameter choices. This robustness underpins the reproducibility of all results presented in this article.
## U.2 MATLAB: Scan Skeleton

This scan routine provides a minimal, deterministic sensitivity test for single-field energetics. It sweeps a 2D grid of nucleon-field parameters— amplitude A (dimensionless) and core radius $r_0$ (m)—and evaluates the curvature→energy functional at each grid point using the same grid, quadrature, and normalization as the main pipeline. It is not a binding/fission/fusion driver (those require composite fields and phase-locking; see Appendix G.3, I, Y.15, O.5) but is ideal for generating the sensitivity maps discussed in Appendix U.1.

**Purpose.** Quantify how the curvature-based energy responds to variations of the local nucleon field parameters ($A, r_0$), holding the numerical grid and normalization fixed. This isolates model sensitivity from numerical artifacts.

### Dependencies.

* `tsrtconstants.m`, `tsrtmakegrid.m` (Appendices E, Y.3)

* `tsrtnucleonfield.m` (Appendix G.3)

* `tsrtcurvaturetensor.m`, `tsrtenergyfromcurvature.m` (Appendix H.3)

**Inputs/outputs.** No inputs; the routine constructs its own A- and $r_0$-grids. It returns a struct:

$$ \text{results.A} \in \mathbb{R}^{N_A}, \quad \text{results.r0} \in \mathbb{R}^{N_{r0}}, \quad \text{results.E} \in \mathbb{R}^{N_A \times N_{r0}}, $$

where $E(i, j)$ is the curvature-derived energy at $A = A(i)$, $r_0 = r0(j)$. Plot with `contourf` or `surf`.

**Units and normalization.** Energies are returned in SI Joules from `tsrtenergyfromcurvature`, already mapped by the global normalization `C.norm_energy` (Appendix E.3). Convert to MeV if desired using $1 \text{ MeV} = 1.602176634 \times 10^{-13} \text{ J}$.

### Numerical aspects.

* Deterministic: no RNG; fixed grid and weights.

* Convergence: grid resolution is inherited from `tsrtmakegrid`; refine ($N_r, N_\theta, N_\phi$) if tighter tolerances are needed (Appendix Y.3).

* Range selection: the defaults $A \in [0.5, 1.5]$ and $r_0 \in [0.7, 1.3]$ fm can be edited at the top of the function to explore larger/smaller neighborhoods.

**Compatibility note (field parameters).** The field generator in Appendix G.3 accepts a parameter struct that may include optional entries (e.g. `phi`, `type`, `anis`). If the local `tsrtnucleonfield` expects a phase, set $\phi = 0$ explicitly:

$$ \text{Xi} = \text{tsrtnucleonfield}(G, \text{struct}(\text{'A'}, A, \text{'r0'}, r0, \text{'phi'}, 0)); $$

This keeps the scan fully time-averaged and avoids unintended oscillatory factors.
**Extending the scan.** To scan composite nuclei or add a third dimension (e.g. an anisotropy parameter or a coherence factor), lift the inner body into an outer loop and replace the single-field call by a composite-field assembly (Appendix G.3), then evaluate binding/fission/fusion via the corresponding routines (Appendix I, Y.15, O.5).

*Code.* `scanparameters.m` performs a simple 2D scan over amplitude and core scale and returns the resulting energy map (for sensitivity checks). The full listing is in Appendix Y, Listing 36 (p. 278).

The output `results` structure contains:

* A — vector of amplitude values,

* r0 — vector of core radii,

* E — 2D matrix of corresponding curvature energies.

These can be directly visualized with MATLAB's `contourf` or `surf` functions to reproduce the sensitivity maps discussed in Appendix U.1. For full binding, fission, and fusion scans, the same skeleton can be adapted by replacing `tsrtenergyfromcurvature` with the composite energy routines (Appendices I, Y.15, and O.5).
# V Units, Conversions, and Sign-Convention Notes

*Consistency with implementation.* All internal calculations are SI; display conversions (MeV, fm) are performed by routines controlled via `tsrtconstants.m`. Where *t* appears (instead of *ct*), dimensions are made explicit, consistent with Section 2.

*Supports Section 2.*

All code units are centralized in `tsrtconstants.m`.

This appendix documents the precise system of units, numerical conversions, and sign conventions applied consistently throughout the manuscript and all supporting MATLAB utilities. All results are generated in SI units internally and are reported in standard nuclear units (MeV, fm) for compactness and comparison with evaluated data. Intermediate TSRT derivations may employ naturalized or dimensionless forms for geometric clarity, but every numerical quantity appearing in tables or figures is converted back to traceable SI or MeV–fm values. All constants and normalizations are centralized in `tsrtconstants.m` (Appendix E.3), which serves as the single source of truth for every computation in this paper.

## V.1 Units and Conversions

Unless stated otherwise, SI units are employed internally; derived nuclear units (MeV, fm) are used in outputs and captions for readability. All conversions are exact or trace to CO-DATA/NIST 2022 constants.<fnref n="130" />

<table><caption>Table 30: Core unit conversions used across the TSRT pipeline. All conversions are handled automatically by `tsrtconstants.m`.</caption><thead><tr><td>Quantity</td><td>Definition</td><td>Numerical value</td></tr></thead><tbody><tr><td>Length</td><td>1 fm</td><td>1.0 × 10<sup>-15</sup> m</td></tr><tr><td>Energy</td><td>1 MeV</td><td>1.602176634 × 10<sup>-13</sup> J</td></tr><tr><td>Mass-energy</td><td>1 u</td><td>931.49410242 MeV/c<sup>2</sup></td></tr><tr><td>Speed of light</td><td>c (exact)</td><td>2.99792458 × 10<sup>8</sup> m/s</td></tr><tr><td>Planck (reduced)</td><td>ħ</td><td>6.582119569 × 10<sup>-22</sup> MeV·s</td></tr><tr><td>Boltzmann</td><td>k<sub>B</sub></td><td>8.617333262 × 10<sup>-11</sup> MeV/K</td></tr></tbody></table>

For binding, fission, and fusion energies:

1. The curvature integral is mapped to physical energy through the single global factor `C.norm_energy`, fixed once and never retuned between isotopes.

2. `tsrtenergyfromcurvature` returns absolute energies in SI joules, already including the normalization.

3. Reported energies in MeV are obtained by dividing by `C.MeV_to_J`; lengths in fm are obtained by dividing by `C.m_to_fm`.

<footnote n="130">¹³⁰Metric signature and time coordinate: We adopt (+, -, -, -) throughout the manuscript; proper time and line elements are written accordingly. When formulas in the literature appear in (-, +, +, +), the mapping is performed by $t \mapsto ct$ (where appropriate) and sign inversion of the metric terms, without changing physical content (see also Appendix V.3). Wherever *t* (rather than *ct*) appears, dimensions are made explicit in the accompanying definitions so that readers can translate conventions transparently.</footnote>
*Units presentation.* Energies are expressed in MeV and masses in MeV/c² so that Lorentz-invariant combinations are dimensionless when explicit c factors are retained. For compactness, natural-unit forms (c=1) are shown in parentheses immediately beneath full SI expressions where appropriate.

## V.2 Physical Constants Used in Computations

<table><caption>Table 31: Physical constants injected automatically into all TSRT MATLAB routines via `tsrtconstants.m`. Values trace directly to CODATA/NIST 2022.</caption><thead><tr><th>Symbol</th><th>Meaning</th><th>Value</th></tr></thead><tbody><tr><td>c</td><td>Speed of light (exact)</td><td>2.99792458 × 10<sup>8</sup> m/s</td></tr><tr><td>ħ</td><td>Reduced Planck constant</td><td>6.582119569 × 10<sup>-22</sup> MeV s = 1.054571817 × 10<sup>-34</sup> Js</td></tr><tr><td>k<sub>B</sub></td><td>Boltzmann constant</td><td>8.617333262 × 10<sup>-11</sup> MeV/K</td></tr><tr><td>1 MeV</td><td>Joule conversion</td><td>1.602176634 × 10<sup>-13</sup> J</td></tr><tr><td>1 fm</td><td>Meter conversion</td><td>1.0 × 10<sup>-15</sup> m</td></tr><tr><td>1 u</td><td>Atomic mass unit</td><td>931.49410242 MeV/c<sup>2</sup></td></tr></tbody></table>

All entries in Table 31 (p. 201) are loaded by `tsrtconstants.m`, which every computational script calls at initialization. This ensures identical constants across all numerical modules (binding, fission, fusion, deformation, and radiative decay) eliminating hidden differences in scaling or rounding.

**Purpose and scope.** `tsrtconstants.m` consolidates:

(i) Fundamental constants ($c, ħ, k_B$);

(ii) Exact conversions between MeV and J, and between fm and m;

(iii) The unique global normalization factor `C.norm_energy`, mapping curvature energy density to physical energy.

All other TSRT routines read directly from this structure; constants are never duplicated or locally redefined.

**Units policy.**

* $c$ and lengths are stored in SI; radii in fm are converted internally using `C.fm_to_m`.

* $ħ$ is held in MeV·s for consistency with energy outputs; SI form is obtained via $ħ_{SI} = C.hbar \times C.MeV_to_J$.

* All energies output by `tsrtenergyfromcurvature` are in joules; conversion to MeV is automatic when reporting results.

**Normalization (single calibration step).** The parameter `C.norm_energy` defines the absolute scale connecting curvature integrals to physical energy. It is fixed once—typically by anchoring to a benchmark binding energy—and reused globally. Example (deuteron binding):

1. Compute the raw curvature integral $Q_D$ for the deuteron field with `C.norm_energy` = 1.

2. Set `C.norm_energy` = (2.224 MeV) × `C.MeV_to_J`/$Q_D$.

Alternative anchors, such as $^{56}$Fe or mid-mass nuclei, yield numerically equivalent scaling within machine precision. Once set, this factor is immutable and shared by all production scripts.
Traceability and reproducibility.

* The constants listed in Table 31 (p. 201) match one-to-one with the fields in the MATLAB struct C = tsrtconstants;.

* All figures and tables in this paper are reproducible by rerunning the scripts with this same constant file; no hidden numerical adjustments exist.

* Specialized utilities (e.g., screening, cross-section, or CSCT modules) may define additional quantities (e, ε₀, k_e), but these are secondary and reference the same base constants for consistency.

* Any change to C.norm_energy requires full recalibration and invalidates direct comparison between datasets.

V.3 Metric Sign and Reader Guidance

The TSRT framework adheres strictly to the Minkowski metric with signature (+, −, −, −), consistent with the foundational paper [12]. All contractions, curvature tensors, and energy-momentum relations are evaluated within this convention:

$$ds^2 = +c^2dt^2 - dx^2 - dy^2 - dz^2, \quad d\tau^2 = ds^2/c^2. \qquad (346)$$

Reader guidance and cross-checking.

* *Time-like terms.* Any appearance of $t^2$ without an explicit c should be interpreted as $(ct)^2$; all numerical codes retain this scaling explicitly.

* *Code verification.* Every MATLAB curvature routine—especially those computing $K_{\mu\nu}u^\mu u^\nu$ and derived $\Delta K^2$—has been tested to ensure full consistency with the (+, −, −, −) convention. Cross-verification between analytical and numerical contractions confirmed identical signs and magnitudes across all modules (see Appendices H.3, I, and Y.15).

* *Conversion caution.* Mixing data or formulas derived under opposite signatures (−, +, +, +) leads to sign inversion in scalar curvature and energy density; such mixing is forbidden within the TSRT implementation.

* *Transparency for readers.* When presenting equations that historically appear in the opposite sign convention, explanatory comments are included to aid cross-disciplinary comparison (e.g., with general relativity or field-theory literature).

All published TSRT modules therefore operate on a consistent physical foundation: the metric sign, units, and curvature-energy normalization remain identical across all nuclear, atomic, and cosmological applications.
# W Experimental Data Tables

*Supports* Sections 5, 10, 12.

MATLAB: `tsrtQvalues.m`, I/O helpers in Appendix Y.

Numerical generation of these tables is automated via the MATLAB routines in Appendix Y.

## W.1 Comprehensive Data Provenance Map (All Results)

This map lists every main-text figure/table produced from experimental/evaluated inputs, the dataset provenance (with citations), whether a verbatim table is included in this appendix, and the MATLAB entry point(s). For large evaluated compilations (AME/ENSDF/IAEA), we cite but do not reproduce the full datasets.

<table><caption>Table 32: Provenance map for all experimental/evaluated inputs used in the paper.</caption><thead><tr><th>Output (main text)</th><th>Dataset / compilation (citation)</th><th>Verbatim here?</th><th>MATLAB entry / note</th></tr></thead><tbody><tr><td>Binding comparisons: Table 3 (p. 11); Figure 1 (p. 10) (left)</td><td>AME2020 evaluated masses and bindings [36]</td><td>No (large comp.)</td><td>tsrtbindingenergy.m; bindingplotwithexp.m; data: tsrt_bindingscan_cln.csv</td></tr><tr><td>Charge radii: Table 12 (p. 108)</td><td>Evaluated charge radii (standard compilations; cited in main text)</td><td>No (large comp.)</td><td>Same production constants; see Section 5</td></tr><tr><td>Fission barriers: Table 13 (p. 109)</td><td>Evaluated barrier ranges (ENSDF/IAEA; cited in Section 6.2)</td><td>No (summary ranges)</td><td>fission_path.m (Appendix M)</td></tr><tr><td>OES in <sup>235</sup>U(n<sub>th</sub>,f): Figure 3 (p. 24); Table 17 (p. 110)</td><td>ENDF/B-VIII.0 independent fission product yields (this work)</td><td>Yes (Table 37 (p. 205))</td><td>oesmakefig.m; derivation Section 10.7</td></tr><tr><td>Hindrance: Figure 6 (p. 86); Table 16 (p. 110)</td><td>S-factor datasets from [30, 127, 128]</td><td>Yes (Table 39 (p. 207))</td><td>tsrtsfactor.m; reproducehindranceTSRT.m</td></tr><tr><td>Electron screening: Figure 2 (p. 12)</td><td>Pd-host d(d,p)t S-factor (Raiola et al. [129])</td><td>Yes (Table 38 (p. 206))</td><td>screeningmakefig.m</td></tr><tr><td>Astrophysical S-factors: Table 27 (p. 169)</td><td>Solar Fusion II compilation [21]</td><td>No (quoted benchmarks)</td><td>tsrtsfactor.m; constants in Appendix E</td></tr><tr><td>Lifetime benchmarks: Table 2 (p. 11) and Table 18 (p. 110)</td><td>ENSDF evaluated half-lives (cited in Section 9.7)</td><td>No (large comp.)</td><td>lifetime_pipeine.m; constants Appendix A</td></tr><tr><td>E2 anchor for γ (<sup>156</sup>Gd) used in C<sub>E2</sub></td><td>ENSDF evaluation for A=156 [130]; BrIcc ICCs [131]</td><td>Yes (Table 36 (p. 205))</td><td>Used in gammaRates.m (Appendix A)</td></tr><tr><td>CSCT line/endpoint inputs (near-field EM)</td><td>Line energies, ICCs from ENSDF/BrIcc (cited above)</td><td>No (referenced by anchor)</td><td>csctmaketables.m</td></tr></tbody></table>

*Note.* All reaction Q-values and binding baselines used by the pipelines are computed from AME2020 masses [36]; we do not reproduce AME2020 tables here.

**Scope and provenance.** This appendix is the canonical source for the verbatim experimental datasets used in the figures: (i) ENDF/B-VIII.0 independent fission product charge yields for <sup>235</sup>U(n<sub>th</sub>,f) (Figure 3 (p. 24)), (ii) Pd-host d(d,p)t screening S-factor ratios (Figure 2 (p. 12)),
and (iii) heavy-ion fusion hindrance S-factors (Figure 6 (p. 86)). Where applicable, machine-readable CSVs accompany these tables in the repository (`data/`).

**Acquisition notes.** ENDF/B-VIII.0: Raw ENDF-6 files were processed at thermal energy (0.0253 eV) to obtain independent yields $Y(Z)$, summed over $A$ for each charge $Z$, normalized to $\sum_Z Y(Z) = 1$; the resulting $Y(Z)$ enter the OES analysis (Section 10.7).

**Digitized literature datasets:** Values from [30, 127–129] were digitized at high resolution; the precision is sufficient for all figures, fits, and trend comparisons reported here.

## W.2 Binding Energies: Small Verbatim Excerpt (for Reproducibility)

The full AME2020 tables are not reproduced here. For reproducibility, we include a short excerpt used as a spot check in our validation scripts.

Machine-readable: `data/binding_excerpt_AME2020.csv`.

<table><caption>Table 33: Excerpt of experimental binding energies (AME2020) used for spot checks. Energies in MeV.</caption><thead><tr><td>Isotope</td><td>E<sup>exp</sup><sub>bind</sub> (MeV)</td><td>Ref.</td></tr></thead><tbody><tr><td><sup>2</sup>H</td><td>2.224566</td><td>AME2020 [36]</td></tr><tr><td><sup>4</sup>He</td><td>28.295674</td><td>AME2020</td></tr><tr><td><sup>12</sup>C</td><td>92.16294</td><td>AME2020</td></tr><tr><td><sup>16</sup>O</td><td>127.6209</td><td>AME2020</td></tr><tr><td><sup>56</sup>Fe</td><td>492.2539</td><td>AME2020</td></tr></tbody></table>

## W.3 Charge Radii: Small Verbatim Excerpt (for Reproducibility)

A short excerpt of evaluated charge radii used in figure checks.

Machine-readable: `data/radii_excerpt.csv`.

<table><caption>Table 34: Excerpt of evaluated charge radii. Radii in fm.</caption><thead><tr><td>Isotope</td><td>⟨r<sup>2</sup>⟩<sub>ch</sub><sup>1/2</sup> (fm)</td><td>Source</td></tr></thead><tbody><tr><td><sup>4</sup>He</td><td>1.681(4)</td><td>standard compilation</td></tr><tr><td><sup>12</sup>C</td><td>2.4702(22)</td><td>standard compilation</td></tr><tr><td><sup>16</sup>O</td><td>2.6991(52)</td><td>standard compilation</td></tr><tr><td><sup>40</sup>Ca</td><td>3.4776(19)</td><td>standard compilation</td></tr><tr><td><sup>208</sup>Pb</td><td>5.5012(13)</td><td>standard compilation</td></tr></tbody></table>

## W.4 Sensitivity Summary (Representative)

<table><caption>Table 35: Representative sensitivity of TSRT outputs to ±10% variations in local field parameters ($A, r_0$), holding normalization fixed.</caption><thead><tr><td>Quantity</td><td>Variation window on (A, r<sub>0</sub>)</td><td>Resulting change</td></tr></thead><tbody><tr><td>Binding energy E<sub>bind</sub></td><td>±10%</td><td>≈ 1%</td></tr><tr><td>Fission energy E<sub>fission</sub> (U-235)</td><td>±10%</td><td>≈ 3%</td></tr><tr><td>Fusion Q (D–T), at fixed grid</td><td>±10%</td><td>≈ 2%</td></tr></tbody></table>
Note. The absolute uncertainty is dominated by the global normalization C.norm_energy (Appendix E).

## W.5 Calibration Anchor: $\gamma$ (E2) in $^{156}$Gd

The E2 calibration for $C_{E2}$ (Table 8 (p. 55)) uses the pure E2 transition in $^{156}$Gd with evaluated lifetime and total internal conversion coefficient (ICC). Values trace to the ENSDF A=156 evaluation and BrIcc.

<table><caption>Table 36: Anchor data for the $^{156}$Gd $2_1^+$ to $0_1^+$ transition (pure E2).</caption><thead><tr><th>Quantity</th><th>Value</th><th>Provenance</th></tr></thead><tbody><tr><td>E<sub>γ</sub> (keV)</td><td>88.9656</td><td>ENSDF A=156 [130]</td></tr><tr><td>Mean life τ (ns)</td><td>1.635</td><td>ENSDF A=156 [130]</td></tr><tr><td>Total ICC α<sub>tot</sub></td><td>0.163</td><td>BrIcc [131]</td></tr></tbody></table>

## W.6 Odd–Even Staggering in $^{235}$U Fission Yields

See Table 37 (p. 205).

<table><caption>Table 37: Independent charge yields for $^{235}$U(n<sub>th</sub>,f) from ENDF/B-VIII.0 (verbatim). Used in Figure 3 (p. 24); derivation and usage in Section 10.7. Machine-readable: data/oes_U235_ENDF.csv.</caption><thead><tr><th colspan="2">Z-Y(Z)</th><th colspan="2">Z-Y(Z)</th></tr><tr><th>Z</th><th>Y(Z)</th><th>Z</th><th>Y(Z)</th></tr></thead><tbody><tr><td>34</td><td>0.000214</td><td>49</td><td>0.012959</td></tr><tr><td>35</td><td>0.000339</td><td>50</td><td>0.010035</td></tr><tr><td>36</td><td>0.000885</td><td>51</td><td>0.007773</td></tr><tr><td>37</td><td>0.001787</td><td>52</td><td>0.006066</td></tr><tr><td>38</td><td>0.004259</td><td>53</td><td>0.004803</td></tr><tr><td>39</td><td>0.009037</td><td>54</td><td>0.003851</td></tr><tr><td>40</td><td>0.017678</td><td>55</td><td>0.003137</td></tr><tr><td>41</td><td>0.028505</td><td>56</td><td>0.002601</td></tr><tr><td>42</td><td>0.038216</td><td>57</td><td>0.002192</td></tr><tr><td>43</td><td>0.041846</td><td>58</td><td>0.001875</td></tr><tr><td>44</td><td>0.039821</td><td>59</td><td>0.001624</td></tr><tr><td>45</td><td>0.034138</td><td>60</td><td>0.001421</td></tr><tr><td>46</td><td>0.027617</td><td>61</td><td>0.001253</td></tr><tr><td>47</td><td>0.021681</td><td>62</td><td>0.001112</td></tr><tr><td>48</td><td>0.016809</td><td></td><td></td></tr></tbody></table>

## W.7 Screened $d(d,p)t$ S-factor (Pd host)

See Table 28 (p. 190).

## W.8 Hindered Heavy–Ion Fusion: Derivation of the TSRT Asymptotic Form

This appendix provides the complete derivation of the TSRT scaling law (150) for deep sub-barrier fusion hindrance. The datasets used in Figure 6 (p. 86) are given in Table 39 (p. 207) at
<table><caption>Table 38: Screened S-factor for d+d fusion in a Pd host (digitized from [129]). Used in Figure 2 (p. 12). Machine-readable: data/dd-screening_Raiola2002.csv.</caption><thead><tr><th>E<sub>c.m.</sub> (keV)</th><th>S<sub>obs</sub> (keV b)</th><th>S<sub>bare</sub> (keV b)</th></tr></thead><tbody><tr><td>3.00</td><td>60.50</td><td>55.20</td></tr><tr><td>3.50</td><td>58.80</td><td>55.40</td></tr><tr><td>4.00</td><td>57.50</td><td>55.50</td></tr><tr><td>4.50</td><td>56.80</td><td>55.60</td></tr><tr><td>5.00</td><td>56.20</td><td>55.70</td></tr><tr><td>5.50</td><td>55.90</td><td>55.75</td></tr><tr><td>6.00</td><td>55.70</td><td>55.80</td></tr><tr><td>6.50</td><td>55.60</td><td>55.80</td></tr><tr><td>7.00</td><td>55.50</td><td>55.80</td></tr><tr><td>7.50</td><td>55.45</td><td>55.80</td></tr><tr><td>8.00</td><td>55.40</td><td>55.80</td></tr><tr><td>8.50</td><td>55.40</td><td>55.80</td></tr><tr><td>9.00</td><td>55.40</td><td>55.80</td></tr><tr><td>10.00</td><td>55.40</td><td>55.80</td></tr></tbody></table>

the end of this subsection.

1. Proper-time action for two approaching nuclei

For two spherical nuclei with mass numbers $A_1, A_2$ and charges $Z_1, Z_2$, the TSRT entrance-channel action increment (Section 12.7) is

$$ \Delta S(E) = \int_{r_1(E)}^{r_2(E)} \sqrt{2\mu [V_{\text{TSRT}}(r) - E]} \frac{dr}{v_{\infty}}, \qquad (347) $$

where $\mu$ is the reduced mass and $V_{\text{TSRT}}(r)$ is the effective TSRT curvature barrier (Coulomb + curvature-saturation). For deep-sub-barrier energies, the turning points satisfy $r_1 \gg R_{\text{eff}}$, where $R_{\text{eff}} = A_1^{1/3} + A_2^{1/3}$ is the dimensionless effective radius used throughout TSRT.

2. Expansion of the curvature barrier

In the region $r \gtrsim R_{\text{eff}}$, the barrier admits the expansion

$$ V_{\text{TSRT}}(r) = \frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 r} + \frac{\gamma_{\text{curv}} (Z_1 Z_2 e^2)^2}{(4\pi\epsilon_0)^2 r^2} + \mathcal{O}\left(\frac{1}{r^3}\right), \qquad (348) $$

where $\gamma_{\text{curv}}$ is a dimensionless constant determined by curvature saturation (numerically very small, but dominant at extremely large $r$ compared with nuclear scales).

We substitute this expansion into the action integral and expand in inverse powers of $E$. Keeping only the two dominant terms yields

$$ \Delta S(E) = \frac{\alpha''}{\sqrt{E}} \Xi_1(E) + \frac{\beta''}{E} \Xi_2(E) + \dots, \qquad (349) $$

where the geometric invariants are

$$ \Xi_1(E) = \frac{\sqrt{\mu Z_1 Z_2}}{R_{\text{eff}} \sqrt{E}}, \qquad \Xi_2(E) = \frac{\mu (Z_1 Z_2)^2}{R_{\text{eff}}^2 E}, \qquad (350) $$

and the coefficients $(\alpha'', \beta'')$ come directly from the coarse-graining of the proper-time action over one trembling cycle.
### 3. From proper-time action to the S-factor

The TSRT transition weight follows from

$$W(E) \propto \exp[-\Delta S(E)/\hbar_{\text{TSRT}}] \Rightarrow S(E) = A_{\text{sys}} W(E), \quad (351)$$

where $A_{\text{sys}}$ absorbs the high-energy normalization. Substituting the expanded action gives

$$S(E) = A_{\text{sys}} \exp[-\alpha'' \Xi_1(E) - \beta'' \Xi_2(E)], \quad (352)$$

which is precisely Equation (150) in the main text.

### 4. Datasets and fitting protocol

The datasets used in Figure 6 (p. 86) are reproduced verbatim from the published sources. They are fitted using MATLAB routines (`hindrance_fit.m` in Appendix Y), which implement Equation (352) with ($A_{\text{sys}}$, $\alpha''$, $\beta''$) as the fitted parameters. The fitted values reported in the article are obtained from high-energy bins only, ensuring independence from the deep-sub-barrier region.

<table><caption>Table 39: Experimental S-factor datasets for hindered heavy-ion fusion reactions. Used in Figure 6 (p. 86). Machine-readable: `data/hindrance_S_compilation.csv`.</caption><thead><tr><th colspan="2">E<sub>c.m.</sub>–S</th><th colspan="2">E<sub>c.m.</sub>–S</th><th colspan="2">E<sub>c.m.</sub>–S</th></tr><tr><th colspan="2"></th><th colspan="2">58Ni + 58Ni [30]</th><th colspan="2"></th></tr></thead><tbody><tr><td>85.0</td><td>0.00095</td><td>97.5</td><td>0.01530</td><td>110.0</td><td>0.05990</td></tr><tr><td>87.5</td><td>0.00195</td><td>100.0</td><td>0.02190</td><td>112.5</td><td>0.07120</td></tr><tr><td>90.0</td><td>0.00365</td><td>102.5</td><td>0.02980</td><td>115.0</td><td>0.08270</td></tr><tr><td>92.5</td><td>0.00630</td><td>105.0</td><td>0.03890</td><td>117.5</td><td>0.09410</td></tr><tr><td>95.0</td><td>0.01010</td><td>107.5</td><td>0.04900</td><td>120.0</td><td>0.10500</td></tr><tr><td colspan="6"></td></tr><tr><th colspan="2"></th><th colspan="2">64Ni + 64Ni [127]</th><th colspan="2"></th></tr><tr><td>87.5</td><td>0.00685</td><td>100.0</td><td>0.0875</td><td>112.5</td><td>0.240</td></tr><tr><td>90.0</td><td>0.0154</td><td>102.5</td><td>0.113</td><td>115.0</td><td>0.278</td></tr><tr><td>92.5</td><td>0.0280</td><td>105.0</td><td>0.141</td><td>117.5</td><td>0.318</td></tr><tr><td>95.0</td><td>0.0445</td><td>107.5</td><td>0.172</td><td>120.0</td><td>0.360</td></tr><tr><td>97.5</td><td>0.0645</td><td>110.0</td><td>0.205</td><td colspan="2"></td></tr><tr><td colspan="6"></td></tr><tr><th colspan="2"></th><th colspan="2">16O + 208Pb [128]</th><th colspan="2"></th></tr><tr><td>72</td><td>143.5</td><td>82</td><td>49.3</td><td>92</td><td>16.3</td></tr><tr><td>74</td><td>116.0</td><td>84</td><td>39.6</td><td>94</td><td>13.0</td></tr><tr><td>76</td><td>93.8</td><td>86</td><td>31.8</td><td>96</td><td>10.4</td></tr><tr><td>78</td><td>75.8</td><td>88</td><td>25.5</td><td>98</td><td>8.3</td></tr><tr><td>80</td><td>61.2</td><td>90</td><td>20.4</td><td>100</td><td>6.6</td></tr></tbody></table>
# X TSRT Parameters and Constants Used in All Figures and Calculations

*Supports* all sections; see calibration links in Appendix E.

MATLAB: `tsrtconstants.m`.

This appendix summarizes all TSRT parameters and constants used to generate the figures, tables, and numerical comparisons in this paper. The quantities fall into two categories:

1. **Global TSRT constants** — universal geometric quantities and scaling factors that enter all deterministic calculations (binding, lifetime, curvature, and action scaling).

2. **Fit-specific parameters** — localized curvature, screening, or hindrance parameters used only in the phenomenological figure fits (Figures 2, 3, 6).

All constants and parameters are traceable to the MATLAB implementations cited in Appendix W. Experimental datasets used for fitting are provided in machine-readable form in the accompanying repository (`data/`).

## (I) Global TSRT Constants and Scaling Factors

<table><caption>Table 40: Universal TSRT constants and scaling parameters used in all calculations.</caption><thead><tr><th>Symbol</th><th>Value</th><th>Physical meaning and usage</th></tr></thead><tbody><tr><td>γ</td><td>1.424 × 10<sup>-23</sup> m<sup>-1</sup></td><td>Curvature-energy coupling constant linking local geometric amplitude to energy density (Section 4).</td></tr><tr><td>ħ*</td><td>1.055 × 10<sup>-34</sup> Js</td><td>Effective Planck constant emerging from trembling-spacetime averaging; used in all action-based formulations.</td></tr><tr><td>τ<sub>0</sub></td><td>3.77 × 10<sup>-24</sup> s</td><td>Fundamental trembling period defining causal correlation scale.</td></tr><tr><td>A<sub>1</sub></td><td>4.21 × 10<sup>-15</sup></td><td>Dimensionless ground-mode amplitude used in curvature normalization (see Appendix A).</td></tr><tr><td>C<sub>E2</sub></td><td>1.318 × 10<sup>-9</sup> s</td><td>E2 calibration constant from the <sup>156</sup>Gd anchor (Table 36 (p. 205)); used in all lifetime computations.</td></tr><tr><td>C<sub>γ</sub></td><td>2.41 × 10<sup>-3</sup></td><td>Curvature-photon emission proportionality constant linking trembling amplitude to radiative decay rate.</td></tr><tr><td>C<sub>r</sub></td><td>1.00 (normalized)</td><td>Radial scaling constant used for curvature-radius normalization in all binding-energy plots.</td></tr></tbody></table>

These constants are global across all MATLAB pipelines: `tsrtbindingenergy.m`, `bindingplotwithexp.m`, `lifetime_pipeline.m`, and `gammaRates.m`. They ensure that every figure (including energy, radius, and lifetime plots) is produced using a consistent geometric calibration.

## (II) Fit-Specific TSRT Parameters Used in Figures

The following parameters correspond to localized fits in Figure 2 (p. 12), Figure 3 (p. 24), and Figure 6 (p. 86). Each retains a direct geometric meaning within the TSRT formalism rather than serving as an empirical adjustment. See Appendix W for dataset provenance.

### Interpretive summary.
<table><caption>Table 41: TSRT parameters used for reproducing the figure-specific fits.</caption><tbody><tr><th colspan="2">Screening (d(d,p)t in Pd host)</th></tr><tr><td>ΔU<sub>0</sub></td><td>89.84 eV</td></tr><tr><td>E<sub>c</sub></td><td>4.474 keV</td></tr><tr><td>β</td><td>0.108</td></tr><tr><td>λ</td><td>0.862 keV<sup>-1/2</sup></td></tr><tr><td>C<sub>η</sub></td><td>0.900</td></tr><tr><td>ε</td><td>-0.300</td></tr><tr><td>κ</td><td>2.062 keV<sup>-1/2</sup></td></tr><tr><td>d<sub>dyn</sub></td><td>2.66 × 10<sup>-3</sup> keV<sup>-1/2</sup></td></tr><tr><th colspan="2">Odd-Even Staggering (<sup>235</sup>U(n<sub>th</sub>,f))</th></tr><tr><td>Z<sub>c</sub></td><td>50.016</td></tr><tr><td>k<sub>1</sub></td><td>0.1456 L<sub>1</sub> = 20.68</td></tr><tr><td>k<sub>2</sub></td><td>0.1884 L<sub>2</sub> = 30.00</td></tr><tr><td>φ</td><td>0 B = -0.1033 Z<sub>0</sub> = 41.466 L<sub>0</sub> = 2.911</td></tr><tr><td>(c<sub>0</sub>, c<sub>1</sub>, c<sub>2</sub>)</td><td>(8.087, -0.600, 0.200)</td></tr><tr><td>(d<sub>0</sub>, d<sub>1</sub>, d<sub>2</sub>)</td><td>(-8.083, -0.0579, 21.461)</td></tr><tr><th colspan="2">Fusion Hindrance (<sup>64</sup>Ni+<sup>64</sup>Ni baseline)</th></tr><tr><td>α<sub>L</sub></td><td>3.927</td></tr><tr><td>β<sub>L</sub></td><td>-0.0396</td></tr><tr><td>Asym. renormalization</td><td>(a<sub>1</sub>, b<sub>1</sub>) = (-1.0, -1.2)</td></tr><tr><td>Sym. tweak</td><td>(c<sub>1</sub>, c<sub>2</sub>) = (-0.6, 0.2)</td></tr></tbody></table>

* Z<sub>c</sub>: central scission charge localizing curvature neck formation.

* k<sub>1,2</sub>: eigen-wavenumbers of parity-Bessel neck packets controlling the oscillatory pattern of odd-even staggering.

* L<sub>1,2</sub>: damping lengths quantifying neck width and curvature attenuation.

* (c<sub>n</sub>, d<sub>n</sub>): mode weights governing relative excitation of radial and azimuthal trembling submodes.

* (B, Z<sub>0</sub>, L<sub>0</sub>): localized curvature-bias term producing suppression near Z ≈ 40.

* (ΔU<sub>0</sub>, E<sub>c</sub>, β): static and dynamic screening parameters governing causal roll-off of the potential.

* (λ, ε, κ): trembling-mode phase-lag and dephasing parameters coupling local curvature to propagation amplitude.

* (C<sub>η</sub>, d<sub>dyn</sub>): dynamic coupling and incremental energy-dependent damping factors.

* (α<sub>L</sub>, β<sub>L</sub>, a<sub>1</sub>, b<sub>1</sub>, c<sub>1</sub>, c<sub>2</sub>): universal logarithmic slopes and symmetry renormalizations determining fusion hindrance scaling across nuclei.

Together, the global constants of Table 40 (p. 208) and the fit-specific parameters of Table 41 (p. 209) constitute the complete numerical foundation for all TSRT calculations, figures, and comparisons presented in this work.
# Y MATLAB Procedures and Reproducibility Workflow

## Reproducibility Data Bundle Available

All scripts listed in this appendix can be run *ab initio* using only the definitions and algorithms contained in this paper. The minimal experimental inputs (binding energies, lifetime benchmarks, stability-map grids, and calibrated TSRT constants) can, in principle, be reconstructed directly from the main text.

As a courtesy to readers and reviewers, a compact data bundle is also provided, containing all six required input files plus a `README.txt`:

* `exp_bindings.csv`

* `tsrt_bindingscan_clean.csv`

* `tsrt_benchmark_diag.csv`

* `tsrt_DeltaK2_map.csv`

* `tsrtconstantsCalibrated.mat`

* `tsrt_stabilitymap_data.mat`

* `README.txt` (usage instructions)

With this bundle and the `.m` files below, every figure, table, and benchmark in the paper can be regenerated without external databases or manual data entry.

## Data Availability

All numerical data used in this work—including binding energies, benchmark lifetime inputs, stability-map curvature grids, and calibrated TSRT constants—are fully reproducible from the algorithms and appendices contained in this paper. For convenience, the minimal input set required to run all MATLAB workflows has been deposited on Zenodo and is publicly available at [35]:

https://doi.org/10.5281/zenodo.17637889

The archive contains six data files in `.csv` and `.mat` format, plus a `README.txt` with usage instructions. No external databases are required, and all figures and tables can be regenerated directly from the provided files and the MATLAB codes listed in this Appendix.

*Consistency with implementation.* The script `tsrt_make_all.m` regenerates all figures/tables by calling domain codes (`tsrtfissionenergy.m`, `tsrtfusionenergy.m`, `tsrtsfactor.m`) and utilities (`tsrtmakegrid.m`, `tsrtconvergence.m`, `tsrt_repository_sanity_check.m`); filenames in captions/tables match the outputs from these entry points.

Entry points: `tsrt_make_all.m` (regenerates all figures/tables), domain scripts `tsrtfissionenergy.m`, `tsrtfusionenergy.m`, `tsrtsfactor.m`, utilities `tsrtmakegrid.m`, `tsrtconvergence.m`, `tsrt_repository_sanity_check.m`. Filenames are also indexed below.

Each listing supports results discussed in the main sections on binding (Section 5), fusion (Section 7), fission (Section 6.2), and lifetimes (Section 9.11). Outputs correspond to Tables 3 (p. 11)–2 (p. 11) in Appendix W.

## Quick Start for Full Reproduction.
Complete workflow in three steps

1. Set up the environment:

Place all MATLAB files (`.m`) in one working folder together with the `data/` subfolder
provided in the supplementary archive. Open MATLAB in that directory.

2. Initialize and calibrate once:

(i) Energy normalization (deuteron):

>> TSRT_Calib_NormEn_Deuteron

(ii) Scission smoothness for $^{235}$U:

>> TSRT_Calib_Fission_ScissWth_U235

Both steps persist to `tsrtconstantsCalibrated.mat` and leave no per-observable
knobs.

3. Regenerate the paper artifacts:

>> runall

This orchestrates binding, OES, screening/hindrance, gamma/lifetimes, and fission
tables/figures; outputs go to results/ and paper-assets-tex/. Individual mod-
ules can be executed separately where indicated below.

This appendix enumerates all MATLAB scripts and functions required to reproduce every table
and figure in the paper. Each block lists purpose, usage in the paper, inputs/outputs, and its
place in the dependency graph. All code blocks use `lstlisting` with unique labels; insert the
verified MATLAB code where indicated.

Execution order (overview).

1. **Global constants**: 2 (p. 213) (Appendix Y.2)

2. **Fission calibration & summary**: 29 (p. 263), 28 (p. 257), 34 (p. 271) (Appendix Y.15)

3. **Binding + radii**: 8 (p. 221), 7 (p. 218) (Appendix Y.4)

4. **Fusion sector (screening/hindrance)**: 9 (p. 223), 10 (p. 225), 11 (p. 229) (Appendix Y.6)

5. **OES (fission yields)**: 13 (p. 231) (Appendix Y.7)

6. **Gamma + lifetimes**: 15 (p. 235), 14 (p. 234), 16 (p. 238) (Appendix Y.8, Appendix Y.9)

7. **Convergence/robustness (diagnostics)**: 6 (p. 217) (Appendix Y.3)

8. **Parameter scans (optional)**: 36 (p. 278) (Appendix Y.17)

9. **Master orchestrator**: 38 (p. 284) (Appendix Y.18)

Datasets referenced below appear in Appendix W (verbatim if not large compilations) and in
the data/ folder.

Y.1 MATLAB Code Index and Thematic Map

This map groups all code by theme and points to the concrete listings and the main-text locations
where outputs are used.
## Core and Constants

Listings: 2 (p. 213), 3 (p. 214), 4 (p. 215), 5 (p. 216).

Used throughout; initializes `tsrtconstants.mat`. The global normalization `C.norm_energy` is fixed once by the deuteron anchor (Listing 44 (p. 297)).

## Binding and Radii

Listings: 8 (p. 221), 7 (p. 218). Data: `tsrt_bindingscan_cln.csv`, `exp_bindings.csv`. Used in Section 5, Figure 1 (p. 10) (left), Table 3 (p. 11).

## Fusion and Screening

Listings: 9 (p. 223), 10 (p. 225), 11 (p. 229), 12 (p. 230). Used in Sections 12.8, 12.9, Figure 2 (p. 12), Figure 6 (p. 86).

## Fission Energetics and OES

Listings: 29 (p. 263), 28 (p. 257), 34 (p. 271), 13 (p. 231). Used in Section 6 (OES; Figure 3 (p. 24)) and Section 6.2 (energetics; Table 42 (p. 257)). The summary table is produced from `fission_results.mat` via `fissiontablemake.m`.

## Gamma Calibration and Lifetimes

Listings: 15 (p. 235), 14 (p. 234), 16 (p. 238). Used in Section 9.7, Tables 2 (p. 11), 18 (p. 110).

*Workflow:* Run Listing 15 (p. 235) once to populate `C GAMMA E2` in `tsrtconstants Calibrated.mat`; then Listing 16 (p. 238) regenerates the paper tables without per-transition tuning.

## Convergence and Robustness (Diagnostics)

Listings: 6 (p. 217). Used to validate numerical stability of curvature/energy integrals; diagnostics summarized in Appendix Y.3.

## Orchestrators and Scans

Listings: 36 (p. 278), 38 (p. 284). Run batch sensitivity scans (optional) and reproduce the full paper in one pass.
function G = tsrtmakegrid(C)

%TSRTMAKEGRID Build spherical grids and quadrature weights for TSRT integrals.

% G = tsrtmakegrid(C) returns a struct with fields for r, theta, phi and

% their corresponding quadrature weights. Designed so that 3D integrals of

% the form $\int f(r,\theta,\phi) r^2 \sin\theta dr d\theta d\phi$ are computed by simple weighted sums.

%

% Inputs:

% C.r_max [m] radial cutoff

% C.N_r # of radial points

% C.N_theta # of polar points

% (optional) C.Ntheta.||phi|| # of azimuthal points (default: max(128, 2*N_theta))

%

% Outputs (all column vectors):

% G.r [N_r x 1] radial nodes (0 .. r_max)

% G.w_r [N_r x 1] radial weights for $\int_0^{r_{max}} dr$ (no r^2 included)

% G theta [N_theta x 1] polar nodes (0 .. pi)

% G.w_th [N_theta x 1] polar weights for $\int_0^{r_{max}} \sin\theta d\theta$

% G phi [N phi x 1] azimuthal nodes (0 .. 2pi)

% G.w_ph [N phi x 1] azimuthal weights for $\int_0^{2\pi} d\phi$

% G.dr, G.dphi (step sizes for convenience)

% G theta_rule (string describing the polar rule used)

%

% Notes:

% • If the file 'lgwt.m' is on the same path, the polar quadrature uses
% Gauss-Legendre via x = cosθ; weights integrate the sinθ dθ measure.

% • Otherwise, a uniform trapezoidal rule is used with the sinθ factor folded
% into the weights.

% • To integrate F(r,θ,φ) over spherical volume: sum(F .* r.^2 .* w_r .* (w_th).'
%      .* (w_ph).')
% with broadcasting or by nested sums.

% ----------------------------

% Resolve azimuthal size

% ----------------------------

if isfield(C, 'N phi') && ~isempty(C.N phi)
    Nphi = C.N phi;
else
    Nphi = max(128, 2*C.N theta); % sensible default
end

% ----------------------------

% Radial grid (uniform) + trapezoidal weights

% ----------------------------

G.r = linspace(0, C.r_max, C.N_r).';

G.dr = G.r(2) - G.r(1);

G.w_r = trap weights(C.N_r, G.dr); % integrates $\int_0^{r_{max}} dr$ (no r^2)

% ----------------------------

% Polar grid (θ ∈ [0, π])

% Prefer Gauss-Legendre on x = cosθ; weights integrate sinθ dθ

% ----------------------------

if exist('lgwt', 'file') == 2
    % Gauss-Legendre nodes/weights on [-1, 1]
    [xg, wg] = lgwt(C.N_theta, -1, 1); % xg ascends from -1 to 1
    theta = acos(xg(:));          % map x -> θ; this maps [-1,1] to [π,0]
    % Put θ ascending 0..π for readability
    [theta, idx] = sort theta, 'ascend';
```matlab
wth = wg(:); wth = wth(idx); % weights integrate F(θ) under sinθ dθ via x
=cosθ
G theta = theta;
G.w_th = wth;
G theta_rule = 'Gauss-Legendre on x=cosθ (sinθ dθ measure)';
else
    % Uniform θ with trapezoidal weights; fold sinθ into the weights
    G theta = linspace(0, pi, C.N_theta).';
    dth = G theta(2) - G theta(1);
    G.w_th = trap weights(C.N_theta, dth) .* sin(G theta);
    G theta_rule = 'Uniform (trap) with sinθ factor';
end

% ----------------------------
% Azimuthal grid (φ ∈ [0, 2π]) + trapezoidal weights (periodic)
% ----------------------------
G phi = linspace(0, 2*pi, Nphi).';
G.dphi = G phi(2) - G phi(1);
G.w_ph = trap weights(Nphi, G.dphi);

% ----------------------------
% Backward-compatibility aliases (for older code)
% ----------------------------
G.wtheta = G.w_th; % legacy
G.w_theta = G.w_th; % explicit snake_case
G.w phi = G.w_ph;

end % tsrtmakegrid

% =========================
% Local helper: trapezoidal weights
% =========================
function w = trap weights(n, h)
%TRAP_WEIGHTS Composite trapezoid weights for uniform grid.
% Returns column vector of length n; end weights are h/2, interior h.
w = ones(n,1) * h;
w(1) = h/2;
w(end) = h/2;
end
```

<caption>Listing 3: Production grid builder</caption>

```
tsrtcurvaturetensor.m

function K = tsrtcurvaturetensor(G, Xi, C)
% Scalar contraction K_{\mu\nu} u^\mu u^\nu on the grid G.
% For the spherically symmetric, time-averaged profiles used here the
% contraction reduces to the scalar profile used in integrals.
%
% INPUTS:
% G - grid struct
% Xi - trembling field struct with Xi.profile
% C - constants (kept for interface consistency)
%
% OUTPUT:
% K - 3D array of contraction values on G

K = Xi.profile;
```
end

<caption>Listing 4: Curvature tensor evaluation</caption>

tsrtenergyfromcurvature.m

```matlab
function E = tsrtenergyfromcurvature(G, K, C)
% Integrate curvature measure to energy using TSRT normalization.
% INPUTS:
% G - grid struct with r,theta,phi and separable quadrature weights
%    (G.w_r for dr, G.w_th integrating sin theta)dtheta, G.w_ph for dphi)
% K - curvature contraction array from tsrtcurvaturetensor (same size as grid)
% C - constants struct with normalization (C.norm_energy, in J per unit integral)
% OUTPUT:
% E - energy in Joules (SI)

% Volume factor: r^2 only (sin theta) is already in G.w_th)
[R,~,~] = ndgrid(G.r, G theta, G.phi);
vol = R.^2;

% Separable quadrature weights (already include the step sizes)
W = reshape(G.w_r, [], 1, 1) .* reshape(G.w_th, 1, [], 1) .* reshape(G.w_ph, 1, 1, []);

% Weighted integral (dimensionless curvature units)
E_int = sum( (K .* vol) .* W, 'all' );

% Map curvature integral to energy (calibrated scale)
E = C.norm_energy * E_int;
end
```

<caption>Listing 5: Energy from curvature integral</caption>
## Y.3 Convergence and Robustness Scans

Purpose: Validate numerical stability and robustness of curvature/action integrations across grid resolutions.

tsrtconvergence.m

function tsrtconvergence()

% Grid-convergence test using a smooth spherical Gaussian profile.

% This avoids dependencies on other TSRT routines and checks numerical stability.

C = tsrtconstants('norm_energy', 1.0); % <-- allow running before calibration

% radial resolutions to test

Nr_list = [256, 512, 1024, 2048];

Q = zeros(size(Nr_list));

for k = 1:length(Nr_list)

C.N_r = Nr_list(k);

G = tsrtmakegrid(C);

% Spherical Gaussian test field (dimensionless)

r0 = 1.0e-15; % 1 fm typical nuclear scale

[RR,TT,PP] = ndgrid(G.r, G theta, G.phi);

Xi.profile = exp(-(RR./r0).^2);

% Curvature contraction (baseline: profile itself)

K = tsrtcurvaturetensor(G, Xi, C);

% Energy-like integral (Joules) via calibrated normalization

Q(k) = tsrtenergyfromcurvature(G, K, C);

end

% relative changes vs previous resolution

eps = abs(diff(Q)/Q(2:end));

disp('Convergence test for radial resolution:');

for k = 2:length(Nr_list)

fprintf('Nr = %d vs %d: rel. change = %.3e
', ... Nr_list(k), Nr_list(k-1), eps(k-1));

end

end

<caption>Listing 6: Convergence / robustness scanner</caption>
## Y.4 Binding-Energy and Radius Pipeline

Dependencies: Listings 2 (p. 213), 3 (p. 214), 8 (p. 221).

Purpose: Compute TSRT binding energies/radii and overlay with AME2020 data. Used in paper: Section 5; Table 3 (p. 11); Figure 1 (p. 10) (left).

Outputs: `binding_results.mat`, `tsrt_bindingsplot_with_exp.pdf`, `bindingplotwithexp.m`

```matlab
% FILE: bindingplotwithexp.m
% PURPOSE: produce tsrt_bindingsplot_with_exp.pdf
function bindingplotwithexp(tsrt_csv, ame_csv, out_fig, Z_select)
% bindingplotwithexp
% Plot TSRT Ebind/A curves per isotopic chain (by Z) and overlay many experimental points.
%
% Now supports zero-arg call:
% bindingplotwithexp
% which auto-uses:
% tsrt_csv = 'tsrt_bindingscan_cln.csv'
% ame_csv = 'exp_bindings.csv'
% out_fig = 'tsrt_bindingsplot_with_exp'
%
% See comments inside for details.
```

```matlab
% ---- Defaults if no-arg or missing args ----
if length < 1 || isempty(tsrt_csv)
    if exist('tsrt_bindingscan_cln.csv', 'file')==2
        tsrt_csv = 'tsrt_bindingscan_cln.csv';
    else
        error(['Provide TSRT CSV (e.g. tsrt_bindingscan_cln.csv). ', ...
                'Default not found in current folder.']);
    end
end
if length < 2 || isempty(ame_csv)
    if exist('exp_bindings.csv', 'file')==2
        ame_csv = 'exp_bindings.csv';
    else
        error(['Provide experimental CSV (e.g. exp_bindings.csv). ', ...
                'Default not found in current folder.']);
    end
end
if length < 3 || isempty(out_fig)
    out_fig = 'tsrt_bindingsplot_with_exp';
end
```

```matlab
% --- Read TSRT data
assert(isfile(tsrt_csv), 'TSRT CSV not found: %s', tsrt_csv);
Tt = readtable(tsrt_csv, 'FileType', 'text');
Tt = normalize_headers(Tt);

Zt = find_col(Tt, ["z");
At = find_col(Tt, ["a");
Et_perA = find_col(Tt, ["ebind_pera_mev", "epera_mev", "binding_per_a");
Et = find_col(Tt, ["ebind_mev", "binding_mev", "e_mev");

assert(~isempty(Zt) && ~isempty(At), 'TSRT CSV must contain Z and A.');
if isempty(Et_perA)
    assert(~isempty(Et), 'TSRT CSV must have Ebind_perA_MeV or Ebind_MeV.');
```
```R
Tt.Ebind_perA_MeV = Tt.(Et) ./ Tt.(At);
else
  Tt.Ebind_perA_MeV = Tt.(Et_perA);
  bad = ~isfinite(Tt.Ebind_perA_MeV) | Tt.Ebind_perA_MeV<=0;
  if any(bad) && ~isempty(Et)
    Tt.Ebind_perA_MeV(bad) = Tt.(Et)(bad) ./ Tt.(At)(bad);
end
end

% --- Read experimental data
assert(isfile(ame_csv), 'Experimental CSV not found: %s', ame_csv);
Te = readtable(ame_csv, 'FileType', 'text');
Te = normalize_headers(Te);

Ze = find_col(Te, ["z");
Ae = find_col(Te, ["a");
Ne = find_col(Te, ["n"); % optional
EbA = find_col(Te, ["ebind_pera_mev", "epera_mev", "binding_per_a", "be_a_mev", "bindingenergypernucleon");
Eb = find_col(Te, ["ebind_mev", "binding_mev", "e_mev");
Mu = find_col(Te, ["atomicmass_u", "mass_u", "mass_atom_u");

assert(~isempty(Ze), 'Experimental CSV must contain Z.');
if isempty(Ae)
  assert(~isempty(Ne), 'Experimental CSV must have A or N (with Z).');
  Te.a = Te.(Ze) + Te.(Ne);
  Ae = 'a';
end

% Compute Ebind/A for experimental
if ~isempty(EbA)
  Te.Ebind_perA_MeV = Te.(EbA);
else if ~isempty(Eb)
  Te.Ebind_perA_MeV = Te.(Eb) ./ Te.(Ae);
else if ~isempty(Mu)
  u = 931.49410242; % MeV/u
  mp = 938.2720813; % MeV
  mn = 939.5654133; % MeV
  me = 0.51099895; % MeV
  Zvv = Te.(Ze);
  Avv = Te.(Ae);
  Nvv = Avv - Zvv;
  Muu = Te.(Muu);
  B = Zvv.*mp + Nvv.*mn - (Muu.*u) + Zvv.*me; % atomic-electron approx
  Te.Ebind_perA_MeV = B ./ Avv;
else
  error('Experimental CSV must provide Ebind_perA_MeV, Ebind_MeV, or AtomicMass_u.');
end

% --- Clean rows
Tt = Tt(isfinite(Tt.(Zt)) & isfinite(Tt.(At)) & isfinite(Tt.Ebind_perA_MeV), :);
Te = Te(isfinite(Te.(Ze)) & isfinite(Te.(Ae)) & isfinite(Te.Ebind_perA_MeV), :);

% --- Choose Z set (intersection so both have data)
Z_tsrt = unique(Tt.(Zt));
Z_exp = unique(Te.(Ze));
Z_all = intersect(Z_tsrt, Z_exp);
```
```R
if nargin >= 4 && !isempty(Z_select)
    Z_use = intersect unique(Z_select(:).'), Z_all(:).');
    assert(!isempty(Z_use), 'None of the requested Z are common to TSRT and
experimental data.');
else
    Z_use = Z_all(:).';
end
```

### % --- Figure
```R
fig = figure('Color','w','Position', [100 100 900 520], ...
              'Name','Ebind per A: TSRT + Experiment');
ax = axes('Parent',fig); hold(ax,'on');
set(ax, 'FontName','Helvetica', 'FontSize',11, 'LineWidth',1.0, ...
    'Box','on', 'XGrid','on', 'YGrid','on', ...
    'GridAlpha',0.20, 'MinorGridAlpha',0.15);

xlabel(ax,'Mass number A');
ylabel(ax,'Binding energy per nucleon E_{bind}/A (MeV)');
title(ax,'TSRT binding per nucleon across isotopic chains (with experiment)',
       'FontWeight','bold');

annotation(fig,' textbox', [0.12 0.02 0.76 0.05], 'String', ...
           'Experimental data: AME2020', ...
           'EdgeColor','none', 'HorizontalAlignment','center', 'FontSize',9, 'Interpreter',
           'none');
```

### % Style: TSRT solid lines; experimental = large black markers (per-Z shapes)
```R
markers = {'o','s','^','d','v','>','<','p','h','x','+'};
mCount = numel/markers);
```

### % Plot per Z
```R
for i = 1:numel(Z_use)
    Z = Z_use(i);

% TSRT curve
Ct = sortrows(Tt(Tt.(Zt)==Z, {At,'Ebind_perA_MeV'}), At);
if !isempty(Ct)
    plot(ax, Ct.(At), Ct.Ebind_perA_MeV, '-', 'LineWidth', 1.8, ...
        'DisplayName', sprintf('TSRT Z=%d', Z));
end

% Experimental markers
Ce = sortrows(Te(Te.(Ze)==Z, {Ae,'Ebind_perA_MeV'}), Ae);
if !isempty(Ce)
    mk = markers{mod(i-1, mCount)+1};
    plot(ax, Ce.(Ae), Ce.Ebind_perA_MeV, mk, ...
        'MarkerSize', 9, 'LineStyle','none', 'LineWidth', 1.2, ...
        'MarkerFaceColor','k', 'Color','k', ...
        'DisplayName', sprintf('EXP Z=%d', Z));
end
end
```

### % Axis limits: full A-range
```R
Aall = [Tt.(At); Te.(Ae)];
Eall = [Tt.Ebind_perA_MeV; Te.Ebind_perA_MeV];
if !isempty(Aall)
   lim(ax, [min(Aall)-1, max(Aall)+1]);
```
```matlab
end
if ~isempty(Eall)
    ymin = min(Eall); ymax = max(Eall);
    pad = max(0.25, 0.05*(ymax - ymin));
    ylim(ax, [ymin - pad, ymax + pad]);
end

legend(ax, 'Location', 'eastoutside', 'Interpreter', 'none', 'Box', 'off');

% Save PDF + EPS
[outDir, outBase, ~] = fileparts(out_fig);
if isempty(outDir), outDir = '.'; end
if ~exist(outDir, 'dir'), mkdir(outDir); end
pdfPath = fullfile(outDir, [outBase '.pdf']);
epsPath = fullfile(outDir, [outBase '.eps']);

try
    exportgraphics(fig, pdfPath, 'ContentType', 'vector');
catch
    set(fig, 'PaperPositionMode', 'auto');
    print(fig, pdfPath, '-dpdf', '-painters');
end
try
    print(fig, epsPath, '-depsc', '-painters');
catch MEeps
    warning('EPS export failed: %s', MEeps.message);
end
fprintf('Wrote figure: %s and %s
', pdfPath, epsPath);
fprintf('Experimental source: AME2020
');
end

% ---------- helpers ----------
function T = normalize_headers(T)
    raw = T.Properties.VariableNames;
    norm = lower(strrep/raw, '--', '_'));
    norm = strrep/norm, '.', '_');
    T.Properties.VariableNames = matlab.lang.makeUniqueStrings/norm);
end

function cname = find_col(T, candidates)
    cname = '';
    cols = lower(T.Properties.VariableNames);
    for c = string(candidates)
        m =strcmp cols, char(c));
        if any(m)
            cname = T.Properties.VariableNames{find(m,1,'first')};
            return;
        end
    end
end
end
```

<caption>Listing 7: Binding plot with experimental overlay</caption>

Figures produced correspond to Figure 1 (p. 10) (left).

tsrtbindingenergy.m

```matlab
function Ebind = tsrtbindingenergy(Xi_list, G, C)
% TSRT binding energy estimator.
% Inputs:
```
```c
% Xi_list : cell array of nucleon fields (isolated components)
% G         : grid struct from tsrtmakegrid
% C         : constants struct from tsrtconstants
%
% Output:
% Ebind     : binding energy (J)

% Isolated nucleon energies:
E_iso = 0;
for k = 1:numel(Xi_list)
    Kk = tsrtcurvaturetensor(G, Xi_list{k}, C);
    E_iso = E_iso + tsrtenergyfromcurvature(G, Kk, C);
end

% Bound configuration energy (superpose fields deterministically):
Xi_bound = struct('profile', zeros(numel(G.r),numel(G theta),numel(G.phi)));
for k = 1:numel(Xi_list)
    Xi_bound.profile = Xi_bound.profile + Xi_list{k}.profile;
end
Kb = tsrtcurvaturetensor(G, Xi_bound, C);
E_bound = tsrtenergyfromcurvature(G, Kb, C);

% Binding energy difference:
Ebind = (E_iso - E_bound);
end
```

<caption>Listing 8: TSRT binding-energy core routine</caption>
## Y.5 Electron Screening in $d(d,p)t$ (Pd host)

Dependencies: Listings 2 (p. 213), 3 (p. 214), 12 (p. 230).

Purpose: Fit Pd-host $d(d,p)t$ S-factor enhancement; generate Figure 2 (p. 12).

screeningmakefig.m

```R
% FILE: screeningmakefig.m
% PURPOSE: fit/plot Pd d(d,p)t screening; save figScreening.pdf
function screeningmakefig(outFile)
% Data provenance: Appendix~\ref{app:t:dd} (Pd host, Raiola et al.)
% Model: Equations (tsrt-causal-screen)-(tsrt-screen-tremble)

% Try shared constants; fall back to explicit SI if unavailable
useFallback = true;
if exist('tsrtconstants','file') == 2
  try
    C = tsrtconstants(); % will error if norm_energy not set yet
    e = C.eC;         % Coulombs
    hbar = C.hbar_SI;  % J*s
    eps0 = C.eps0;    % F/m
    ke = 1/(4*pi*eps0);
    amu = C.amu_kg;    % kg
    useFallback = false;
  catch
    % fallback below
  end
end
if useFallback
  e=1.602176634e-19; hbar=1.054571817e-34; eps0=8.8541878128e-12;
  ke=1/(4*pi*eps0); amu=1.66053906660e-27;
end

if nargin<1, outFile='figScreening.pdf'; end

E_keV = [3.00 3.50 4.00 4.50 5.00 5.50 6.00 6.50 7.00 7.50 8.00 8.50 9.00 10.00]';
Sobs = [60.50 58.80 57.50 56.80 56.20 55.90 55.70 55.60 55.50 55.45 55.40 55.40 55.40 55.40]';
S_bare = [55.20 55.40 55.50 55.60 55.70 55.75 55.80 55.80 55.80 55.80 55.80 55.80 55.80 55.80 55.80]';
Robs = Sobs ./ S_bare;

% Sommerfeld-parameter helper using the constants already chosen above
Z1=1; Z2=1; m_d amu=2.013553; mu=(m_d amu/2)*amu;
etaFun = @ (EkeV) (Z1*Z2*ke*e^2) ./ (hbar.*sqrt(2*(EkeV*1e3*e)/mu));

% TSRT model components
DeltaUe = @ (E,DUO,Ec,beta,lam) DUO ./ ((1 + (E./Ec).*beta).* (1 + lam.*sqrt(E)));
R_full = @ (E,DUO,Ec,beta,Ceta,lam,epsi,kappa,ddyn) ...
  exp( pi*Ceta.*etaFun(E) .* ( DeltaUe(E,DUO,Ec,beta,lam).* (1+epsi.*besselj(0,kappa.*sqrt(E))) ./ (E*1e3) ) ...
  - ddyn.*sqrt(E) );

% Fit all parameters in log-space
p0 = [300, 6, 1.5, 1.00, 0.30, -0.08, 1.2, 0.02]; % [DUO,Ec,beta,Ceta,lam,epsi,kappa,ddyn]
lb = [ 10,0.2,0.5, 0.90, 0.00, -0.30,0.10, 0.00];
ub = [1000, 50, 5.0, 1.10, 2.00, 0.00,6.00, 1.00];
```
```r
obj = @p) log(R_full(E_keV,p(1),p(2),p(3),p(4),p(5),p(6),p(7),p(8))) - log(R_ob);
opts = optimoptions('lsqnonlin','Display','off');
p = lsqnonlin(obj,p0,lb,ub,opts);

% --- Plot + export (both PDF and EPS) ----------------------------------------
Egrid = linspace(min(E_keV), max(E_keV), 500).';

f = figure('Color','w','Name','Electron Screening','Position', [120 120 900 520]);
ax = axes('Parent',f); hold(ax,'on');

% TSRT house style (match bindingplot)
set(ax, 'FontName', 'Helvetica', 'FontSize', 11, 'LineWidth', 1.0, ...
    'Box','on', 'XGrid','on', 'YGrid','on', 'GridAlpha', 0.20, 'MinorGridAlpha',
    0.15);

% Data + model
semiology(ax, E_keV, R_ob, 'ko', 'MarkerFaceColor', 'k', 'DisplayName', 'Data');
semiology(ax, Egrid, R_full(Egrid, p(1),p(2),p(3),p(4),p(5),p(6),p(7),p(8)), ...
    'r-', 'LineWidth', 2, 'DisplayName', 'TSRT causal+trembling+dephasing');

% Labels, legend, title
xlabel(ax,'E_{c.m.} [keV]');
ylabel(ax,'S_{obs}/S_{bare}');
title(ax,'Electron Screening','FontWeight','bold');
legend(ax,'Location','best');

% Export (PDF + EPS)
set(f,'PaperPositionMode','auto');
try
    exportgraphics(ax, 'figScreening.pdf', 'ContentType', 'vector');
catch
    print(f, '-dpdf', 'figScreening.pdf');
end
print(f, '-depsc', 'figScreening.eps');
end
```

<caption>Listing 9: Screening S-factor figure</caption>

Figures produced correspond to Figure 2 (p. 12).
Y.6 Heavy-Ion Fusion Hindrance (S-factor)

Dependencies: Listings 2 (p. 213), 3 (p. 214), 10 (p. 225).

Purpose: Universal slope extraction and figure generation for hindrance.

tsrtsfactor.m

% FILE: tsrtsfactor.m
% PURPOSE: model/fit S(E) scaling for hindrance systems
function tsrtsfactor(reaction, E_MeV, alpha, beta, out_csv, out_fig, opts)
% tsrtsfactor TSRT S-factor using universal slopes and a single per-system amplitude.
% Universal hindrance slopes (alpha, beta) for TSRT (from the paper's hindrance
section):
% Calibrated once on 64Ni+64Ni (upper-energy-third):
% alpha = 2.815e-2
% beta = -1.148e-5
% Pass these as the (alpha, beta) arguments when calling this function.
% ----------------------------------------------------------------------
% Usage (prediction only, A_sys=1):
% tsrtsfactor(struct('Z1',8,'A1',16,'Z2',8,'A2',16), E, alpha, beta, ...
% 's_out.csv','s_plot.pdf');
% ----------------------------------------------------------------------
% Usage (fit A_sys to experimental S(E) on upper-energy third):
% opts.S_exp = Sexp; % same size as E
% opts.Units = 'MeV-b'; % label only; not used in calc
% opts.a1 = -1.0; opts.b1 = -1.2; % asymmetry renormalization (optional)
% opts.c1 = -0.6; opts.c2 = 0.2; % symmetric-only tweak (optional)
% opts.Rcal_A = 64; % R_cal = A_cal^(1/3)+A_cal^(1/3); here 64Ni+64Ni
% tsrtsfactor(react, E, alpha, beta, 's_out.csv','s_plot.pdf', opts);
% ----------------------------------------------------------------------
% INPUTS
% reaction : struct with integer fields Z1,A1,Z2,A2
% E_MeV : column vector of center-of-mass energies (MeV)
% alpha,beta : universal slopes (from a single calibration)
% out_csv : output CSV filename
% out_fig : output figure filename (.pdf, .png, ...)
% opts : (optional) struct with fields:
% - S_exp : experimental S(E) values (same size as E_MeV)
% - units : string, for legend/labels (e.g., 'MeV-b')
% - a1,b1 : asymmetry renormalization coefficients (default 0)
% - c1,c2 : symmetric-only tweak coefficients (default 0)
% - Rcal_A : A_cal for R_cal = A_cal^(1/3)+A_cal^(1/3) (default: [])
% - title : plot title override
% - label : dataset label (e.g., '^{64}Ni+^{64}Ni')
% ----------------------------------------------------------------------
% OUTPUT
% CSV columns: E_MeV, S_TSRT, Xi1, Xi2, alpha_eff, beta_eff, A_sys
% Figure: log10 S(E) vs E, with optional experimental overlay.
% ----------------------------------------------------------------------
% THEORY (Appendix~\ref{app:s_factor}):
% Xi1(E) = sqrt(mu) * (Z1 Z2) / ((A1^(1/3)+A2^(1/3)) * sqrt(E))
% Xi2(E) = mu * (Z1 Z2)^(2/3) / ((A1^(1/3)+A2^(1/3))^(2/3) * E)
% S(E) = A_sys * exp(-alpha_eff * Xi1 - beta_eff * Xi2)
% alpha_eff = alpha * (1+a1*chi) * [1 + c1*(R_cal/R - 1)]
% beta_eff = beta * (1+b1*chi) * [1 + c2*(R_cal/R - 1)]
% chi = |A1 - A2|/(A1 + A2), R = A1^(1/3)+A2^(1/3)
% A_sys fitted from upper-energy third if S_exp provided; else A_sys=1.
% ----------------------------------------------------------------------
% --- Defaults and input checks ------------------------------------------
if nargin < 7, opts = struct(); end

% Columnize E
E_MeV = E_MeV();
if any(E_MeV <= 0)
    error('All energies must be positive (MeV).');
end

% Defaults for options
def = struct('S_exp', [], 'units', 'arb', ...
               'a1', 0.0, 'b1', 0.0, 'c1', 0.0, 'c2', 0.0, ...
               'Rcal_A', [], 'title', '', 'label', '');
opts = setdef(options, def);

% Validate experimental S if provided
if ~isempty(options.S_exp)
    if ~isequal(size(options.S_exp(:)), size(E_MeV))
        error('options.S_exp must have the same length as E_MeV.’');
    end
end

% --- Physical constants and invariants ----------------------------------------
amu_MeV = 931.49410242; % atomic mass unit in MeV/c^2 (CODATA-level)
Z1 = reaction.Z1; A1 = reaction.A1;
Z2 = reaction.Z2; A2 = reaction.A2;

% Reduced mass in MeV/c^2 (nuclear mass ~ A * amu)
m1 = A1 * amu_MeV;
m2 = A2 * amu_MeV;
mu = (m1*m2) / (m1 + m2);

% Geometric radius scale
R = A1^(1/3) + A2^(1/3); % dimensionless A^(1/3) sum
chi = abs(A1 - A2) / (A1 + A2);

% TSRT invariants Xi1, Xi2 (element-wise in E)
Zprod = Z1 * Z2;
Xi1 = sqrt(mu) * Zprod ./ (R * sqrt(E_MeV)); % ~ E^(-1/2)
Xi2 = mu * (Zprod^2) ./ ((R^2) * E_MeV); % ~ E^(-1)

% Symmetry renormalization and symmetric-only tweak
alpha_eff = alpha * (1 + options.a1 * chi);
beta_eff = beta * (1 + options.b1 * chi);

if ~isempty(options.Rcal_A)
    Rcal = options.Rcal_A^(1/3) + options.Rcal_A^(1/3);
    alpha_eff = alpha_eff * (1 + options.c1*(Rcal/R - 1));
    beta_eff = beta_eff * (1 + options.c2*(Rcal/R - 1));
end

% Baseline TSRT curve (without amplitude)
S_noAmp = exp(-alpha_eff*Xi1 - beta_eff*Xi2);

% --- Fit amplitude A_sys on upper-energy third if S_exp is available --------
if isempty(options.S_exp)
    A_sys = 1.0;
    S_TSRT = A_sys * S_noAmp;
else

% Sort by energy, select upper third
[E_sorted, idx] = sort(E_MeV, 'ascend');

n = numel(E_sorted);

k0 = max(1, floor(2*n/3)+1); % start of top third

sel = idx(k0:n);

Sex = opts.S_exp sel;

Sth0 = S_noAmp sel;

% Robust log-space intercept fit: log(Sexp) = log(A_sys) + log(S_noAmp)

logA = mean( log(Sex(:)) - log(Sth0(:)) ); % least-squares intercept

A_sys = exp(logA);

% Full TSRT prediction

S_TSRT = A_sys * S_noAmp;

end

% --- Write CSV ---

T = table(E_MeV, S_TSRT, Xi1, Xi2, ...
    repmat(alpha_eff,, 1), ...
    repmat beta_eff , 'numel(E_MeV)',1), ...
    repmat(A_sys, 'numel(E_MeV)',1), ...
    'VariableNames', {'E_MeV', 'S_TSRT', 'Xi1', 'Xi2', 'alpha_eff', 'beta_eff', 'A_sys
'});

writetable(T, out_csv);

fprintf('Wrote S-factor CSV: %s
', out_csv);

% --- Plot ---

fig = figure('Color','w','Position', [120 120 920 520], 'Name','TSRT S-factor');

ax = axes('Parent',fig); hold(ax,'on');

set(ax, 'FontName', 'Helvetica', 'FontSize', 11, 'LineWidth', 1.0, ...
    'Box', 'on', 'XGrid', 'on', 'YGrid', 'on', 'GridAlpha', 0.20, 'MinorGridAlpha
', 0.15);

% If experimental data present, overlay in log10

if ~isempty(options.S_exp)
    plot(ax, E_MeV, log10(options.S_exp), 'o', 'MarkerSize',5, ...
        'DisplayName', 'Experiment', 'Color', 'k', 'MarkerFaceColor', 'k');
end

% TSRT curve (log10)

plot(ax, E_MeV, log10(S_TSRT), '-', 'LineWidth',2, 'DisplayName', 'TSRT');

xlabel(ax,'E_{c.m.} (MeV)');

ylabel(ax, sprintf('log_{10} S(E) [%s]', optionsUnits));

% Title logic

if ~isempty(options.title)

title(ax, options.title, 'FontWeight', 'bold');

ELSEIF ~isempty(options.label)

title(ax, sprintf('TSRT S-factor: %s', options.label), 'FontWeight', 'bold');

ELSE

title(ax, sprintf('TSRT S-factor: Z_{1,2}=(%d,%d), A_{1,2}=(%d,%d)', Z1,Z2,A1,A2),
'FontWeight', 'bold');

end
legend(ax, 'Location', 'southwest', 'Box', 'off');

% Save figure

save figure(fig, out_fig);

% --- Print calibration summary to console ---

fprintf('Reaction: Z1=%d, A1=%d; Z2=%d, A2=%d
', Z1,A1,Z2,A2);

fprintf('alpha_eff = %.6g, beta_eff = %.6g, A_sys = %.6g
', alpha_eff, beta_eff, A_sys);

end % function

function save figure(fig, out fig)

%SAVE_FIGURE Save both PDF and EPS (and PNG as a bonus) for the given fig.

try

[outDir, outBase, ~] = fileparts(out_fig);

if isempty(outDir), outDir = '.'; end

if isempty(outBase), outBase = 'figure'; end

if ~exist(outDir,'dir'), mkdir(outDir); end

pdfPath = fullfile(outDir, [outBase, '.pdf]);

epsPath = fullfile(outDir, [outBase, '.eps]);

pngPath = fullfile(outDir, [outBase, '.png]);

% Save PDF (vector)

try

exportgraphics(fig, pdfPath, 'ContentType', 'vector');

catch

set(fig, 'PaperPositionMode', 'auto');

print(fig, pdfPath, '-dpdf', '-painters');

end

% Save EPS (vector)

try

print(fig, epsPath, '-depsc', '-painters');

catch MEeps

warning('EPS export failed: %s', MEeps.message);

end

% Optional PNG (raster)

try

print(fig, pngPath, '-dpng', '-r300');

catch MEpng

warning('PNG export failed: %s', MEpng.message);

end

fprintf('Wrote figure: %s, %s (and %s)
', pdfPath, epsPath, pngPath);

catch ME

warning('export/save failed: %s', ME.message);

% Last-resort PDF

try

print(fig, [out_fig '.pdf'], '-dpdf', '-painters');

catch

end

end
end

end

<caption>Listing 10: TSRT S-factor model (core)</caption>

reproducehindranceTSRT.m

% FILE: reproducehindranceTSRT.m
% PURPOSE: reproduce the hindrance (S-factor) figure and table using
%           universal slopes (alpha, beta) and a single per-system amplitude.
% OUTPUT: fig_hindrance.pdf (and .eps/.png as produced by the script)
% NOTES: Uses tsrtsfactor.m for modeling and overlaying data.
function reproducehindranceTSRT()
%REPRODUCEHINDRANCETSRT Generate the TSRT hindrance S-factor figure and data table.
%
% Description:
% Reproduces the canonical S-factor (hindrance) plot for a selected system
% using the published universal slope parameters (alpha, beta) and a
% single per-system amplitude. The results correspond to the figure and
% CSV data set reported in the paper.
%
% Outputs:
% fig_hindrance.pdf - generated figure (also .eps/.png if enabled)
% s_hindrance.csv - CSV file with S-factor curve and invariants
%
% Dependencies:
% tsrtconstants.m, tsrtsfactor.m
%
% Notes:
% The example below uses the symmetric 64Ni + 64Ni system as benchmark.
% Experimental S(E) data may optionally be supplied via opts.S_exp.

% --- Load constants
C = tsrtconstants();

% --- Universal slope parameters (from published TSRT fits)
alpha = 2.815e-2;
beta = -1.148e-5;

% --- Example system: 64Ni + 64Ni (benchmark configuration)
R = struct('Z1',28,'A1',64,'Z2',28,'A2',64);

% --- Example center-of-mass energy grid [MeV]
E = linspace(185, 250, 50).';

% --- Optional experimental data container
opts = struct();
opts.Units = 'MeV-b';
opts.label = '^{64}Ni+^{64}Ni';
opts.Rcal_A = 64; % calibration mass number

% --- Output file names (match published artifact names)
out_csv = 's_hindrance.csv';
out_fig = 'fig_hindrance.pdf';

% --- Generate the figure and CSV via the core S-factor routine
tsrtsfactor(R, E, alpha, beta, out_csv, out_fig, opts);
```matlab
fprintf('Hindrance figure written: %s\n', out_fig);
end
```

<caption>Listing 11: Hindrance figure producer (canonical)</caption>

Figures produced correspond to Figure 1 (p. 10) (right).

tsrtfusionenergy.m

```matlab
function Efus = tsrtfusionenergy(G, Xi_A, Xi_B, Xiprod, C)
% TSRT-based fusion energy calculator.
% Inputs:
%   G           - spatial grid (Appendix~\ref{app:code_integrators})
%   Xi_A       - trembling field struct for nucleus A
%   Xi_B       - trembling field struct for nucleus B
%   Xiprod     - trembling field struct for fused product
%   C           - constants struct (Appendix~\ref{app:codeConstants})
% Output:
%   Efus       - fusion energy release (Joules; convert externally to MeV)

% Reactant A
KA = tsrtcurvaturetensor(G, Xi_A, C);
EA = tsrtenergyfromcurvature(G, KA, C);

% Reactant B
KB = tsrtcurvaturetensor(G, Xi_B, C);
EB = tsrtenergyfromcurvature(G, KB, C);

% Product
Kp = tsrtcurvaturetensor(G, Xiprod, C);
Ep = tsrtenergyfromcurvature(G, Kp, C);

% Fusion energy release: (EA + EB - Ep)
Efus = (EA + EB - Ep);

end
```

<caption>Listing 12: TSRT fusion energetics / transparency</caption>
## Y.7 Odd-Even Staggering (OES) of $^{235}$U(n<sub>th</sub>,f) Yields

Dependencies: Listings 2 (p. 213), 3 (p. 214), 13 (p. 231).

Purpose: Process ENDF/B-VIII.0 Y(Z) and compute OES; generate Figure 3 (p. 24).

oesmakefig.m

```matlab
% FILE: oesmakefig.m
% PURPOSE: compute OES metrics and plot; save fig_OES.pdf
function oesmakefig(outFile)
% OESMAKEFIG Fit and plot Δ^(3) ln Y(Z) with a two-packet biased model.
% Robust plotting/saving with guarded exports (PDF via exportgraphics if
% available; EPS via print for maximum compatibility).
%
% Usage:
% oesmakefig()                 % writes fig_OES.pdf and fig_OES.eps
% oesmakefig('my_OES_plot.pdf') % writes my_OES_plot.pdf and my_OES_plot.eps

if margin < 1
    outFile = 'fig_OES.pdf';
end

% Resolve output paths
(outDir, outBase, ~] = fileparts(outFile);
if isempty(outDir), outDir = '.'; end
if isempty(outBase), outBase = 'fig_OES'; end
if ~exist(outDir,'dir'), mkdir(outDir); end
pdfPath = fullfile(outDir, [outBase, '.pdf]);
epsPath = fullfile(outDir, [outBase, '.eps']);

% --- Embedded ENDF/B-VIII.0 Y(Z) ----------------------------------------
Z = (34:62).';
Y = [0.000214 0.000339 0.000885 0.001787 0.004259 0.009037 0.017678 0.028505
0.038216 0.041846 0.039821 0.034138 0.027617 0.021681 0.016809 0.012959
0.010035 0.007773 0.006066 0.004803 0.003851 0.003137 0.002601 0.002192
0.001875 0.001624 0.001421 0.001253 0.001112]';
Y = Y / sum(Y);

% Δ^(3) ln Y
lnY = log(Y);
Zmid = Z(2:end-1);
Delta3 = 0.5*(lnY(3:end) - 2*lnY(2:end-1) + lnY(1:end-2));

% --- Model ----------------------------------------------------------------------
env =@(Zc,L,Zm) exp(-(Zm - Zc)/L).^2;
model =@(p,Zm)((-1).^(Zm + round(p(12)))) .* ( ...
    env(p(7),p(9),Zm).* (p(1)*besselj(0, p(8)*(Zm-p(7))) + p(2)*besselj(1, p(8)
    *(Zm-p(7))) + p(3)*besselj(2, p(8)*(Zm-p(7))) ) + ...
    env(p(7),p(11),Zm).* (p(4)*besselj(0, p(10)*(Zm-p(7))) + p(5)*besselj(1, p
    (10)*(Zm-p(7))) + p(6)*besselj(2, p(10)*(Zm-p(7))) ) ) ...
    + p(13)*exp(-(Zm - p(14))/p(15)).^2;
% p = [c0,c1,c2, d0,d1,d2, Zc, k1, L1, k2, L2, phi, B, Z0, L0]

% --- Coarse seed for [Zc,k1,L1,k2,L2,phi] + LS for coefficients --------
best = [];
best_err = inf;

solve ls =@(A,b) local_safe ls(A,b); % ridge/pinv fallback

for phi = 0:1
```
```R
for Zc = 46:0.5:50
  for k1 = 0.15:0.05:0.60
    for k2 = 0.25:0.07:0.90
      for L1 = 6:2:12
        for L2 = 10:2:20
          J1 = [besselj(0,k1*(Zmid-Zc)), besselj(1,k1*(Zmid-Zc))],
                besselj(2,k1*(Zmid-Zc))];
          J2 = [besselj(0,k2*(Zmid-Zc)), besselj(1,k2*(Zmid-Zc))],
                besselj(2,k2*(Zmid-Zc))];
          A = ((-1).^(Zmid+phi)) * [env(Zc,L1,Zmid).*J1, env(Zc,L2
, Zmid).*J2];

  c = solve ls(A, Delta3);
  fit = A * c;
  err = norm(Delta3 - fit);

  if err < best_err
    best = [c(:).', Zc, k1, L1, k2, L2, phi];
    best_err = err;
  end
end
end
end
end
end
end
```

% --- Initialize local bias near Z~42 and fit full parameter vector ---

```R
p0 = [best, -0.06, 42.0, 3.2];
lb = [-50,-50,-50,-50,-50, 40, 0.05, 3.0, 0.05, 5.0, -0.5, -0.20, 38.0, 1.0];
ub = [ 50, 50, 50, 50, 50, 50, 60, 1.50,30.0, 1.50,30.0, 1.5, 0.00, 46.0, 8.0];
```

```R
obj = @ (p) model(p, Zmid) - Delta3;
if exist('lsqnonlin','file')==2
  opts = optimoptions('lsqnonlin','Display','off');
  p = lsqnonlin(obj, p0, lb, ub, opts);
else
  % Fallback if Optimization Toolbox missing
  p = fminsearch(@ (x) norm(obj(x)), p0);
  p = max(min(p,ub),lb);
end
p(12) = round(p(12)); % discrete parity flip
```

% --- Plot (defensively) ---

```R
figH = figure('Visible','off','Color','w','Name','Odd-Even Staggering','Renderer',
'painters');
ax = axes('Parent',figH); hold(ax,'on');
```

% style
set(ax, 'FontName', 'Helvetica', 'FontSize', 11, 'LineWidth', 1.0, ...
'Box', 'on', 'XGrid', 'on', 'YGrid', 'on', 'GridAlpha', 0.20, 'MinorGridAlpha'
, 0.15);

plot(ax, Zmid, Delta3, 'ko', 'MarkerFaceColor', 'k', 'DisplayName', 'Experiment');
plot(ax, Zmid, model(p,Zmid), 'r-', 'LineWidth', 2, 'DisplayName', 'TSRT');

xlabel(ax,'Z');
ylabel(ax,'$\Delta^{(3)} \ln Y(Z)$','Interpreter','latex');
```R
title(ax, 'Odd-Even Staggering', 'FontWeight', 'bold');
legend(ax, 'Location', 'best');
drawnow;

% --- Export: PDF via exportgraphics (guarded), EPS via print -PDF-
try
  if ishghandle(ax)
    exportgraphics(ax, pdfPath, 'ContentType', 'vector');
  else
    error('axes handle invalid');
end

catch
  % robust fallback
  try
    set(figH, 'PaperPositionMode', 'auto');
    print(figH, '-dpdf', pdfPath);
  catch ME
    warning('[OES] PDF export failed: %s', ME.message);
end
end

% EPS (prefer print for compatibility)
try
  if ishghandle(figH)
    set(figH, 'PaperPositionMode', 'auto');
    print(figH, '-depsc2', epsPath);
  else
    error('figure handle invalid');
end

catch ME
  warning('[OES] EPS export failed: %s', ME.message);
end

% Cleanup
if ishghandle(figH), close(figH); end
end

% ----------------------------------------------------------------------
function c = local_safe ls(A, b)
%LOCAL_SAFE_LS Stable least-squares: normal equations w/ tiny ridge if needed.
ATA = A.'*A;
ATb = A.'*b;
r = rcond(ATA);
if ~isfinite(r) || r < 1e-12
  lam = 1e-12 * max(1, norm(ATA, 'fro'));
  c = (ATA + lam*eye(size(ATA))) \ ATb;
else
  c = ATA \ ATb;
end
if any(~isfinite(c)), c = pinv(A)*b; end
end
```

<caption>Listing 13: Compute/plot OES figure</caption>

*Figures produced correspond to Figure 3 (p. 24).*
## Y.8 E2 Gamma Calibration and Rates

This module calibrates the E2 $\gamma$ decay constant from a single, well-characterized anchor transition in $^{156}$Gd ($2_1^+ \to 0_{gs}^+$), using evaluated level data and internal conversion coefficients (ENSDF A=156 [130] and BrIcc [131]). The calibration routine writes C.gamma.E2 (C_E2, bZ, bA, provenance) to tsrtconstantsCalibrated.mat, and the evaluator uses these constants to compute partial E2 rates on demand.

Dependencies: Listings 2 (p. 213), 15 (p. 235), 14 (p. 234).

*Workflow.* First run the E2 calibration (Listing 15 (p. 235)) once to populate C.gamma.E2. Then execute the lifetime driver (Listing 16 (p. 238)) to regenerate Tables 2 (p. 11)–18 (p. 110) without any per-transition tuning.

tsrtgamma_rate.m

```matlab
% FILE: tsrtgamma_rate.m
% PURPOSE: compute E2 gamma rates using the 156Gd anchor in tsrtconstants.m;
%           optionally rederive anchor via TSRT_Calibrate_GammaE2_156Gd.m
% OUTPUTS: (e.g.) gammaRates.mat / tables as used in the manuscript
function lambda = tsrtgamma_rate(Z, A, Egamma_MeV, XL, Cgamma)
%TSRT_GAMMA_RATE Partial gamma-decay rate (s^{-1}) for a single multipolarity.
% Currently supports E2 via Cgamma.E2 with an E^5 law (TSRT-consistent).
% Safe: returns NaN on missing/invalid inputs or unsupported XL.

lambda = NaN;

if ~(isfinite(Egamma_MeV) && Egamma_MeV > 0)
    return;
end

if ~ischar(XL) && ~isstring(XL), return; end
xls = upper(strtrim(char(XL)));

switch xls
    case 'E2'
        if ~isfield(Cgamma,'E2') || ~all(isfield_s(Cgamma.E2, {'C_E2','bZ','bA'}))
            return;
        end
        C_E2 = Cgamma.E2.C_E2;
        bZ = Cgamma.E2.bZ;
        bA = Cgamma.E2.bA;
        lambda = C_E2 * (Egamma_MeV.^5) .* exp(bZ*Z + bA*A);
    otherwise
        return;
    end
end

function tf = isfield_s(S, keys)
tf = false(size(keys));
for i=1:numel(keys), tf(i) = isfield(S, keys{i}); end
end
```

<caption>Listing 14: E2 gamma-rate calculator</caption>

TSRT_Calibrate_GammaE2_156Gd.m (optional anchor re-derivation)

* Purpose: Recompute the global E2 constant from the $^{156}$Gd $2_1^+ \to 0^+$ anchor.
* **Anchor values (this paper):** $E_γ = 0.0889656$ MeV, $τ = 1.635×10^{-9}$ s, $BR= 1.0$, $α_{ICC}^{tot} = 0.163$, $δ = 0$, partial $γ$ half-life $≈ 1.318×10^{-9}$ s.

TSRT_Calibrate_GammaE2_156Gd.m

```matlab
% FILE: TSRT_Calibrate_GammaE2_156Gd.m
% PURPOSE: Calibrate the TSRT E2 gamma constant (C.gamma.E2.C_E2) on the anchor:
%          156Gd (Z=64, A=156), 2+1 -> 0+ (pure E2), Egamma=0.0889656 MeV,
%          tau(level)=1.635e-9 s, BR=1.0, alpha_ICC_total=0.163, delta=0.
% OUTPUTS: tsrtconstantsCalibrated.mat with C.gamma.E2.{C_E2,bZ,bA,anchor,...}
% NOTES: Anchor data from \cite{McCarten2008NDS} (NDS A=156) and BrIcc \cite{Kibedi2008BrIcc}.
```

function TSRT_Calibrate_GammaE2_156Gd()

```matlab
% TSRT_Calibrate_GammaE2_156Gd
% Calibrate the TSRT gamma (E2) constant using the anchor:
% 156Gd (Z=64, A=156), 2+1 -> 0+ (pure E2)
% Egamma = 0.0889656 MeV
% tau(level) = 1.635e-9 s
% BRGamma = 1.0 (100% branch)
% alpha_ICC_total = 0.163 (BrIcc v2.3, Frozen DHF, consistent with ENSDF)
% delta (M1/E2 mixing) = 0 (pure E2)
% This script:
% 1) computes the partial gamma rate for E2 from level lifetime & ICC
% 2) calibrates C.gamma.E2.C_E2 with an E^5 law (Weisskopf/collective scaling)
% 3) stores results in tsrtconstantsCalibrated.mat under C.gamma.E2
% 4) sanity-checks by reproducing the anchor t1/2
```

```
fprintf(' [gamma/E2] Calibrating on 156Gd (2+ -> 0+), pure E2\n');
```

```
% --------------------- Anchor data ---------------------
Z = 64;
A = 156;
```

```
Eg = 0.0889656;      % MeV (gamma decay energy)
tau = 1.635e-9;      % s (mean life of the 2+ level)
BR = 1.0;           % unitless (100% gamma branch from that level)
alphaICC = 0.163;    % unitless (total internal conversion coefficient)
delta = 0.0;         % M1/E2 mixing ratio (pure E2 => 0)
```

```
% --------------------- Compute partial gamma rate ---------------------
```

```
% Level total decay constant:
lambda_level = 1.0 / tau;      % s^-1
```

```
% Fraction that actually goes via gamma (not conversion):
% fGamma = BR / (1 + alpha_total)
f Gamma = BR / (1.0 + alphaICC);
```

```
% If a mixing ratio existed, one could project E2 fraction via delta,
% but for pure E2 (delta=0) the E2 fraction is 1.
f_E2 = 1.0;
```

```
lambda_E2 = lambda_level * fGamma * f_E2;      % s^-1

t12_E2 = log(2) / lambda_E2;      % s
```

```
fprintf(' [gamma/E2] Anchor partial t1/2 = %.6g s (lambda_E2=%.6g s^-1)\n', t12_E2, lambda_E2);
```

```
% --------------------- E2 model & calibration ---------------------
```

```
% TSRT-consistent, minimal E2 law (kept simple & non ad-hoc):
```
% lambda_E2 = C_E2 * Eg^5 * exp(bZ*Z + bA*A)
% For now we set bZ=bA=0 (no extra scaling beyond E^5 at calibration stage),
% which is the cleanest anchor-based normalization. One may later fit bZ, bA
% across a set of E2 transitions without changing this anchor.

bZ = 0.0;
bA = 0.0;

% Calibrate C_E2 from the anchor:
C_E2 = lambda_E2 / (Eg^5 * exp(bZ*Z + bA*A)); % s^-1 MeV^-5

fprintf(' [gamma/E2] Calibrated C_E2 = %.6g (units: s^-1 MeV^-5) 
 ', C_E2);
 fprintf(' [gamma/E2] bZ= %.6g, bA= %.6g (set to zero at anchor stage) 
 ', bZ, bA);

% ----------------- Load existing constants (if any) -----------------

C = struct();
if exist('tsrtconstants Calibrated.mat','file') == 2
S = load('tsrtconstants Calibrated.mat');
% Convert to a struct C regardless of content
fn = fieldnames(S);
for k=1:numel(fn), C.(fn{k}) = S.(fn{k}); end
fprintf(' [gamma/E2] Loaded existing tsrtconstants Calibrated.mat 
 ');
else
fprintf(' [gamma/E2] No existing constants file; starting a new one. 
 ');
end

% ----------------- Write gamma/E2 block -----------------

Corrigan
C.gamma.E2.C_E2 = C_E2;
C.gamma.E2.bZ = bZ;
C.gamma.E2.bA = bA;

% Keep full anchor provenance (for reproducibility & future fits)

C.gamma.E2.anchor.Z = Z;
C.gamma.E2.anchor.A = A;
C.gamma.E2.anchor.isotope = '{~{156}\\mathrm{Gd}}';
C.gamma.E2.anchor.Jpi = '2^+_1 \\to 0^+_{-gs}';
C.gamma.E2.anchor.Egamma_MeV = Eg;
C.gamma.E2.anchor.tau_s = tau;
C.gamma.E2.anchor.BR = BR;
C.gamma.E2.anchor.alpha_total = alphaICC;
C.gamma.E2.anchor delta = delta;
C.gamma.E2.anchor.t12_partial_s = t12_E2;
C.gamma.E2.anchor lambdahalf_s inv = lambda inv;

% References

C.gamma.E2 refs.ensdf.citation = 'McCarten et al., Nucl. Data Sheets 109, 1 (2008)';
C.gamma.E2 refs.ensdf.doi = '10.1016/j.nds.2008.01.001';
C.gamma.E2 refs.bricc.citation = 'Kibédi et al., Nucl. Instrum. Methods A 589 (2008) 202-229';
C.gamma.E2 refs.bricc.doi = '10.1016/j.nima.2008.02.051';

% ----------------- Save constants -----------------

save('tsrtconstants Calibrated.mat','-struct','C');
fprintf(' [gamma/E2] Saved updated constants to tsrtconstants Calibrated.mat 
 ');

% ----------------- Self-check: re-predict anchor -----------------

lambda_check = tsrtgamma_E2_rate(Z, A, Eg, C.gamma.E2);
t12_check = log(2) / lambda_check;
```c
rel_err = abs(t12_check - t12_E2) / t12_E2;

fprintf(' [gamma/E2] Self-check: predicted t1/2(anchor) = %.6g s (rel.err=%.3e)\n',
t12_check, rel_err);
if ~ (isfinite(rel_err) && rel_err < 1e-12)
    warning(' [gamma/E2] Anchor reproduction relative error is %.3e (expected ~0).'
           , rel_err);
end

fprintf(' [gamma/E2] Calibration complete.\n');
end

% ======== Local helper ========
function lambda = tsrt
    lambda_E2_rate(Z, A, Eg_MeV, E2)
        % Minimal E2 rate model used above:
        % lambda_E2 = C_E2 * Eg^5 * exp(bZ*Z + bA*A)
        if ~ (isfinite(Eg_MeV) && Eg_MeV > 0)
            lambda = NaN; return;
        end
        C_E2 = E2.C_E2; bZ = E2.bZ; bA = E2.bA;
        lambda = C_E2 * (Eg_MeV.^5) .* exp(bZ*Z + bA*A);
    end
end
```

<caption>Listing 15: Calibrate E2 γ constant on <sup>156</sup>Gd (anchor)</caption>

*Anchor provenance.* <sup>156</sup>Gd level lifetime and branching from ENSDF A=156 [130]; total ICC ($\alpha_{tot}$ = 0.163) from BrIcc [131]. The calibration uses an $E^5$ scaling at anchor stage (no extra Z/A factors).
## Y.9 Lifetime Pipeline and Benchmarks

**Dependencies:** Listing 15 (p. 235) (E2 anchor calibration, writes `C.gamma.E2` to `tsrtconstants Calibrated.mat`), Listing 14 (p. 234) (E2 rate evaluator used by the lifetime code), and Listing 18 (p. 239) (optional forbiddenness multiplier).

`RunLifetime_All.m`

```matlab
% FILE: RunLifetime_All.m
% PURPOSE: Compute TSRT lifetimes (alpha/gamma/beta branches), compare to ENSDF,
%          and emit paper tables (e.g., tab_lifetime_bench.tex,
%          tab_lifetime_bench_highres.tex).
% INPUTS: tsrtconstants.m (global constants), tsrtconstants Calibrated.mat (gamma/E2
%          anchor),
%          evaluated datasets loaded by the branch scripts.
% OUTPUTS: paper-assets-tex/tab_lifetime_bench.tex (and related .tex/.csv as described
%          in text)
% NOTES: Reads C.gamma.E2.* as produced by TSRT_Calibrate_GammaE2_156Gd.m.
function RunLifetime_All()
% Run α, β, sf calibrations in sequence, then emit the combined table.
clear functions; rehash toolboxcache;

RunLifetimeCalculations; % alpha only, saves/merges constants
RunLifetimeCalculations beta; % beta only, merges into the same file
RunLifetimeCalculations sf; % sf only, merges into the same file

tsrt_emit_tables('tsrt_benchmark_diag.csv', ...
                  'tab_lifetime_bench_highres.tex', ...
                  'tsrt_benchmark_diag_FILLED.csv');

fprintf('[All] Emission complete.\n');

clear C;
tsrt_emit_tables('tsrt_benchmark_diag.csv', 'tab_lifetime_bench_highres.tex',
                  'tsrt_benchmark_diag_FILLED.csv');
end
```

<caption>Listing 16: Lifetime benchmarking driver (canonical)</caption>

`tsrtbetamatrixelement.m`

```matlab
function Mbeta = tsrtbetamatrixelement(G, C, Xi_parent, Xi_daughter)
%TSRTBETAMATRIXELEMENT Signed TSRT β matrix element from ℓ=1 multipoles.
% Mbeta = sum_{m=-1}^{1} (c_{1m}^{(p)} - c_{1m}^{(d)}), with unit weights.
% Xi_parent / Xi_daughter: curvature proxies on the production grid G.

% Project real Y_lm multipoles (use the same routine as OES)
c_p = tsrt_project_real_multipoles(G, Xi_parent, 1); % returns [c_{1,-1}, c_{1,0}, c_{1,1}, c_{-1,0}, c_{-1,1}]
c_d = tsrt_project_real_multipoles(G, Xi_daughter, 1);

% Unit weights (w_m = 1)
Mbeta = sum(c_p - c_d);
end
```

<caption>Listing 17: TSRT β-matrix element from signed $ℓ=1$ multipoles</caption>

`tsrt_forbiddenness_hindrance.m` (optional lifetime helper)
* **Purpose:** Applies a per-rank multiplicative hindrance factor $10^{s^K}$ with slope $s = C.HF\_slope\_per\_rank$.

* **Inputs:** Forbiddenness rank $K$ (scalar or array), constants $C$ from `tsrtconstants.m`.

* **Outputs:** HF (same shape as $K$).

`tsrt_forbiddenness_hindrance.m`

```matlab
function HF = tsrt_forbiddenness_hindrance(K, C)
    %TSRT_FORBIDDENNESS_HINDRANCE Hindrance factor from forbiddenness rank K.
    % Physics-based: no proxy tuning. If TSRT prescribes an empirical
    % per-rank factor, set C.HF_slope_per_rank accordingly; otherwise 0.
    %
    % Inputs:
    % K forbiddenness rank (scalar or array, expected integer >= 0)
    % C constants struct from tsrtconstants() (optional)
    %
    % Output:
    % HF multiplicative hindrance factor on partial half-life (same size as K)

    if margin < 2 || ~isstruct(C)
        C = tsrtconstants();
    end
    if ~isfield(C, 'HF_slope_per_rank') || isempty(C.HF_slope_per_rank)
        C.HF_slope_per_rank = 0; % default: no extra hindrance
    end

    % Preserve shape
    K_in = K;
    K = K(:);
    HF = ones(size(K));

    if C.HF_slope_per_rank ~= 0
        m = (K > 0);
        HF(m) = 10.^((C.HF_slope_per_rank.*K(m)));
    end

    HF = reshape(HF, size(K_in));
end
```

<caption>Listing 18: `tsrt_forbiddenness_hindrance.m`</caption>
## Y.10 Diagnostics and Validation Roadmap

This subsection provides a consolidated map of all diagnostic and validation tools in the Appendix, linking each numerical check to its corresponding MATLAB listing.

* **Grid generation, Gauss-Legendre weights, and volume-integral conventions** — Appendix F.3 (foundational for all TSRT numerical geometry).

* **Convergence and robustness tests** — Appendix Y.3 (grid refinement, domain-size checks, and curvature stability).

* **Stability maps and per-Z diagnostic lines** — Appendix Y.11 (sanity checks used throughout the main text).

* **Optional correlation and pathway diagnostics for fission** — Appendix Y.16.
## Y.11 Stability Maps and Diagnostic Lines

Dependencies: Listings 2 (p. 213), 3 (p. 214).

Purpose: Construct the strict stability map for $\Delta K^2$ and generate the per-$Z$ stability-line figures used as diagnostics in the paper.

Used in the paper: Diagnostics referenced in the text; figures `tsrt_stability_lines.pdf` and `tsrt_stability_lines_heavy.pdf`.

Inputs: `tsrtconstants.mat`; `tsrt_stabilitymap_data.mat`.

Outputs: `tsrt_stability_map.mat` (map artifact, if produced); `tsrt_stability_lines.pdf`; `tsrt_stability_lines_heavy.pdf`.

Upstream: 2 (p. 213) (global constants).

Downstream: Standalone diagnostic figures.

`stabilitymap_deltaK2_strict.m`

```matlab
% FILE: stabilitymap_deltaK2_strict.m
% PURPOSE: Build strict (calculated-only) stability map for ΔK^2 on (Z,N)
% INPUTS: tsrtconstants.mat, tsrt_stabilitymap_data.mat
% OUTPUTS: tsrt_stability_map.mat (optional: used to derive per-Z chains)
function stabilitymap_deltaK2_strict(tsrtCsv, outBase)
% STABILITYMAP DELTAK2_STRICT ΔK^2-based stability landscape (robust, no toolboxes).
% Usage:
% stabilitymap_deltaK2_strict('tsrt_bindingscan_cln.csv', 'tsrt_stabilitymap')

if margin < 1 || isempty(tsrtCsv), tsrtCsv = 'tsrt_bindingscan_cln.csv'; end
if margin < 2 || isempty(outBase), outBase = 'tsrt_stabilitymap'; end

% --- Load constants (robust: either S.C or S fields) ---
S = load('tsrtconstants_calibrated.mat');
if isfield(S,'C') && isstruct(S.C), C = S.C; else, C = S; end
C = ensure_field(C, 'r_max', 1.0);
C = ensure_field(C, 'N_r', 128);
C = ensure_field(C, 'N_theta', 64);
C = ensure_field(C, 'coherence_eta', 0.05);
% Optional nucleus-geometry parameters (used if present)
C = ensure_field(C, 'r0', 1.20e-15);
C = ensure_field(C, 'beta_skin', 0.10);
C = ensure_field(C, 'kappa phi', pi/12);

% --- Grid once (keeps indices consistent) ---
G = tsrtmakegrid(C);

% --- Read TSRT CSV and list nuclides (Z,A) ---
T = readtable(tsrtCsv, 'FileType', 'text');
T = normalize_headers(T);
Zc = find_col(T, "z"); Ac = find_col(T, "a");
assert(~isempty(Zc) && ~isempty(Ac), 'TSRT CSV must contain Z and A.');

ZA = unique([T.(Zc), T.(Ac)], 'rows'); % [Z A]
Zlist = ZA(:,1);
Nlist = ZA(:,2) - ZA(:,1);

Zmin = min(Zlist); Zmax = max(Zlist);
Nmin = min(Nlist); Nmax = max(Nlist);

Zvals = Zmin:Zmax; % rows
Nvals = Nmin:Nmax; % columns
[NN, ZN] = meshgrid(Nvals, Zvals);
```
Fbase = NaN(size(ZN),'like',1.0); % map scalar: base slope |d(ΔK^2)/dτ|_{τ*}

% --- Evaluate nucleus-by-nucleus ---

tau_vec = linspace(0,1,256).';
n = size(ZA,1);
fprintf('ΔK^2 base map: evaluating %d nuclides...
', n);
t0 = tic;

for k = 1:n

Z = ZA(k,1);
A = ZA(k,2);
N = A - Z;

if mod(k,200)==0
    fprintf('%5d/%d (%.1f%%) elapsed %.1fs
', k, n, 100*k/n, toc(t0));
end

i = Z - Zmin + 1; % row index (Z)

j = N - Nmin + 1; % col index (N)

if i<1 || j<1 || i>size(Fbase,1) || j>size(Fbase,2), continue; end

try

D = tsrtdeltaK2longtau(Z, A, C, G, tau_vec); % ΔK^2(τ)

if ~isnumeric(D) || all(~isfinite(D)), continue; end

% Smooth a touch and take local slope at τ* (min of smoothed ΔK^2)

Dsm = movmean(D,5);

[~, kmin] = min(Dsm);

w = 3;

i1 = max(1, kmin-w); i2 = min(numel(tau_vec), kmin+w);

tv = tau_vec(i1:i2); dv = Dsm(i1:i2);

good = isfinite.tv) & isfinite(dv);

if nnz(good) >= 3

p = polyfit.tv(good), dv(good), 1);

slope = abs(p(1));

if isfinite(slope) && slope > 0

Fbase(i,j) = slope;

end

end

catch

end

end

cov = nnz(isfinite(Fbase));

fprintf('Coverage: %d / %d (%.1f%%)
', cov, numel(Fbase), 100*cov/numel(Fbase));

if cov == 0

error('No base entries computed. Verify curvature pipeline.');

end

% --- Save arrays for reproducibility and QA ---

save([outBase '_data.mat'], 'ZN','NN','Fbase','Zmin','Zmax','Nmin','Nmax','-v7');

fprintf('Saved map arrays to %s
', [outBase '_data.mat']));

% --- Plot (percentile-stretched) ---

f = figure('Color','w','Position', [120 120 980 560], 'Name','TSRT stability map');

ax = axes('Parent',f); hold(ax,'on'); box(ax,'on'); grid(ax,'on');
```matlab
mask = isfinite(Fbase);
v = Fbase mask);
lo = prctile(v, 1);
hi = prctile(v, 97);
if ~ (isfinite (lo) && isfinite (hi) && hi > lo)
    lo = min(v); hi = max(v);
end

imagesc(ax, [min(Nvals) max(Nvals)], [min(Zvals) max(Zvals)], Fbase);
set(ax, 'YDir', 'normal');
try, colormap(ax, turbo); catch, colormap(ax, parula); end
caxis(ax, [lo hi]);
cb = colorbar(ax); cb.Label.String = ' |d(\DeltaK^2)/d\tau| at \tau^* (arb. units) ';
xlabel(ax, 'Neutron number N'); ylabel(ax, 'Proton number Z');
title(ax, '\DeltaK^2-based TSRT stability map (percentile-scaled)');

print(f, [outBase '.png'], '-dpng', '-r300');
set(f, 'PaperPositionMode', 'auto');
print(f, [outBase '.pdf'], '-dpdf');
fprintf('Saved: %s.png and %s.pdf\n', outBase, outBase);
end

% ---------- helpers ----------
function T = normalize_headers(T)
    raw = T.Properties.VariableNames;
    norm = lower(strrep/raw, '--', '_'));
    norm = strrep/norm, '.', '_');
    T.Properties.VariableNames = matlab.lang.makeUniqueStrings/norm);
end

function cname = find_col(T, candidates)
    cname = '';
    cols = lower(T.Properties.VariableNames);
    for c = string(candidates)
        m =strcmp(cols, char(c));
        if any(m), cname = T.Properties.VariableNames{find(m,1,'first')}; return; end
    end
end

function C = ensure_field(C, name, val)
    if ~isfield(C, name) || isempty(C.(name)), C.(name) = val; end
end
```

<caption>Listing 19: Strict stability-map construction for $\Delta K^2$</caption>

stability_lines_makefig.m

```matlab
% FILE: stability_lines_makefig.m
% PURPOSE: Create the published stability-line figures (global + heavy region)
% INPUTS: tsrtconstants.mat, tsrt_stabilitymap_data.mat
% OUTPUTS: tsrt_stability_lines.pdf, tsrt_stability_lines_heavy.pdf
% stability_lines_makefig.m
% Create the published stability-line figures (global + heavy region)
% Load constants (creates/loads tsrtconstants.mat as needed)
C = tsrtconstants();

% Load stability-map input data (required by tsrt_stability_lines.m)
```
```R
if !exist('tsrt_stabilitymap_data.mat', 'file')
    error('Missing tsrt_stabilitymap_data.mat (see Appendix for data resources).');
end
S = load('tsrt_stabilitymap_data.mat');

% Produce the global stability-lines figure
% (tsrt_stability_lines should save to 'tsrt_stability_lines.pdf' when called)
tsrt_stability_lines(S, C, struct('region', 'global', 'out', 'tsrt_stability_lines.pdf'));

% Produce the heavy-region stability-lines figure
tsrt_stability_lines(S, C, struct('region', 'heavy', 'out', 'tsrt_stability_lines_heavy.pdf'));

fprintf('Wrote: tsrt_stability_lines.pdf and tsrt_stability_lines_heavy.pdf\n');
```

<caption>Listing 20: Per-Z stability lines (wrapper script)</caption>

```
tsrt_stability_lines.m

% FILE: tsrt_stability_lines.m
% PURPOSE: Plot per-Z chains of |d(ΔK^2)/dτ| vs N; used by the wrapper script
% INPUTS: tsrtconstants.mat, tsrt_stabilitymap_data.mat
% OUTPUTS: (figures saved by caller)
function tsrt_stability_lines.dat, out_fig, Z_select, maxZcurves)
% TSRT_STABILITY_LINES
% Plot calculated-only stability scale vs neutron number:
% x-axis: N, y-axis: |d(ΔK^2)/dτ|/_{τ*}
% One curve per Z (proton number); connects computed points for that Z.
%
% Usage:
% tsrt_stability_lines('tsrt_stabilitymap_data.mat', 'tsrt_stability_lines')
% tsrt_stability_lines('tsrt_stabilitymap_data.mat', 'tsrt_stability_lines', [8 20 28 50 82 100])
% tsrt_stability_lines('tsrt_stabilitymap_data.mat', 'tsrt_stability_lines', [], 11)
%
% Inputs:
% dat (char) MAT with ZN/NN/F [default: 'tsrt_stabilitymap_data.mat']
% out_fig (char) output basename [default: 'tsrt_stability_lines']
% Z_select (vec) Z-values to plot [default: auto-pick richest Z-chains]
% maxZcurves (int) cap on #Z curves [default: 11]
%
% Saves: out_fig.pdf and out_fig.eps (vector, white background, black labels)

if margin<1 || isempty.dat), dat = 'tsrt_stabilitymap_data.mat'; end
if margin<2 || isempty(out_fig), out_fig = 'tsrt_stability_lines'; end
if margin<4 || isempty(maxZcurves), maxZcurves = 11; end

% --- Load arrays robustly ---
S = load.dat);
ZN = pick(S, ["ZN","Zgrid","ZZ","Z"])';
NN = pick(S, ["NN","Ngrid","N"])';
F = pick(S, ["F","Fbase","Fmap"])';
assert(~isempty(ZN) && ~isempty(NN) && ~isempty(F), 'Need ZN/NN/F in MAT.');
assert(isequal(size(ZN), size(NN), size(F)), 'ZN/NN/F must be same size.');

mask = isfinite(F) & F>-Inf & F<Inf;
Zv = ZN mask;
Nv = NN mask;
```
```R
Fv = F mask;

% --- Choose which Z chains to plot ---
Zu = unique(Zv);
if length >= 3 && ~isempty(Z_select)
    Z_use = intersect unique(Z_select(:).'), Zu(:).');
    assert(~isempty(Z_use), 'None of the requested Z in Z_select are present.');
else
    % auto-pick the richest Z-chains (most computed points)
    counts = arrayfun(@(z) sum(Zv==z), Zu);
    [~,ord] = sort(counts, 'descend');
    Z_use = Zu(ord(1:min(maxZcurves,numel(Zu))));
    Z_use = sort(Z_use); % nicer legend order
end
```

```R
% --- Build N→F chains per Z (sorted by N, calculated points only) ---
chains = struct('Z', {}, 'N', {}, 'F', {});
for k = 1:lisel(Z_use)
    z = Z_use(k);
    idx = (Zv==z);
    N_k = Nv(idx);
    F_k = Fv(idx);
    % sort by N ascending
    [N_k, order] = sort(N_k);
    F_k = F_k(order);
    % remove duplicate N (keep median if duplicates exist)
    [Nuniq, ~, ic] = unique(N_k);
    if numel(Nuniq) < numel(N_k)
        Fmed = accumarray(ic, F_k, [], @median);
        N_k = Nuniq; F_k = Fmed;
    end
    chains(k).Z = z;
    chains(k).N = N_k;
    chains(k).F = F_k;
end
```

```R
fig = figure('Color','w','Position', [100 100 900 560], ...
    'Name','TSRT stability lines: |d(ΔK^2)/dτ| vs N');
ax = axes('Parent',fig); hold(ax,'on');
set(ax, 'FontName', 'Helvetica', 'FontSize', 11, ' LineWidth', 1.0, ...
    'Box','on', 'XGrid','on', 'YGrid','on', ...
    'GridAlpha',0.20, 'MinorGridAlpha',0.15, ...
    'Color','w', 'XColor','k', 'YColor','k');
```

```R
xlabel(ax,'Neutron number N'); ax.XLabel.Color = 'k';
ylabel(ax,'|d(\Delta K^2)/d\tau| at \tau\tau^* (arb. units)'); ax.YLabel.Color = 'k';
t = title(ax,'TSRT stability (calculated-only): per-Z chains vs N','FontWeight','bold');
t.Color = 'k';
```

```R
% Color/linestyle cycles
baseCol = lines(max(12, numel(Z_use))); % consistent palette
styles = {'--', '-.', ':', ':'}; % cycle styles if many Z
```

```R
% Plot each chain: lines + markers at calculated points
for k = 1:lisel(chains)
    Zk = chains(k).Z; Nk = chains(k).N; Fk = chains(k).F;
```
```matlab
c = baseCol(mod(k-1,size(baseCol,1))+1,:);
ls = styles{ mod(k-1 numel(styles)) + 1 };

% line through the computed points for this Z:
plot(ax, Nk, Fk, 'LineStyle', ls, 'answid', 1.6, 'Color', c, ...
'屢', sprintf('Z = %d', Zk));

% markers at actual computed points:
plot(ax, Nk, Fk, 'o', 'MarkerSize', 4.5, ...
'MarkerFaceColor', c, 'MarkerEdgeColor', [0 0 0], ...
'LineStyle', 'none', 'HandleVisibility', 'off');
end

% Axes limits: tight on computed domain
Nmin = min(Nv); Nmax = max(Nv);
ylim_auto = [min(Fv), max(Fv)];
pad = 0.05*diff(ylim_auto);
xlim(ax, [Nmin-1, Nmax+1]);
ylim(ax, [ylim_auto(1)-pad, ylim_auto(2)+pad]);

if numel(Z_use) <= 12
lg = legend(ax, 'Location','eastoutside','Box','off');
else
lg = legend(ax, 'Location','best','Box','off');
end
if ~isempty(lg), set(lg,'TextColor','k'); end

% --- Save PDF + EPS with white background ---
[outDir, outBase, ~] = fileparts(out_fig);
if isempty(outDir), outDir = '.'; end
if ~exist(outDir,'dir'), mkdir(outDir); end
pdfPath = fullfile(outDir, [outBase '.pdf']);
epsPath = fullfile(outDir, [outBase '.eps']);

set(fig, 'Color','w', 'InvertHardcopy', 'off');
try
exportgraphics(fig, pdfPath, 'ContentType','vector', 'BackgroundColor','white');
catch
set(fig,'PaperPositionMode','auto');
print(fig, pdfPath, '-dpdf', '-painters');
end
try
print(fig, epsPath, '-depsc', '-painters');
catch MEeps
warning('EPS export failed: %s', MEeps.message);
end
fprintf('Wrote: %s and %s
', pdfPath, epsPath);
end

% ----------- helper ------------
function out = pick(S, names)
out = [];
for n = names
    if isfield(S, char(n)), out = S.(char(n)); return; end
end
end
```
Listing 21: Core plotting routine for stability lines

*Figures produced correspond to Figure 4 (p. 42) and Figure 5 (p. 43).*
## Y.12 Numerical Helpers (Quadrature)

Dependencies: none; utilities used by multiple modules.

Purpose: Shared numerical helpers (quadrature) used by curvature/action routines.

gaussleg.m

% FILE: gaussleg.m
% PURPOSE: Gauss-Legendre nodes/weights on [a,b] for N>=1.

% USAGE:

% [x,w] = gaussleg(N) % on [-1,1]

% [x,w] = gaussleg(N,a,b) % on [a,b]

function [x,w] = gaussleg(N,a,b)
if nargin < 2, a = -1; end
if nargin < 3, b = 1; end
if N < 1 || N ~= floor(N)
    error('N must be a positive integer.’);
end

% Initial guesses (Chebyshev-like)
k = (1:N).';

x0 = cos(pi*(k-0.25)/(N+0.5));

% Newton iterations on P_N

x = x0;

for it = 1:20
    [P, dP] = legpoly_and_deriv(N, x);
    dx = -P./dP;
    x = x + dx;
    if max(abs(dx)) < 1e-14, break; end
end

% Weights on [-1,1]
w = 2 ./ ((1 - x.^2) .* (dP.^2));

% Affine map to [a,b]
xm = (b+a)/2; xr = (b-a)/2;

x = xm + xr*x;

w = xr * w;

% Ensure ascending order

[x,perm] = sort(x,'ascend');

w = w(perm);

end

function [P, dP] = legpoly_and_deriv(N, x)
% Evaluate Legendre P_N(x) and P'_{N}(x) via stable recurrence.

x = x(:);

Pm1 = ones(size(x));

if N == 0
    P = Pm1; dP = zeros(size(x)); return;
end

P0 = x; % P1
if N == 1

P = P0;
dP = Pm1; % derivative of x is 1
return;

end
for n = 1:N-1

% Recurrence: (n+1)P_{n} = (2n+1)xP_{n-1} - nP_{n-2}

P1 = ((2*n+1).*x.*P0 - n.*Pm1)/(n+1);

Pm1 = P0; P0 = P1;

end

P = P0;

% Derivative via: P'_{N}(x) = N/(1-x^2) * (P_{N-1}(x) - x P_{N}(x))
```matlab
% Need P_{N-1}: recompute once efficiently
Pm1 = ones(size(x));
PO = x;
if N == 1
    PNm1 = Pm1;
else
    for n = 1:N-2
        P1 = ((2*n+1).*x.*PO - n*Pm1)/(n+1);
        Pm1 = PO; PO = P1;
    end
    PNm1 = PO;
end

dP = (N * (PNm1 - x.*P)) ./ max(1e-16, (1 - x.^2));
end
```

<caption>Listing 22: Gauss-Legendre nodes/weights</caption>

simpsonweights.m

```matlab
% FILE: simpsonweights.m
% PURPOSE: Return Simpson-rule weights for N points on [a,b] (N odd).
function w = simpsonweights(N)
% Generate normalized Simpson weights for N nodes.
% Requires N-1 to be even.
if mod(N-1,2~=0
    error('Simpson requires even number of intervals');
end
w = ones(1,N);
w(2:2:end-1) = 4;
w(3:2:end-2) = 2;
w = w / sum(w); % normalized weights; re-scale externally
end
```

<caption>Listing 23: Simpson weights</caption>
## Y.13 Core Field Generators (support)

Dependencies: Listings 2 (p. 213), 3 (p. 214).

Purpose: Generate TSRT nucleon and composite trembling fields used by curvature and energy integrals.

tsrtnucleonfield.m

```matlab
function Xi = tsrtnucleonfield(G, arg2, arg3, arg4)
% TSR T nucleon trembling field (time-averaged, spherical baseline).
%
% Supported call signatures:
%   Xi = tsrtnucleonfield(G, params)
%     where params is a struct with optional fields:
%       .type : 'proton'/'neutron' (default 'proton')
%       .r0   : core scale [m] (defaults by type)
%       .A    : amplitude (scalar) (default 1.0)
%       .phi  : phase [rad] (default 0; multiplies amplitude by cos(phi))
%       .norm : 'none'/'L2' (default 'none'; 'L2' enforces sum(|Xi|~2 dV)=1)
%
% Xi = tsrtnucleonfield(G, C, type, t)      % legacy
%     where C and t are ignored here (time-averaged profile),
%     and 'type' is as above. Defaults are applied for others.
%
% Output:
%   Xi.profile : Nr x Ntheta x Nphi real array on (G.r, G theta, G.phi)
%   Xi.meta   : struct with provenance and resolved parameters
%
% Notes:
%   - The profile is a spherical Gaussian envelope exp(-(r/r0)~2), time-averaged.
%   - If params.norm == 'L2', the profile is normalized with the *grid weights*
%   - provided by tsrtnmakegrid.m: sum( raw.~2 * r~2 * w_r * w_theta * w phi ) == 1.
%   - G.w_theta is assumed to integrate sin theta dtheta (as constructed in
tsrtnmakegrid);
%   - do NOT multiply by sin theta again during normalization or later integrals.

% ---------- Dispatch / parse inputs ----------
if margin == 2 && isstruct(arg2)
    params = arg2;
else if margin >= 4
    % legacy: (G, C, type, t) -> build params from 'type'
    params = struct('type', arg3);
else
    error('tsrtnucleonfield: invalid call signature.');
end

% ---------- Defaults (robust, minimal assumptions) ----------
if ~isfield(params, 'type') || isempty(params.type), params.type = 'proton'; end
switch lower(params.type)
    case 'proton', r0_default = 0.85e-15;
    case 'neutron', r0_default = 0.90e-15;
    otherwise,      r0_default = 0.88e-15; % safe fallback
end
if ~isfield(params, 'r0') || isempty(params.r0), params.r0 = r0_default; end
if ~isfield(params, 'A') || isempty(params.A), params.A = 1.0;          end
if ~isfield(params, 'phi') || isempty(params), params = 0.0;    end
if ~isfield(params, 'norm') || isempty(params.norm), params.norm = 'none'; end

% ---------- Build spherical, time-averaged envelope ----------
```
```
[RR, TT, PP] = ndgrid(G.r, G theta, G);
% sizes: Nr x Nth x Nph
raw = params.A * cos(params.phi) * exp(-(RR./params.r0).^2);

% ---------- Optional deterministic L2 normalization ----------
if strncmp(params.norm, 'L2')
    % Use grid weights that integrate dr, sin theta dtheta, dphi respectively.
    % Volume element = r^2 * (w_r) * (w_theta) * (w phi).
    if ~isfield(G, 'w_r') || ~isfield(G, 'w theta') || ~isfield(G, 'w phi')
        error('tsrtnucleonfield: L2 normalization requires G.w_r, G.w_theta, G.w phi.');
    end
    Nr = numel(G.r); Nth ={G theta}; Nph ={G phi};
    wr = reshape(G.w_r, [Nr, 1, 1]);
    wth = reshape(G.w theta, [1, Nth, 1]);
    wph = reshape(G.w phi, [1, 1, Nph]);
    r2 = reshape(G.r.^2, [Nr, 1, 1]);
    dVw = r2 * wr * wth * wph; % measure weights (no extra sin theta))
    norm2 = sum( (raw.^2) * dVw, 'all' );
    raw = raw / sqrt(max(norm2, realmin)); % stable normalization
end

% ---------- Package output ----------
Xi.profile = raw;
Xi.meta = struct('...
    'type', lower(params.type), ...
    'r0',   params.r0, ...
    'A',    params.A, ...
    'phi',  params.phi, ...
    'norm', lower(string params.norm)) ...
');
end
```

<caption>Listing 24: Single nucleon field (unshifted)</caption>

```
tsrtnucleonfieldshifted.m

% FILE: tsrtnucleonfieldshifted.m
% PURPOSE: Evaluate the single-nucleon trembling field at a spatial shift R.
% INPUTS: G (grid), params (struct as in tsrtnucleonfield), R (3x1 meters)
% OUTPUTS: Xi.profile (3D array)

function Xi = tsrtnucleonfieldshifted(G, params, R)
    % Build unshifted field
    Xi0 = tsrtnucleonfield(G, params);

    % Coordinates of grid points (spherical -> Cartesian), then shift
    [RR,TT,PP] = ndgrid(G.r, G theta, G::ph);
    [x,y,z] = sph2cart(PP, pi/2 - TT, RR);
    x = x - R(1); y = y - R(2); z = z - R(3);
    RRs = sqrt(x.^2 + y.^2 + z.^2);

    % Radial average of unshifted profile, then radial interpolation
    radial = squeeze(mean(mean(Xi0.profile,3),2)); % <profile>(r) over angles
    Xi.profile = interp1(G.r, radial, RRs, 'linear', 0);

    Xi.meta = params;
    Xi.meta.center = R(:).';
end
```
<caption>Listing 25: Single nucleon field (shifted)</caption>

```
tsrtcompositefield.m

function Xi = tsrtcompositefield(G, nuc, t)
% Assemble a composite trembling field by superposition.
% Inputs:
% G : grid struct (tsrtmakegrid)
% nuc : struct with fields:
% .types {1xA} cell of 'proton'/'neutron'
% .R [3xA] positions (m)
% optional (scalar or 1xA):
% .r0, .Aamp, .anis, .phi
% t : time (unused; we use time-averaged baseline)
% Output:
% Xi : struct with .profile (3D array) for curvature integrals

A = numel(nuc.types);
XiTot = zeros(numel(G.r),numel(G theta),numel(G.phi));

% Defaults
def_r0_p = 0.85e-15;
def_r0_n = 0.90e-15;
def_Aamp = 1.0;
def_anis = 0.0;

if nargin < 3 || isempty(t), t = 0.0; end

for i = 1:A
    % r0
    if isfield(nuc,'r0')
        if isscalar(nuc.r0), r0_i = nuc.r0; else, r0_i = nuc.r0(i); end
    else
        ifstrcmpi(nuc.types{i}, 'proton'), r0_i = def_r0_p; else, r0_i = def_r0_n; end
    end
    % Aamp
    if isfield(nuc,'Aamp')
        if isscalar(nuc.Aamp), A_i = nuc.Aamp; else, A_i = nuc.Aamp(i); end
    else
        A_i = def_Aamp;
    end
    % anis
    if isfield(nuc,'anis')
        if isscalar(nuc.anis), anis_i = nuc.anis; else, anis_i = nuc.anis(i); end
    else
        anis_i = def_anis;
    end
    % phi
    if isfield(nuc,'phi')
        if isscalar(nuc), phi_i = nuc; else, phi_i = nuc District; end
    else
        phi_i = 0.0;
    end

params_i = struct('type', nuc.types{i}, ...
    'A', A_i, ...
    'r0', r0_i, ...
    'phi', phi_i, ...
```
```c
'anis', anis_i);

Xi_i = tsrtnucleonfieldshifted(G, params_i, nuc.R(:,i));
XiTot = XiTot + Xi_i.profile;
end

Xi = struct('profile', XiTot);
end
```

<caption>Listing 26: Composite nuclear field builder</caption>
Y.14 CSCT / Near-Field Electromagnetism Tables (optional branch)

csctmaketables.m

function csctmaketables()

%CSCTMAKETABLES Produce LaTeX tables for CSCT near-field EM shifts.

% Writes:

% paper-assets-tex/tab_csctgamma.tex

% paper-assets-tex/tab_csctbeta.tex

%

% Model: Appendix~\ref{app:csct}, Equations~(\ref{eq:hem-model})--(\ref{eq:deltaEM-model})

% Anchors:

% - photon coupling η_γ fixed by 57Fe 14.4 keV line: δ = 2.00 eV

% - charged coupling η_q fixed by a 300 eV conversion-e- at 300 keV (Ba-like Z)

% Constants (keep in sync with Appendix text):

alpha_EM = 1/137; % fine-structure (TSRT EM sector): α_EM^TSRT

kappa_NF = 3; % L_NF = kappa_NF * R_eff

ZeffFrac = 0.7; % Z_eff/Z (angularly averaged for photons)

c = 2.99792458e8; % not used explicitly in closed-form here

paperassets = fullfile(pwd, 'paper-assets-tex');

if ~exist(paperassets, 'dir'), mkdir(paperassets); end

% ---- Helper inline functions --------------------------------------------

R_eff =@(A) 1.2e-15.*A.^(1/3); % meters (not explicitly used in ratio form)

L_NF =@(A) kappa_NF.*R_eff(A); % meters
% Photon: δ = η_γ * E_src * (α_EM * Z_eff / R_eff) * (L_NF/c) -> scales as α_EM * Z_eff * kappa_NF

% Charged: δ = η_q * q * (α_EM * Z_eff / R_eff) * L_NF -> scales as α_EM * Z_eff * kappa_NF

% We absorb geometry into η_γ, η_q via anchors (Eq. usage).

% ---- Calibrate photon coupling η_γ on 57Fe 14.4 keV line ----------------------

Z_Fe = 26; A_Fe = 57;

E_Fe_eV = 14.4e3; % 14.4 keV

delta_Fe_eV = 2.00; % anchor shift

Zeff_Fe = ZeffFrac * Z_Fe;

% δ/E = η_γ * α_EM * ZeffFrac * Z * kappa_NF * (unit factors absorbed)

eta = (delta_Fe_eV / E_Fe_eV) / (alpha_EM * ZeffFrac * Z_Fe * kappa_NF);

% ---- Calibrate charged coupling η_q on 300 eV at 300 keV (Ba-like Z≈56) -----

Z_Ba = 56; A_Ba = 137;

E_ref_eV = 300e3; delta_ref_eV = 300;

eta_q = (delta_ref_eV / E_ref_eV) / (alpha_EM * ZeffFrac * Z_Ba * kappa_NF);

% ---- Define gamma-line set (extend as needed) ---------------------------------

gamma = struct('label', {'~{57}Fe (14.4 keV)', '~{137}Ba (661.7 keV)'}, ... 
'Z', [26, 56], ... 
'A', [57, 137], ... 
'E_eV', [14.4e3, 661.7e3] ... 
);

% ---- Define charged-ejectile set (extend as needed) ----------------------------

beta = struct('label', {'conv e^- (0.2-0.5 MeV mid)', 'β endpoint (few MeV mid)'}, ... 
'Z', [56, 56], ... % Ba-like environment for illustration 
'A', [137, 137], ...
```matlab
for i = 1:numel(B.label)
    fprintf fid, '%s & %.3f & %.2f & %.2f \\\\n', ...
    latexesc(B.label{i}), B.E_eV(i)/1e6, delta_keV(i), frac_1e4(i));
end
fprintf fid, '\\bottomrule\end{tabular}\n');
end

function s = latexesc(x)
s = string(x);
s = strrep(s, '\', '\textbackslash{}');
s = strrep(s, '_', '\_');
s = strrep(s, '%', '\%');
s = strrep(s, '&', '\&');
s = strrep(s, '#', '\#');
s = strrep(s, '{', '\{');
s = strrep(s, '}', '\}');
end
```

<caption>Listing 27: Make CSCT tables from the closed-form model (canonical)</caption>

*Figures and tables produced correspond to supplementary near-field diagnostics (Appendix K).*
## Y.15 Fission Energetics and Tables

**Dependencies:** Listings 2 (p. 213), 44 (p. 297), 29 (p. 263), 30 (p. 266), 31 (p. 268), 32 (p. 270), 28 (p. 257), 34 (p. 271), 4 (p. 215), 5 (p. 216), 26 (p. 252).

**Workflow.** Run the one-time energy normalization in Listing 44 (p. 297). Then calibrate the scission smoothness for $^{235}$U with Listing 29 (p. 263) (to set C.fission.scission_w_fm). Finally, regenerate the fission summary with Listing 28 (p. 257). No per-observable tuning is performed.

**Calibration step (single amplitude rescale).** The daughters are rescaled by a single global amplitude factor $s$, determined by a bracketing+bisection solve inside eval_Efiss so that $E_1+E_2 \rightarrow E_0 - E_{fiss}^{\text{target}}$ with $E_{fiss}^{\text{target}} = 170.0$ MeV for thermal-neutron-induced $^{235}$U(n,f). This is a one-time scission calibration (no per-observable tuning).

Table 42 (p. 257) shows the TSRT fission summary for thermal-neutron-induced $^{235}$U(n,f). Table 42 (p. 257) is generated by Listing 28 (p. 257) using the one-time calibrations from Listings 44 (p. 297) (energy scale) and 29 (p. 263) (scission smoothness); no per-observable tuning is applied.

**Reproducibility note.** The numeric value in Table 42 (p. 257) is generated deterministically by Listing 28 (p. 257) using the fixed energy normalization (Listing 44 (p. 297)) and the scission-width calibration (Listing 29 (p. 263)).

Table 42: TSRT fission summary for thermal-neutron-induced $^{235}$U(n,f). Values are produced by the workflow of Appendix Y.15 using the calibrated energy normalization (Listing 44 (p. 297)) and scission width (Listing 29 (p. 263)). No per-observable tuning is applied. This table documents the absolute energy calibration verified by the workflow; observable comparisons appear in Table 4 (p. 12).

<table><thead><tr><td>System</td><td>E<sub>fiss</sub> (MeV)</td></tr></thead><tbody><tr><td><sup>235</sup>U(n,f)</td><td>170.028</td></tr></tbody></table>

Generated directly by Listing 28 (p. 257), which computes the absolute fission energy using the globally calibrated C.norm_energy (from Listing 44 (p. 297)) and the fixed scission width (Listing 29 (p. 263)).

**Purpose:** Compute TSRT fission energetics (barriers, action integrals, and derived quantities) and assemble summary tables used throughout the main text and figures.

**Used in the paper:** Main text Section 6.2, as well as the composite fission-energy and hindrance figures.

**Inputs:** tsrtconstants.m (with calibrated C.norm_energy and C.fission.scission_w_fm), curvature-field builders (tsrt_buildVO235.m, scission builders), and the core curvature/energy evaluators.

**Outputs:** fission_results.mat (struct array F) and paper-assets-tex/tab_fission_summary.tex (or ./tab_fission_summary.tex if the folder is absent).

**Reproducibility note.** The value in Table 42 (p. 257) is generated by Listing 28 using the fixed energy normalization (Listing 44) and the single scission-width calibration (Listing 29); no per-observable tuning is applied.

make_fission_summary_now.m

% FILE: make_fission_summary_now.m

% PURPOSE: Build the U-235 fission summary (no legacy edits, no tsrtmakegrid call).

% OUTPUTS: fission_results.mat (struct array F), and LaTeX table (Listing~\ref{lst:fissiontablemake} if being used)

% NOTES: Uses calibrated C.norm_energy (Listing~\ref{lst:TSRT_Calib_NormEn_Deuteron})
```c
% and calibrated C.fission.scission_w_fm (Listing~\refp{lst:TSRT_Calib_Fission_ScissWth_U235}).
function make_fission_summary_now()
% MAKE_FISSION_SUMMARY_NOW
% End-to-end fission summary WITHOUT changing legacy files.
% Strategy:
% 1) Build Simpson grid (no tsrtmakegrid).
% 2) Build U-235 parent (tsrt_buildVO parent_U235), compute E0.
% 3) Try tuned scission builder; sweep targetfraction' downward to avoid overlap
% energy blow-up.
% 4) If still Efiss<=0, fall back to partition builder.
% 5) Emit fission_results.mat (struct array F) and tab_fission_summary.tex.

clc;
if empty(strfind(path, pwd)), addpath(genpathling)); end

% ---- constants (merge calibrated) ----
C = tsrtconstants();
if exist('tsrtconstantsCalibrated.mat', 'file') == 2
S = load('tsrtconstantsCalibrated.mat');
if isstruct(S)
    f = fieldnames(S); for i=1:numel(f), C.(f{i}) = S.(f{i}); end
disp(' [FISSION] Merged tsrtconstantsCalibrated.mat into C.');
end
end
if ~isfield(C,'norm_energy') || ~isscalar(C.norm_energy) || ~isfinite(C.
norm_energy) || C.norm_energy<=0
    error('FISSION:norm_energy','C.norm_energy missing/invalid. Set it in
tsrtconstants.m and re-run.');
end
if ~isfield(C,'fission') || ~isstruct(C.fission), C.fission = struct(); end
if ~isfield(C.fission,'scission_w_fm') || ~isfinite(C.fission.scission_w_fm)
    C.fission.scission_w_fm = 1.00; % fallback
end
if ~isfield(C.fission,'targetfraction') || ~isfinite(C.fission.target_fraction)
    C.fission.target_fraction = 0.98; % start
end

% Simpson spherical grid (with aliases) ----
if ~isfield(C,'N_r')|| empty(C.N_r), C.N_r=96; end
if ~isfield(C,'N_theta')|| empty(C.N_theta), C.N_theta=64; end
if ~isfield(C,'N phi')|| empty(C.N phi), C.N phi=64; end
if ~isfield(C,'r_max')|| empty(C.r_max), C.r_max=12e-15; end

Nr=max(3,C.N_r); Nt=max(3,C.N_theta); Np=max(3,C.N phi);
r=linspace(0,C.r_max,Nr); th=linspace(0,pi,Nt); ph=linspace(0,2*pi,Np);
wr = simpson weights(Nr) * (C.r_max/(Nr-1));
wt = simpson weights(Nt) * (pi/(Nt-1));
wp = simpson weights(Np) * ((2*pi)/(Np-1));
wt = wt.* sin(th); % match with energy integral

G = struct();
G.r=r; G theta=th; G phi=ph;
G.w_r=wr; G.w_theta=wt; G.w phi=wp;
G.wr=wr; G.wrad=wr; G.w_th=wt; G.wtheta=wt; G.wth=wt; G.wph=wp; G.wphi=wp; G.wph=wp;
disp(' [FISSION] Using manual Simpson-weighted grid (robust fallback).');
```
% ---- Parent and E0 ----
assert_file('tsrt_buildVO235.m');
XiVO parent = tsrt_buildVO235(G, C);
K0 = tsrtcurvaturetensor(G, XiVO parent, C);
E0_J = tsrtenergyfromcurvature(G, K0, C);
E0_MeV = E0_J / 1.602176634e-13;

% ---- Try tuned scission builder first, sweeping target
bname = '';
if exist('tsrt_build_daughters_scission_tuned.m', 'file') == 2
bname = 'tsrt_build_daughters_scission_tuned';
else if exist('tsrt_build_daughters_scission.m', 'file') == 2
bname = 'tsrt_build_daughters_scission';
end

Xi_d1 = []; Xi_d2 = []; which_path = '';
Efiss_MeV = NaN;

if ~isempty(bname)
tf_list = unique([C.fission.target
0.96,0.94,0.92,0.90,0.88,0.86,0.85,0.82,0.80]);
for tf = tf_list
[Xd1, Xd2] = try_build_scission(bname, G, C, XiVO parent, E0_J, tf);
[Ef_MeV, ok, Xd1, Xd2] = eval_Efiss(G, C, E0_J, E0_MeV, Xd1, Xd2);
if ok
    Xi_d1 = Xd1; Xi_d2 = Xd2; Efiss_MeV = Ef_MeV;
    which_path = sprintf('%s(target fraction=%2f)', bname, tf);
    break
end
end
end

% ---- Fallback: partition builder if still not positive ----
if ~isfinite(Efiss_MeV) || Efiss_MeV <= 0
    if exist('tsrt_build_daughters_partition.m', 'file') == 2
        z0 = 0.0; w_m = C.fission.scission_w FM * 1e-15;
        [Xd1, Xd2] = tsrt_build_daughters_partition(G, XiVO parent, z0, w_m);
        [Ef_MeV, ok, Xd1, Xd2] = eval_Efiss(G, C, E0_J, E0_MeV, Xd1, Xd2);
        if ok
            Xi_d1 = Xd1; Xi_d2 = Xd2; Efiss_MeV = Ef_MeV;
            which_path = sprintf('partition(z0=0,w=%3ffm)', C.fission.scission_w FM);
        else
            % keep the numeric value even if <=0 so the table shows the truth
            Xi_d1 = Xd1; Xi_d2 = Xd2; Efiss_MeV = Ef_MeV;
            which_path = sprintf('partition(z0=0,w=%3ffm, Efiss<=0)', C.fission.scission_w FM);
        end
    else
        which_path = sprintf('%s(no positive Efiss found)', ternary(~isempty(bname), bname, 'noVO builder'));
    end
end

% ---- Emit results (struct array F) ----
F = struct('nucleus', 'B_TSRT_MeV', 'S_action', ...
'Efiss_MeV', sprintf('%.*g', Efiss_MeV), ...
'note', sprintf('%s, E0=%.*g MeV', which_path, E0_MeV));
F = F();
save('fission_results.mat', 'F');
disp('[FISSION] Saved fission_results.mat with struct array F.');

% ---- Table ----
outA = fullfile('paper-assets-tex','tab_fission_summary.tex');
try
    fissiontablemake(outA);
    fprintf(' [FISSION] Table written: %s
', outA);
    fprintf(' [FISSION] In LaTeX, use: \
    input{\paperassets/tab_fission_summary.
tex}

');
catch ME
    warning('Cannot open %s (%s). Falling back to ./tab_fission_summary.tex', outA
        , ME.message);
    outB = 'tab_fission_summary.tex';
    fissiontablemake(outB);
    fprintf(' [FISSION] Table written: %s
', outB);
    fprintf(' [FISSION] In LaTeX, use: \
    input{tab_fission_summary.tex}

');
end
end

% ---- helpers ----
function [Efiss_MeV, ok, Xi_d1, Xi_d2] = eval_Efiss(G, C, E0_J, E0_MeV, Xi_d1, Xi_d2)
% Target Efiss (MeV)
Efiss_target_MeV = 170.0; % target

% Helper: compute sum and Efiss for given daughters
function [E_sum_MeV, Efiss_MeV] = sum_and_Efiss(X1, X2)
    K1 = tsrtcurvaturetensor(G, X1, C);
    E1_J = tsrtenergyfromcurvature(G, K1, C);
    K2 = tsrtcurvaturetensor(G, X2, C);
    E2_J = tsrtenergyfromcurvature(G, K2, C);
    E_sum_MeV = (E1_J + E2_J) / 1.602176634e-13;
    Efiss_MeV = E0_MeV - E_sum_MeV;
end

% Helper: scale daughters by s (without mutating originals)
function [X1b, X2b] = scaled(X1, X2, s)
    X1b = X1; X2b = X2;
    if isfield(X1b, 'profile')
        X1b.profile = s * X1b.profile;
    elsif isfield(X1b, 'amp')
        X1b.amp = s * X1b.amp;
    else
        error('eval_Efiss: noField', 'Daughter 1 has neither .profile nor .amp to
scale.');
    end
    if isfield(X2b, 'profile')
        X2b.profile = s * X2b.profile;
    elsif isfield(X2b, 'amp')
        X2b.amp = s * X2b.amp;
    else
        error('eval_Efiss: noField', 'Daughter 2 has neither .profile nor .amp to
scale.');
    end
end
% Current energies

[E_sum0, Efiss0] = sum_and_Efiss(Xi_d1, Xi_d2);

% Desired sum

E_target_sum = max(E0_MeV - Efiss_target_MeV, 0.0);
if E_target_sum == 0
    % Trivial corner case
    Efiss_MeV = Efiss0;
    ok = isfinite(Efiss_MeV) && Efiss_MeV > 0;
    fprintf(' [FISSION] Final Efiss = %.3f MeV (target %.1f) \n ', Efiss_MeV,
            Efiss_target_MeV);
    return;
end
if ~ (isfinite(E_sum0) && E_sum0 > 0)
    Efiss_MeV = Efiss0;
    ok = isfinite(Efiss_MeV) && Efiss_MeV > 0;
    fprintf(' [FISSION] Final Efiss = %.3f MeV (target %.1f) | unable to scale (
            E_sum0 invalid) \n ',
            Efiss_MeV, Efiss_target_MeV);
    return;
end

% Initial quadratic guess for scale

s0 = sqrt(E_target_sum / E_sum0);

% Bisection bracket on s (Efiss decreases as s increases -- monotonic)

s_lo = max(s0*1e-3, 1e-6);
s_hi = s0*1e+3;

% Ensure bracket actually brackets the target (Efiss_lo >= target >= Efiss_hi)

[X1b, X2b] = scaled(Xi_d1, Xi_d2, s_lo);
[~, Ef_lo] = sum_and_Efiss(X1b, X2b);
[X1b, X2b] = scaled(Xi_d1, Xi_d2, s_hi);
[~, Ef_hi] = sum_and_Efiss(X1b, X2b);

% If not bracketed, expand until it is or until limits explode

expand_count = 0;
while ~((Ef_lo >= Efiss_target_MeV) && (Ef_hi <= Efiss_target_MeV)) &&
    expand_count < 20
    if Ef_lo < Efiss_target_MeV
        s_lo = s_lo * 0.1;
    end
    if Ef_hi > Efiss_target_MeV
        s_hi = s_hi * 10;
    end
    [X1b, X2b] = scaled(Xi_d1, Xi_d2, s_lo);
    [~, Ef_lo] = sum_and_Efiss(X1b, X2b);
    [X1b, X2b] = scaled(Xi_d1, Xi_d2, s_hi);
    [~, Ef_hi] = sum_and_Efiss(X1b, X2b);
    expand_count = expand_count + 1;
end

% If still not bracketed, fall back to closest we have

if ~((Ef_lo >= Efiss_target_MeV) && (Ef_hi <= Efiss_target_MeV))
    % pick whichever is closer
    if abs(Ef_lo - Efiss_target_MeV) < abs(Ef_hi - Efiss_target_MeV)
        s_star = s_lo;
    end