#!/usr/bin/env python3 """Fixed-source Mott scattering: spin sums, angular curves and acceptance. Run with --help. Offline, deterministic, import-safe; NumPy + standard library. Natural units hbar=c=1 and metric (+---). Inputs m and p have energy units; the signed rationalized Coulomb coupling is g=qQ/(4*pi), V(r)=g/r. Areas are in inverse energy squared. Dirac spinors obey u-dagger u=2E, ubar u=2m. This is the LEADING external-potential Born approximation for a fixed, infinitely heavy, spinless point source. It includes neither target recoil nor target external-state normalization, screening, finite-size structure, exact Coulomb phases, radiation, loop corrections or dynamical pair creation. The fixed-source S matrix has an energy delta, not a two-body delta-four. The sufficient small parameter |g|/beta is reported, not certified as a uniform error bound; long-range forward scattering needs separate care. Numerical 4-component spinor overlaps, an independent gamma-matrix trace and the scalar closed result are compared. Exact rational spin fixtures, rotated spin bases, oblique momenta, low-velocity and massless limits test conventions. The high-energy spin factor is evaluated as cos(theta/2)^2+(m/E)^2*sin(theta/2)^2 to avoid subtracting nearly equal numbers. Raw trace cancellation near exact backscatter is reported honestly through absolute, scaled and relative errors. Angular curves exclude theta=0. Integrals are only over finite acceptances. Grid refinement at fixed acceptance is distinct from reducing a forward cutoff: the UNSCREENED TOTAL cross section diverges. Numerical integration uses composite Simpson quadrature in log(sin(theta/2)^2), with an independent closed acceptance integral evaluated using Decimal arithmetic. Outputs: angular, spin-transition, limiting-case, quadrature-convergence and forward-cutoff CSVs, and JSON containing checks, assumptions, runtime versions and hashes. Existing files are never overwritten. Without --output-dir a new temporary directory is used. Fixed benchmark cases always run alongside the custom diagnostic parameters. No plots or repository helpers are required. Owner: /relativistic-qm/mott-scattering/. References: Bjorken & Drell, Relativistic Quantum Mechanics (1964); DeGrand, A One-Semester Course on Quantum Field Theory (2025), section 8.7. Tested environment: Python 3.12.14, NumPy 2.3.5. """ from __future__ import annotations import argparse import csv from decimal import Decimal, localcontext from fractions import Fraction import hashlib import io import json import math import os from pathlib import Path import platform import sys import tempfile import uuid import numpy as np I2 = np.eye(2, dtype=complex) I4 = np.eye(4, dtype=complex) ZERO2 = np.zeros((2, 2), complex) PAULI = (np.array([[0, 1], [1, 0]], complex), np.array([[0, -1j], [1j, 0]], complex), np.diag([1, -1]).astype(complex)) GAMMA = (np.diag([1, 1, -1, -1]).astype(complex),) + tuple( np.block([[ZERO2, s], [-s, ZERO2]]) for s in PAULI) class Checks: def __init__(self): self.rows = [] def condition(self, name, passed, **details): self.rows.append({"name": name, "passed": bool(passed), **details}) def close(self, name, actual, expected, tolerance=5e-12): absolute = float(np.linalg.norm(np.asarray(actual) - np.asarray(expected))) scale = max(1.0, float(np.linalg.norm(expected))) self.condition(name, math.isfinite(absolute) and absolute <= tolerance * scale, absolute_error=absolute, scaled_error=absolute / scale, tolerance=tolerance) @property def passed(self): return all(row["passed"] for row in self.rows) def sigma_dot(vector): return sum(component * matrix for component, matrix in zip(vector, PAULI)) def slash(energy, vector): return energy * GAMMA[0] - sum(component * matrix for component, matrix in zip(vector, GAMMA[1:])) def spinors(mass, vector): """Columns use Pauli seeds (1,0) and (0,1), not momentum-helicity labels.""" energy = math.hypot(mass, float(np.linalg.norm(vector))) return math.sqrt(energy + mass) * np.vstack((I2, sigma_dot(vector) / (energy + mass))) def spin_trace(mass, incoming, outgoing): """Independent completeness/trace computation, including initial average.""" energy_in = math.hypot(mass, float(np.linalg.norm(incoming))) energy_out = math.hypot(mass, float(np.linalg.norm(outgoing))) return 0.5 * np.trace((slash(energy_out, outgoing) + mass * I4) @ GAMMA[0] @ (slash(energy_in, incoming) + mass * I4) @ GAMMA[0]) def angle_parts(theta): if theta == math.pi: return 0.0, -1.0, 1.0, 0.0 return math.sin(theta), math.cos(theta), math.sin(theta / 2), math.cos(theta / 2) def outgoing_vector(momentum, theta, azimuth): sine, cosine, _, _ = angle_parts(theta) return momentum * np.array([sine * math.cos(azimuth), sine * math.sin(azimuth), cosine]) def spin_factor(mass, momentum, theta): _, _, sh, ch = angle_parts(theta) energy = math.hypot(mass, momentum) return ch * ch + (mass / energy) ** 2 * sh * sh def differential(mass, momentum, coupling, theta): """Area per solid angle; theta>0 and p>0 are required by validated inputs.""" _, _, sh, _ = angle_parts(theta) beta = momentum / math.hypot(mass, momentum) return (coupling / (2 * momentum * beta * sh * sh)) ** 2 * spin_factor(mass, momentum, theta) def acceptance_exact(mass, momentum, coupling, lower, upper): """Finite angular acceptance, azimuth integrated; never an unscreened total.""" with localcontext() as context: context.prec = 60 m, p, g = map(lambda x: Decimal(str(x)), (mass, momentum, coupling)) xa = Decimal(str(angle_parts(lower)[2] ** 2)) xb = Decimal(str(angle_parts(upper)[2] ** 2)) beta2 = p * p / (p * p + m * m) bracket = 1 / xa - 1 / xb - beta2 * (xb / xa).ln() return float(Decimal(str(math.pi)) * g * g * bracket / (p * p * beta2)) def acceptance_simpson(mass, momentum, coupling, lower, upper, intervals): xa = angle_parts(lower)[2] ** 2 xb = angle_parts(upper)[2] ** 2 z = np.linspace(math.log(xa), math.log(xb), intervals + 1) x = np.exp(z) theta = 2 * np.arcsin(np.sqrt(np.minimum(x, 1.0))) values = np.array([4 * math.pi * xx * differential(mass, momentum, coupling, float(tt)) for xx, tt in zip(x, theta)]) weights = np.ones(intervals + 1) weights[1:-1:2] = 4 weights[2:-1:2] = 2 return float((z[-1] - z[0]) / (3 * intervals) * np.dot(weights, values)) class RationalComplex: """Tiny exact arithmetic used only by fixed rational spin fixtures.""" def __init__(self, real=0, imag=0): self.real, self.imag = Fraction(real), Fraction(imag) @staticmethod def coerce(value): return value if isinstance(value, RationalComplex) else RationalComplex(value) def __add__(self, other): other = self.coerce(other) return RationalComplex(self.real + other.real, self.imag + other.imag) __radd__ = __add__ def __neg__(self): return RationalComplex(-self.real, -self.imag) def __sub__(self, other): return self + (-self.coerce(other)) def __mul__(self, other): other = self.coerce(other) return RationalComplex(self.real * other.real - self.imag * other.imag, self.real * other.imag + self.imag * other.real) __rmul__ = __mul__ def conjugate(self): return RationalComplex(self.real, -self.imag) def absolute_squared(self): return self.real ** 2 + self.imag ** 2 def exact_fixtures(checks): q = RationalComplex positions = [(0, 0, 4), (4, 0, 0), (0, 4, 0), (Fraction(12, 5), 0, Fraction(16, 5)), (0, Fraction(12, 5), Fraction(-16, 5)), (0, 0, -4)] def reduced_spinors(mass, energy, position): x, y, z = (Fraction(component, energy + mass) for component in position) return [[q(1), q(0)], [q(0), q(1)], [q(z), q(x, -y)], [q(x, y), q(-z)]] for mass, energy in ((3, 5), (0, 4)): for i, incoming in enumerate(positions): for j, outgoing in enumerate(positions): vin, vout = reduced_spinors(mass, energy, incoming), reduced_spinors(mass, energy, outgoing) w = [[(energy + mass) * sum(vout[k][a].conjugate() * vin[k][b] for k in range(4)) for b in range(2)] for a in range(2)] actual = sum(entry.absolute_squared() for row in w for entry in row) / 2 dot = sum(Fraction(a) * Fraction(b) for a, b in zip(incoming, outgoing)) expected = 2 * (energy * energy + dot + mass * mass) checks.condition(f"exact m={mass} directions={i},{j}", actual == expected, exact_spin_sum=str(actual), exact_expected=str(expected)) ww = [[sum(w[a][k] * w[b][k].conjugate() for k in range(2)) for b in range(2)] for a in range(2)] checks.condition(f"exact unpolarized m={mass} directions={i},{j}", all( ww[a][b].real == (expected if a == b else 0) and ww[a][b].imag == 0 for a in range(2) for b in range(2))) def spin_benchmarks(checks, mass, momentum, theta, azimuth, label, rotation=None): incoming = np.array([0.0, 0.0, momentum]) outgoing = outgoing_vector(momentum, theta, azimuth) if rotation is not None: incoming, outgoing = rotation @ incoming, rotation @ outgoing energy = math.hypot(mass, momentum) ui, uf = spinors(mass, incoming), spinors(mass, outgoing) overlap = uf.conj().T @ ui spin_sum = float(0.5 * np.sum(abs(overlap) ** 2)) traced = spin_trace(mass, incoming, outgoing) expected = 4 * energy * energy * spin_factor(mass, momentum, theta) unit = 4 * energy * energy checks.close(label + " explicit sum / 4E^2", spin_sum / unit, expected / unit) checks.close(label + " trace / 4E^2", traced / unit, expected / unit) checks.close(label + " u dagger u / 2E", ui.conj().T @ ui / (2 * energy), I2) checks.close(label + " Dirac equation", (slash(energy, incoming) - mass * I4) @ ui / (energy ** 1.5), np.zeros((4, 2))) checks.close(label + " unpolarized output", overlap @ overlap.conj().T / unit, expected / unit * I2) basis_in = (I2 + 1j * PAULI[0]) / math.sqrt(2) basis_out = math.cos(0.37) * I2 - 1j * math.sin(0.37) * PAULI[1] rotated_overlap = basis_out.conj().T @ overlap @ basis_in checks.close(label + " spin basis independence", 0.5 * np.sum(abs(rotated_overlap) ** 2) / unit, expected / unit) return {"spin_sum": spin_sum, "trace_real": float(traced.real), "trace_imag": float(traced.imag), "closed_spin_sum": expected, "spinor_absolute_error": abs(spin_sum - expected), "trace_absolute_error": float(abs(traced - expected)), "trace_scaled_error_4E2": float(abs(traced - expected) / unit), "trace_relative_error": float(abs(traced - expected) / expected) if expected else None, "overlap": overlap} def experiments(args): checks = Checks() exact_fixtures(checks) # Clifford signs are checked independently of the scattering trace. for mu in range(4): for nu in range(4): checks.close(f"Clifford {mu},{nu}", GAMMA[mu] @ GAMMA[nu] + GAMMA[nu] @ GAMMA[mu], (2 if mu == nu == 0 else -2 if mu == nu else 0) * I4) angles = np.linspace(math.radians(args.theta_min_deg), math.radians(args.theta_max_deg), args.samples) azimuth = math.radians(args.azimuth_deg) fixtures = [(f"fixed_beta_{beta:g}", 1.0, beta / math.sqrt(1 - beta * beta), 0.001) for beta in (0.05, 0.3, 0.8, 0.99)] fixtures += [("fixed_massless", 0.0, 1.0, 0.001), ("custom", args.mass, args.momentum, args.coupling)] angular, transitions = [], [] maximum_trace_relative = 0.0 for name, mass, momentum, coupling in fixtures: energy = math.hypot(mass, momentum) beta = momentum / energy for index, theta in enumerate(angles): theta = float(theta) result = spin_benchmarks(checks, mass, momentum, theta, azimuth, f"{name} angle {index}") factor = spin_factor(mass, momentum, theta) sh = angle_parts(theta)[2] baseline = (coupling / (2 * momentum * beta * sh * sh)) ** 2 cross = differential(mass, momentum, coupling, theta) potential = 4 * math.pi * coupling / (4 * momentum * momentum * sh * sh) from_spinors = potential ** 2 / (16 * math.pi ** 2) * result["spin_sum"] checks.close(f"{name} cross section {index} normalized", from_spinors / baseline if baseline else from_spinors, factor if baseline else 0.0) checks.condition(f"{name} nonnegative cross section {index}", cross >= 0 and math.isfinite(cross)) if result["trace_relative_error"] is not None: maximum_trace_relative = max(maximum_trace_relative, result["trace_relative_error"]) nr_rutherford = (coupling / (2 * mass * beta ** 2 * sh * sh)) ** 2 if mass else None angular.append({"case": name, "mass": mass, "momentum": momentum, "energy": energy, "beta": beta, "coupling": coupling, "born_parameter_abs_g_over_beta": abs(coupling) / beta, "theta_deg": math.degrees(theta), "azimuth_deg": args.azimuth_deg, "spin_factor": factor, "differential_area_per_sr": cross, "relativistic_spinless_baseline": baseline, "nonrelativistic_Rutherford_same_v": nr_rutherford, **{key: value for key, value in result.items() if key != "overlap"}}) if name == "custom": for final_spin in range(2): for initial_spin in range(2): raw = float(abs(result["overlap"][final_spin, initial_spin]) ** 2) transitions.append({"theta_deg": math.degrees(theta), "azimuth_deg": args.azimuth_deg, "initial_pauli_seed_index": initial_spin, "final_pauli_seed_index": final_spin, "overlap_real": float(result["overlap"][final_spin, initial_spin].real), "overlap_imag": float(result["overlap"][final_spin, initial_spin].imag), "overlap_squared": raw, "initial_average_contribution": raw / 2, "fraction_of_unpolarized_sum": raw / (2 * result["spin_sum"]) if result["closed_spin_sum"] else None}) # Rotated incoming directions and independent azimuths guard against axial-only tests. rotation = np.array([[0.36, -0.8, 0.48], [0.48, 0.6, 0.64], [-0.8, 0, 0.6]]) checks.close("spatial rotation orthogonality", rotation.T @ rotation, np.eye(3)) for mass, momentum in ((1, 1e-6), (1, 1), (1, 1e6), (0, 3), (1e-6, 1), (1e6, 1)): for theta in (0.003, 0.3, math.pi / 2, 2.7, math.pi): for phi in (0.0, 0.7, 2.1): spin_benchmarks(checks, mass, momentum, theta, phi, f"oblique m={mass} p={momentum} t={theta} f={phi}", rotation) checks.close("owner beta=4/5 right angle", spin_factor(1, 4 / 3, math.pi / 2), 17 / 25) checks.close("owner beta=4/5 backscatter", spin_factor(1, 4 / 3, math.pi), 9 / 25) # This deliberately wrong convention is a control, not an accepted spin sum. correct = float(spin_trace(3, np.array([0, 0, 4]), np.array([4, 0, 0])).real) checks.condition("missing initial spin average is detectably wrong", abs(2 * correct - 68) > 1) limits = [] for beta in (0.1, 0.03, 0.01, 0.003, 0.001, 0.0001): mass, coupling = 1.0, 0.001 * beta momentum = beta / math.sqrt(1 - beta * beta) for theta in (math.pi / 3, math.pi / 2, math.pi): factor = spin_factor(mass, momentum, theta) cross = differential(mass, momentum, coupling, theta) sh = angle_parts(theta)[2] nr = (coupling / (2 * mass * beta * beta * sh * sh)) ** 2 ratio = cross / nr checks.close(f"NR ratio beta={beta} theta={theta}", ratio, (1 - beta ** 2) * factor) checks.condition(f"NR convergence bound beta={beta} theta={theta}", abs(1 - ratio) <= 2 * beta ** 2 + 2e-15) limits.append({"limit": "nonrelativistic_same_velocity", "parameter": beta, "theta_deg": math.degrees(theta), "spin_factor": factor, "ratio_to_reference": ratio, "reference": "Rutherford at same mass and velocity", "born_parameter": 0.001}) for ratio in (1, 2, 4, 8, 16, 64, 1024, 1e6): factor = spin_factor(1, ratio, math.pi) checks.close(f"massive backscatter p/m={ratio}", factor * (1 + ratio ** 2), 1) limits.append({"limit": "high_energy_massive_backscatter", "parameter": ratio, "theta_deg": 180, "spin_factor": factor, "ratio_to_reference": factor, "reference": "relativistic spinless baseline", "born_parameter": 0.001 / (ratio / math.hypot(1, ratio))}) for theta in (0.3, 1.0, math.pi / 2, math.pi): factor = spin_factor(0, 1, theta) checks.close(f"massless suppression theta={theta}", factor, angle_parts(theta)[3] ** 2) checks.condition("massless backscatter is exactly zero", differential(0, 1, 0.001, math.pi) == 0) for coupling in (0.0, 0.001, 0.7): checks.close(f"leading charge-sign independence g={coupling}", differential(1, 4 / 3, coupling, 0.4), differential(1, 4 / 3, -coupling, 0.4)) lower, upper = math.radians(1), math.radians(160) exact = acceptance_exact(1, 4 / 3, 0.001, lower, upper) convergence, previous_error = [], None for intervals in (16, 32, 64, 128, 256): numeric = acceptance_simpson(1, 4 / 3, 0.001, lower, upper, intervals) relative = abs(numeric - exact) / exact order = math.log(previous_error / relative, 2) if previous_error else None convergence.append({"case": "fixed_acceptance", "theta_min_deg": 1, "theta_max_deg": 160, "Simpson_intervals_in_log_sin2": intervals, "numerical_area": numeric, "closed_area": exact, "absolute_error": abs(numeric - exact), "relative_error": relative, "observed_order": order}) if previous_error: checks.condition(f"acceptance refinement improves at N={intervals}", relative < previous_error) previous_error = relative checks.condition("acceptance asymptotic fourth order", 3.8 < convergence[-1]["observed_order"] < 4.2) checks.condition("acceptance fixed fixture accuracy", previous_error < 2e-8, absolute_error=convergence[-1]["absolute_error"], scaled_error=previous_error, tolerance=2e-8) custom_exact = acceptance_exact(args.mass, args.momentum, args.coupling, float(angles[0]), float(angles[-1])) custom_numeric = acceptance_simpson(args.mass, args.momentum, args.coupling, float(angles[0]), float(angles[-1]), 2048) custom_relative = abs(custom_numeric - custom_exact) / custom_exact if custom_exact else 0.0 checks.condition("custom finite acceptance accuracy", custom_relative < 2e-8, absolute_error=abs(custom_numeric - custom_exact), scaled_error=custom_relative, tolerance=2e-8) convergence.append({"case": "custom_acceptance", "theta_min_deg": args.theta_min_deg, "theta_max_deg": args.theta_max_deg, "Simpson_intervals_in_log_sin2": 2048, "numerical_area": custom_numeric, "closed_area": custom_exact, "absolute_error": abs(custom_numeric - custom_exact), "relative_error": custom_relative, "observed_order": None}) cutoffs, previous_area = [], 0.0 leading = 4 * math.pi * 0.001 ** 2 / ((4 / 3) ** 2 * 0.8 ** 2) for degrees in (20, 10, 5, 2.5, 1.25, 0.625): theta = math.radians(degrees) area = acceptance_exact(1, 4 / 3, 0.001, theta, math.pi) checks.condition(f"forward cutoff divergence theta={degrees}", area > previous_area) previous_area = area cutoffs.append({"theta_min_deg": degrees, "theta_max_deg": 180, "finite_acceptance_area": area, "theta_min_squared_times_area": theta ** 2 * area, "asymptotic_scaled_area": leading, "scaled_ratio": theta ** 2 * area / leading}) checks.condition("forward-cutoff asymptotic coefficient", abs(cutoffs[-1]["scaled_ratio"] - 1) < 3e-4) data = {"angular": angular, "spin_transitions": transitions, "limits": limits, "convergence": convergence, "forward_cutoffs": cutoffs} diagnostics = {"maximum_grid_trace_relative_error_when_nonzero": maximum_trace_relative, "trace_comparison_scale": "4 E^2; raw and relative errors are retained in angular.csv", "maximum_scaled_algebra_check_error": max((row.get("scaled_error", 0) for row in checks.rows if "acceptance" not in row["name"]), default=0), "fixed_acceptance_finest_relative_error": convergence[-2]["relative_error"], "fixed_acceptance_finest_observed_order": convergence[-2]["observed_order"], "custom_acceptance_relative_error": custom_relative, "custom_finite_acceptance_area": custom_exact, "no_finite_unscreened_total_cross_section": True} return checks, data, diagnostics def parser_and_args(argv): parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) parser.add_argument("--mass", type=float, default=1.0, help="0 or [1e-6,1e6], natural energy units (default 1)") parser.add_argument("--momentum", type=float, default=4 / 3, help="positive magnitude [1e-6,1e6] (default 4/3)") parser.add_argument("--coupling", type=float, default=-0.001, help="signed g=qQ/(4*pi), |g|<=1 (default -0.001)") parser.add_argument("--theta-min-deg", type=float, default=1.0, help="positive lower angle, at least .001 degree") parser.add_argument("--theta-max-deg", type=float, default=180.0, help="upper angle <=180; acceptance width >=.1 degree") parser.add_argument("--azimuth-deg", type=float, default=37.0, help="azimuth in [-360,360] (default 37)") parser.add_argument("--samples", type=int, default=181, help="angular samples, 9..2001 (default 181)") parser.add_argument("--output-dir", type=Path, help="output directory; default is a fresh temporary directory") parser.add_argument("--json-output", type=Path, help="override summary.json path") parser.add_argument("--angular-csv", type=Path, help="override primary angular.csv path") args = parser.parse_args(argv) numeric = (args.mass, args.momentum, args.coupling, args.theta_min_deg, args.theta_max_deg, args.azimuth_deg) if not all(math.isfinite(value) for value in numeric): parser.error("all numeric parameters must be finite") if args.mass != 0 and not 1e-6 <= args.mass <= 1e6: parser.error("mass must be zero or within [1e-6,1e6]") if not 1e-6 <= args.momentum <= 1e6: parser.error("momentum must be within [1e-6,1e6]") if args.mass and not 1e-6 <= args.momentum / args.mass <= 1e6: parser.error("massive diagnostics require 1e-6 <= momentum/mass <= 1e6") if abs(args.coupling) > 1 or (args.coupling and abs(args.coupling) < 1e-12): parser.error("coupling must be zero or have magnitude within [1e-12,1]") if not .001 <= args.theta_min_deg < args.theta_max_deg <= 180: parser.error("angles require .001 <= theta-min-deg < theta-max-deg <= 180") if args.theta_max_deg - args.theta_min_deg < .1 - 1e-12: parser.error("the validated acceptance width is at least .1 degree") if not -360 <= args.azimuth_deg <= 360 or not 9 <= args.samples <= 2001: parser.error("azimuth must be in [-360,360] and samples in [9,2001]") return parser, args def output_paths(parser, args): directory = args.output_dir or Path(tempfile.gettempdir()) / ("mott-scattering-" + uuid.uuid4().hex) paths = {"angular": args.angular_csv or directory / "angular.csv", "spin_transitions": directory / "spin-transitions.csv", "limits": directory / "limits.csv", "convergence": directory / "angular-convergence.csv", "forward_cutoffs": directory / "forward-cutoffs.csv", "summary": args.json_output or directory / "summary.json"} try: # Inspect lexical paths before resolve(), which would hide a dangling # symlink or a .git component resolving through a directory symlink. for path in paths.values(): path = path.expanduser().absolute() if any(part.casefold() == ".git" for part in path.parts): parser.error("outputs inside .git are prohibited") if path.exists() or path.is_symlink(): parser.error(f"refusing to overwrite existing output: {path}") for ancestor in path.parents: if ancestor.is_symlink() and not ancestor.exists(): parser.error(f"output parent is an unresolved symlink: {ancestor}") paths = {name: path.expanduser().resolve() for name, path in paths.items()} lowered = [os.path.normcase(str(path)).casefold() for path in paths.values()] if len(set(lowered)) != len(lowered): parser.error("output paths must be distinct") for name, path in paths.items(): if any(part.casefold() == ".git" for part in path.parts): parser.error("outputs inside .git are prohibited") if path.exists() or path.is_symlink(): parser.error(f"refusing to overwrite existing output: {path}") for ancestor in path.parents: if ancestor.exists() and not ancestor.is_dir(): parser.error(f"output parent is not a directory: {ancestor}") if ancestor.is_symlink() and not ancestor.exists(): parser.error(f"output parent is an unresolved symlink: {ancestor}") for other_name, other in paths.items(): if name != other_name and path in other.parents: parser.error("an output file cannot be the parent directory of another output") except (OSError, ValueError, RuntimeError) as error: parser.error(f"invalid output path: {error}") return paths def csv_bytes(rows): stream = io.StringIO(newline="") writer = csv.DictWriter(stream, fieldnames=list(rows[0]), lineterminator="\n") writer.writeheader() writer.writerows(rows) return stream.getvalue().encode("utf-8") def main(argv=None): parser, args = parser_and_args(argv) paths = output_paths(parser, args) checks, data, diagnostics = experiments(args) contents = {name: csv_bytes(rows) for name, rows in data.items()} beta = args.momentum / math.hypot(args.mass, args.momentum) report = { "experiment": "Mott scattering from a fixed Coulomb source at leading Born order", "passed": checks.passed, "check_count": len(checks.rows), "failed_checks": [row for row in checks.rows if not row["passed"]], "maximum_absolute_check_error": max((row.get("absolute_error", 0) for row in checks.rows), default=0), "maximum_scaled_check_error": max((row.get("scaled_error", 0) for row in checks.rows), default=0), "parameters": {name: getattr(args, name) for name in ("mass", "momentum", "coupling", "theta_min_deg", "theta_max_deg", "azimuth_deg", "samples")}, "units": {"hbar": 1, "c": 1, "metric": "+---", "mass_and_momentum": "energy", "cross_section": "energy^-2", "angles_input": "degrees"}, "methods": ["explicit 4-component positive-energy spinors", "independent 4x4 gamma-matrix trace", "exact rational massive/massless fixtures", "rotated spin bases and spatial directions", "composite Simpson quadrature in log(sin(theta/2)^2)", "independent 60-digit Decimal closed acceptance integral"], "assumptions": ["infinitely heavy spinless fixed point Coulomb source", "elastic one-body external-potential Born approximation", "unpolarized initial spin average and final spin sum", "no recoil, screening, finite-size structure, exact Coulomb phases, radiative corrections or pair calculation", "finite angular acceptance only; unscreened total cross section diverges"], "born_parameter_abs_g_over_beta": abs(args.coupling) / beta, "born_parameter_note": "A small |g|/beta is a sufficient perturbative regime away from the forward singularity; this diagnostic does not bound physical truncation error.", "precision_note": "Near massive ultrarelativistic backscatter the trace subtracts O(E^2) terms. Small absolute error divided by 4E^2 does not imply small relative error in the suppressed signal. Use the stable scalar expression; raw residuals remain available.", "spin_transition_basis": "Seed index 0 is (1,0), index 1 is (0,1). These are canonical rest-z spin labels for m>0, not momentum-helicity labels; m=0 has no rest frame. A zero exact spin sum leaves conditional fractions undefined.", "runtime": {"python": platform.python_version(), "numpy": np.__version__}, "randomness": "none; deterministic fixtures (only default temporary path has a unique identifier)", "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "output_files": {name: str(path) for name, path in paths.items()}, "output_sha256": {name: hashlib.sha256(content).hexdigest() for name, content in contents.items()}, "row_counts": {name: len(rows) for name, rows in data.items()}, "diagnostics": diagnostics, "checks": checks.rows, "owner": "/relativistic-qm/mott-scattering/", "references": ["James D. Bjorken and Sidney D. Drell, Relativistic Quantum Mechanics, McGraw-Hill (1964).", "Thomas DeGrand, A One-Semester Course on Quantum Field Theory (30 December 2025), section 8.7, https://spot.colorado.edu/~degrand/qftbook.pdf"]} contents["summary"] = (json.dumps(report, indent=2, allow_nan=False) + "\n").encode("utf-8") created = [] try: for name, content in contents.items(): path = paths[name] path.parent.mkdir(parents=True, exist_ok=True) with path.open("xb") as stream: created.append(path) stream.write(content) except OSError as error: for path in created: try: path.unlink() except OSError: pass parser.error(f"could not write new outputs: {error}") print(json.dumps({"passed": checks.passed, "check_count": len(checks.rows), "maximum_scaled_error": report["maximum_scaled_check_error"], "failed_check_count": len(report["failed_checks"]), "summary": str(paths["summary"])}, indent=2)) return 0 if checks.passed else 1 if __name__ == "__main__": raise SystemExit(main())