symplectic integrator hamiltonian chaos particle tracking

# Symplectic Integrator Dynamics, Phase Space Geometry Preservation, and Chaotic Trajectory Kinetics in Particle Beams

## 1. Introduction: Hamilton's Equations and Canonical Structure

Classical mechanics is rooted in Hamilton's equations:

(dq_i)/(dt) = (∂ H)/(∂ p_i), (dp_i)/(dt) = -(∂ H)/(∂ q_i)

where (q_i, p_i) are canonical coordinates and momenta, and H(q, p, t) is the Hamiltonian.

These equations are not merely differential equations—they are structure-preserving mathematical transformations that conserve the fundamental geometric property of phase space: symplecticity.

## 2. Symplectic Structure and Liouville's Theorem

The phase space of a Hamiltonian system has a canonical symplectic form:

Omega = ∑_i dq_i wedge dp_i

This 2-form encodes the fact that Hamiltonian flows are volume-preserving in phase space. Formally, if Phi_t is the flow map at time t, then:

Phi_t^* Omega = Omega

This is Liouville's theorem: the volume of any region in phase space is invariant under Hamiltonian evolution.

### Symplectic Matrices

A linear transformation is symplectic if:

M^T J M = J

where J is the symplectic matrix:

J = begin(pmatrix) 0 & I -I & 0 end(pmatrix)

and I is the identity. Symplectic transformations preserve the symplectic form and thus phase space volume.

## 3. Standard Integrators: Runge-Kutta and Its Failure for Long-Term Dynamics

The classical 4th-order Runge-Kutta method:

begin(align)
k_1 &= h f(t_n, x_n)
k_2 &= h f(t_n + h/2, x_n + k_1/2)
k_3 &= h f(t_n + h/2, x_n + k_2/2)
k_4 &= h f(t_n + h, x_n + k_3)
x_n+1 &= x_n + (1)/(6)(k_1 + 2k_2 + 2k_3 + k_4)
end(align)

is highly accurate for short times (Delta E ~ h^5 per step), but systematically violates symplecticity:

M_RK4^T J M_RK4 ≠ J

Over long trajectories, even minute violations accumulate. The phase space volume drifts by ~1-10% over 1000 steps, and more critically, energy errors accumulate systematically (not randomly cancelled).

For long-term predictions (chaotic beams, ion trap dynamics over microseconds), this is catastrophic.

## 4. Störmer-Verlet Symplectic Integrator

The Störmer-Verlet (or leapfrog) algorithm is explicitly symplectic:

begin(align)
p_n+1/2 &= p_n - (h)/(2) (∂ H)/(∂ q)(q_n)
q_n+1 &= q_n + h (∂ H)/(∂ p)(p_n+1/2)
p_n+1 &= p_n+1/2 - (h)/(2) (∂ H)/(∂ q)(q_n+1)
end(align)

This is a composition of two transformations:
- Position update using half-step momentum
- Momentum update using new position
- Another half-step momentum update

Key property: This is an exact product of two symplectic maps, hence symplectic.

### Error Scaling

The local truncation error is O(h^3) per step, but—crucially—energy conservation is O(h^2) rather than O(h^5) for RK4. The total accumulated error over time T is:

Delta E_total ~ O(h^2)

Not O(h^2 · T/h) = O(hT) as naive analysis would suggest. This is because symplecticity constrains the error structure.

## 5. Higher-Order Symplectic Integrators: Yoshida's Algorithm

Yoshida's 4th-order symplectic integrator uses a composition of Störmer-Verlet steps with carefully chosen coefficients:

Phi = Phi_c_1 h ° Phi_c_2 h ° Phi_c_3 h ° Phi_c_4 h ° Phi_c_3 h ° Phi_c_2 h ° Phi_c_1 h

where the coefficients satisfy:

c_1 = c_7 = (1)/(2(2-2^1/3)), c_2 = c_6 = (1)/(2(2-2^1/3)), c_4 = 1 - 2(c_1 + c_2)

The Yoshida integrator:
- Preserves symplecticity exactly
- Has 4th-order accuracy: local error O(h^5), energy error O(h^4)
- Requires 7 function evaluations per step (vs. 4 for RK4), but superior long-term stability justifies the cost

## 6. Chaotic Dynamics in Hamiltonian Systems: The Lyapunov Exponent

A dynamical system is chaotic if nearby trajectories diverge exponentially:

delta r(t) ~ delta r_0 e^lambda t

The Lyapunov exponent lambda > 0 quantifies this divergence rate:

lambda = lim_t → ∞ (1)/(t) ln((|delta r(t)|)/(|delta r_0|))

For Hamiltonian systems, the phase space structure ensures that if one direction is chaotic (lambda > 0), another direction must be "anti-chaotic" (-lambda), preserving volume.

### Computation via Variational Equations

To compute lambda, integrate the variational equations alongside the original dynamics:

(d)/(dt)delta r = ∇_r vec(v) · delta r

where vec(v) = dot(r) is the velocity field. The Jacobian ∇_r vec(v) evolves the tangent space.

Over time T:

lambda = (1)/(T) ln((|delta r(T)|)/(|delta r(0)|))

## 7. Application to Charged Particle Beams and Ion Traps

### Magnetic Focusing Systems

In a magnetic quadrupole lens, particles undergo complex 3D motion:

m(d^2x)/(dt^2) = q(vec(v) × vec(B))_x + F_edge

For the ideal quadrupole vec(B) = B'(yhat(z) - zhat(y)):

(d^2x)/(dt^2) = (qB')/(m)(v_y z - v_z y)

This is Hamiltonian with H = (1)/(2m)(p_x^2 + p_y^2 + p_z^2) + q Phi(x,y,z).

Over long propagation distances, tracking particle positions requires integrating for 10^3–10^6 steps. Non-symplectic methods accumulate energy errors, causing spurious particle loss or gain.

### Ion Trap Dynamics

In a Paul trap (RF-driven quadrupole), ions oscillate at the RF frequency (~1 MHz) while slowly drifting. The timescale separation demands:
- Accurate harmonic oscillation (requires high-order accuracy)
- Long integration time (requires energy stability)

Symplectic integrators naturally handle both: the 4th-order Yoshida method achieves excellent accuracy even with large timesteps.

## 8. Python Implementation: Symplectic vs. Non-Symplectic Comparison

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, odeint

def hamiltonian_simple(q, p, t=0):
    """
    Simple 1D harmonic oscillator: H = p²/2m + ½mω²q²
    """
    omega = 1.0  # Frequency
    return omega**2 * q  # m = 1

def verlet_symplectic_step(q, p, dqdt_func, dt):
    """
    One Störmer-Verlet symplectic step.
    
    Parameters:
    -----------
    q, p : float
        Position and momentum
    dqdt_func : callable
        Function returning ∂H/∂p (e.g., p/m for quadratic kinetic energy)
    dt : float
        Timestep
    
    Returns:
    --------
    q_new, p_new : float
        Updated position and momentum
    """
    
    # Half-step momentum update
    a = 1.0  # d²q/dt² = -∂²H/∂q² (spring constant = 1)
    p_half = p - 0.5 * dt * a * q
    
    # Full-step position update
    q_new = q + dt * p_half  # dq/dt = ∂H/∂p = p (for m=1)
    
    # Half-step momentum update
    p_new = p_half - 0.5 * dt * a * q_new
    
    return q_new, p_new

def yoshida_symplectic_step(q, p, a_coeff, dt):
    """
    4th-order Yoshida symplectic integrator.
    
    Parameters:
    -----------
    q, p : float
        State
    a_coeff : float
        Spring constant (acceleration = -a * q)
    dt : float
        Timestep
    
    Returns:
    --------
    q_new, p_new : float
        Updated state
    """
    
    c1 = 1.0 / (2.0 * (2.0 - 2.0**(1/3)))
    c2 = 1.0 / (2.0 * (2.0 - 2.0**(1/3)))
    c3 = 1.0 - 2*(c1 + c2)
    c4 = 1.0 - 2*(c1 + c2)
    
    coeffs = [c1, c2, c3, c2, c1]
    
    for c in coeffs:
        # Half-step p update
        p -= 0.5 * dt * c * a_coeff * q
        # Full-step q update
        q += dt * c * p
        # Half-step p update
        p -= 0.5 * dt * c * a_coeff * q
    
    return q, p

def energy_harmonic(q, p, omega=1.0):
    """
    Energy of harmonic oscillator: E = p²/2 + ω²q²/2
    """
    return 0.5 * p**2 + 0.5 * omega**2 * q**2

# Simulation parameters
t_final = 1000  # Long integration time
dt = 0.1

# Initial conditions
q0, p0 = 1.0, 0.0  # Start at maximum displacement

# Time arrays
t_verlet = np.arange(0, t_final, dt)
t_yoshida = np.arange(0, t_final, dt)

# Initialize trajectories
q_verlet = np.zeros_like(t_verlet)
p_verlet = np.zeros_like(t_verlet)
E_verlet = np.zeros_like(t_verlet)

q_yoshida = np.zeros_like(t_yoshida)
p_yoshida = np.zeros_like(t_yoshida)
E_yoshida = np.zeros_like(t_yoshida)

q_verlet[0], p_verlet[0] = q0, p0
q_yoshida[0], p_yoshida[0] = q0, p0
E_verlet[0] = energy_harmonic(q0, p0)
E_yoshida[0] = energy_harmonic(q0, p0)

# Integrate with Verlet
for i in range(1, len(t_verlet)):
    q_verlet[i], p_verlet[i] = verlet_symplectic_step(
        q_verlet[i-1], p_verlet[i-1], lambda x: x, dt
    )
    E_verlet[i] = energy_harmonic(q_verlet[i], p_verlet[i])

# Integrate with Yoshida
for i in range(1, len(t_yoshida)):
    q_yoshida[i], p_yoshida[i] = yoshida_symplectic_step(
        q_yoshida[i-1], p_yoshida[i-1], 1.0, dt
    )
    E_yoshida[i] = energy_harmonic(q_yoshida[i], p_yoshida[i])

# RK4 for comparison (non-symplectic)
def rk4_step(t, q, p, dt):
    def f(t, y):
        q, p = y
        dqdt = p
        dpdt = -q  # Harmonic oscillator
        return [dqdt, dpdt]
    
    sol = solve_ivp(f, [t, t+dt], [q, p], method='RK45', max_step=dt/10)
    return sol.y[0, -1], sol.y[1, -1]

t_rk4 = np.arange(0, min(t_final, 100), dt)  # Limited to t=100 for visualization
q_rk4 = np.zeros_like(t_rk4)
p_rk4 = np.zeros_like(t_rk4)
E_rk4 = np.zeros_like(t_rk4)

q_rk4[0], p_rk4[0] = q0, p0
E_rk4[0] = energy_harmonic(q0, p0)

for i in range(1, len(t_rk4)):
    q_rk4[i], p_rk4[i] = rk4_step(t_rk4[i-1], q_rk4[i-1], p_rk4[i-1], dt)
    E_rk4[i] = energy_harmonic(q_rk4[i], p_rk4[i])

# Plotting
fig, axes = plt.subplots(2, 2, figsize=(13, 10))

# Panel A: Energy conservation (long time)
axes[0, 0].semilogy(t_verlet, np.abs(E_verlet - E_verlet[0]), 'b-', linewidth=2, label='Verlet (symplectic)')
axes[0, 0].semilogy(t_yoshida, np.abs(E_yoshida - E_yoshida[0]), 'g-', linewidth=2, label='Yoshida (4th-order)')
axes[0, 0].semilogy(t_rk4, np.abs(E_rk4 - E_rk4[0]), 'r-', linewidth=2, label='RK4 (non-symplectic)')
axes[0, 0].set_xlabel('Time (arbitrary units)', fontsize=11)
axes[0, 0].set_ylabel('|ΔE| (log scale)', fontsize=11)
axes[0, 0].set_title('Energy Error Accumulation', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=10)
axes[0, 0].grid(True, alpha=0.3, which='both')

# Panel B: Phase space trajectory (early times)
axes[0, 1].plot(q_verlet[:500], p_verlet[:500], 'b-', alpha=0.7, linewidth=1.5, label='Verlet')
axes[0, 1].plot(q_yoshida[:500], p_yoshida[:500], 'g--', alpha=0.7, linewidth=1.5, label='Yoshida')
axes[0, 1].plot(q_rk4[:500], p_rk4[:500], 'r:', alpha=0.7, linewidth=2, label='RK4')
# True circle
theta = np.linspace(0, 2*np.pi, 100)
axes[0, 1].plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3, linewidth=2, label='True orbit')
axes[0, 1].set_xlabel('Position q', fontsize=11)
axes[0, 1].set_ylabel('Momentum p', fontsize=11)
axes[0, 1].set_title('Phase Space Trajectory (t=0 to 50)', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].set_xlim([-1.5, 1.5])
axes[0, 1].set_ylim([-1.5, 1.5])

# Panel C: Position vs time
axes[1, 0].plot(t_verlet, q_verlet, 'b-', alpha=0.7, label='Verlet')
axes[1, 0].plot(t_yoshida, q_yoshida, 'g--', alpha=0.7, label='Yoshida')
axes[1, 0].plot(t_rk4, q_rk4, 'r:', alpha=0.7, linewidth=2, label='RK4')
axes[1, 0].set_xlabel('Time', fontsize=11)
axes[1, 0].set_ylabel('Position q(t)', fontsize=11)
axes[1, 0].set_title('Position Evolution (Symplectic vs Non-Symplectic)', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=10)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_xlim([0, 100])

# Panel D: Error scaling law verification
dt_array = np.array([0.01, 0.02, 0.05, 0.1, 0.2])
error_verlet = []
error_yoshida = []

for dt_test in dt_array:
    n_steps = int(100 / dt_test)
    q_v, p_v = q0, p0
    E_v_init = energy_harmonic(q_v, p_v)
    
    for _ in range(n_steps):
        q_v, p_v = verlet_symplectic_step(q_v, p_v, lambda x: x, dt_test)
    E_v_final = energy_harmonic(q_v, p_v)
    error_verlet.append(np.abs(E_v_final - E_v_init))
    
    # Yoshida
    q_y, p_y = q0, p0
    E_y_init = energy_harmonic(q_y, p_y)
    for _ in range(n_steps):
        q_y, p_y = yoshida_symplectic_step(q_y, p_y, 1.0, dt_test)
    E_y_final = energy_harmonic(q_y, p_y)
    error_yoshida.append(np.abs(E_y_final - E_y_init))

axes[1, 1].loglog(dt_array, error_verlet, 'b-o', linewidth=2, markersize=8, label='Verlet (O(h²))')
axes[1, 1].loglog(dt_array, error_yoshida, 'g-s', linewidth=2, markersize=8, label='Yoshida (O(h⁴))')
axes[1, 1].loglog(dt_array, dt_array**2, 'k--', alpha=0.5, linewidth=1.5, label='O(h²) reference')
axes[1, 1].loglog(dt_array, dt_array**4, 'k:', alpha=0.5, linewidth=1.5, label='O(h⁴) reference')
axes[1, 1].set_xlabel('Timestep h (log scale)', fontsize=11)
axes[1, 1].set_ylabel('Total Energy Error (log scale)', fontsize=11)
axes[1, 1].set_title('Error Scaling: Symplectic Integrators', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('symplectic_integrator_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

print("=== Symplectic Integrator Comparison Complete ===")
print(f"Final energy error (Verlet): {np.abs(E_verlet[-1] - E_verlet[0]):.2e}")
print(f"Final energy error (Yoshida): {np.abs(E_yoshida[-1] - E_yoshida[0]):.2e}")
print(f"Verlet integrator is symplectic: energy error scales as O(h²)")
print(f"Yoshida integrator is symplectic: energy error scales as O(h⁴)")

## 9. Chaos Detection: Lyapunov Exponent Calculation

Computing Lyapunov exponents reveals chaotic dynamics:

lambda = (1)/(T)ln|(delta r(T))/(delta r_0)|

For an ion in a Paul trap undergoing resonant heating, lambda > 0 signals chaotic motion and energy diffusion.

## 10. Phase Space Volume Conservation

Unlike non-symplectic methods, symplectic integrators maintain:

∫∫ dq dp = constant

This is verified numerically by tracking an ensemble of particles: their phase space area remains constant.

## 11. Applications to Plasma Particle-in-Cell Simulations

In PIC simulations with 10^6–10^9 particles, symplectic methods prevent:
- Spurious heating (non-physical energy gain)
- Numerical instabilities
- Long-term validity of temperature and pressure calculations

## 12. Limitations and Higher-Order Methods

Even symplectic integrators have limitations:
- Reversibility: True time-reversibility only in exact arithmetic
- Long-term error: Still accumulates, though at controlled rate
- Stiff systems: Large timestep may violate accuracy for fast oscillations

Forest-Ruth 4th-order and Omelyan-Folland-Kusnezov 6th-order symplectic methods achieve higher accuracy at increased computational cost.

## 13. Molecular Dynamics and Statistical Mechanics

Symplectic integrators are standard in MD simulations for thermostats (Nosé-Hoover, Lennard-Jones). Preserving phase space volume ensures proper ensemble statistics (NVE, NVT).

## 14. Comparison Table: Integration Methods

## 15. Future Directions: Adaptive and Geometric Integration

Emerging techniques:
- Adaptive symplectic integrators: Variable timestep maintaining symplecticity
- Splitting methods: Decompose Hamiltonian into integrable parts
- Machine learning-augmented: Neural networks predict next state while preserving symplecticity constraints