Free Particle and Dipole Interaction

This example demonstrates a free charged particle responding to an oscillating dipole field. The dipole (charges at \(z = \pm 0.5\) nm) starts oscillating at \(t=0\). The particle, located 1500 nm away, experiences:

1. Initial motion: Static electric field pushes particle downward (\(-z\)).

2. Propagation delay: Oscillating radiation travels at speed of light \(c\). Delay time \(= d/c \approx 0.5\) fs.

3. Oscillations begin: After delay, particle responds to time-varying electromagnetic wave with sinusoidal motion.

Import necessary libraries

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
from scipy.constants import c, e, m_e  # type: ignore

from pycharge import dipole_source, free_particle_source, simulate

jax.config.update("jax_enable_x64", True)

Define sources and simulation timeline

# Oscillating dipole starting at t=0 with charges at z=±0.5 nm
initial_moment = [0.0, 0.0, 1e-9]
q_dipole = e * 50
omega_0 = 1e15 * 2 * jnp.pi  # 1 PHz
m_dipole = m_e
dipole_origin = [0.0, 0.0, 0.0]

dipole = dipole_source(initial_moment, omega_0, dipole_origin, q_dipole, m_dipole)

# Free particle at rest, 1500 nm away along x
particle_position = [15e-7, 0.0, 0]
q_particle = e
m_particle = m_e


def position_0_fn(t):
    """Initial position of the free particle."""
    return jnp.array(particle_position)


free_particle = free_particle_source(position_0_fn, q_particle, m_particle)

# Calculate light propagation delay
distance = jnp.linalg.norm(jnp.array(particle_position) - jnp.array(dipole_origin))
light_travel_time = distance / c

print(f"Distance from dipole to particle: {distance * 1e9:.1f} nm")
print(f"Light travel time: {light_travel_time * 1e15:.3f} fs")

# Simulation parameters
num_steps = 20_000
dt = 1e-18  # 1 attosecond
ts = jnp.arange(num_steps) * dt

Distance from dipole to particle: 1500.0 nm
Light travel time: 5.003 fs

Run the simulation

sim_fn = jax.jit(simulate([dipole, free_particle], ts))
states = sim_fn()

Timestep 0
Timestep 100
Timestep 200
Timestep 300
Timestep 400
Timestep 500
Timestep 600
Timestep 700
Timestep 800
Timestep 900
Timestep 1000
Timestep 1100
Timestep 1200
Timestep 1300
Timestep 1400
Timestep 1500
Timestep 1600
Timestep 1700
Timestep 1800
Timestep 1900
Timestep 2000
Timestep 2100
Timestep 2200
Timestep 2300
Timestep 2400
Timestep 2500
Timestep 2600
Timestep 2700
Timestep 2800
Timestep 2900
Timestep 3000
Timestep 3100
Timestep 3200
Timestep 3300
Timestep 3400
Timestep 3500
Timestep 3600
Timestep 3700
Timestep 3800
Timestep 3900
Timestep 4000
Timestep 4100
Timestep 4200
Timestep 4300
Timestep 4400
Timestep 4500
Timestep 4600
Timestep 4700
Timestep 4800
Timestep 4900
Timestep 5000
Timestep 5100
Timestep 5200
Timestep 5300
Timestep 5400
Timestep 5500
Timestep 5600
Timestep 5700
Timestep 5800
Timestep 5900
Timestep 6000
Timestep 6100
Timestep 6200
Timestep 6300
Timestep 6400
Timestep 6500
Timestep 6600
Timestep 6700
Timestep 6800
Timestep 6900
Timestep 7000
Timestep 7100
Timestep 7200
Timestep 7300
Timestep 7400
Timestep 7500
Timestep 7600
Timestep 7700
Timestep 7800
Timestep 7900
Timestep 8000
Timestep 8100
Timestep 8200
Timestep 8300
Timestep 8400
Timestep 8500
Timestep 8600
Timestep 8700
Timestep 8800
Timestep 8900
Timestep 9000
Timestep 9100
Timestep 9200
Timestep 9300
Timestep 9400
Timestep 9500
Timestep 9600
Timestep 9700
Timestep 9800
Timestep 9900
Timestep 10000
Timestep 10100
Timestep 10200
Timestep 10300
Timestep 10400
Timestep 10500
Timestep 10600
Timestep 10700
Timestep 10800
Timestep 10900
Timestep 11000
Timestep 11100
Timestep 11200
Timestep 11300
Timestep 11400
Timestep 11500
Timestep 11600
Timestep 11700
Timestep 11800
Timestep 11900
Timestep 12000
Timestep 12100
Timestep 12200
Timestep 12300
Timestep 12400
Timestep 12500
Timestep 12600
Timestep 12700
Timestep 12800
Timestep 12900
Timestep 13000
Timestep 13100
Timestep 13200
Timestep 13300
Timestep 13400
Timestep 13500
Timestep 13600
Timestep 13700
Timestep 13800
Timestep 13900
Timestep 14000
Timestep 14100
Timestep 14200
Timestep 14300
Timestep 14400
Timestep 14500
Timestep 14600
Timestep 14700
Timestep 14800
Timestep 14900
Timestep 15000
Timestep 15100
Timestep 15200
Timestep 15300
Timestep 15400
Timestep 15500
Timestep 15600
Timestep 15700
Timestep 15800
Timestep 15900
Timestep 16000
Timestep 16100
Timestep 16200
Timestep 16300
Timestep 16400
Timestep 16500
Timestep 16600
Timestep 16700
Timestep 16800
Timestep 16900
Timestep 17000
Timestep 17100
Timestep 17200
Timestep 17300
Timestep 17400
Timestep 17500
Timestep 17600
Timestep 17700
Timestep 17800
Timestep 17900
Timestep 18000
Timestep 18100
Timestep 18200
Timestep 18300
Timestep 18400
Timestep 18500
Timestep 18600
Timestep 18700
Timestep 18800
Timestep 18900
Timestep 19000
Timestep 19100
Timestep 19200
Timestep 19300
Timestep 19400
Timestep 19500
Timestep 19600
Timestep 19700
Timestep 19800
Timestep 19900

Extract and convert trajectories

# Extract z-positions for all charges (in nm)
d1_z = states[0][:, 0, 0, 2] * 1e9  # Dipole charge 1 (-50e)
d2_z = states[0][:, 1, 0, 2] * 1e9  # Dipole charge 2 (+50e)
p_z = states[1][:, 0, 0, 2] * 1e9  # Free particle (+1e)

Plot z-positions for all charges

Dipole charges: Oscillate at \(\omega_0 = 1\) PHz along \(z\)-axis around \(z = \pm 0.5\) nm.

Free particle: Initially drifts downward due to static field. After \(\sim 0.5\) fs (light travel time), sinusoidal oscillations begin as the electromagnetic wave arrives. Red dashed line marks wave arrival.

def plot_z_position(ax, times, z_pos, color, title, vline_time=None):
    """Plot z-position with optional vertical line marker."""
    ax.plot(times, z_pos, color=color, linewidth=1.5)

    # Add vertical line if specified
    if vline_time is not None:
        ax.axvline(vline_time, color="red", linestyle="--", linewidth=1.5, alpha=0.7, label="Light arrival")
        ax.legend(loc="upper right", fontsize=9)

    ax.set_ylabel("z-position (nm)")
    ax.set_xlabel("Time (fs)")
    ax.set_title(title)
    ax.grid(True, which="both", ls=":")


fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Plot z-position for each charge
ts_fs = ts * 1e15  # Time in femtoseconds

# Create titles with charge and position information


def charge_title(label: str, q: float, pos: list) -> str:
    return f"{label} (q={q / e:.0f}e)\n(x={pos[0] * 1e9:.1f} nm, y={pos[1] * 1e9:.1f} nm)"


d1_title = charge_title("Dipole Charge 1", -q_dipole, dipole_origin)
d2_title = charge_title("Dipole Charge 2", q_dipole, dipole_origin)
p_title = charge_title("Free Particle", q_particle, particle_position)

plot_z_position(axes[0], ts_fs, d1_z, "C0", d1_title)
plot_z_position(axes[1], ts_fs, d2_z, "C1", d2_title)
plot_z_position(axes[2], ts_fs, p_z, "C2", p_title, vline_time=light_travel_time * 1e15)

plt.tight_layout()
plt.show()