#!/usr/bin/env python3 """Calibrate the quartic double well and its one-instanton splitting. The convention used on the companion page is H(g) = -(g/2) d^2/dq^2 + q^2 (1-q)^2/(2g), g > 0. In the centered coordinate x=q-1/2 this is H(g) = -(g/2) d^2/dx^2 + (x^2-1/4)^2/(2g). The program first generates the local ground-state perturbation series by an exact rational Bender--Wu recurrence. Its factorial-stripped ratios estimate the neutral return action 2*S_I=1/3, not the parity-changing one-instanton action S_I=1/6. It then computes the lowest even--odd doublet in two independent, parity-resolved representations: 1. a centered harmonic-oscillator Galerkin basis, refined by increasing the nested basis cutoff; 2. a uniform coordinate grid, reduced into even and odd Jacobi matrices, whose lowest eigenvalues are found by Sturm bisection and Richardson extrapolated after step halving. For the ground doublet, the one-instanton prediction is Delta E_0 ~ 2 exp[-1/(6g)] / sqrt(pi*g) * (1 - 71 g/12 - 6299 g^2/288 + ...). The calibrated double-precision interval is 0.0125 <= g <= 0.05. Refinement shifts and cross-representation differences are empirical diagnostics, not certified enclosures. Requirements: CPython 3.9+ and NumPy 1.24+. """ from __future__ import annotations import argparse import csv from dataclasses import dataclass from fractions import Fraction import math from pathlib import Path import platform from typing import Iterable import numpy as np PI = math.pi INSTANTON_ACTION = 1.0 / 6.0 MIN_CALIBRATED_G = 0.0125 MAX_CALIBRATED_G = 0.05 FIGURE_G_VALUES = ( 0.0125, 0.015, 0.0175, 0.02, 0.025, 0.03, 0.04, 0.05, ) def require(condition: bool, message: str) -> None: """Raise an optimization-safe exception when an audit fails.""" if not condition: raise RuntimeError(message) def potential(x: np.ndarray | float, g: float) -> np.ndarray | float: """Return the centered quartic potential.""" return (np.asarray(x) ** 2 - 0.25) ** 2 / (2.0 * g) @dataclass(frozen=True) class ParityPair: """Lowest even and odd energies in one numerical representation.""" even: float odd: float @property def splitting(self) -> float: return self.odd - self.even @dataclass(frozen=True) class GalerkinResult: """Refined oscillator-basis result and cutoff diagnostics.""" pair: ParityPair coarse_pair: ParityPair energy_shift: float splitting_shift: float @dataclass(frozen=True) class GridResult: """Richardson-refined grid result and step diagnostics.""" pair: ParityPair fine_pair: ParityPair energy_shift: float splitting_shift: float sturm_bracket: float @dataclass(frozen=True) class CalibrationRow: """All independently computed data at one coupling.""" g: float galerkin: GalerkinResult grid: GridResult @property def cross_energy_residual(self) -> float: return max( abs(self.galerkin.pair.even - self.grid.pair.even), abs(self.galerkin.pair.odd - self.grid.pair.odd), ) @property def cross_splitting_relative(self) -> float: scale = max( self.galerkin.pair.splitting, self.grid.pair.splitting, np.finfo(float).tiny, ) return abs( self.galerkin.pair.splitting - self.grid.pair.splitting ) / scale @dataclass(frozen=True) class Settings: """Numerical cutoffs for one run mode.""" galerkin_coarse: int galerkin_fine: int grid_coarse_cells: int domain: float perturbative_order: int def oscillator_hamiltonian(g: float, size: int) -> np.ndarray: """Return the centered harmonic-oscillator Galerkin matrix.""" require(size >= 20 and size % 2 == 0, "Galerkin size must be even and >= 20") length = math.sqrt(g) work_size = size + 4 coordinate = np.zeros((work_size, work_size)) for n in range(work_size - 1): value = math.sqrt((n + 1.0) / 2.0) coordinate[n, n + 1] = value coordinate[n + 1, n] = value coordinate_squared = coordinate @ coordinate coordinate_fourth = coordinate_squared @ coordinate_squared minus_second_derivative = np.diag( np.arange(work_size, dtype=float) + 0.5 ) for n in range(work_size - 2): value = -0.5 * math.sqrt((n + 1.0) * (n + 2.0)) minus_second_derivative[n, n + 2] = value minus_second_derivative[n + 2, n] = value identity = np.eye(work_size) kinetic = ( g * minus_second_derivative / (2.0 * length * length) ) quartic = ( length**4 * coordinate_fourth - 0.5 * length**2 * coordinate_squared + 0.0625 * identity ) / (2.0 * g) matrix = kinetic + quartic return matrix[:size, :size] def galerkin_pair(g: float, size: int) -> ParityPair: """Diagonalize the even and odd oscillator-basis blocks.""" matrix = oscillator_hamiltonian(g, size) even_indices = np.arange(0, size, 2) odd_indices = np.arange(1, size, 2) even_block = matrix[np.ix_(even_indices, even_indices)] odd_block = matrix[np.ix_(odd_indices, odd_indices)] even = float(np.linalg.eigvalsh(even_block)[0]) odd = float(np.linalg.eigvalsh(odd_block)[0]) return ParityPair(even, odd) def refined_galerkin( g: float, coarse_size: int, fine_size: int, ) -> GalerkinResult: """Refine nested Rayleigh--Ritz cutoffs.""" require(fine_size > coarse_size, "fine Galerkin size must be larger") coarse = galerkin_pair(g, coarse_size) fine = galerkin_pair(g, fine_size) energy_shift = max( abs(fine.even - coarse.even), abs(fine.odd - coarse.odd), ) splitting_shift = abs(fine.splitting - coarse.splitting) return GalerkinResult(fine, coarse, energy_shift, splitting_shift) def grid_parity_blocks( g: float, domain: float, cells: int, ) -> tuple[ tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray], ]: """Return even and odd tridiagonal blocks of the centered grid.""" require(domain > 1.0, "coordinate domain must contain both wells") require(cells >= 200, "coordinate grid is too coarse") step = domain / cells kinetic_diagonal = g / step**2 kinetic_off_diagonal = -g / (2.0 * step**2) even_x = step * np.arange(cells, dtype=float) even_diagonal = kinetic_diagonal + potential(even_x, g) even_off = np.full(cells - 1, kinetic_off_diagonal) # This sqrt(2) is the orthonormal parity reduction of the full-line # coupling between x=0 and the symmetric pair x=+-h. even_off[0] *= math.sqrt(2.0) odd_x = step * np.arange(1, cells, dtype=float) odd_diagonal = kinetic_diagonal + potential(odd_x, g) odd_off = np.full(cells - 2, kinetic_off_diagonal) return ( (np.asarray(even_diagonal), even_off), (np.asarray(odd_diagonal), odd_off), ) def sturm_count( diagonal: np.ndarray, off_diagonal: np.ndarray, value: float, ) -> int: """Count eigenvalues strictly below value by an LDL Sturm sequence.""" scale = max(1.0, float(np.max(np.abs(diagonal)))) pivot_floor = 8.0 * np.finfo(float).eps * scale pivot = float(diagonal[0] - value) if abs(pivot) < pivot_floor: pivot = -pivot_floor count = int(pivot < 0.0) for index in range(1, diagonal.size): pivot = float( diagonal[index] - value - off_diagonal[index - 1] ** 2 / pivot ) if abs(pivot) < pivot_floor: pivot = -pivot_floor count += int(pivot < 0.0) return count def lowest_tridiagonal_eigenvalue( diagonal: np.ndarray, off_diagonal: np.ndarray, ) -> tuple[float, float]: """Return the lowest Jacobi eigenvalue and final bisection width.""" require( diagonal.size == off_diagonal.size + 1, "invalid tridiagonal dimensions", ) radii = np.zeros_like(diagonal) radii[:-1] += np.abs(off_diagonal) radii[1:] += np.abs(off_diagonal) scale = max(1.0, float(np.max(np.abs(diagonal)))) lower = float(np.min(diagonal - radii) - 8.0 * np.finfo(float).eps * scale) upper = float(np.max(diagonal + radii) + 8.0 * np.finfo(float).eps * scale) for _ in range(96): midpoint = 0.5 * (lower + upper) if midpoint == lower or midpoint == upper: break if sturm_count(diagonal, off_diagonal, midpoint) == 0: lower = midpoint else: upper = midpoint return 0.5 * (lower + upper), upper - lower def grid_pair(g: float, domain: float, cells: int) -> tuple[ParityPair, float]: """Compute the parity pair on one uniform full-line grid.""" even_block, odd_block = grid_parity_blocks(g, domain, cells) even, even_width = lowest_tridiagonal_eigenvalue(*even_block) odd, odd_width = lowest_tridiagonal_eigenvalue(*odd_block) return ParityPair(even, odd), max(even_width, odd_width) def refined_grid( g: float, domain: float, coarse_cells: int, ) -> GridResult: """Cancel the leading second-order grid error by step halving.""" coarse, coarse_bracket = grid_pair(g, domain, coarse_cells) fine, fine_bracket = grid_pair(g, domain, 2 * coarse_cells) extrapolated = ParityPair( (4.0 * fine.even - coarse.even) / 3.0, (4.0 * fine.odd - coarse.odd) / 3.0, ) energy_shift = max( abs(extrapolated.even - fine.even), abs(extrapolated.odd - fine.odd), ) splitting_shift = abs(extrapolated.splitting - fine.splitting) return GridResult( extrapolated, fine, energy_shift, splitting_shift, max(coarse_bracket, fine_bracket), ) def leading_splitting(g: float) -> float: """Return the leading ground-doublet one-instanton splitting.""" return ( 2.0 * math.exp(-INSTANTON_ACTION / g) / math.sqrt(PI * g) ) def scaled_instanton_ratio(splitting: float, g: float) -> float: """Remove the leading instanton exponential and determinant.""" return splitting / leading_splitting(g) def fluctuation_ratio(g: float, order: int) -> float: """Return the one-instanton fluctuation factor through order g^order.""" require(order in (0, 1, 2), "implemented fluctuation order is 0, 1, or 2") value = 1.0 if order >= 1: value -= 71.0 * g / 12.0 if order >= 2: value -= 6299.0 * g * g / 288.0 return value def local_ground_perturbation(order: int) -> list[Fraction]: """Generate exact coefficients of the common local ground series. Set q=sqrt(g)*y near the left minimum and lambda=sqrt(g). With psi = exp(-y^2/2) sum_n lambda^n P_n(y), E = sum_n lambda^n epsilon_n, and P_n(0)=0 for n>0, coefficient matching in H = H_0 - lambda*y^3 + lambda^2*y^4/2 gives a descending polynomial recurrence. Reflection forces every odd epsilon_n to vanish exactly. The returned entry k is epsilon_(2k), i.e. the coefficient of g^k. """ require(order >= 4, "perturbative order must be at least four") polynomials: list[list[Fraction]] = [[Fraction(1)]] lambda_energies: list[Fraction] = [Fraction(1, 2)] for n in range(1, 2 * order + 1): maximum_degree = 3 * n coefficients = [ Fraction(0) for _ in range(maximum_degree + 3) ] for m in range(maximum_degree, 0, -1): source = Fraction(0) previous = polynomials[n - 1] if m >= 3 and m - 3 < len(previous): source += previous[m - 3] if n >= 2: previous_two = polynomials[n - 2] if m >= 4 and m - 4 < len(previous_two): source -= Fraction(1, 2) * previous_two[m - 4] for j in range(1, n): lower_polynomial = polynomials[n - j] if m < len(lower_polynomial): source += lambda_energies[j] * lower_polynomial[m] coefficients[m] = ( source + Fraction((m + 2) * (m + 1), 2) * coefficients[m + 2] ) / m # The constant term of # (-P_n''/2 + y P_n') = sum_j epsilon_j P_(n-j) + ... # is -P_{n,2}=epsilon_n under the intermediate normalization. lambda_energies.append(-coefficients[2]) polynomials.append(coefficients[: maximum_degree + 1]) require( all(value == 0 for value in lambda_energies[1::2]), "reflection audit failed in the perturbative recurrence", ) return [lambda_energies[2 * k] for k in range(order + 1)] def neutral_action_estimator( coefficients: list[Fraction], order: int, ) -> Fraction: """Return A_k=k E_(k-1)/E_k for the local perturbative series.""" require(1 <= order < len(coefficients), "invalid estimator order") require(coefficients[order] != 0, "zero perturbative denominator") return Fraction(order) * coefficients[order - 1] / coefficients[order] def audit_perturbation(coefficients: list[Fraction]) -> None: """Check low orders and the independently inferred neutral action.""" expected = [ Fraction(1, 2), Fraction(-1), Fraction(-9, 2), Fraction(-89, 2), Fraction(-5013, 8), ] require( coefficients[:5] == expected, "local perturbative recurrence audit failed", ) estimate = float( neutral_action_estimator(coefficients, len(coefficients) - 1) ) require( 0.0 < estimate < 1.0 / 3.0 and abs(estimate - 1.0 / 3.0) < 5.0e-3, "neutral-action large-order audit failed", ) def calibrate(g: float, settings: Settings) -> CalibrationRow: """Compute both independently refined representations.""" galerkin = refined_galerkin( g, settings.galerkin_coarse, settings.galerkin_fine, ) grid = refined_grid(g, settings.domain, settings.grid_coarse_cells) return CalibrationRow(g, galerkin, grid) def audit_rows(rows: Iterable[CalibrationRow]) -> None: """Reject failed parity, refinement, or cross-method diagnostics.""" for row in rows: require(row.galerkin.pair.splitting > 0.0, "Galerkin parity order failed") require(row.grid.pair.splitting > 0.0, "grid parity order failed") require( row.galerkin.energy_shift < 2.0e-10, "Galerkin cutoff audit failed", ) require( row.grid.energy_shift < 4.0e-6, "coordinate-grid step audit failed", ) require( row.cross_energy_residual < 2.0e-6, "Galerkin/grid energy audit failed", ) require( row.cross_splitting_relative < 1.0e-2, "Galerkin/grid splitting audit failed", ) CSV_FIELDS = ( "g", "E_even_galerkin", "E_odd_galerkin", "splitting_galerkin", "E_even_grid", "E_odd_grid", "splitting_grid", "scaled_ratio_galerkin", "scaled_ratio_grid", "ratio_leading", "ratio_through_g", "ratio_through_g2", "galerkin_energy_shift", "galerkin_splitting_shift", "grid_energy_shift", "grid_splitting_shift", "cross_energy_residual", "cross_splitting_relative", ) def row_values(row: CalibrationRow) -> tuple[float, ...]: """Serialize one calibration row in the declared CSV order.""" g = row.g return ( g, row.galerkin.pair.even, row.galerkin.pair.odd, row.galerkin.pair.splitting, row.grid.pair.even, row.grid.pair.odd, row.grid.pair.splitting, scaled_instanton_ratio(row.galerkin.pair.splitting, g), scaled_instanton_ratio(row.grid.pair.splitting, g), fluctuation_ratio(g, 0), fluctuation_ratio(g, 1), fluctuation_ratio(g, 2), row.galerkin.energy_shift, row.galerkin.splitting_shift, row.grid.energy_shift, row.grid.splitting_shift, row.cross_energy_residual, row.cross_splitting_relative, ) def write_rows(path: Path, rows: Iterable[CalibrationRow]) -> None: """Write the complete numerical audit table.""" with path.open("w", newline="", encoding="utf-8") as handle: writer = csv.writer(handle) writer.writerow(CSV_FIELDS) for row in rows: writer.writerow(row_values(row)) def turning_points(scaled_energy: float) -> tuple[float, float, float, float]: """Return the four real roots of q^2(1-q)^2=scaled_energy.""" require(0.0 < scaled_energy < 1.0 / 16.0, "energy must lie below barrier") root_energy = math.sqrt(scaled_energy) inner_root = math.sqrt(1.0 - 4.0 * root_energy) outer_root = math.sqrt(1.0 + 4.0 * root_energy) return ( 0.5 * (1.0 - outer_root), 0.5 * (1.0 - inner_root), 0.5 * (1.0 + inner_root), 0.5 * (1.0 + outer_root), ) def write_figure_data( prefix: Path, rows: list[CalibrationRow], ) -> tuple[Path, Path, Path]: """Export every numerical input used by the standalone page figure.""" splitting_path = Path(f"{prefix}-splitting.csv") potential_path = Path(f"{prefix}-potential.csv") turning_points_path = Path(f"{prefix}-turning-points.csv") write_rows(splitting_path, rows) q_values = np.linspace(-0.25, 1.25, 121) with potential_path.open("w", newline="", encoding="utf-8") as handle: writer = csv.writer(handle) writer.writerow(["q", "U"]) for q in q_values: writer.writerow( [float(q), float(0.5 * q * q * (1.0 - q) ** 2)] ) well_energy = 0.01 roots = turning_points(2.0 * well_energy) with turning_points_path.open("w", newline="", encoding="utf-8") as handle: writer = csv.writer(handle) writer.writerow(["name", "q", "gE"]) for index, root in enumerate(roots, start=1): writer.writerow([f"q_{index}", root, well_energy]) return splitting_path, potential_path, turning_points_path def print_report( benchmark: CalibrationRow, rows: list[CalibrationRow], settings: Settings, perturbative_coefficients: list[Fraction], ) -> None: """Print the benchmark, refinement audits, and instanton scaling.""" print("Symmetric quartic double-well calibration") print("H = -(g/2) d^2/dq^2 + q^2(1-q)^2/(2g)") print( f"Python {platform.python_version()}; NumPy {np.__version__}; " f"calibrated {MIN_CALIBRATED_G:g} <= g <= {MAX_CALIBRATED_G:g}" ) print( "cutoffs: Galerkin " f"{settings.galerkin_coarse}->{settings.galerkin_fine}; " f"grid {settings.grid_coarse_cells}->{2 * settings.grid_coarse_cells} " f"cells on [0,{settings.domain:g}]" ) print() print("Exact local perturbation recurrence") print("E_0^(0)(g) = sum_k e_k g^k (common to both parities)") for k, coefficient in enumerate(perturbative_coefficients[:5]): print(f" e_{k} = {coefficient}") print("Neutral two-event action estimator A_k = k e_(k-1)/e_k") estimator_orders = [8, 12, 16, 20] estimator_orders.extend( order for order in (24, 28, 32) if order < len(perturbative_coefficients) ) print(f"{'k':>8s} {'A_k':>16s} {'1/3 - A_k':>16s}") for order in estimator_orders: estimate = float( neutral_action_estimator(perturbative_coefficients, order) ) print( f"{order:8d} {estimate:16.12f} " f"{1.0 / 3.0 - estimate:16.3e}" ) print( "The limit is 1/3=2*S_I; local perturbation theory cannot see " "the parity-changing S_I=1/6 sector directly." ) print() print(f"Detailed benchmark at g={benchmark.g:g}") print(f"{'method':>12s} {'E_even':>18s} {'E_odd':>18s} {'Delta E':>16s}") print( f"{'Galerkin':>12s} {benchmark.galerkin.pair.even:18.12f} " f"{benchmark.galerkin.pair.odd:18.12f} " f"{benchmark.galerkin.pair.splitting:16.9e}" ) print( f"{'grid/Sturm':>12s} {benchmark.grid.pair.even:18.12f} " f"{benchmark.grid.pair.odd:18.12f} " f"{benchmark.grid.pair.splitting:16.9e}" ) print() print("Benchmark diagnostics") print( f" Galerkin maximum cutoff shift: " f"{benchmark.galerkin.energy_shift:.3e}" ) print( f" Galerkin splitting shift: " f"{benchmark.galerkin.splitting_shift:.3e}" ) print( f" grid maximum Richardson shift: " f"{benchmark.grid.energy_shift:.3e}" ) print( f" grid splitting Richardson shift: " f"{benchmark.grid.splitting_shift:.3e}" ) print( f" cross-method energy residual: " f"{benchmark.cross_energy_residual:.3e}" ) print( f" cross-method splitting residual: " f"{benchmark.cross_splitting_relative:.3e}" ) print( f" Sturm bisection stagnation width (not an error bound): " f"{benchmark.grid.sturm_bracket:.3e}" ) print() print("Ground-doublet instanton scaling") print( f"{'g':>8s} {'Delta E':>14s} {'scaled R':>12s} " f"{'1-71g/12':>12s} {'through g^2':>12s} {'cross rel.':>12s}" ) for row in rows: print( f"{row.g:8.4f} {row.galerkin.pair.splitting:14.7e} " f"{scaled_instanton_ratio(row.galerkin.pair.splitting, row.g):12.8f} " f"{fluctuation_ratio(row.g, 1):12.8f} " f"{fluctuation_ratio(row.g, 2):12.8f} " f"{row.cross_splitting_relative:12.2e}" ) print() print("Caveat: all shifts are empirical floating-point diagnostics.") print( "The scaled exact splitting contains higher one-instanton " "fluctuations and exponentially smaller odd-instanton sectors." ) def calibrated_g(text: str) -> float: """Parse a coupling inside the double-precision calibration interval.""" try: value = float(text) except ValueError as error: raise argparse.ArgumentTypeError("g must be a real number") from error if not MIN_CALIBRATED_G <= value <= MAX_CALIBRATED_G: raise argparse.ArgumentTypeError( f"g must lie in {MIN_CALIBRATED_G:g} <= g <= " f"{MAX_CALIBRATED_G:g}; outside this interval use higher " "precision or retune both discretizations" ) return value def parse_arguments() -> argparse.Namespace: """Parse command-line options.""" parser = argparse.ArgumentParser(description=__doc__) parser.add_argument( "--g", type=calibrated_g, default=0.025, help=( "detailed benchmark coupling in the calibrated interval " "0.0125 <= g <= 0.05" ), ) parser.add_argument( "--mode", choices=("default", "high"), default="default", help=( "refinement profile: default, or high for larger Galerkin " "cutoffs, a finer coordinate grid, and more recurrence orders" ), ) parser.add_argument( "--csv", metavar="PATH", help="write the full coupling scan and diagnostics to PATH", ) parser.add_argument( "--figure-data", metavar="PREFIX", help=( "write PREFIX-splitting.csv, PREFIX-potential.csv, and " "PREFIX-turning-points.csv" ), ) return parser.parse_args() def main() -> None: """Run the two-representation double-well calibration.""" arguments = parse_arguments() settings = ( Settings(88, 120, 1200, 1.5, 32) if arguments.mode == "high" else Settings(64, 88, 700, 1.5, 20) ) perturbative_coefficients = local_ground_perturbation( settings.perturbative_order ) audit_perturbation(perturbative_coefficients) requested = sorted(set(FIGURE_G_VALUES) | {float(arguments.g)}) rows = [calibrate(g, settings) for g in requested] audit_rows(rows) row_by_g = {row.g: row for row in rows} benchmark = row_by_g[float(arguments.g)] figure_rows = [row_by_g[g] for g in FIGURE_G_VALUES] print_report( benchmark, figure_rows, settings, perturbative_coefficients, ) if arguments.csv: csv_path = Path(arguments.csv) write_rows(csv_path, rows) print(f"\nWrote {csv_path}") if arguments.figure_data: paths = write_figure_data(Path(arguments.figure_data), figure_rows) print("\nWrote figure data:") for path in paths: print(f" {path}") if __name__ == "__main__": main()