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:

$$\frac{dq_i}{dt} = \frac{\partial H}{\partial p_i}, \quad \frac{dp_i}{dt} = -\frac{\partial H}{\partial 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 = \sum_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 + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{align}$$

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

$$M_{ ext{RK4}}^T J M_{ ext{RK4}} eq 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 - \frac{h}{2} \frac{\partial H}{\partial q}(q_n) \\ q_{n+1} &= q_n + h \frac{\partial H}{\partial p}(p_{n+1/2}) \\ p_{n+1} &= p_{n+1/2} - \frac{h}{2} \frac{\partial H}{\partial 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_{ ext{total}} \sim O(h^2)$$

Not $O(h^2 \cdot 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} \circ \Phi_{c_2 h} \circ \Phi_{c_3 h} \circ \Phi_{c_4 h} \circ \Phi_{c_3 h} \circ \Phi_{c_2 h} \circ \Phi_{c_1 h}$$

where the coefficients satisfy:

$$c_1 = c_7 = \frac{1}{2(2-2^{1/3})}, \quad c_2 = c_6 = \frac{1}{2(2-2^{1/3})}, \quad 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) \sim \delta r_0 e^{\lambda t}$$

The Lyapunov exponent $\lambda > 0$ quantifies this divergence rate:

$$\lambda = \lim_{t o \infty} \frac{1}{t} \ln\left(\frac{|\delta r(t)|}{|\delta r_0|} ight)$$

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:

$$\frac{d}{dt}\delta r = abla_r \vec{v} \cdot \delta r$$

where $\vec{v} = \dot{r}$ is the velocity field. The Jacobian $
abla_r \vec{v}$ evolves the tangent space.

Over time $T$:

$$\lambda = \frac{1}{T} \ln\left(\frac{|\delta r(T)|}{|\delta r(0)|} ight)$$

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

### Magnetic Focusing Systems

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

$$m\frac{d^2x}{dt^2} = q(\vec{v} imes \vec{B})_x + F_{ ext{edge}}$$

For the ideal quadrupole $\vec{B} = B'(y\hat{z} - z\hat{y})$:

$$\frac{d^2x}{dt^2} = \frac{qB'}{m}\left(v_y z - v_z y ight)$$

This is Hamiltonian with $H = \frac{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 = \frac{1}{T}\ln\left|\frac{\delta r(T)}{\delta r_0} ight|$$

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:

$$\int\int dq \, dp = ext{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