The Classical Helium Atom in Two Dimensions

---

Introduction

In this notebook, we explore the classical two-dimensional (2D) helium atom. The model is a simplified three-body problem consisting of a stationary nucleus of charge $+2$ at the origin and two electrons of charge $-1$. We assume each electron has a mass of 1 in these units, and the electrostatic force between them is given by: $$ F = \frac{1}{r^2}, $$ where $r$ is the distance between the two charges.

The equations of motion for the two electrons ($\mathbf{r}_1$ and $\mathbf{r}_2$) in vector form are:

$$ \mathbf{a}_1 = -2 \cdot \frac{\mathbf{r}_1}{\|\mathbf{r}_1\|^3} \;+\; \frac{\mathbf{r}_1 - \mathbf{r}_2}{\|\mathbf{r}_1 - \mathbf{r}_2\|^3}, $$ $$ \mathbf{a}_2 = -2 \cdot \frac{\mathbf{r}_2}{\|\mathbf{r}_2\|^3} \;+\; \frac{\mathbf{r}_2 - \mathbf{r}_1}{\|\mathbf{r}_2 - \mathbf{r}_1\|^3}. $$

Here,

• $\mathbf{r}_1$ and $\mathbf{r}_2$ are the position vectors of the electrons relative to the nucleus,
• $\mathbf{a}_1$ and $\mathbf{a}_2$ are their respective accelerations,
• The factor of $-2$ in each expression comes from the nucleus (of charge $+2$) attracting each electron,
• The remaining terms represent the inter-electron repulsion.

Because one electron can sometimes approach the nucleus very closely, the acceleration can become very large, requiring very small time steps for stable numerical integration. Rather than use a single, uniformly small time step, we employ an adaptive time-stepping method (specifically, an embedded Runge-Kutta 4(5) scheme, RK45) that refines the time step automatically based on an error tolerance.

The goal of the simulations is to investigate:

1. Trajectories of the two electrons under various initial conditions.
2. Autoionization, where one electron can escape to large distances.
3. Conservation of Angular Momentum, to evaluate numerical accuracy.
4. Chaotic behavior, by studying how small changes in initial conditions can drastically change the resulting trajectories and ionization times.

---

Methods

Overview of Numerical Methods

1. Equations of Motion
   We model the helium atom by integrating the coupled ordinary differential equations (ODEs) for the two electrons: $$ \frac{d\mathbf{r}_i}{dt} = \mathbf{v}_i, \quad \frac{d\mathbf{v}_i}{dt} = \mathbf{a}_i, $$ for $i = 1,2$. The function helium_derivs in the code returns the time derivatives of the state vector $$ y = [x_1,\, y_1,\, v_{x_1},\, v_{y_1},\, x_2,\, y_2,\, v_{x_2},\, v_{y_2}]. $$

2. Adaptive Step Size (RK45)
   We use an embedded Runge-Kutta 4(5) method, which provides both a 4th- and 5th-order estimate of the solution at each step. By comparing these two estimates, we can estimate the local truncation error and adapt the time step accordingly. The implementation is in the RK45 function. The user sets tolerances rtol (relative tolerance) and atol (absolute tolerance), and the algorithm refines or coarsens the time step to keep the error below these thresholds.

3. Implementation Details

   • solve_ode_RK45: A driver function that repeatedly calls RK45, storing the state and time in arrays for later analysis.
   • simulate_until_ionization: A helper function that runs the simulation until a maximum time (tf) or until one electron’s distance from the origin exceeds 6. This latter condition is taken to indicate ionization.

4. Verification and Reproducibility

   • We monitor angular momentum to assess the accuracy of our method. Because we ignore recoil of the nucleus, the total angular momentum should be conserved in theory. Observing small variations in total angular momentum over time indicates the level of numerical precision.
   • To reproduce the experiments, one can simply run the provided cells in sequence. All parameters (initial conditions, step sizes, tolerances, and maximum times) are clearly indicated in the code.

---

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('dark_background')

def RK45(dt, f, t, y, rtol, atol, args):
    # atol = 1e-10
    # rtol = 1e-14

    c_params = np.array([0, 1/5, 3/10, 4/5, 8/9, 1, 1])
    a_params = np.array([[0, 0, 0, 0, 0, 0],
                          [1/5, 0, 0, 0, 0, 0],
                          [3/40, 9/40, 0, 0, 0, 0],
                          [44/45, -56/15, 32/9, 0, 0, 0],
                          [19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0],
                          [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0],
                          [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]])
    b_params = np.array([35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0])
    b_star_params = np.array([5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40])

    k = np.zeros((7, len(y)))

    k[0] = f(t, y, *args) # first step
    for i in range(1, 7): # iterate over the rest of the steps (1 to 6)
        y_sum = np.zeros_like(y) # this will be our sum of y's over each k step
        for j in range(i): # sum over all previous k's
            y_sum += a_params[i, j] * k[j] 
        k[i] = f(t + c_params[i] * dt, y + dt * y_sum, *args) # calculate the next k
    y_new = y + dt * np.dot(b_params, k) # calculate the new y
    y_star = y + dt * np.dot(b_star_params, k) # y_star is used to estimate the error

    scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol # scale is used to calculate the error
    error = np.linalg.norm((y_new - y_star) / scale) # error is described as a norm in the book
    # print(f'dt: ', dt)
    # print(f'error: ', error)
    # delta = np.abs(y_new - y_star) <= scale
    if error < 1:
        return y_new # if the error is small enough, we return the new y
    else:
        return RK45(dt/2, f, t, y,rtol,atol, args) # if the error is too large, we half the step size and try again

def solve_ode_RK45(f, tspan, y0, method=RK45, args=(), **options):
    t0 = tspan[0]
    tf = tspan[1]
    y = [y0]
    t = [t0]
    dt = options.get('first_step')
    rtol = options.get('rtol', 1e-6)
    atol = options.get('atol', 1e-6)
    
    
    while t[-1] < tf:
        y_new = method(dt, f, t[-1], y[-1], rtol, atol, args)
        y.append(y_new)
        t.append(t[-1] + dt)
    return np.array(t), np.array(y)

# Define the derivative function for the helium atom.
def helium_derivs(t, y, *args):
    # State vector y = [x1, y1, vx1, vy1, x2, y2, vx2, vy2]
    x1, y1, vx1, vy1, x2, y2, vx2, vy2 = y
    r1 = np.array([x1, y1])
    r2 = np.array([x2, y2])
    # Avoid division by zero:
    norm_r1 = np.linalg.norm(r1) if np.linalg.norm(r1) > 1e-8 else 1e-8
    norm_r2 = np.linalg.norm(r2) if np.linalg.norm(r2) > 1e-8 else 1e-8
    r12 = r1 - r2
    norm_r12 = np.linalg.norm(r12) if np.linalg.norm(r12) > 1e-8 else 1e-8
    # Compute accelerations:
    a1 = -2 * r1 / (norm_r1**3) + r12 / (norm_r12**3)
    # Note: (r2 - r1)/|r12|^3 = - (r1 - r2)/|r12|^3
    a2 = -2 * r2 / (norm_r2**3) - r12 / (norm_r12**3)
    return np.array([vx1, vy1, a1[0], a1[1], vx2, vy2, a2[0], a2[1]])

Analysis

1. Trajectories Under Different Initial Conditions

Initial Conditions (Part 1)

• $r_1 = (2, 0)$, $r_2 = (-1, 0)$, $v_1 = (0, 0.95)$, $v_2 = (0, -1)$

Using solve_ode_RK45, we simulated up to $t = 100$ with an initial step $\Delta t = 0.001$. Below is the code:

# Problem from the book:

"""1. For simplicity, we restrict our atom to two dimensions. Simulate the classical helium atom. 
Choose units such that the electron mass is one and the other constants are absorbed into the unit of 
charge so that the force between two electrons is $$ \left|F\right| = \frac{1}{r^2}$$
Choose the initial value of the time step to be $\Delta t = 0.001$. For example, try these initial conditions:
$$r_1 = (2,0), r_2 = (−1,0),v_1 = (0,0.95),\text{ and }v_2 = (0,−1).$$"""

# Part 1: Simulation for initial conditions
# r1 = (2,0), r2 = (-1,0), v1 = (0,0.95), v2 = (0,-1)
y0 = np.array([2, 0, 0, 0.95, -1, 0, 0, -1])
tspan = [0, 100]  # simulate for 100 time units
dt = 0.001
rtol = 1e-12
atol = 1e-14

t, y = solve_ode_RK45(helium_derivs, tspan, y0, args=(), first_step=dt, rtol=rtol, atol=atol)

# Plot trajectories of electron 1 and electron 2.
plt.figure(figsize=(8,6))
plt.plot(y[:,0], y[:,1], label='Electron 1')
plt.plot(y[:,4], y[:,5], label='Electron 2')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Trajectories of Electrons (Initial Conditions, Part 1)')
plt.legend()
plt.grid(True)
plt.show()

<>:3: SyntaxWarning: invalid escape sequence '\l'
<>:3: SyntaxWarning: invalid escape sequence '\l'
/tmp/ipykernel_16693/2460938740.py:3: SyntaxWarning: invalid escape sequence '\l'
  """1. For simplicity, we restrict our atom to two dimensions. Simulate the classical helium atom.

From the plot:

• Electron 1 (teal) orbits at larger radii.
• Electron 2 (yellow) remains in a tighter orbit.

No immediate ionization occurs; the outer electron remains bound despite having a large orbit.

---

2. Braiding Orbits and Autoionization

Braiding Orbits (Part 2)

• $r_1 = (1.4, 0)$, $r_2 = (-1, 0)$, $v_1 = (0, 0.86)$, $v_2 = (0, -1)$

This set of initial conditions yields braiding orbits, in which the electrons appear to weave around each other without either escaping.

Autoionization with a Slight Perturbation

If we slightly modify $v_1$ from 0.86 to 0.87, one electron escapes after some time.

# Problem from the book:

""" 2. Most initial conditions result in unstable orbits in which one electron eventually leaves the atom (autoionization). 
The initial conditions $$r_1 = (1.4,0), r_2 = (−1,0),v_1 = (0,0.86),\text{ and }v_2 = (0,−1)$$ gives “braiding” orbits. 
Make small changes in this initial condition to observe autoionization."""

# Part 2: Investigate braiding orbits
# Base braiding initial conditions:
# r1 = (1.4, 0), r2 = (-1, 0), v1 = (0, 0.86), v2 = (0, -1)
y0_braiding = np.array([1.4, 0, 0, 0.86, -1, 0, 0, -1])
tspan = [0, 100] # simulate for 100 time units

t_b, y_b = solve_ode_RK45(helium_derivs, tspan, y0_braiding, args=(), first_step=dt, rtol=rtol, atol=atol)

# Now perturb the initial condition slightly (e.g. change v1_y from 0.86 to 0.87)
y0_perturbed = np.array([1.4, 0, 0, 0.87, -1, 0, 0, -1])
t_a, y_a = solve_ode_RK45(helium_derivs, tspan, y0_perturbed, args=(), first_step=dt, rtol=rtol, atol=atol)

# Create a single plot with two subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

# Plot for braiding orbits
ax1.plot(y_b[:,0], y_b[:,1], label='Electron 1 (braiding)')
ax1.plot(y_b[:,4], y_b[:,5], label='Electron 2 (braiding)')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_title('Braiding Orbits (Part 2)')
ax1.legend()
ax1.grid(True)

# Plot for perturbed orbits
ax2.plot(y_a[:,0], y_a[:,1], label='Electron 1 (perturbed)')
ax2.plot(y_a[:,4], y_a[:,5], label='Electron 2 (perturbed)')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('Autoionization with a Slight Perturbation (Part 2)')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.show()

The two images demonstrate the high sensitivity to initial conditions: small changes in velocity can lead to drastically different outcomes (bound orbits vs. one electron eventually escaping).

---

3. Angular Momentum Conservation

I tested two different initial conditions (IC1 and IC2) to investigate how well angular momentum is conserved. Below are sample plots of electron trajectories along with time evolution of total angular momentum:

# Problem from the book:

"""3. The classical helium atom is capable of very complex orbits (see Figure 5.8). 
Investigate the motion for the initial condition 
$$r_1 = (3,0), r_2 = (1,0),v_1 = (0,0.4),\text{ and }v_2 = (0,−1).$$ 
Does the motion conserve the total angular momentum? 
Also try $$r_1 = (2.5,0), r_2 = (1,0),v_1 = (0,0.4),\text{ and }v_2 = (0,−1).$$"""


# Helper function to compute total angular momentum (in 2D, L = x*v_y - y*v_x)
def total_angular_momentum(state):
    # state: [x1, y1, vx1, vy1, x2, y2, vx2, vy2]
    L1 = state[0]*state[3] - state[1]*state[2]
    L2 = state[4]*state[7] - state[5]*state[6]
    return L1 + L2

# Initial condition set 1: r1 = (3,0), r2 = (1,0), v1 = (0,0.4), v2 = (0,-1)
y0_1 = np.array([3, 0, 0, 0.4, 1, 0, 0, -1])
tspan = [0, 50]
t1, y1 = solve_ode_RK45(helium_derivs, tspan, y0_1, args=(), first_step=dt, rtol=rtol, atol=atol)
L1_values = [total_angular_momentum(state) for state in y1]

# Initial condition set 2: r1 = (2.5,0), r2 = (1,0), v1 = (0,0.4), v2 = (0,-1)
y0_2 = np.array([2.5, 0, 0, 0.4, 1, 0, 0, -1])
t2, y2 = solve_ode_RK45(helium_derivs, tspan, y0_2, args=(), first_step=dt, rtol=rtol, atol=atol)
L2_values = [total_angular_momentum(state) for state in y2]

# Create a single plot with two subplots
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
(ax3, ax4), (ax1, ax2) = axes

# Plot for initial condition set 1
ax1.plot(t1, L1_values, label='Total Angular Momentum (IC1)')
ax1.set_xlabel('Time')
ax1.set_ylabel('Angular Momentum')
ax1.set_title('Angular Momentum Conservation (IC1)')
ax1.legend()
ax1.grid(True)

# Plot for initial condition set 2
ax2.plot(t2, L2_values, label='Total Angular Momentum (IC2)')
ax2.set_xlabel('Time')
ax2.set_ylabel('Angular Momentum')
ax2.set_title('Angular Momentum Conservation (IC2)')
ax2.legend()
ax2.grid(True)

ax3.plot(y1[:,0], y1[:,1], label='Electron 1')
ax3.plot(y1[:,4], y1[:,5], label='Electron 2')
ax3.set_xlabel('x')
ax3.set_ylabel('y')
ax3.set_title('IC1: Trajectories of Electrons')
ax3.legend()
ax3.grid(True)

# Plot for perturbed orbits
ax4.plot(y2[:,0], y2[:,1], label='Electron 1')
ax4.plot(y2[:,4], y2[:,5], label='Electron 2')
ax4.set_xlabel('x')
ax4.set_ylabel('y')
ax4.set_title('IC2: Trajectories of Electrons')
ax4.legend()
ax4.grid(True)

plt.tight_layout()
plt.show()

Although the total angular momentum varies by a minuscule amount, it remains relatively close to a constant value over time. The variation is on the order of $10^{-13}$, while machine precision is about $10^{-15}$. Small deviations are due to the adaptive time-stepping and floating-point arithmetic. I would state with reasonable confidence that this momentum conservation is about as good as it can get. Thus, our numerical scheme is stable and reasonably accurate for the problem at hand.

---

4. Chaotic Ionization Times

Finally, we varied the initial $v_{1x}$ from 0.6 to 1.3 in increments of 0.02, measuring how long it takes before one electron's distance from the nucleus exceeds 6. The result is the following Ionization Time vs. $v_{1x}$ curve:

# Problem from the book:

"""4. Choose the initial condition $$r_1 = (2,0), r_2 = (−1,0),\text{ and }v_2 = (0,−1).$$ 
Then vary the initial value of $v_1$ from $(0.6,0)$ to $(1.3,0)$ in steps of $\Delta v = 0.02$. 
For each set of initial conditions, calculate the time it takes for autoionization. 
Assume that ionization occurs when either electron exceeds a distance of six from the nucleus. 
Run each simulation for a maximum time of 2000. Plot the ionization time versus $v_{1x}$. 
Repeat for a smaller interval of ∆v centered about one of the longer ionization times. 
These calculations require much computer resources. 
Do the two plots look similar? 
If so, such behavior is called “self-similar” and is characteristic of chaotic systems and the geometry of fractals."""

def simulate_until_ionization(y0, dt, rtol, atol, tf=2000):
    """
    Simulate the system using RK45 until either the simulation time reaches tf or 
    one electron's squared distance from the nucleus exceeds 36 (i.e. distance > 6).
    """
    t_current = 0.0
    state = y0.copy()
    threshold2 = 6.0**2  # 6^2 = 36
    while t_current < tf:
        state = RK45(dt, helium_derivs, t_current, state, rtol, atol, args=())
        t_current += dt
        if np.dot(state[0:2], state[0:2]) > threshold2 or np.dot(state[4:6], state[4:6]) > threshold2:
            return t_current
    return tf

<>:3: SyntaxWarning: invalid escape sequence '\D'
<>:3: SyntaxWarning: invalid escape sequence '\D'
/tmp/ipykernel_16693/4137021622.py:3: SyntaxWarning: invalid escape sequence '\D'
  """4. Choose the initial condition $$r_1 = (2,0), r_2 = (−1,0),\text{ and }v_2 = (0,−1).$$

# Sweep v1_x from 0.6 to 1.3 (step 0.02) using a list comprehension
v1x_values = np.arange(0.6, 1.3 + 0.02, 0.02)
ionization_times = [
    simulate_until_ionization(np.array([2, 0, v1x, 0, -1, 0, 0, -1]), dt, rtol, atol, tf=2000)
    for v1x in v1x_values
]

plt.figure(figsize=(8,6))
plt.plot(v1x_values, ionization_times, marker='o')
plt.xlabel('$v_{1x}$')
plt.ylabel('Ionization Time')
plt.title('Ionization Time vs. $v_{1x}$ (Full Range)')
plt.grid(True)
plt.show()

We observe:

• Rapid fluctuations in ionization time for slightly different initial velocities.
• For certain “windows” of $v_{1x}$, the electron remains bound for much longer.

To investigate self-similarity, we zoom in around one of these peaks:

# Zoom in around the region with longer ionization times.
max_index = np.argmax(ionization_times)
v1x_max = v1x_values[max_index]
# Choose a zoom window around v1x_max (e.g. ±0.05) with a finer step (0.002)
v1x_zoom = np.arange(v1x_max - 0.05, v1x_max + 0.05, 0.005)
ionization_times_zoom = [
    simulate_until_ionization(np.array([2, 0, v1x, 0, -1, 0, 0, -1]), dt, rtol, atol, tf=2000)
    for v1x in v1x_zoom
]

plt.figure(figsize=(8,6))
plt.plot(v1x_zoom, ionization_times_zoom, marker='o')
plt.xlabel('$v_{1x}$ (zoomed)')
plt.ylabel('Ionization Time')
plt.title('Ionization Time vs. $v_{1x}$ (Zoomed In)')
plt.grid(True)
plt.show()

The zoomed-in plot shows a pattern of spiky fluctuations resembling the behavior seen on the full range, indicative of fractal-like or self-similar structure. This is a hallmark of chaotic systems: small changes in initial conditions produce large changes in outcomes (in this case, the time before ionization).

---

Interpretation

1. Summary of Results

   • Braiding Orbits: Certain initial conditions lead to stable, complex paths where the electrons orbit each other and the nucleus without either electron escaping.
   • Autoionization: Many initial conditions eventually lead to one electron escaping; the system is inherently unstable in most cases.
   • Sensitivity to Initial Conditions: Tiny changes in velocity (e.g., from 0.86 to 0.87) can drastically change whether the electron escapes and when it does so, showcasing chaotic dynamics.
   • Angular Momentum: While not perfectly conserved (due to numerical approximation), it remains reasonably stable, indicating the solver is functioning as intended.

2. Insights

   • The 2D helium atom provides a rich and chaotic system with orbits that can be braided or can lead to ionization.
   • Long ionization times in certain intervals of initial velocity hint at quasi-stable orbits. These intervals often repeat in a self-similar pattern.

3. Conclusions and Future Directions

   • Chaotic Nature: The system exhibits hallmark features of chaos, such as sensitive dependence on initial conditions and fractal-like structure in the ionization times.
   • Computational Approaches: An adaptive time-step method (RK45) is crucial for capturing rapid changes when electrons approach the nucleus without using prohibitively small uniform time steps.
   • Extensions: A future project might include
     • investigating the 3D helium atom,
     • adding relativistic corrections,
     • or allowing the nucleus to recoil (instead of treating it as fixed).