#!/usr/bin/env python3 """Scalar mode Green functions: source normalization and boundary conditions. Run with --help. Python standard library + NumPy only; deterministic/offline. Natural units hbar=c=1. At ONE fixed spatial momentum p and m>0, set omega=sqrt(p^2+m^2), L=d_t^2+omega^2, and normalize L*g=delta(t): g_R=theta(t)*sin(omega*t)/omega, g_A=-theta(-t)*sin(omega*t)/omega, g_F=i*exp(-i*omega*abs(t))/(2*omega). All are continuous at zero and have first-derivative jump +1. The conventional correlator D_F=-i*g_F instead has source -i*delta. It is not a unit-source classical response kernel. The code does not evaluate a pointwise delta. Convolve each kernel directly with the smooth compact source J(t)=exp(1-1/(1-(t/T)^2)) for |t|0, Im Omega<0. It inverts L-i*epsilon, not L; the off-source L residual is i*epsilon*g_epsilon. Epsilon refinement is separate from numerical quadrature. Fixed-momentum time support is NOT a test of spacetime microcausality. There is no inverse spatial Fourier transform, light-cone plot, finite spatial box, frequency-cutoff integration, or particle-probability interpretation here. Feynman response before the source means different boundary data, not a retarded signal or backwards communication. These are free massive kernels. Four CSVs export kernels, compact-source responses, numerical refinement, and regulated kernels. JSON records conventions, checks, versions and SHA-256. Only new output files are written; defaults use a fresh temporary directory. Owner: /relativistic-qm/klein-gordon-propagators/. References: Schwartz, Quantum Field Theory and the Standard Model (2014); Tong, Quantum Field Theory lectures, sections 2.6--2.7 (2006). Tested environment: Python 3.12.14, NumPy 2.3.5. """ from __future__ import annotations import argparse import cmath import csv from fractions import Fraction import hashlib import json import math from pathlib import Path import platform import sys import tempfile import numpy as np KINDS = ("retarded", "advanced", "feynman") 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=3e-12): absolute = float(np.max(np.abs(np.asarray(actual) - np.asarray(expected)))) scale = max(1.0, float(np.max(np.abs(expected)))) self.condition(name, 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 kernel(kind, time, omega): t = np.asarray(time, dtype=float) if kind == "retarded": return np.where(t >= 0, t * np.sinc(omega * t / np.pi), 0).astype(complex) if kind == "advanced": return np.where(t <= 0, -t * np.sinc(omega * t / np.pi), 0).astype(complex) if kind == "feynman": return 1j * np.exp(-1j * omega * np.abs(t)) / (2 * omega) raise ValueError("unknown kernel kind") def kernel_derivative(kind, time, omega): t = np.asarray(time, dtype=float) if np.any(t == 0): raise ValueError("the derivative at the source needs a specified one-sided limit") if kind == "retarded": return np.where(t > 0, np.cos(omega * t), 0).astype(complex) if kind == "advanced": return np.where(t < 0, -np.cos(omega * t), 0).astype(complex) if kind == "feynman": return .5 * np.sign(t) * np.exp(-1j * omega * np.abs(t)) raise ValueError("unknown kernel kind") def source(time, half_width): t = np.asarray(time, dtype=float) scaled = t / half_width inside = np.abs(scaled) < 1 result = np.zeros_like(t) result[inside] = np.exp(1 - 1 / (1 - scaled[inside] ** 2)) return result def simpson(function, left, right, panels): if panels < 2 or panels % 2: raise ValueError("Simpson panels must be a positive even number >=2") if left == right: return 0j x = np.linspace(left, right, panels + 1) values = np.asarray(function(x), dtype=complex) return (right - left) / (3 * panels) * (values[0] + values[-1] + 4 * np.sum(values[1:-1:2]) + 2 * np.sum(values[2:-1:2])) def convolve(kind, times, omega, half_width, panels): result = [] for t in np.asarray(times): split = min(half_width, max(-half_width, float(t))) integrand = lambda s: kernel(kind, t - s, omega) * source(s, half_width) result.append(simpson(integrand, -half_width, split, panels) + simpson(integrand, split, half_width, panels)) return np.asarray(result) def source_fourier(omega, half_width, panels, sign=-1): return simpson(lambda s: np.exp(sign * 1j * omega * s) * source(s, half_width), -half_width, half_width, panels) def rk4_forced(omega, half_width, horizon, steps, backward=False): """Independent time-domain ODE; no Green kernel is called here.""" times = np.linspace(horizon, -horizon, steps + 1) if backward else np.linspace(-horizon, horizon, steps + 1) dt = float(times[1] - times[0]) state = np.zeros(2, dtype=complex) states = [state.copy()] def rhs(t, y): return np.array([y[1], source(t, half_width) - omega ** 2 * y[0]], complex) for t in times[:-1]: k1 = rhs(t, state) k2 = rhs(t + dt / 2, state + dt * k1 / 2) k3 = rhs(t + dt / 2, state + dt * k2 / 2) k4 = rhs(t + dt, state + dt * k3) state = state + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6 states.append(state.copy()) states = np.asarray(states) return (times[::-1], states[::-1]) if backward else (times, states) def regulated_feynman(time, omega, epsilon): frequency = cmath.sqrt(omega ** 2 - 1j * epsilon) return 1j * np.exp(-1j * frequency * np.abs(np.asarray(time))) / (2 * frequency) def kernel_checks(checks): rows = [] for kind, left, right in (("retarded", Fraction(0), Fraction(1)), ("advanced", Fraction(-1), Fraction(0)), ("feynman", Fraction(-1, 2), Fraction(1, 2))): checks.condition("exact derivative jump " + kind, right - left == 1, arithmetic="Fraction") checks.condition("exact correlator source factor", -1j * (Fraction(1, 2) - Fraction(-1, 2)) == -1j) for omega in (.4, 1.3, 3.0): times = np.array([-2.1, -.73, -.1, 0, .1, .73, 2.1]) r, a, f = (kernel(kind, times, omega) for kind in KINDS) checks.close(f"omega={omega} R-A sine solution", r - a, np.sin(omega * times) / omega) checks.close(f"omega={omega} F-R homogeneous coefficient", f - r, 1j * np.exp(1j * omega * times) / (2 * omega)) checks.close(f"omega={omega} F-A homogeneous coefficient", f - a, 1j * np.exp(-1j * omega * times) / (2 * omega)) checks.close(f"omega={omega} exact retarded past support", r[times < 0], np.zeros(sum(times < 0))) checks.close(f"omega={omega} exact advanced future support", a[times > 0], np.zeros(sum(times > 0))) checks.condition(f"omega={omega} Feynman nonzero before impulse", np.all(abs(f[times < 0]) > 0)) for kind in KINDS: previous = None for dt in (.04, .02, .01, .005): point = .73 if kind != "advanced" else -.73 residual = (kernel(kind, point + dt, omega) - 2 * kernel(kind, point, omega) + kernel(kind, point - dt, omega)) / dt ** 2 + omega ** 2 * kernel(kind, point, omega) error = abs(complex(residual)) order = None if previous is None else math.log2(previous / error) rows.append({"study": "off_source_centered_derivative", "kind": kind, "omega": omega, "resolution": dt, "error": error, "observed_order": order}) previous = error checks.condition(f"omega={omega} {kind} off-source derivative second order", abs(order - 2) < .02, observed_order=order) for half_interval in (.15, .7): integral = simpson(lambda t: kernel(kind, t, omega), -half_interval, 0, 512) + simpson(lambda t: kernel(kind, t, omega), 0, half_interval, 512) weak_source = kernel_derivative(kind, half_interval, omega) - kernel_derivative(kind, -half_interval, omega) + omega ** 2 * integral checks.close(f"omega={omega} {kind} integrated unit delta width={half_interval}", weak_source, 1, tolerance=2e-11) # One-sided slopes approach the jump; F has O(dt), R/A O(dt^2). dt = 1e-5 jump = (kernel(kind, dt, omega) + kernel(kind, -dt, omega) - 2 * kernel(kind, 0, omega)) / dt checks.close(f"omega={omega} {kind} numerical one-sided derivative jump", jump, 1, tolerance=2e-5) return rows def source_studies(checks): omega, half_width, horizon = 1.3, 2.0, 6.0 rows = [] output_times = np.linspace(-horizon, horizon, 41) reference = {kind: convolve(kind, output_times, omega, half_width, 2048) for kind in KINDS} shell_minus = source_fourier(omega, half_width, 2048, -1) shell_plus = source_fourier(omega, half_width, 2048, 1) checks.close("source F-R homogeneous response", reference["feynman"] - reference["retarded"], 1j * np.exp(1j * omega * output_times) * shell_minus / (2 * omega), tolerance=2e-12) checks.close("source F-A homogeneous response", reference["feynman"] - reference["advanced"], 1j * np.exp(-1j * omega * output_times) * shell_plus / (2 * omega), tolerance=2e-12) checks.close("retarded response vanishes before compact source", reference["retarded"][output_times < -half_width], 0) checks.close("advanced response vanishes after compact source", reference["advanced"][output_times > half_width], 0) checks.condition("Feynman response nonzero before compact source", np.max(abs(reference["feynman"][output_times < -half_width])) > .1) for kind in KINDS: evaluation = np.array([-.7, .3, .9]) exact = convolve(kind, evaluation, omega, half_width, 4096) previous = None for panels in (16, 32, 64, 128, 256, 512): approximation = convolve(kind, evaluation, omega, half_width, panels) error = float(np.max(abs(approximation - exact))) order = None if previous is None else math.log2(previous / error) rows.append({"study": "source_quadrature", "kind": kind, "omega": omega, "resolution": panels, "error": error, "observed_order": order}) previous = error checks.condition(kind + " independent quadrature convergence", error < 2e-10, final_error=error) checks.condition(kind + " asymptotic Simpson fourth order", abs(order - 4) < .03, observed_order=order) for backward, kind in ((False, "retarded"), (True, "advanced")): previous = None for steps in (120, 240, 480, 960): _, states = rk4_forced(omega, half_width, horizon, steps, backward) sampled = states[::steps // (len(output_times) - 1), 0] error = float(np.max(abs(sampled - reference[kind]))) order = None if previous is None else math.log2(previous / error) rows.append({"study": "independent_RK4_time_step", "kind": kind, "omega": omega, "resolution": 2 * horizon / steps, "error": error, "observed_order": order}) previous = error checks.condition(kind + " independent forced ODE accuracy", error < 1e-7, final_error=error) checks.condition(kind + " independent RK4 fourth order", abs(order - 4) < .2, observed_order=order) if not backward: f_from_ode = sampled + 1j * np.exp(1j * omega * output_times) * shell_minus / (2 * omega) checks.close("Feynman convolution versus ODE plus homogeneous boundary data", reference["feynman"], f_from_ode, tolerance=1e-7) return rows def regulator_study(checks): rows = [] for omega in (.7, 1.3, 3): times = np.array([-2.0, -.5, 0, .5, 2.0]) target = kernel("feynman", times, omega) previous = None for epsilon in (.1, .025, .00625, .0015625): frequency = cmath.sqrt(omega ** 2 - 1j * epsilon) regulated = regulated_feynman(times, omega, epsilon) checks.condition(f"epsilon={epsilon} omega={omega} decaying pole branch", frequency.real > 0 and frequency.imag < 0) regularized_residual = (-frequency ** 2 + omega ** 2 - 1j * epsilon) * regulated checks.close(f"epsilon={epsilon} omega={omega} regulated homogeneous operator", regularized_residual[times != 0], np.zeros(sum(times != 0))) unregularized_residual = (-frequency ** 2 + omega ** 2) * regulated checks.close(f"epsilon={epsilon} omega={omega} unregularized residual is i epsilon G", unregularized_residual[times != 0], 1j * epsilon * regulated[times != 0]) half_interval = .5 integral = simpson(lambda t: regulated_feynman(t, omega, epsilon), -half_interval, 0, 512) + simpson(lambda t: regulated_feynman(t, omega, epsilon), 0, half_interval, 512) weak_source = np.exp(-1j * frequency * half_interval) + (omega ** 2 - 1j * epsilon) * integral checks.close(f"epsilon={epsilon} omega={omega} regulated integrated unit source", weak_source, 1, tolerance=2e-11) error = float(np.max(abs(regulated - target))) order = None if previous is None else math.log(previous / error, 4) rows.append({"study": "finite_Feynman_epsilon", "omega": omega, "resolution": epsilon, "error_on_fixed_time_window": error, "observed_order": order, "max_unregularized_off_source_residual": float(np.max(abs(unregularized_residual[times != 0])))}) previous = error checks.condition(f"omega={omega} epsilon approaches zero linearly", abs(order - 1) < .04, observed_order=order) return rows def write_csv(path, rows): fields = list(dict.fromkeys(key for row in rows for key in row)) with path.open("x", encoding="utf8", newline="") as stream: writer = csv.DictWriter(stream, fieldnames=fields) writer.writeheader() writer.writerows(rows) def parser(): p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) p.add_argument("--output-dir", type=Path, help="new output files; default fresh temporary directory") p.add_argument("--json-output", type=Path, help="override summary JSON path") p.add_argument("--kernels-csv", type=Path, help="override sampled-kernel CSV path") p.add_argument("--mass", type=float, default=1.0, help="positive mass .1..4 (default 1)") p.add_argument("--momentum", type=float, default=.75, help="spatial-momentum magnitude 0..5 (default .75)") p.add_argument("--source-width", type=float, default=2.0, help="source support half-width T, .25..4 (default 2)") p.add_argument("--horizon", type=float, default=6.0, help="time interval [-horizon,horizon], >T and <=20 (default 6)") p.add_argument("--samples", type=int, default=121, help="odd output sample count 21..401 (default 121)") p.add_argument("--panels", type=int, default=512, help="even Simpson panels per smooth segment, 32..4096 (default 512)") p.add_argument("--ode-steps", type=int, default=960, help="RK4 steps, multiple of samples-1 and <=32768 (default 960)") p.add_argument("--epsilon", type=float, default=.05, help="positive regulator <=.5, units of energy squared (default .05)") return p def main(argv=None): p = parser() args = p.parse_args(argv) bounds = {"mass": (.1, 4), "momentum": (0, 5), "source_width": (.25, 4), "horizon": (.5, 20), "epsilon": (1e-6, .5)} for name, (low, high) in bounds.items(): value = getattr(args, name) if not math.isfinite(value) or not low <= value <= high: p.error(f"{name.replace('_', '-')} must be finite and in {low}..{high}") if args.horizon <= args.source_width: p.error("horizon must exceed the source half-width") if not 21 <= args.samples <= 401 or args.samples % 2 != 1: p.error("samples must be odd and in 21..401") if not 32 <= args.panels <= 4096 or args.panels % 2: p.error("panels must be even and in 32..4096") if not 40 <= args.ode_steps <= 32768 or args.ode_steps % (args.samples - 1): p.error("ode-steps must be 40..32768 and a multiple of samples-1") omega = math.hypot(args.mass, args.momentum) dt = 2 * args.horizon / args.ode_steps if omega * dt > .2 or args.source_width / dt < 20: p.error("increase ode-steps: require omega*dt<=.2 and at least20 steps per source half-width") output = args.output_dir or Path(tempfile.mkdtemp(prefix="scalar-propagators-")) paths = {"kernels": args.kernels_csv or output / "kernels.csv", "responses": output / "responses.csv", "refinement": output / "refinement.csv", "regularization": output / "regularization.csv", "summary": args.json_output or output / "summary.json"} if len({path.resolve() for path in paths.values()}) != len(paths): p.error("output paths must be distinct") for path in paths.values(): if path.exists() or path.is_symlink(): p.error(f"existing outputs are never overwritten: {path}") if ".git" in [part.lower() for part in path.resolve().parts]: p.error(f"output cannot be inside .git: {path}") for path in paths.values(): path.parent.mkdir(parents=True, exist_ok=True) checks = Checks() refinement = kernel_checks(checks) + source_studies(checks) regulator_checks = regulator_study(checks) times = np.linspace(-args.horizon, args.horizon, args.samples) kernels = {kind: kernel(kind, times, omega) for kind in KINDS} responses = {kind: convolve(kind, times, omega, args.source_width, args.panels) for kind in KINDS} shell = source_fourier(omega, args.source_width, args.panels) ode_values = {} for backward, kind in ((False, "retarded"), (True, "advanced")): _, states = rk4_forced(omega, args.source_width, args.horizon, args.ode_steps, backward) ode_values[kind] = states[::args.ode_steps // (args.samples - 1), 0] ode_values["feynman"] = ode_values["retarded"] + 1j * np.exp(1j * omega * times) * shell / (2 * omega) kernel_rows, response_rows, regularized_rows = [], [], [] regulated = regulated_feynman(times, omega, args.epsilon) for index, time in enumerate(times): for kind in KINDS: g, response, independent = kernels[kind][index], responses[kind][index], ode_values[kind][index] kernel_rows.append({"kind": kind, "time": float(time), "omega": omega, "real": float(g.real), "imag": float(g.imag), "source_jump": 1}) response_rows.append({"kind": kind, "time": float(time), "source": float(source(time, args.source_width)), "quadrature_real": float(response.real), "quadrature_imag": float(response.imag), "independent_ODE_real": float(independent.real), "independent_ODE_imag": float(independent.imag), "absolute_difference": float(abs(response - independent)), "ODE_boundary_data": "retarded ODE plus homogeneous Feynman data" if kind == "feynman" else "zero final data" if kind == "advanced" else "zero initial data"}) g = regulated[index] regularized_rows.append({"time": float(time), "omega": omega, "epsilon": args.epsilon, "real": float(g.real), "imag": float(g.imag), "difference_from_unregulated_F": float(abs(g - kernels["feynman"][index])), "unregularized_off_source_residual_real": float((1j * args.epsilon * g).real) if time != 0 else None, "unregularized_off_source_residual_imag": float((1j * args.epsilon * g).imag) if time != 0 else None}) for name, rows in (("kernels", kernel_rows), ("responses", response_rows), ("refinement", refinement + regulator_checks), ("regularization", regularized_rows)): write_csv(paths[name], rows) report = {"program": Path(__file__).name, "passed": checks.passed, "check_count": len(checks.rows), "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "environment": {"python": platform.python_version(), "numpy": np.__version__, "platform": platform.platform()}, "randomness": "none", "parameters": {key: getattr(args, key) for key in ("mass", "momentum", "source_width", "horizon", "samples", "panels", "ode_steps", "epsilon")}, "omega": omega, "checks": checks.rows, "configured_response_errors": {kind: float(np.max(abs(responses[kind] - ode_values[kind]))) for kind in KINDS}, "interpretation": ["All exported g kernels use unit delta source; D_F=-i*g_F has source -i delta.", "A temporal Fourier-mode experiment does not test full spacetime microcausality.", "The source is smooth and exactly compact in time; retarded/advanced ODE data vanish outside its support.", "Feynman response differs from retarded response by homogeneous boundary data, not a causal signal.", "Fixed fixture checks and configured output diagnostics are distinct; passing fixtures does not certify arbitrary quadrature accuracy.", "Quadrature panels and ODE time step are refined independently; output sample spacing is only for plotting.", "Finite epsilon inverts L-i epsilon; its nonzero unregularized off-source residual is expected.", "No particle probability, spatial Fourier cutoff, finite box, or interacting-background calculation is performed."], "output_files": {key: str(path) for key, path in paths.items()}, "output_sha256": {key: hashlib.sha256(path.read_bytes()).hexdigest() for key, path in paths.items() if key != "summary"}} with paths["summary"].open("x", encoding="utf8") as stream: json.dump(report, stream, indent=2, allow_nan=False) stream.write("\n") print(json.dumps({"passed": checks.passed, "checks": len(checks.rows), "failed": [row for row in checks.rows if not row["passed"]], "summary_json": str(paths["summary"])}, indent=2)) return 0 if checks.passed else 1 if __name__ == "__main__": sys.exit(main())