"""Generate deterministic Lecture 27 finite-size and collapse records. The calculation uses the exact even-parity Ising-ring pair gap Delta_+(L, gamma) / J = 4 sqrt(1 + gamma^2 - 2 gamma cos(pi/L)). Only the Python standard library is required. The quartic common-curve collapse score is a descriptive residual for deterministic data, not a statistical covariance estimate. """ from __future__ import annotations import json import math from pathlib import Path LENGTHS = [4, 6, 8, 10, 12] REFERENCE_X = [-2.0, -1.0, 0.0, 1.0, 2.0] CANDIDATES = [ {"name": "Ising", "gamma_c": 1.0, "nu": 1.0, "z": 1.0}, {"name": "shifted-critical-point", "gamma_c": 0.98, "nu": 1.0, "z": 1.0}, {"name": "changed-exponents", "gamma_c": 1.0, "nu": 0.8, "z": 1.1}, ] def gap_over_j(length: int, gamma: float) -> float: return 4.0 * math.sqrt( 1.0 + gamma * gamma - 2.0 * gamma * math.cos(math.pi / length) ) def linear_fit(xs: list[float], ys: list[float]) -> tuple[float, float]: x_bar = sum(xs) / len(xs) y_bar = sum(ys) / len(ys) denominator = sum((x - x_bar) ** 2 for x in xs) slope = sum((x - x_bar) * (y - y_bar) for x, y in zip(xs, ys)) / denominator return y_bar - slope * x_bar, slope def solve_linear(matrix: list[list[float]], vector: list[float]) -> list[float]: size = len(vector) augmented = [row[:] + [value] for row, value in zip(matrix, vector)] for column in range(size): pivot = max(range(column, size), key=lambda row: abs(augmented[row][column])) if abs(augmented[pivot][column]) < 1.0e-14: raise ValueError("singular normal-equation matrix") augmented[column], augmented[pivot] = augmented[pivot], augmented[column] scale = augmented[column][column] augmented[column] = [value / scale for value in augmented[column]] for row in range(size): if row == column: continue factor = augmented[row][column] augmented[row] = [ left - factor * right for left, right in zip(augmented[row], augmented[column]) ] return [augmented[row][-1] for row in range(size)] def polynomial_fit(xs: list[float], ys: list[float], degree: int) -> list[float]: matrix = [ [sum(x ** (row + column) for x in xs) for column in range(degree + 1)] for row in range(degree + 1) ] vector = [sum(y * x**row for x, y in zip(xs, ys)) for row in range(degree + 1)] return solve_linear(matrix, vector) def polynomial_value(coefficients: list[float], x: float) -> float: return sum(coefficient * x**power for power, coefficient in enumerate(coefficients)) def main() -> None: output = Path("artifacts/validation/lec29_scaling/scaling.json") output.parent.mkdir(parents=True, exist_ok=True) records = [] for length in LENGTHS: row = [] for reference_x in REFERENCE_X: gamma = 1.0 + reference_x / length gap = gap_over_j(length, gamma) row.append( { "reference_x": reference_x, "gamma": gamma, "gap_over_J": gap, "L_gap_over_J": length * gap, } ) records.append({"L": length, "points": row}) critical_fits = [] for minimum_length in [4, 6, 8]: selected = [length for length in LENGTHS if length >= minimum_length] log_lengths = [math.log(length) for length in selected] log_gaps = [math.log(gap_over_j(length, 1.0)) for length in selected] intercept, slope = linear_fit(log_lengths, log_gaps) critical_fits.append( { "minimum_L": minimum_length, "included_L": selected, "intercept": intercept, "amplitude": math.exp(intercept), "slope": slope, "z_estimate": -slope, } ) exact_limit_residuals = [] for minimum_length in [4, 6, 8]: squared = [] for record in records: if record["L"] < minimum_length: continue for point in record["points"]: scaling_limit = 4.0 * math.sqrt( point["reference_x"] ** 2 + math.pi**2 ) squared.append((point["L_gap_over_J"] - scaling_limit) ** 2) exact_limit_residuals.append( { "minimum_L": minimum_length, "rms_to_exact_scaling_limit": math.sqrt(sum(squared) / len(squared)), } ) candidate_scores = [] for minimum_length in [4, 6, 8]: window = {"minimum_L": minimum_length, "candidates": []} for candidate in CANDIDATES: rescaled_x = [] rescaled_y = [] for record in records: length = record["L"] if length < minimum_length: continue for point in record["points"]: reduced = ( point["gamma"] - candidate["gamma_c"] ) / candidate["gamma_c"] rescaled_x.append(reduced * length ** (1.0 / candidate["nu"])) rescaled_y.append( length ** candidate["z"] * point["gap_over_J"] ) coefficients = polynomial_fit(rescaled_x, rescaled_y, degree=4) residuals = [ y - polynomial_value(coefficients, x) for x, y in zip(rescaled_x, rescaled_y) ] window["candidates"].append( { **candidate, "quartic_coefficients_low_to_high": coefficients, "rms_common_curve_residual": math.sqrt( sum(value * value for value in residuals) / len(residuals) ), "point_count": len(residuals), } ) candidate_scores.append(window) thermal_cases = [] for temperature_over_j in [0.05, 0.20]: gap = 0.20 thermal_cases.append( { "gamma": 0.9, "bulk_gap_over_J": gap, "kBT_over_J": temperature_over_j, "J_tau_T_over_hbar": 1.0 / temperature_over_j, "J_xi_tau_over_hbar": 1.0 / gap, "fermion_occupation": 1.0 / (math.exp(gap / temperature_over_j) + 1.0), "fan_edge_gamma": [ 1.0 - temperature_over_j / 2.0, 1.0 + temperature_over_j / 2.0, ], } ) payload = { "format_version": 1, "model": "periodic transverse-field Ising ring, even-parity pair gap", "formula": "Delta_+/J = 4 sqrt(1 + gamma^2 - 2 gamma cos(pi/L))", "lengths": LENGTHS, "reference_x": REFERENCE_X, "records": records, "critical_log_fits": critical_fits, "exact_limit_residuals": exact_limit_residuals, "collapse_metric": { "description": ( "ordinary least-squares quartic common curve in each candidate's " "rescaled coordinates; RMS is descriptive for exact data" ), "candidate_scores": candidate_scores, "statistical_boundary": ( "No sampling noise is present, so these scores do not define a " "parameter covariance matrix or confidence interval." ), }, "thermal_cases": thermal_cases, } output.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8") print(output) if __name__ == "__main__": main()