Mass Spectrum and Hessian Analysis

Mass Spectrum and Hessian Analysis#

What’s in this notebook? This notebook demonstrates how to compute the scalar mass spectrum at flux vacua using the Hessian of the scalar potential. It covers the four Hessian computation modes available in JAXVacua, the mass matrix after canonical normalisation, eigenvalue analysis, and mass hierarchies.

In this notebook, you will learn:

How to compute the Hessian of the scalar potential via model.hessian(mode=...)

The four computation modes: None (autodiff), "SUSY", "SUGRA", "real"

How to obtain the canonically normalised mass matrix via model.mass_matrix

How to analyse eigenvalue spectra and identify light/heavy directions

How to verify consistency between Hessian implementations

Prerequisites: NB05: Finding flux vacua.

Related: NB05: Finding flux vacua for the vacuum data used below.

(Created: March 2026)

Outline#

Imports and setup
Model and test points
The Hessian: four computation modes
Stress test: random points
Mass matrix and eigenvalue spectrum
Eigenvalue distributions across an ensemble
No-scale vs full SUGRA Hessian
Appendix: SUGRA Hessian validation
Take-aways
Further reading

Imports and setup#

import numpy as np
import time
import warnings
warnings.filterwarnings('ignore')

import jax
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)
from jax import vmap
import time

import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams.update({"figure.dpi": 200, "font.size": 10, "axes.labelsize": 11,
                      "figure.figsize": (7, 3), "legend.frameon": False,
                      "xtick.direction": "in", "ytick.direction": "in"})

import jaxvacua as jvc

Model and test points#

We set up a two-modulus LCS model and prepare two types of test points:

A generic point — random moduli and fluxes where \(D_AW \neq 0\)
A SUSY vacuum — a true solution to \(D_AW = 0\) found by Newton’s method

# Model
model = jvc.FluxVacuaFinder(h12=2, model_ID=1, maximum_degree=2, limit="LCS", model_type="KS")

# ── Generic point ──
rng = np.random.default_rng(42)
z_gen  = jnp.array(rng.uniform(-0.3, 0.3, 2) + 1j * rng.uniform(2, 4, 2))
zc_gen = jnp.conj(z_gen)
tau_gen  = complex(rng.uniform(-0.3, 0.3), rng.uniform(2, 5))
tauc_gen = jnp.conj(tau_gen)
fl_gen = jnp.array(rng.integers(-5, 6, 2 * model.n_fluxes).astype(float))

DW_gen = model.DW(z_gen, zc_gen, tau_gen, tauc_gen, fl_gen)
print(f"Generic point: |DW| = {float(jnp.sum(jnp.abs(DW_gen))):.2f}")

# ── SUSY vacuum ──
sampler = jvc.data_sampler(model, flux_bounds=[-8, 8], moduli_bounds=[1, 4], dilaton_bounds=[1, 4])
rns_key = jvc.PRNGSequence(123)
moduli_all, tau_all, fluxes_all, res_all = model.sample_SUSY_flux_vacua(
    N=20, sampler=sampler, rns_key=rns_key,
    max_iters=100, tol=1e-12, max_tadpole=model.D3_tadpole,
    mode="ISD", vmap_dim=20, print_progress=False,
)

# Pick the vacuum with smallest residual
best = int(jnp.argmin(res_all))
z_susy  = moduli_all[best]
zc_susy = jnp.conj(z_susy)
tau_susy  = tau_all[best]
tauc_susy = jnp.conj(tau_susy)
fl_susy = fluxes_all[best]

DW_susy = model.DW(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy)
W0 = model.superpotential_gauge_invariant(z_susy, tau_susy, fl_susy)
print(f"SUSY vacuum:   |DW| = {float(jnp.sum(jnp.abs(DW_susy))):.2e},  |W0| = {float(jnp.abs(W0)):.4f}")
print(f"  z  = {z_susy}")
print(f"  τ  = {tau_susy}")

The Hessian: four computation modes#

The FluxEFT class provides four modes for computing the Hessian via model.hessian(mode=...):

Mode	Method	Valid at	Approach
`None`	`_hessian_general`	Any point	`jacrev(dV)` — exact autodiff
`"SUSY"`	`_hessian_SUSY`	\(D_AW = 0\) only	Analytical SUGRA formula (simplified at SUSY)
`"SUGRA"`	`_hessian_SUGRA`	Any point	Explicit SUGRA building blocks (\(\Gamma\), \(R\), \(D_IW\), \(K\), \(W\))
`"real"`	`ddV_x`	Any point	Hessian in real coordinates \((\mathrm{Re}\,z, \mathrm{Im}\,z, \ldots)\)

All modes should agree at SUSY points. At generic points, "SUSY" mode will differ.

Comparison at the generic point#

# Adapted from hessian_SUGRA_analysis.ipynb, Section 12

modes = [None, "SUSY", "SUGRA", "real"]

print("Generic point (|DW| >> 0):")
print(f"|DW| = {float(jnp.sum(jnp.abs(DW_gen))):.2f}\n")
print(f"{'Mode':<10} {'Time (ms)':>10} {'Max |err| vs None':>20} {'Rel err':>12} {'Eigenvalues (real part)':>30}")
print("-" * 90)

H_ref = None
for mode in modes:
    # Warm up JIT
    _ = model.hessian(z_gen, zc_gen, tau_gen, tauc_gen, fl_gen, noscale=True, mode=mode)
    
    t0 = time.perf_counter()
    H = model.hessian(z_gen, zc_gen, tau_gen, tauc_gen, fl_gen, noscale=True, mode=mode)
    jax.block_until_ready(H)
    dt = (time.perf_counter() - t0) * 1000
    
    eigs = jnp.sort(jnp.linalg.eigvals(H).real)
    
    if mode == "real":
        eigs = eigs/2
    if mode is None:
        H_ref = H
        err_str = "—"
        rel_str = "—"
    else:
        if mode == "real":
            # Real mode has different shape — compare eigenvalues instead
            eigs_ref = jnp.sort(jnp.linalg.eigvals(H_ref).real)
            err = float(jnp.max(jnp.abs(eigs_ref - eigs)))
            scale = float(jnp.max(jnp.abs(eigs_ref)))
        else:
            err = float(jnp.max(jnp.abs(H_ref - H)))
            scale = float(jnp.max(jnp.abs(H_ref)))
        err_str = f"{err:.2e}"
        rel_str = f"{err/scale:.2e}" if scale > 0 else "—"
    
    eig_str = ", ".join(f"{float(e):.2f}" for e in eigs[:4])
    label = str(mode) if mode else "None (autodiff)"
    print(f"{label:<16} {dt:8.1f}ms   {err_str:>18}   {rel_str:>10}   [{eig_str}]")

The "SUGRA" mode shows a relative error of \(\sim 10^{-8}\) against the autodiff reference (mode=None). This is not a bug — it reflects the accumulated floating-point error from chaining many intermediate autodiff operations: the SUGRA formula builds the Hessian from third derivatives of the Kähler potential (\(\partial_j\partial_k\partial_{\bar l} K\)), Christoffel symbols, the Riemann tensor, and a 9-term product-rule expansion. Each differentiation step loses \(\sim 1\)–\(2\) digits of precision.

In contrast, mode=None computes jacrev(dV) in a single second-order autodiff pass, reaching machine precision (\(\sim 10^{-15}\)).

For practical purposes, \(10^{-8}\) relative error is negligible — it is well below the sensitivity of any physical observable. Use mode="SUGRA" when you need access to individual SUGRA building blocks (\(D_IW\), \(\Gamma^i_{jk}\), \(R_{i\bar j k\bar l}\)); use mode=None when maximum numerical precision is required.

Comparison at the SUSY vacuum#

print("SUSY vacuum (|DW| ≈ 0):")
print(f"|DW| = {float(jnp.sum(jnp.abs(DW_susy))):.2e}\n")
print(f"{'Mode':<10} {'Time (ms)':>10} {'Max |err| vs None':>20} {'Rel err':>12}")
print("-" * 65)

H_ref_susy = None
for mode in modes:
    _ = model.hessian(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy, noscale=True, mode=mode)
    
    t0 = time.perf_counter()
    H = model.hessian(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy, noscale=True, mode=mode)
    jax.block_until_ready(H)
    dt = (time.perf_counter() - t0) * 1000
    
    if mode is None:
        H_ref_susy = H
        print(f"{'None (autodiff)':<16} {dt:8.1f}ms   {'—':>18}   {'—':>10}")
    else:
        if mode == "real":
            eigs_ref = jnp.sort(jnp.linalg.eigvals(H_ref_susy).real)
            eigs = jnp.sort(jnp.linalg.eigvals(H).real)/2
            err = float(jnp.max(jnp.abs(eigs_ref - eigs)))
            scale = float(jnp.max(jnp.abs(eigs_ref)))
        else:
            err = float(jnp.max(jnp.abs(H_ref_susy - H)))
            scale = float(jnp.max(jnp.abs(H_ref_susy)))
        rel = err/scale if scale > 0 else 0
        print(f"{str(mode):<16} {dt:8.1f}ms   {err:>18.2e}   {rel:>10.2e}")

print("\n→ All modes agree at SUSY (as expected).")

Stress test: random points#

Verify hessian(mode="SUGRA") matches hessian(mode=None) across 50 random field-space points. This checks the SUGRA Hessian implementation away from special points, where cancellations at SUSY vacua would make the comparison less discriminating.

rng = np.random.default_rng(42)
n_tests = 50
rel_errors = []

for i in range(n_tests):
    re_z = rng.uniform(-0.5, 0.5, model.h12)
    im_z = rng.uniform(2.0, 5.0, model.h12)
    z_i = jnp.array(re_z + 1j * im_z)
    zc_i = jnp.conj(z_i)
    tau_i = complex(rng.uniform(-0.4, 0.4), rng.uniform(2, 6))
    tauc_i = jnp.conj(tau_i)
    fl_i = jnp.array(rng.integers(-8, 9, 2 * model.n_fluxes).astype(float))
    ns = bool(rng.choice([True, False]))
    
    H_ref = model._hessian_general(z_i, zc_i, tau_i, tauc_i, fl_i, noscale=ns)
    H_sug = model._hessian_SUGRA(z_i, zc_i, tau_i, tauc_i, fl_i, noscale=ns)
    
    err = float(jnp.max(jnp.abs(H_ref - H_sug)))
    scale = float(jnp.max(jnp.abs(H_ref)))
    rel_errors.append(err / scale if scale > 0 else 0)

rel_errors = np.array(rel_errors)
print(f"Tested {n_tests} random points:")
print(f"  Max relative error: {rel_errors.max():.2e}")
print(f"  Mean relative error: {rel_errors.mean():.2e}")
print(f"  All < 1e-10: {np.all(rel_errors < 1e-10)} ✓")

Mass matrix and eigenvalue spectrum#

The mass matrix is the Hessian after canonical normalisation by the Kähler metric:

\[ M^2_{A\bar B} = K^{A\bar C}\, V_{C\bar B} \]

The eigenvalues of \(M^2\) give the physical scalar masses-squared. At a SUSY Minkowski vacuum (\(V=0\), \(DW=0\)), we expect:

All eigenvalues \(\geq 0\) (positive semi-definite)
A mass hierarchy controlled by the flux choice

# Mass matrix at the SUSY vacuum
M2 = model.mass_matrix(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy, mode=None, noscale=True)
eigenvalues = jnp.sort(jnp.linalg.eigvals(M2).real)

print("Mass matrix eigenvalues at SUSY vacuum:")
for i, ev in enumerate(eigenvalues):
    print(f"  m²_{i} = {float(ev):+.4e}")

print(f"\nAll ≥ 0: {bool(jnp.all(eigenvalues > -1e-8))} ✓")
print(f"Hierarchy: m²_max / m²_min = {float(eigenvalues[-1] / eigenvalues[0]):.1f}")

Eigenvalue distributions across an ensemble#

Compute the mass spectrum at every vacuum in our sample and visualise the distributions.

# Compute mass eigenvalues for all vacua
all_eigenvalues = []
all_W0 = []

for i in range(len(moduli_all)):
    z_i = moduli_all[i]
    tau_i = tau_all[i]
    fl_i = fluxes_all[i]
    
    M2_i = model.mass_matrix(z_i, jnp.conj(z_i), tau_i, jnp.conj(tau_i), fl_i,
                              mode="SUSY", noscale=True)
    eigs_i = jnp.sort(jnp.linalg.eigvals(M2_i).real)
    all_eigenvalues.append(np.array(eigs_i))
    all_W0.append(float(jnp.abs(model.superpotential_gauge_invariant(z_i, tau_i, fl_i))))

all_eigenvalues = np.array(all_eigenvalues)
all_W0 = np.array(all_W0)

n_fields = all_eigenvalues.shape[1]
print(f"Computed mass spectra for {len(all_eigenvalues)} vacua ({n_fields} fields each)")

fig, axes = plt.subplots(2, 2, figsize=(8.0, 6.0), constrained_layout=True)
axes = axes.ravel()

colors = plt.cm.viridis(np.linspace(0.2, 0.9, n_fields))
positive_eigs = all_eigenvalues[np.isfinite(all_eigenvalues) & (all_eigenvalues > 0)]
log_all = np.log10(positive_eigs)
lo, hi = np.nanpercentile(log_all, [2, 98])
bins = np.linspace(lo, hi, 28)

# (a) Histogram of eigenvalues by field index, using common finite bins.
for j in range(n_fields):
    ev_j = all_eigenvalues[:, j]
    ev_j = ev_j[np.isfinite(ev_j) & (ev_j > 0)]
    if len(ev_j) > 0:
        axes[0].hist(np.log10(ev_j), bins=bins, alpha=0.5, color=colors[j],
                     label=f"$m^2_{j}$", edgecolor="white", linewidth=0.3)
axes[0].set_xlabel(r"$\log_{10}(m^2)$")
axes[0].set_ylabel("Count")
axes[0].legend(fontsize=8, ncol=2)
axes[0].text(0.05, 0.92, "(a)", transform=axes[0].transAxes, fontsize=12, fontweight="bold")

# (b) Mass hierarchy: lightest vs heaviest, with percentile limits.
m2_min = np.min(all_eigenvalues, axis=1)
m2_max = np.max(all_eigenvalues, axis=1)
valid = (m2_min > 0) & np.isfinite(m2_min) & np.isfinite(m2_max)
x = np.log10(m2_min[valid])
y = np.log10(m2_max[valid])
W_abs = np.abs(all_W0[valid])
sc = axes[1].scatter(x, y, s=12, alpha=0.65, c=W_abs, cmap="viridis",
                     edgecolors="none", rasterized=True)
xlo, xhi = np.nanpercentile(x, [2, 98])
ylo, yhi = np.nanpercentile(y, [2, 98])
lo_diag, hi_diag = min(xlo, ylo), max(xhi, yhi)
axes[1].plot([lo_diag, hi_diag], [lo_diag, hi_diag], 'k--', lw=0.6, alpha=0.35)
axes[1].set_xlim(xlo, xhi)
axes[1].set_ylim(ylo, yhi)
axes[1].set_xlabel(r"$\log_{10}(m^2_{\mathrm{min}})$")
axes[1].set_ylabel(r"$\log_{10}(m^2_{\mathrm{max}})$")
cb = fig.colorbar(sc, ax=axes[1], shrink=0.9, pad=0.02)
cb.set_label(r"$|W_0|$", fontsize=10)
cb.ax.tick_params(labelsize=8)
axes[1].text(0.05, 0.92, "(b)", transform=axes[1].transAxes, fontsize=12, fontweight="bold")

# (c) Hierarchy ratio versus |W0|.
hierarchy = m2_max[valid] / m2_min[valid]
axes[2].scatter(W_abs, hierarchy, s=12, alpha=0.65, color="#7c3aed",
                edgecolors="none", rasterized=True)
axes[2].set_xscale("log")
axes[2].set_yscale("log")
axes[2].set_xlabel(r"$|W_0|$")
axes[2].set_ylabel(r"$m^2_{\mathrm{max}}/m^2_{\mathrm{min}}$")
axes[2].text(0.05, 0.92, "(c)", transform=axes[2].transAxes, fontsize=12, fontweight="bold")

# (d) Lightest eigenvalue distribution.
axes[3].hist(np.log10(m2_min[valid]), bins=24, color="#16a34a", alpha=0.7,
             edgecolor="white", linewidth=0.3)
axes[3].set_xlabel(r"$\log_{10}(m^2_{\mathrm{min}})$")
axes[3].set_ylabel("Count")
axes[3].text(0.05, 0.92, "(d)", transform=axes[3].transAxes, fontsize=12, fontweight="bold")

plt.show()

No-scale vs full SUGRA Hessian#

The no-scale potential has \(\lambda = 0\) (Kähler moduli F-terms cancel \(-3|W|^2\)), while the full SUGRA potential has \(\lambda = 3\). We compare the eigenvalue spectra.

# At a SUSY point we can use the specialised Hessian formula.
# This avoids compiling the generic autodiff Hessian twice for this illustrative comparison.
H_noscale = model.hessian(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy, noscale=True, mode="SUSY")
H_full    = model.hessian(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy, noscale=False, mode="SUSY")

eigs_ns = jnp.sort(jnp.linalg.eigvals(H_noscale).real)
eigs_full = jnp.sort(jnp.linalg.eigvals(H_full).real)

print(f"{'':>12} {'no-scale (λ=0)':>20} {'full SUGRA (λ=3)':>20}")
print("-" * 55)
for i in range(len(eigs_ns)):
    print(f"  eig_{i}:   {float(eigs_ns[i]):>+18.4e}   {float(eigs_full[i]):>+18.4e}")

print(f"\nDifference from -3|W|² correction: max|Δeig| = {float(jnp.max(jnp.abs(eigs_ns - eigs_full))):.4e}")

Appendix: SUGRA Hessian validation#

This appendix cross-validates the SUGRA Hessian implementation at many random points and tabulates all four mode values against each other.

Stress test: random moduli, dilaton, and flux choices#

Verify model._hessian_SUGRA matches _hessian_general across many randomly chosen field-space points and flux vectors. This is the generic-point SUGRA formula, so it should match the autodiff reference for both no-scale and full SUGRA potentials. The specialised mode="SUSY" formula is checked separately at SUSY points.

rng = np.random.default_rng(42)
n_tests = 50
max_rel_errors = []
max_abs_errors = []

print(f"Running {n_tests} random test points (h12={model.h12})...\n")
print(f"{'#':>3} {'Im(z1)':>7} {'Im(z2)':>7} {'Im(τ)':>7} {'|DW|':>10} {'noscale':>8} "
      f"{'|err|':>10} {'rel err':>10} {'ok?':>4}")
print("-" * 82)

for i in range(n_tests):
    # Random moduli: Re(z) in [-0.5, 0.5], Im(z) in [2, 5].
    re_z = rng.uniform(-0.5, 0.5, model.h12)
    im_z = rng.uniform(2.0, 5.0, model.h12)
    z_i = jnp.array(re_z + 1j * im_z)
    zc_i = jnp.conj(z_i)

    # Random axio-dilaton: Re(tau) in [-0.5, 0.5], Im(tau) in [1, 8].
    tau_i = complex(rng.uniform(-0.5, 0.5) + 1j * rng.uniform(1.0, 8.0))
    tauc_i = jnp.conj(tau_i)

    # Random integer fluxes in [-5, 5].
    fl_i = jnp.array(rng.integers(-5, 6, size=2 * model.n_fluxes), dtype=float)
    if jnp.all(fl_i == 0):
        fl_i = fl_i.at[0].set(1.0)

    ns_i = bool(rng.choice([True, False]))

    H_ref_i = model._hessian_general(z_i, zc_i, tau_i, tauc_i, fl_i, noscale=ns_i)
    H_sugra_i = model._hessian_SUGRA(z_i, zc_i, tau_i, tauc_i, fl_i, noscale=ns_i)

    abs_err = float(jnp.max(jnp.abs(H_ref_i - H_sugra_i)))
    scale_i = float(jnp.max(jnp.abs(H_ref_i)))
    rel_err = abs_err / scale_i if scale_i > 0 else 0.0

    max_rel_errors.append(rel_err)
    max_abs_errors.append(abs_err)

    DW_i = model.DW(z_i, zc_i, tau_i, tauc_i, fl_i)
    dw_norm = float(jnp.sum(jnp.abs(DW_i)))
    ok = "yes" if rel_err < 1e-10 else "no"

    if i < 10 or rel_err > 1e-10:  # print first 10 + any failures
        print(f"{i:>3} {float(im_z[0]):>7.2f} {float(im_z[1]):>7.2f} "
              f"{float(jnp.imag(tau_i)):>7.2f} {dw_norm:>10.2e} "
              f"{'True' if ns_i else 'False':>8} "
              f"{abs_err:>10.2e} {rel_err:>10.2e} {ok:>4}")

max_rel = max(max_rel_errors)
max_abs = max(max_abs_errors)
n_pass = sum(1 for r in max_rel_errors if r < 1e-10)

print("-" * 82)
print(f"\nResults: {n_pass}/{len(max_rel_errors)} passed (rel err < 1e-10)")
print(f"  Max relative error: {max_rel:.2e}")
print(f"  Max absolute error: {max_abs:.2e}")
print(f"  Median rel error:   {np.median(max_rel_errors):.2e}")

if n_pass == len(max_rel_errors):
    print(f"\nAll {n_pass} random test points match to machine precision.")

Comparison of all Hessian implementations#

The FluxEFT class provides four modes for computing the Hessian via model.hessian(mode=...):

Mode	Method	Valid at	Approach
`None`	`_hessian_general`	Generic points	`jacrev(dV)` — exact autodiff
`"SUSY"`	`_hessian_SUSY`	\(D_AW = 0\) only	Analytical SUGRA formula (simplified at SUSY)
`"SUGRA"`	`_hessian_SUGRA`	Generic points	Explicit SUGRA building blocks (\(\Gamma\), \(R\), \(D_IW\), \(K\), \(W\))
`"real"`	`ddV_x`	Generic points	Real-field Hessian via `jacrev` of real-field potential

Below we compare all four at generic and SUSY points.

Interpreting the Hessian comparison#

The modes below do not all have the same domain of validity or coordinate basis. mode=None is the autodiff reference at generic points. mode="SUGRA" should agree with it and is backed by tests for the Christoffel symbols, Riemann tensor and SUGRA Hessian blocks. mode="SUSY" uses formulae simplified by \(D_IW=0\), so it is expected to disagree away from a SUSY point. mode="real" is a real-coordinate Hessian; the table compares eigenvalues rather than raw matrix entries because the basis differs.

modes_all = [None, "SUSY", "SUGRA", "real"]


def _mode_label(mode):
    return str(mode) if mode is not None else "None"


def _mode_error_against_reference(H, H_ref, mode, scale):
    if mode == "real":
        # The real-coordinate Hessian lives in a different basis.  Compare the
        # physical eigenvalues, using the same normalisation as above.
        eig_ref = jnp.sort(jnp.linalg.eigvalsh(H_ref).real)
        eig_real = jnp.sort(jnp.linalg.eigvalsh(H).real) / 2
        err = float(jnp.max(jnp.abs(eig_ref - eig_real)))
    else:
        err = float(jnp.max(jnp.abs(H_ref - H)))
    rel = err / scale if scale > 0 else 0.0
    return err, rel


# -- Generic point from Section 2 ------------------------------------------------
print("=" * 75)
print("A. GENERIC POINT (|DW| >> 0)")
print("=" * 75)
print(f"|DW| = {float(jnp.sum(jnp.abs(DW_gen))):.2f}\n")

H_all = {}
timings = {}
for mode in modes_all:
    _ = model.hessian(z_gen, zc_gen, tau_gen, tauc_gen, fl_gen, noscale=True, mode=mode)
    t0 = time.perf_counter()
    H_all[mode] = model.hessian(z_gen, zc_gen, tau_gen, tauc_gen, fl_gen, noscale=True, mode=mode)
    jax.block_until_ready(H_all[mode])
    timings[mode] = time.perf_counter() - t0

H_ref = H_all[None]
scale = float(jnp.max(jnp.abs(H_ref)))

print(f"{'Mode':>8} {'Shape':>12} {'|err| vs None':>15} {'Rel err':>12} {'Time (ms)':>10}")
print("-" * 65)
for mode in modes_all:
    H = H_all[mode]
    err, rel = _mode_error_against_reference(H, H_ref, mode, scale)
    ms = timings[mode] * 1000
    print(f"{_mode_label(mode):>8} {str(H.shape):>12} {err:>15.2e} {rel:>12.2e} {ms:>10.1f}")

# -- Also test noscale=False ------------------------------------------------------
print("\nnoscale=False (lambda=3):")
H_ref_f = model.hessian(z_gen, zc_gen, tau_gen, tauc_gen, fl_gen, noscale=False, mode=None)
scale_f = float(jnp.max(jnp.abs(H_ref_f)))
for mode in modes_all:
    H_f = model.hessian(z_gen, zc_gen, tau_gen, tauc_gen, fl_gen, noscale=False, mode=mode)
    _, rel = _mode_error_against_reference(H_f, H_ref_f, mode, scale_f)
    print(f"  {_mode_label(mode):>8}: rel err = {rel:.2e}")

# -- SUSY point from Section 2 ----------------------------------------------------
print("\n" + "=" * 75)
print("B. SUSY POINT (|DW| approx 0)")
print("=" * 75)
print(f"|DW| = {float(jnp.sum(jnp.abs(DW_susy))):.2e}\n")

H_ref_s = model.hessian(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy, noscale=True, mode=None)
scale_s = float(jnp.max(jnp.abs(H_ref_s)))

print(f"{'Mode':>8} {'Rel err vs None':>18} {'Status':>8}")
print("-" * 40)
for mode in modes_all:
    H_s = model.hessian(z_susy, zc_susy, tau_susy, tauc_susy, fl_susy, noscale=True, mode=mode)
    _, rel = _mode_error_against_reference(H_s, H_ref_s, mode, scale_s)
    status = "yes" if rel < 1e-10 else ("near" if rel < 1e-5 else "no")
    print(f"{_mode_label(mode):>8} {rel:>18.2e} {status:>8}")

# -- Stress test: all modes at 20 random points ----------------------------------
print("\n" + "=" * 75)
print("C. STRESS TEST: 20 random points, all modes")
print("=" * 75)

rng2 = np.random.default_rng(123)
results = {m: [] for m in modes_all}

for i in range(20):
    z_ = jnp.array(rng2.uniform(-0.3, 0.3, 2) + 1j * rng2.uniform(2.5, 4.5, 2))
    tau_ = complex(rng2.uniform(-0.4, 0.4) + 1j * rng2.uniform(2, 7))
    fl_ = jnp.array(rng2.integers(-4, 5, 12), dtype=float)
    if jnp.all(fl_ == 0):
        fl_ = fl_.at[0].set(1.0)
    zc_ = jnp.conj(z_)
    tauc_ = jnp.conj(tau_)
    ns_ = bool(rng2.choice([True, False]))

    H_ref_ = model.hessian(z_, zc_, tau_, tauc_, fl_, noscale=ns_, mode=None)
    scale_ = float(jnp.max(jnp.abs(H_ref_)))

    for mode in modes_all:
        H_ = model.hessian(z_, zc_, tau_, tauc_, fl_, noscale=ns_, mode=mode)
        _, rel_ = _mode_error_against_reference(H_, H_ref_, mode, scale_)
        results[mode].append(rel_)

print(f"\n{'Mode':>8} {'Max rel err':>12} {'Median':>12} {'All < 1e-10':>12}")
print("-" * 50)
for mode in modes_all:
    errs = results[mode]
    mx = max(errs)
    med = np.median(errs)
    ok = all(e < 1e-10 for e in errs)
    print(f"{_mode_label(mode):>8} {mx:>12.2e} {med:>12.2e} {'yes' if ok else 'no':>12}")

print("\nNote: mode='SUSY' is only exact at DW=0; at generic points it has O(1) error.")

Take-aways#

Hessian mode	Speed	Accuracy / valid regime	Use case
`None` (autodiff)	Baseline	Reference implementation at generic points	Default; always correct
`"SUGRA"`	Comparable after compilation	Matches the generic reference when evaluated in the same complex basis	Explicit SUGRA building-block diagnostics
`"SUSY"`	Fastest at SUSY loci	Valid only when \(D_AW=0\)	Large-scale mass computations at SUSY vacua
`"real"`	Comparable after compilation	Same physics in real coordinates; compare spectra, not raw complex-basis entries	Real-coordinate solvers and diagnostics

Use model.mass_matrix(mode="SUSY") for efficient mass computation at SUSY vacua.
Use model.hessian(mode=None) for a conservative generic-point reference.
Use mode="SUGRA" when you want to inspect the explicit Christoffel/Riemann building blocks.
Tachyonic directions (\(m^2 < 0\)) signal instability — expected at non-SUSY critical points.

Mass Spectrum and Hessian Analysis

Contents

Mass Spectrum and Hessian Analysis#

Outline#

Imports and setup#

Model and test points#

The Hessian: four computation modes#

Comparison at the generic point#

Comparison at the SUSY vacuum#

Stress test: random points#

Mass matrix and eigenvalue spectrum#

Eigenvalue distributions across an ensemble#

No-scale vs full SUGRA Hessian#

Appendix: SUGRA Hessian validation#

Stress test: random moduli, dilaton, and flux choices#

Comparison of all Hessian implementations#

Interpreting the Hessian comparison#

Take-aways#

Further reading#