Theory primer

This page distills the physics encoded in the LW Integrator and mirrors the notation used in the internal design notes that accompany the project. It connects the covariant Liénard–Wiechert formalism to the concrete data structures exposed in core/trajectory_integrator.py and the validation studies under examples/validation.

Note that throughout the codebase Gaussian units are used. The unorthodox choice of amu-millimeter-nanosecond units is implemented to avoid numerical overflows.

Retarded fields

The solver models every source particle as a point charge whose fields are sampled at the observer’s retarded time. Starting from Jackson’s form of the Liénard–Wiechert fields, the magnetic field is obtained from the electric field via a cross product, while the electric field splits into a velocity term and an acceleration term:

\[\mathbf{B} = \bigl[\mathbf{n} \times \mathbf{E}\bigr]_{\text{ret}},\]

\[\mathbf{E} = e\left[\frac{\mathbf{n} - \boldsymbol{\beta}}{\gamma^{2}\,\kappa^{3} R^{2}}\right]_{\text{ret}} + \frac{e}{c} \left[ \frac{\mathbf{n} \times \bigl((\mathbf{n} - \boldsymbol{\beta}) \times \dot{\boldsymbol{\beta}}\bigr)}{\kappa^{3} R} \right]_{\text{ret}},\]

where \(\kappa = 1 - \boldsymbol{\beta} \cdot \mathbf{n}\), \(R\) is the retarded source–observer separation, \(\boldsymbol{\beta} = \mathbf{v}/c\), and \(\gamma = (1-\beta^{2})^{-1/2}\). Each quantity is evaluated at the retarded time \(t - R/c\). The implementation samples these terms inside core.trajectory_integrator.retarded_integrator(), looping over all available source trajectories.

A key limit for the benchmark problems is the near head-on configuration where \(\mathbf{n}\) aligns with \(\boldsymbol{\beta}\). Neglecting transverse components yields

(1)\[\mathbf{E}_{\parallel} \approx e\,\frac{1-\beta}{(1+\beta) R^{2}}\,\mathbf{n},\]

which explains the asymptotic growth of the longitudinal field as \(\beta \rightarrow 1\). This is the regime probed by the aperture-loss studies and the recoil-reduction scenarios in examples/validation.

Covariant potentials

Instead of tracking fields directly, the integrator evolves the covariant potential \(A^{\alpha}\) for each source trajectory. Using proper time \(\tau\) as the integration variable, the retarded potential reads

(2)\[A^{\alpha}(x) = \left.\frac{e\, V^{\alpha}(\tau)}{V(\tau) \cdot [x - r(\tau)]}\right|_{\tau = \tau_{0}},\]

with \(V^{\alpha} = \{c\gamma, \gamma \mathbf{u}\}\) the four-velocity, \(r^{\alpha}(\tau)\) the source worldline, and \(\tau_{0}\) obtained from light-cone constraint \([x - r(\tau_{0})]^{2} = 0\). The denominator reduces to \(\gamma c R \kappa\), linking the potential back to the geometry used in (1).

Conjugate momentum and equations of motion

Each observer particle carries a conjugate four-momentum

(3)\[\mathcal{P}^{\alpha} = m V^{\alpha} + \frac{e}{c} A^{\alpha},\]

where \(m\) and \(e\) are the observer mass and charge. Differentiating \(\mathcal{P}^{\alpha}\) with respect to proper time leads to the mixed-field force law used inside the stepping kernel:

(4)\[\frac{d\mathcal{P}^{\alpha}}{d\tau} = \frac{e}{c} V_{\beta} \, \partial^{\alpha} A^{\beta}.\]

Expanding \(\partial^{\alpha} A^{\beta}\) in terms of \(V^{\alpha}\), \(R^{\alpha}\), \(\dot{V}^{\alpha}\), and \(\kappa\) yields the component-wise form implemented in core.trajectory_integrator._update_conjugate_momentum. The spatial components couple velocity, acceleration, and retarded distance, ensuring that head-on image-charge interactions reproduce the steep gradients reported in the reference study.

Position updates follow directly from the Hamiltonian identity

(5)\[\frac{d x^{\alpha}}{d\tau} = \frac{1}{m}\left( \mathcal{P}^{\alpha} - \frac{e}{c} A^{\alpha} \right),\]

which the solver evaluates after each momentum update to keep particle states in sync. Proper-time stepping avoids runaway behaviour at high \(\gamma\) while keeping the integration scheme close to the historical reference implementation.

Relativistic position updates in coordinate time

While the covariant formulation uses proper time \(d\tau\), the numerical implementation steps forward in coordinate time with interval \(h = \Delta t\). The spatial position update relates proper-time and coordinate-time derivatives:

\[\Delta \mathbf{x} = \mathbf{v} \, \Delta t = \frac{\mathbf{P}_{\text{kinetic}}}{\gamma m} \, h,\]

where \(\mathbf{P}_{\text{kinetic}} = \mathcal{P} - (e/c)\mathbf{A}\) is the kinetic (mechanical) momentum. The crucial \(1/\gamma\) factor ensures that velocity \(\mathbf{v} = \mathbf{P}_{\text{kinetic}}/(\gamma m)\) remains subluminal even as momentum grows with \(\gamma\).

The corresponding velocity is then computed from the coordinate-time displacement:

\[\boldsymbol{\beta} = \frac{\mathbf{v}}{c} = \frac{\Delta \mathbf{x}}{c \, h}.\]

Note that this formula does not include a \(\gamma\) factor in the denominator—the time dilation is already accounted for in the position update.

Self-consistency iterations

For ultra-relativistic particles (\(\gamma \gg 1\)), forces depend strongly on \(\gamma\) through the retarded field geometry (via \(\kappa\) and field Lorentz contraction). This creates a circular dependency:

\[\gamma \rightarrow \text{forces} \rightarrow \mathcal{P} \rightarrow \gamma.\]

The integrator resolves this through self-consistency iterations at each timestep. Within iteration \(n\):

1. Use \(\gamma_{n-1}\) from the previous iteration to compute retarded forces

2. Update conjugate momentum \(\mathcal{P}_{n}\) from those forces

3. Compute positions using the same \(\gamma_{n-1}\): \(\Delta \mathbf{x} = (\mathbf{P}_{\text{kinetic}}/(\gamma_{n-1} m)) h\)

4. Compute velocity: \(\boldsymbol{\beta}_{n} = \Delta \mathbf{x}/(c h)\)

5. Derive two independent estimates of \(\gamma_{n}\):

   • From energy: \(\gamma_{\text{E}} = (\mathcal{P}^{0} - e\Phi)/(mc)\)

   • From velocity: \(\gamma_{\text{V}} = 1/\sqrt{1 - \beta^{2}}\)

6. Check convergence: \(|\gamma_{\text{E}} - \gamma_{\text{V}}|/\gamma_{\text{E}} < \epsilon\)

If not converged, iteration \(n+1\) uses \(\gamma_{n} = \gamma_{\text{E}}\) and repeats. Typical tolerance \(\epsilon = 10^{-6}\) achieves convergence within 1–3 iterations even after large energy jumps.

The key to stable convergence is using a consistent \(\gamma\) throughout each iteration for both force calculation and position updates, ensuring that the velocity extracted from \(\Delta \mathbf{x}\) corresponds physically to the momentum computed from those forces.

Implementation details are in core.self_consistency.SelfConsistencyConfig and core.equations.retarded_equations_of_motion().

Radiation pressure and reaction

The validation notebooks explore scenarios where residual fields act on a test particle once a conducting surface or driving bunch is withdrawn. Two secondary forces are monitored to confirm that their contribution is negligible for the reported configurations:

• Radiation pressure. Using Jackson’s scaling, the momentum transfer to an

  observer with area \(a_{T}\) receiving power \(P_{R}\) across solid angle \(\Omega\) is \(\dot{P}_{\text{RP}} = (P_{R}/c)\,(a_{T}/\Omega R^{2})\). For the millimetre-to-micron geometries in this repository, this quantity is orders of magnitude smaller than the Lorentz force recovered from (1).

• Radiation reaction. Passive Liénard radiation diagnostics are tracked separately from any optional self-force. The provisional power_matched_damping mode removes the computed radiated energy by scaling mechanical momentum magnitude; it is an energy-bookkeeping approximation, not a LAD model. The experimental medina_lad mode applies Medina’s reduced-order Lorentz–Abraham–Dirac force to mechanical momentum:

  \[\mathbf{F}_{\text{rad}} = \frac{2}{3}\frac{e^{2}}{m c^{3}}\left[\frac{d\gamma}{dt}\,\mathbf{F}_{\text{ext}} - \frac{\gamma^{3}}{c^{2}} (\mathbf{F}_{\text{ext}} \cdot \mathbf{a})\, \mathbf{v}\right].\]

  This mode is opt-in and currently validated only against controlled prescribed-field cases. Longitudinal acceleration should show the expected near-cancellation of the Medina terms, while transverse bending gives the synchrotron-style recoil mostly opposite the particle velocity. Conducting boundary cases still require dedicated convergence checks before medina_lad should be treated as physics evidence.

COLD_START gating mechanism

The COLD_START startup mode suppresses retarded forces during the initial phase of a simulation until particles have traveled far enough from their origin for light from external sources to have causally reached them. This ensures physical causality is respected and avoids applying forces based on incomplete or unphysical retarded histories.

Gating threshold formula

The threshold distance a particle must travel before forces are applied is:

\[d_{\text{threshold}} = \frac{\beta \cdot R}{1 - \boldsymbol{\beta} \cdot \mathbf{n}},\]

where:

• \(\beta = |\boldsymbol{\beta}|\) is the particle speed (in units of c)

• \(R\) is the current distance from particle to external source

• \(\mathbf{n}\) is the unit vector pointing from source to particle

• \(\boldsymbol{\beta} \cdot \mathbf{n}\) is the velocity component along the separation

Physical interpretation

The threshold represents the distance the particle must travel before light from the source’s initial position could reach it. The formula accounts for the relative closing speed between light and particle:

• Approaching (\(\boldsymbol{\beta} \cdot \mathbf{n} < 0\)): particles and light meet quickly, threshold < \(\beta R\)

  Example: \(\boldsymbol{\beta} \cdot \mathbf{n} = -\beta\) (head-on) → threshold = \(\beta R / (1 + \beta)\)

• Perpendicular (\(\boldsymbol{\beta} \cdot \mathbf{n} = 0\)): light travels full distance, threshold = \(\beta R\)

• Receding (\(\boldsymbol{\beta} \cdot \mathbf{n} > 0\)): light takes longer to catch up, threshold > \(\beta R\)

• Receding at c (\(\boldsymbol{\beta} \cdot \mathbf{n} \rightarrow 1\)): light never catches particle, threshold → ∞ (forces never applied)

Dynamic threshold calculation

The threshold is recalculated every integration step using current values of \(R\) and \(\boldsymbol{\beta}\). This ensures causality is respected as geometry evolves:

1. Distance \(R\) changes as particles move and images reposition

2. Velocity \(\boldsymbol{\beta}\) updates as particles accelerate

3. Threshold automatically decreases as particle approaches sources

4. Forces “turn on” dynamically when travel distance exceeds threshold

Two-stage implementation

For computational efficiency, the code uses a two-stage check:

Stage 1: Early conservative check (performance optimization)

• Uses estimated maximum \(R\) from external particle bounds

• Conservative threshold: \(\beta R_{\max} / (1 + \beta)\) (assumes head-on approach)

• If travel distance < estimated threshold, skip expensive retarded distance calculations

Stage 2: Precise check (accurate gating)

• Uses actual retarded distance \(R\) to each external source

• Per-source thresholds: \(\beta R / (1 - \boldsymbol{\beta} \cdot \mathbf{n})\)

• Forces applied when travel distance ≥ threshold

Velocity regime behavior

The \(\beta\) factor in the numerator is critical for non-relativistic particles:

Regime	\(\beta\)	Threshold (approaching)	Physical meaning
Non-relativistic	0.01	\(0.01 R / 1.01 \approx 0.01 R\)	Forces apply almost immediately
Low velocity	0.1	\(0.1 R / 1.1 \approx 0.09 R\)	Early force application
Moderate	0.5	\(0.5 R / 1.5 \approx 0.33 R\)	Forces apply at 1/3 distance
Relativistic	0.9	\(0.9 R / 1.9 \approx 0.47 R\)	Near half-distance
Ultra-relativistic	→ 1	→ \(R / 2\)	Approach halfway limit (never exceed)

Without the \(\beta\) factor, low-velocity particles would have forces suppressed for distances far exceeding physical interaction regions, producing incorrect physics.

Implementation

The gating mechanism is implemented in core/equations.py:

• _compute_gating_threshold() computes per-source thresholds

• _should_apply_external_forces() performs the Stage 2 check

• Early conservative check embedded in retarded_equations_of_motion()

For conducting walls, the distance \(R\) is computed to image charges (not the wall itself), ensuring correct handling of virtual source positions.

Bridging back to the code

The mathematical relationships above surface in the codebase as follows:

• core.trajectory_integrator.IntegratorConfig captures the physical parameters (\(\Delta\tau\), aperture radius, wall position) implied by the analytical terms.

• core.trajectory_integrator.generate_conducting_image() and core.trajectory_integrator.generate_switching_image() encode the boundary conditions assumed when taking the head-on limit to model conducting apertures and switching walls.

• The reference notebooks under examples/validation/ document historical comparisons and exploratory studies related to (1).

• The notebooks under legacy/ are retained for historical investigations, but they are not part of the modern retarded-field workflows.

For deeper derivations and experimental context, see https://doi.org/10.1016/j.nima.2024.169988.