Stability

stability

Stability analysis using Hessian computation via JAX autodiff.

Implements:

• Tachyonic mode detection
• Linear response theory
• Fluctuation spectrum
• Phase transition indicators
• Tunneling rate estimates

All second derivatives computed automatically via JAX.

compute_mass_matrix(phi, m, lambda_)

Compute effective mass matrix \(\partial^2 V/\partial\phi\,\partial\phi^*\) (Hessian of potential).

Convention: \(m_\text{eff}^2 = \text{sign}(m)\,m^2 + \lambda|\phi|^2\).

• \(m > 0\): quartic potential \(V = m^2/2\,|\phi|^2 + \lambda/4\,|\phi|^4\). Curvature at origin is \(+m^2\) (stable vacuum at \(\phi=0\)).
• \(m < 0\): double-well potential with tachyonic mass. Curvature at origin is \(-m^2\) (unstable vacuum, SSB).

This tells us if the field configuration is stable (all positive) or has tachyonic modes (negative values).

Parameters:

Name	Type	Description	Default
phi	Array	Scalar field (N, N, N) complex	required
m	float	Signed mass parameter. Positive → stable; negative → tachyonic at origin.	required
lambda_	float	Self-coupling	required

Returns:

Type	Description
Array	Effective mass squared at each point (N, N, N) real

Source code in jaxlatt/observables/stability.py

@jit
def compute_mass_matrix(
    phi: Array,
    m: float,
    lambda_: float,
) -> Array:
    r"""
    Compute effective mass matrix $\partial^2 V/\partial\phi\,\partial\phi^*$ (Hessian of potential).

    Convention: $m_\text{eff}^2 = \text{sign}(m)\,m^2 + \lambda|\phi|^2$.

    - $m > 0$: quartic potential $V = m^2/2\,|\phi|^2 + \lambda/4\,|\phi|^4$.
      Curvature at origin is $+m^2$ (stable vacuum at $\phi=0$).
    - $m < 0$: double-well potential with tachyonic mass.
      Curvature at origin is $-m^2$ (unstable vacuum, SSB).

    This tells us if the field configuration is stable (all positive)
    or has tachyonic modes (negative values).

    Args:
        phi: Scalar field (N, N, N) complex
        m: Signed mass parameter. Positive → stable; negative → tachyonic at origin.
        lambda_: Self-coupling

    Returns:
        Effective mass squared at each point (N, N, N) real
    """
    phi_sq = jnp.abs(phi) ** 2
    # sign(m) * m² + λ|φ|²: when m < 0 this gives negative curvature at φ=0
    # (tachyonic / symmetry-breaking). m * abs(m) = sign(m) * m² for all real m.
    m_eff_sq = m * jnp.abs(m) + lambda_ * phi_sq
    return m_eff_sq

detect_tachyonic_modes(phi, m, lambda_)

Detect tachyonic (unstable) modes where m_eff² < 0.

Parameters:

Name	Type	Description	Default
phi	Array	Scalar field configuration.	required
m	float	Mass parameter.	required
lambda_	float	Self-coupling.	required

Returns:

Type	Description
dict[str, float]	Dictionary with counts and fractions of unstable sites, extrema of
dict[str, float]	effective mass-squared, and a boolean instability flag.

Source code in jaxlatt/observables/stability.py

def detect_tachyonic_modes(
    phi: Array,
    m: float,
    lambda_: float,
) -> dict[str, float]:
    """
    Detect tachyonic (unstable) modes where m_eff² < 0.

    Args:
        phi: Scalar field configuration.
        m: Mass parameter.
        lambda_: Self-coupling.

    Returns:
        Dictionary with counts and fractions of unstable sites, extrema of
        effective mass-squared, and a boolean instability flag.
    """
    m_eff_sq = compute_mass_matrix(phi, m, lambda_)

    # Find tachyonic regions
    is_tachyonic = m_eff_sq < 0
    n_tachyonic = jnp.sum(is_tachyonic)
    fraction_tachyonic = jnp.mean(is_tachyonic)

    # Strongest instability
    min_m_sq = jnp.min(m_eff_sq)
    max_m_sq = jnp.max(m_eff_sq)

    return {
        "n_tachyonic_sites": int(n_tachyonic.item()),
        "fraction_tachyonic": float(fraction_tachyonic.item()),
        "min_mass_squared": float(min_m_sq.item()),
        "max_mass_squared": float(max_m_sq.item()),
        "has_instability": bool(min_m_sq < 0),
    }

make_energy_hessian_function(m, lambda_, dx)

Create function to compute full Hessian \(\partial^2 H/\partial\phi\,\partial\phi\) via autodiff.

This includes both potential and gradient terms:

\[H = \frac{1}{2}\int |\nabla\phi|^2 + V(\phi)\]

The Hessian operator is: \(-\nabla^2 + \partial^2 V/\partial\phi^2\)

Parameters:

Name	Type	Description	Default
m	float	Mass	required
lambda_	float	Self-coupling	required
dx	float	Lattice spacing	required

Returns:

Type	Description
	Function \(\phi \to\) Hessian operator applied to \(\phi\)

Source code in jaxlatt/observables/stability.py

def make_energy_hessian_function(
    m: float,
    lambda_: float,
    dx: float,
):
    r"""
    Create function to compute full Hessian $\partial^2 H/\partial\phi\,\partial\phi$ via autodiff.

    This includes both potential and gradient terms:

    $$H = \frac{1}{2}\int |\nabla\phi|^2 + V(\phi)$$

    The Hessian operator is: $-\nabla^2 + \partial^2 V/\partial\phi^2$

    Args:
        m: Mass
        lambda_: Self-coupling
        dx: Lattice spacing

    Returns:
        Function $\phi \to$ Hessian operator applied to $\phi$
    """

    def energy_functional(phi_flat: Array) -> float:
        r"""Energy as scalar function of flattened $\phi$."""
        # phi_flat.shape[0] is a static Python int at trace time; use pure Python
        # arithmetic so the reshape doesn't create a JAX tracer inside hessian.
        N = round(phi_flat.shape[0] ** (1.0 / 3.0))
        phi = phi_flat.reshape((N, N, N))

        # Gradient energy (finite difference Laplacian)
        laplacian = (
            jnp.roll(phi, -1, axis=0)
            + jnp.roll(phi, 1, axis=0)
            + jnp.roll(phi, -1, axis=1)
            + jnp.roll(phi, 1, axis=1)
            + jnp.roll(phi, -1, axis=2)
            + jnp.roll(phi, 1, axis=2)
            - 6 * phi
        ) / (dx**2)
        grad_energy = -0.5 * jnp.sum(jnp.real(jnp.conj(phi) * laplacian))

        # Potential energy
        pot_energy = jnp.sum(scalar_potential_energy(phi, m, lambda_))

        return (grad_energy + pot_energy) * dx**3

    # Compute Hessian via autodiff
    hess_fn = hessian(energy_functional, holomorphic=False)

    @jit
    def apply_hessian(phi: Array) -> Array:
        """Apply Hessian operator to phi."""
        phi_flat = phi.ravel()
        H = hess_fn(phi_flat)
        # For large systems, H is huge - instead compute H @ phi
        result = jnp.dot(H, phi_flat)
        return result.reshape(phi.shape)

    return apply_hessian

compute_fluctuation_spectrum_1d(phi, m, lambda_, dx)

Compute power spectrum of fluctuations \(\delta\phi\) around background.

\(P(k) = \langle |\delta\phi_k|^2 \rangle\) gives the spectrum of quantum/thermal fluctuations.

Parameters:

Name	Type	Description	Default
phi	Array	Field configuration	required
m	float	Mass	required
lambda_	float	Coupling	required
dx	float	Lattice spacing	required

Returns:

Type	Description
tuple[Array, Array]	(k_values, power_spectrum)

Source code in jaxlatt/observables/stability.py

@jit
def compute_fluctuation_spectrum_1d(
    phi: Array,
    m: float,
    lambda_: float,
    dx: float,
) -> tuple[Array, Array]:
    r"""
    Compute power spectrum of fluctuations $\delta\phi$ around background.

    $P(k) = \langle |\delta\phi_k|^2 \rangle$ gives the spectrum of quantum/thermal fluctuations.

    Args:
        phi: Field configuration
        m: Mass
        lambda_: Coupling
        dx: Lattice spacing

    Returns:
        (k_values, power_spectrum)
    """
    # FFT to k-space
    phi_k = jnp.fft.fftn(phi)

    # Power spectrum
    P_k = jnp.abs(phi_k) ** 2 / phi.size

    # Radial average (for 3D → 1D spectrum)
    # First get k-space coordinates
    kx = jnp.fft.fftfreq(phi.shape[0], d=dx) * 2 * jnp.pi
    ky = jnp.fft.fftfreq(phi.shape[1], d=dx) * 2 * jnp.pi
    kz = jnp.fft.fftfreq(phi.shape[2], d=dx) * 2 * jnp.pi

    KX, KY, KZ = jnp.meshgrid(kx, ky, kz, indexing="ij")
    K = jnp.sqrt(KX**2 + KY**2 + KZ**2)

    # Bin by |k| — vectorised with segment_sum (avoids 19× loop unroll at trace time)
    k_bins = jnp.linspace(0, K.max(), 20)
    k_centers = 0.5 * (k_bins[:-1] + k_bins[1:])

    bin_idx = jnp.searchsorted(k_bins[1:], K.ravel())  # [0, 19]; out-of-range dropped
    P_flat = P_k.ravel()
    P_sum = jax.ops.segment_sum(P_flat, bin_idx, num_segments=19)
    counts = jax.ops.segment_sum(
        jnp.ones(P_flat.size, dtype=P_flat.dtype), bin_idx, num_segments=19
    )
    P_avg = jnp.where(counts > 0, P_sum / counts, 0.0)

    return k_centers, P_avg

estimate_effective_mass_from_spectrum(k_values, power_spectrum, dx)

Estimate effective mass from fluctuation spectrum.

For free theory: \(P(k) \approx 1/(k^2 + m_\mathrm{eff}^2)\)

Fit low-\(k\) behavior to extract \(m_\mathrm{eff}\).

Parameters:

Name	Type	Description	Default
k_values	Array	Momentum values	required
power_spectrum	Array	\(P(k)\)	required
dx	float	Lattice spacing	required

Returns:

Type	Description
float	Effective mass estimate

Source code in jaxlatt/observables/stability.py

def estimate_effective_mass_from_spectrum(
    k_values: Array,
    power_spectrum: Array,
    dx: float,
) -> float:
    r"""
    Estimate effective mass from fluctuation spectrum.

    For free theory: $P(k) \approx 1/(k^2 + m_\mathrm{eff}^2)$

    Fit low-$k$ behavior to extract $m_\mathrm{eff}$.

    Args:
        k_values: Momentum values
        power_spectrum: $P(k)$
        dx: Lattice spacing

    Returns:
        Effective mass estimate
    """
    # Use only low-k modes (k < pi/4L)
    k_max = jnp.pi / (4 * dx)
    mask = k_values < k_max

    P_low = power_spectrum[mask]

    # Fit 1/P ~ k^2 + m^2
    # m^2 ~ 1/P(k->0) - k^2
    if int(jnp.sum(mask)) > 2:
        P_zero = P_low[0]
        m_eff_sq = 1.0 / (P_zero + 1e-10)
        return float(jnp.sqrt(jnp.maximum(m_eff_sq, 0)))
    else:
        return 0.0

compute_linear_response(phi_background, delta_phi, m, lambda_, dx)

Compute linear response: how system responds to perturbation.

\[\delta\phi \to \left.\frac{\partial^2 V}{\partial\phi^2}\right|_\mathrm{bg} \delta\phi\]

This is the linearized equation of motion.

Parameters:

Name	Type	Description	Default
phi_background	Array	Background field	required
delta_phi	Array	Perturbation \(\delta\phi\)	required
m	float	Mass	required
lambda_	float	Coupling	required
dx	float	Lattice spacing	required

Returns:

Type	Description
Array	Linear response \(\delta\phi\)

Source code in jaxlatt/observables/stability.py

@jit
def compute_linear_response(
    phi_background: Array,
    delta_phi: Array,
    m: float,
    lambda_: float,
    dx: float,
) -> Array:
    r"""
    Compute linear response: how system responds to perturbation.

    $$\delta\phi \to \left.\frac{\partial^2 V}{\partial\phi^2}\right|_\mathrm{bg} \delta\phi$$

    This is the linearized equation of motion.

    Args:
        phi_background: Background field
        delta_phi: Perturbation $\delta\phi$
        m: Mass
        lambda_: Coupling
        dx: Lattice spacing

    Returns:
        Linear response $\delta\phi$
    """
    # Effective mass at background
    m_eff_sq = compute_mass_matrix(phi_background, m, lambda_)

    # Laplacian of perturbation
    laplacian_delta = (
        jnp.roll(delta_phi, -1, axis=0)
        + jnp.roll(delta_phi, 1, axis=0)
        + jnp.roll(delta_phi, -1, axis=1)
        + jnp.roll(delta_phi, 1, axis=1)
        + jnp.roll(delta_phi, -1, axis=2)
        + jnp.roll(delta_phi, 1, axis=2)
        - 6 * delta_phi
    ) / (dx**2)

    # Linear response: (-nabla^2 + m_eff^2) delta_phi
    return -laplacian_delta + m_eff_sq * delta_phi

compute_curvature_at_origin(m, lambda_)

Compute \(\partial^2 V/\partial\phi^2|_{\phi=0} = m^2\).

• If \(m^2 < 0\): symmetry broken (\(\phi = 0\) unstable)
• If \(m^2 > 0\): symmetric phase (\(\phi = 0\) stable)

Parameters:

Name	Type	Description	Default
m	float	Mass parameter	required
lambda_	float	Coupling	required

Returns:

Type	Description
float	Curvature at origin (\(= m^2\))

Source code in jaxlatt/observables/stability.py

@jit
def compute_curvature_at_origin(
    m: float,
    lambda_: float,
) -> float:
    r"""
    Compute $\partial^2 V/\partial\phi^2|_{\phi=0} = m^2$.

    - If $m^2 < 0$: symmetry broken ($\phi = 0$ unstable)
    - If $m^2 > 0$: symmetric phase ($\phi = 0$ stable)

    Args:
        m: Mass parameter
        lambda_: Coupling

    Returns:
        Curvature at origin ($= m^2$)
    """
    return m * abs(m)

compute_vev_estimate(m, lambda_)

Estimate vacuum expectation value for symmetry breaking.

• If \(m^2 < 0\): VEV \(\approx \sqrt{-m^2/\lambda}\)
• If \(m^2 > 0\): VEV \(= 0\)

Parameters:

Name	Type	Description	Default
m	float	Mass (can be negative for SSB)	required
lambda_	float	Self-coupling	required

Returns:

Type	Description
float	Estimated \(|\langle\phi\rangle|\)

Source code in jaxlatt/observables/stability.py

def compute_vev_estimate(
    m: float,
    lambda_: float,
) -> float:
    r"""
    Estimate vacuum expectation value for symmetry breaking.

    - If $m^2 < 0$: VEV $\approx \sqrt{-m^2/\lambda}$
    - If $m^2 > 0$: VEV $= 0$

    Args:
        m: Mass (can be negative for SSB)
        lambda_: Self-coupling

    Returns:
        Estimated $|\langle\phi\rangle|$
    """
    if m < 0 and lambda_ > 0:
        return float(jnp.sqrt(m**2 / lambda_))
    else:
        return 0.0

check_symmetry_breaking(phi, m, lambda_)

Check if field has non-zero VEV (symmetry breaking).

Parameters:

Name	Type	Description	Default
phi	Array	Scalar field configuration.	required
m	float	Mass parameter.	required
lambda_	float	Self-coupling.	required

Returns:

Type	Description
dict[str, float]	Dictionary with measured field statistics, expected VEV, curvature at the
dict[str, float]	origin, and a symmetry-breaking indicator.

Source code in jaxlatt/observables/stability.py

def check_symmetry_breaking(
    phi: Array,
    m: float,
    lambda_: float,
) -> dict[str, float]:
    """
    Check if field has non-zero VEV (symmetry breaking).

    Args:
        phi: Scalar field configuration.
        m: Mass parameter.
        lambda_: Self-coupling.

    Returns:
        Dictionary with measured field statistics, expected VEV, curvature at the
        origin, and a symmetry-breaking indicator.
    """
    # Mean field
    phi_mean = jnp.mean(phi)
    phi_rms = jnp.sqrt(jnp.mean(jnp.abs(phi) ** 2))

    # Expected VEV
    vev_expected = compute_vev_estimate(m, lambda_)

    # Curvature
    curvature = compute_curvature_at_origin(m, lambda_)

    return {
        "mean_field_magnitude": float(jnp.abs(phi_mean)),
        "rms_field": float(phi_rms),
        "expected_vev": vev_expected,
        "curvature_at_origin": curvature,
        "is_symmetry_broken": bool(curvature < 0),
        "vev_ratio": float(phi_rms / (vev_expected + 1e-10)),
    }

estimate_bounce_action(phi_false, phi_true, m, lambda_)

Estimate bounce action for tunneling \(\phi_\mathrm{false} \to \phi_\mathrm{true}\).

\[S_\mathrm{bounce} \sim (\Delta V)^{-2}\,(\text{barrier height})^4\]

Tunneling rate: \(\Gamma \sim \exp(-S_\mathrm{bounce})\)

Parameters:

Name	Type	Description	Default
phi_false	float	False vacuum field value	required
phi_true	float	True vacuum field value	required
m	float	Mass	required
lambda_	float	Coupling	required

Returns:

Type	Description
float	Bounce action estimate

Source code in jaxlatt/observables/stability.py

def estimate_bounce_action(
    phi_false: float,
    phi_true: float,
    m: float,
    lambda_: float,
) -> float:
    r"""
    Estimate bounce action for tunneling $\phi_\mathrm{false} \to \phi_\mathrm{true}$.

    $$S_\mathrm{bounce} \sim (\Delta V)^{-2}\,(\text{barrier height})^4$$

    Tunneling rate: $\Gamma \sim \exp(-S_\mathrm{bounce})$

    Args:
        phi_false: False vacuum field value
        phi_true: True vacuum field value
        m: Mass
        lambda_: Coupling

    Returns:
        Bounce action estimate
    """
    # Potential at false and true vacua
    V_false = scalar_potential_energy(phi_false, m, lambda_)
    V_true = scalar_potential_energy(phi_true, m, lambda_)

    # Potential difference
    Delta_V = jnp.abs(V_true - V_false)

    # Barrier estimate (crude approximation)
    # For quartic potential with symmetry breaking
    if m < 0:
        # Barrier at φ = 0
        V_barrier = 0.0
        barrier_height = jnp.abs(V_barrier - V_false)
    else:
        # No barrier
        barrier_height = 0.0

    # Bounce action ~ barrier^4 / ΔV²
    if float(Delta_V) > 1e-10:
        S_bounce = barrier_height**4 / (Delta_V**2 + 1e-10)
    else:
        S_bounce = 1e10  # Very suppressed

    return float(S_bounce)

generate_stability_report(phi, m, lambda_, dx)

Generate comprehensive stability analysis report.

Parameters:

Name	Type	Description	Default
phi	Array	Scalar field configuration.	required
m	float	Mass parameter.	required
lambda_	float	Self-coupling.	required
dx	float	Lattice spacing.	required

Returns:

Type	Description
dict[str, float]	Dictionary combining mass-matrix statistics, tachyonic-mode diagnostics,
dict[str, float]	and symmetry-breaking indicators.

Source code in jaxlatt/observables/stability.py

def generate_stability_report(
    phi: Array,
    m: float,
    lambda_: float,
    dx: float,
) -> dict[str, float]:
    """
    Generate comprehensive stability analysis report.

    Args:
        phi: Scalar field configuration.
        m: Mass parameter.
        lambda_: Self-coupling.
        dx: Lattice spacing.

    Returns:
        Dictionary combining mass-matrix statistics, tachyonic-mode diagnostics,
        and symmetry-breaking indicators.
    """
    # Mass matrix
    m_eff_sq = compute_mass_matrix(phi, m, lambda_)

    # Tachyonic modes
    tach_info = detect_tachyonic_modes(phi, m, lambda_)

    # Symmetry breaking
    ssb_info = check_symmetry_breaking(phi, m, lambda_)

    # Combined report
    report = {
        "min_mass_sq": float(jnp.min(m_eff_sq)),
        "max_mass_sq": float(jnp.max(m_eff_sq)),
        "mean_mass_sq": float(jnp.mean(m_eff_sq)),
        "has_tachyonic_modes": tach_info["has_instability"],
        "fraction_unstable": tach_info["fraction_tachyonic"],
        "field_rms": ssb_info["rms_field"],
        "symmetry_broken": ssb_info["is_symmetry_broken"],
        "curvature_at_origin": ssb_info["curvature_at_origin"],
    }

    return report