================================================================================
4D ENGINE — EQUATION SPECIFICATION v2
================================================================================
Purpose: Define EXACTLY what equations the engine implements.
         Each equation traced to its source document.
         NO additions, NO leanings, NO nudges.

Revision: v2 (2026-03-20) — addresses all 20 audit findings from v1.
Audit report: AUDIT_REPORT_v1.txt

================================================================================
UNITS SYSTEM                                              [Fixes audit #20]
================================================================================

All quantities are DIMENSIONLESS (natural simulation units):
  - c_3D = 1.0 (propagation speed normalized)
  - Base wavelength lambda_0 = 1.0
  - Base frequency f_0 = c / lambda_0 = 1.0
  - Grid spacing h = lambda_0 / PPW (points per wavelength)
  - Time step dt from CFL condition

The c_4D sweep values (1.5 to 1.8) are dimensionless ratios relative to c_3D.


================================================================================
1. THE WAVE EQUATION (4D, with potential term)          [Fixes audit #1, #9]
================================================================================

Source: math_framework B.4 (line 293), approach_B (line 89, 106)

  d²ψ/dt² = c² ∇²ψ - ω₀²ψ(1 + A(x)/A_ref) + J(x,t)                  (1)

where:
  ∇² is the 4D Laplacian (NOT ∇⁴ which is biharmonic)        [Fixes #3]
  ∇²ψ = Σᵢ₌₁⁴ [ψ(xᵢ+h) - 2ψ(xᵢ) + ψ(xᵢ-h)] / h²                   (2)

  ω₀ = base angular frequency = 2πf₀ (restoring force scale)
  A(x) = local amplitude (from eq. 6)
  A_ref = reference amplitude (normalization constant)
  J(x,t) = pulsed source term (from eq. 3)

The potential term ω₀²ψ(1 + A/A_ref) provides:
  - A restoring force on the field (Klein-Gordon mass-like term)
  - Coupling between field dynamics and local structure via A(x)
  - The mechanism by which organized regions behave differently
    from unorganized regions (f + A | t in the Lagrangian)

Source Lagrangian density (approach_B line 106):
  L = ½(∂ψ/∂t)² - ½c²(∇ψ)² - ½ω₀²ψ²(1 + A/A_ref) + J·ψ             (1b)

Equation (1) is the Euler-Lagrange equation of (1b).

FDTD update:
  ψ(t+dt) = 2ψ(t) - ψ(t-dt) + dt²[c²∇²ψ(t) - ω₀²ψ(t)(1+A/A_ref) + J(t)]


================================================================================
2. CFL STABILITY CONDITION (4D)                         [Fixes audit #8]
================================================================================

  dt < h / (c × √d)  where d = number of spatial dimensions = 4

  dt_max = h / (c_max × 2)   since √4 = 2

  For the sweep, c_max = 1.8 (highest sweep value).
  Safety factor: Courant number = 0.9

  dt = 0.9 × h / (c_max × 2) = 0.45 × h / c_max                       (2b)

  For h = 1/8 = 0.125 (8 points per wavelength), c_max = 1.8:
    dt = 0.45 × 0.125 / 1.8 = 0.03125

  Grid spacing h is UNIFORM in all 4 spatial dimensions.

  With potential term: ω₀ effectively increases the "speed" of the system.
  Modified CFL: dt < h / √(4c² + ω₀²h²). For ω₀ = 2π and h = 0.125:
    4c² + ω₀²h² ≈ 4(1.8²) + (6.28)²(0.0156) ≈ 12.96 + 0.616 ≈ 13.58
    dt < 0.125 / √13.58 ≈ 0.034
  So the potential term does NOT significantly tighten the CFL for this h.


================================================================================
3. THE SOURCE TERM: f | t                               [Fixes audit #15]
================================================================================

Source: theory.txt lines 145-147, math_framework B.3

  J(x, t) = A_src(x) × sin(2πf t) × W(t; r(x))                        (3)

  A_src(x) = source amplitude (nonzero only at source positions)

Windowing function with SMOOTH TRANSITIONS:                [Fixes audit #15]

  W(t; r) uses a raised cosine ramp over N_ramp timesteps at each
  on/off transition to avoid Gibbs artifacts from discontinuous sources.

  W_smooth(t; r) = W_rect(t; r) × ramp(t)                               (4)

  where ramp(t) transitions from 0→1 over N_ramp steps at pulse onset
  and from 1→0 over N_ramp steps at pulse offset.

  N_ramp = max(4, round(0.05 × T/dt))  (5% of period, minimum 4 steps)

  The rectangular window is the THEORY's specification (structured silence).
  The ramp is a numerical necessity. The ramp width is documented and its
  effect on results should be checked by varying N_ramp in sensitivity tests.

  r(x) = decoherence ratio (from eq. 7)
  T = 1/f = period
  f = source frequency in dimensionless units


================================================================================
4. SOURCE GEOMETRY: 24-CELL VERTICES                    [Fixes audit #5, #6]
================================================================================

Source: established mathematics (Schläfli 1852), 24cell_research.txt

Form A (cell-center form, 24 vertices on unit sphere, radius = 1.0):
  8 vertices: permutations of (±1, 0, 0, 0)
  16 vertices: all sign combinations of (½, ½, ½, ½)

Form B (edge-midpoint form, 24 vertices at radius √2):
  All permutations of (±1, ±1, 0, 0)

NORMALIZATION: Both forms MUST be normalized to the SAME radius.    [#5]
  Standard: normalize all to unit sphere (radius = 1.0).
  Form B normalized: divide by √2 → permutations of (±1/√2, ±1/√2, 0, 0)

The self-dual mapping between Form A and normalized Form B is an
isoclinic (Clifford) rotation of π/4 (45°), COMPOSED WITH the
normalization scaling. This is a rotation-plus-scaling, not a          [#6]
pure rotation. The 45° angle is exact; the scaling is √2.

PRIMARY TEST CONFIGURATION:                                     [Fixes #19]
  Sources at Form A vertices ONLY (24 sources on unit sphere).
  Observe whether the interference pattern produces structure at
  the dual (Form B) positions WITHOUT sources placed there.

SECONDARY CONFIGURATIONS (for comparison):
  (a) Sources at normalized Form B vertices only
  (b) Sources at all 48 vertices (Form A ∪ normalized Form B)


================================================================================
5. AMPLITUDE COUPLING: f + A | t                        [Fixes audit #2, #13]
================================================================================

Source: theory.txt lines 149-153, math_framework B.5

The amplitude A appears in TWO places (both required):                  [#2]
  (a) The potential term: ω₀²ψ(1 + A(x)/A_ref) in eq. (1)
  (b) The source amplitude: A_src(x) modulated by local structure

For the INITIAL TEST (Phase 1):                                         [#13]
  A(x) = A_base (uniform). No CN coupling.
  The CN(x) coordination number requires structure to already exist.
  Phase 1 runs without CN coupling to establish the base interference.

For PHASE 2 (after structure forms):
  CN(x) is computed from the time-averaged intensity field:
  - Identify local maxima of ⟨|ψ|²⟩ above a threshold
  - For each maximum, count neighbors within a cutoff radius
  - CN(x) = number of neighbor maxima
  - A(x) = A_base / (1 + CN(x) / CN_ref)                               (6)
  - Update CN every N_update timesteps (not every step)
  - N_update = round(T / dt) (once per period)

  CN_ref = 8 (BCC coordination, the vertex coordination of the 24-cell)


================================================================================
6. C_POTENTIAL: POSITION-DEPENDENT DECOHERENCE          [Fixes audit #7, #14]
================================================================================

Source: theory.txt lines 155-171, math_framework B.1

  r(x) = r₀ + α × κ(x)                                                 (7)

where:
  r₀ = base decoherence ratio = 0.3 (from math_framework B.1 optimal)
  α = coupling strength (dimensionless, sweep parameter)

  κ(x) = LOCAL CURVATURE of the energy density profile              [#7]
        = -∇²⟨|ψ(x)|²⟩ / max(|∇²⟨|ψ|²⟩|)                             (8)

  NOTE THE NEGATIVE SIGN: ∇² is negative at energy PEAKS (concave    [v2 fix]
  down) and positive at valleys. The minus sign ensures κ(x) is
  POSITIVE at peaks, so r(x) INCREASES where energy concentrates.
  This matches theory.txt line 161: "the space of decoherence
  INCREASES at the peak of the curve."

  Normalized to [-1, 1]. The theory says "the curve is Lagrangian"
  and "time slows as the bandwidth curve steepens."
  Steepness = curvature = second derivative.

  Theory.txt line 161: "the space of decoherence INCREASES at the peak
  of the curve" — the peak of the curvature, not the peak of the energy.

  ⟨|ψ|²⟩ is time-averaged over the most recent period T.

  CEILING: r(x) ≤ 0.5 (hard clamp)                                 [#14]
  Source: math_framework B.1 line 296: "Collapse at r=0.5:
          active/rest phase symmetry (exact)"

  When r(x) reaches 0.5, energy overflows to the 5th dimension (see §8).


================================================================================
7. BOUNDARY CONDITIONS                                  [Fixes audit #11]
================================================================================

PRIMARY: Absorbing boundaries (Perfectly Matched Layer, PML)

  PML layers of thickness N_pml = 10 grid points on all 8 faces of
  the 4D domain. The PML absorbs outgoing waves without reflection,
  simulating an infinite domain.

  PML conductivity profile: σ(d) = σ_max × (d/N_pml)²
  where d = depth into PML layer (0 at inner face, N_pml at outer face)
  σ_max = 3c / (2h × N_pml) (standard calibration)

  SECONDARY TEST: Periodic boundaries (for comparison).
  Periodic BCs create a 4-torus topology. If results differ significantly
  from PML, the boundary effects are important and must be documented.


================================================================================
8. DIMENSIONAL OVERFLOW: 5D CAPTURE                     [Fixes audit #17a]
================================================================================

Source: theory.txt lines 173-191, user directive 2026-03-20

When r(x) → 0.5 at any grid point, excess energy overflows to a
5th spatial dimension following the same scaling rules as lower dimensions.

Implementation: a THIN 5th dimension (N_5 = 8 grid points) attached
to the 4D grid. When r(x) ≥ 0.5 at position x in 4D:

  - Compute excess: E_overflow = |ψ(x)|² × (r(x) - 0.5) / 0.5
  - Inject E_overflow into the 5D layer at position (x, w=0)
  - The 5D layer propagates at c_5D = 21/8 × c_3D = 2.625
    (Fibonacci pair {8,13}, sum=21)
  - Clamp r(x) back to 0.5 after overflow

The 5D layer has its own wave equation:
  d²ψ₅/dt² = c₅² ∂²ψ₅/∂w² + S_overflow(x, w, t)

This is a 1D wave equation along the w axis at each 4D position,
making it computationally cheap (N_5 = 8 points per 4D position).

The energy in the 5D layer is tracked as "overflow energy" and
reported separately from the 4D field energy.


================================================================================
9. GRID AND DOMAIN SPECIFICATION                        [Fixes audit #12]
================================================================================

Domain: [-L, L]⁴ where L = 3.0 (in dimensionless units)
  This gives room for the unit-sphere sources plus propagation space.

Grid points per dimension:
  Phase 1 (feasibility): N = 32 → 32⁴ = 1,048,576 points
    Memory: 3 arrays × 8 bytes × 1M ≈ 24 MB (comfortable)
    h = 6.0/32 = 0.1875
    PPW = lambda_0/h = 1/0.1875 ≈ 5.3 points per wavelength

  Phase 2 (production): N = 48 → 48⁴ = 5,308,416 points
    Memory: 3 arrays × 8 bytes × 5.3M ≈ 127 MB (feasible)
    h = 6.0/48 = 0.125
    PPW = 8.0 points per wavelength (standard for FDTD accuracy)

  Phase 3 (high-res, if needed): N = 64 → 64⁴ = 16,777,216 points
    Memory: 3 arrays × 8 bytes × 16.8M ≈ 403 MB (feasible on VPS)
    h = 6.0/64 = 0.09375
    PPW = 10.7 (high accuracy)

Source-to-grid interpolation:
  Sources rarely fall exactly on grid nodes.
  Use TRILINEAR INTERPOLATION in 4D (16-point stencil):
  distribute source amplitude to the 16 nearest grid nodes
  weighted by proximity.

Runtime estimate (Phase 1, N=32, 100 periods):
  dt ≈ 0.03 (from CFL with c=1.8)
  Steps per period: T/dt = 1.0/0.03 ≈ 33
  Total steps: 100 × 33 = 3,300
  Time per step: ~1M multiply-adds ≈ 1ms on VPS
  Total: ~3 seconds per run, ~40 seconds for 13-run sweep
  (This is very fast. Phase 2 at N=48 is ~25 seconds per run.)


================================================================================
10. FRAMERATE SWEEP                                     [Fixes audit #16]
================================================================================

The sweep tests which value of c produces the most organized
interference pattern in the 4D domain.

  c_4D values: 1.500, 1.525, 1.550, ..., 1.775, 1.800 (13 runs)

Key candidates:
  1.625 = 13/8 (Fibonacci prediction, derived from pair {5,8})
  1.707 = 1 + √2/2 (HEURISTIC candidate, not derived from theory;     [#16]
           related to 24-cell 45° stagger but connection is ad-hoc;
           included in sweep for completeness, NOT as a prediction)
  1.700 = Steinberg central measurement

For each c_4D value, measure:
  - Total field energy vs time (stability check)
  - Interference pattern coherence (spatial autocorrelation)
  - Symmetry order parameters (adapted Steinhardt Q_l for 4D)
  - Peak-to-noise ratio in intensity field
  - Energy balance (injected vs field vs overflow to 5D)


================================================================================
11. WHAT THE ENGINE DOES NOT INCLUDE                    [Fixes audit #17]
================================================================================

The engine does NOT:
  - Impose phi — phi regulation of framerate is a PREDICTION to be
    tested, not an input. If phi emerges in the interference pattern
    or the optimal framerate, that is an observation.                    [#17b]
  - Add nonlinear field interactions beyond the V(ψ,A) potential
  - Add damping terms (beyond PML boundaries)
  - Add noise terms
  - Nudge the decoherence ratio toward any preferred value
  - Pre-select symmetry outcomes
  - Assume matter/antimatter interpretation

The engine DOES include (per source documents):
  - The potential term ω₀²ψ(1+A/A_ref) from the Lagrangian            [#1]
  - Position-dependent decoherence via curvature κ(x)                   [#7]
  - Dimensional overflow to 5D when r→0.5                              [#17a]
  - Smooth source windowing to avoid Gibbs artifacts                    [#15]

DEFERRED TO LATER PHASES:
  - CN(x) coupling (Phase 2, after structure forms)
  - Full 5D simulation (only thin 5D overflow layer in Phase 1)


================================================================================
12. MEASUREMENTS AND OUTPUTS
================================================================================

Per sweep run, output:
  1. Energy time series: E_field(t), E_injected(t), E_overflow(t)
  2. Final intensity field: ⟨|ψ|²⟩ over last 10 periods
  3. 3D cross-sections at w=0 (midplane of 4th dimension)
  4. Symmetry analysis: adapted Q_l order parameters
  5. Peak positions: local maxima of intensity field
  6. Coherence metric: spatial autocorrelation of intensity
  7. Dual vertex check: intensity at Form B positions (not sources)
  8. Energy partition: fraction in each dual orientation's Voronoi cell

Save raw data as JSON + numpy arrays for post-analysis.


================================================================================
13. AUDIT CHECKLIST v2                                  [Fixes audit #18]
================================================================================

Equation fidelity:
  □ Wave equation matches eq. (1) INCLUDING potential term
  □ Potential term = ω₀²ψ(1+A/A_ref), not simplified
  □ Source term matches eq. (3) with smooth windowing eq. (4)
  □ 4D Laplacian uses ∇² notation, implemented as eq. (2)
  □ C_potential uses CURVATURE κ(x) per eq. (7-8), not raw energy density

Geometry:
  □ Form A vertices correct: (±1,0,0,0) ∪ (±½,±½,±½,±½) = 24 points
  □ Form B vertices correct: perms of (±1,±1,0,0) / √2 = 24 points
  □ Both normalized to SAME radius (unit sphere)
  □ Primary test: Form A sources only

Numerics:
  □ CFL condition: dt = 0.9 × h / (2 × c_max), c_max = 1.8
  □ Grid spacing h uniform in all 4 dimensions
  □ Source-to-grid: 4D trilinear interpolation (16-point stencil)
  □ PML absorbing boundaries, N_pml = 10, σ profile specified
  □ Windowing ramp: N_ramp = max(4, 0.05 × T/dt)

Physics:
  □ No damping beyond PML
  □ No noise
  □ No symmetry forcing
  □ No phi imposed
  □ 5D overflow layer active (N_5 = 8, c_5D = 2.625)
  □ r(x) hard clamped at 0.5
  □ CN coupling disabled for Phase 1 (A = A_base uniform)

Outputs:
  □ Energy balance tracked (injected - field - overflow - boundary)
  □ All 13 sweep values produce output files
  □ 3D cross-sections saved for visualization
  □ Dual vertex intensities recorded

Units:
  □ Dimensionless: c_3D = 1, lambda_0 = 1, f_0 = 1
  □ All sweep values are ratios relative to c_3D

================================================================================