Laser Cooling

Author

Daniel Fischer

Introduction

In general, several mechanisms can exert forces on atoms exposed to laser fields. One example is the optical dipole force, a conservative force arising from the Autler–Townes shift (for details see the page on Rabi Oscillations) at positions with a non-vanishing laser-intensity gradient. Another mechanism is the photon scattering force, which originates from the impulse an atom receives during each photon-absorption and photon-emission event.

This chapter examines how the scattering force can be used experimentally to cool atomic gases to extremely low temperatures and to confine them within small spatial regions. A concise and very readable overview of these mechanisms can also be found in the article Laser Cooling and Trapping of Neutral Atoms by Harold J. Metcalf and Peter van der Straten @doi:10.1002/9783527600441.oe005, which closely parallels the material presented here.

After studying this chapter, you will be able to:

  • Explain the origin of the scattering force and how it depends on detuning and saturation.
  • Derive the velocity-dependent force leading to Doppler cooling.
  • Explain the concept of optical molasses and its relation to velocity damping.
  • Compute the Doppler cooling limit and understand its physical meaning.
  • Describe the principles of magneto-optical trapping and the role of the Zeeman effect.

1. Doppler Cooling

Laser cooling relies on the fact that light can exert forces on atoms. Before examining Doppler cooling in detail, we first analyze the scattering force, which arises when atoms repeatedly absorb and spontaneously emit photons. This force is central to all radiative cooling mechanisms and forms the foundation of optical molasses, magneto-optical traps, and most modern cold-atom techniques.

1.1 The Scattering Force

To understand the scattering force, let us recall that a force is defined as the rate of change of momentum, \[ \vec F = \frac{d\vec p}{dt}. \]

Photons carry momentum. A single photon of wave vector \(\vec{k}\) has momentum \[ \vec p_{\text{ph}} = \hbar\vec{k}, \] directed along the laser propagation direction. Therefore, each photon absorption event gives the atom a momentum kick of \(+\hbar\vec k\). A spontaneous emission event also imparts momentum of magnitude \(\hbar k\), but in a random direction, so its average contribution to the force is zero.

Stimulated emission does not contribute to a net force either, because the emitted photon is identical to the absorbed one and travels in the same direction; an absorption–stimulated-emission cycle therefore results in zero net momentum transfer.

Thus the net force comes entirely from directed absorption events, each of which transfers \(\hbar\vec k\), occurring at the spontaneous scattering rate.

Using the spontaneous scattering rate from Equation 13 derived from the Optical Bloch Equations, the resulting average force is

\[ \vec{F} = (\hbar\vec{k})\Gamma = \hbar\vec{k}\,\frac{\gamma}{2}\, \frac{S_0}{1 + S_0 + \frac{4\delta^2}{\gamma^2}}. \tag{1}\]

This force points along the laser beam direction and is largest when the detuning \(\delta\) is zero.

1.2 Doppler Cooling

The expression in Equation 1 describes atoms at rest. For a moving atom, however, the laser frequency is Doppler-shifted in the atom’s rest frame. The first-order (non-relativistic) Doppler shift is

\[ \Delta\omega_D = -\vec{k}\cdot\vec{v}, \]

valid when \(v \ll c\). Substituting this shifted frequency into the scattering-force expression gives

\[ \vec{F} = \hbar\vec{k}\frac{\gamma}{2} \frac{S_0}{1 + S_0 + \frac{4(\delta - \vec{k}\cdot\vec{v})^2}{\gamma^2}}. \tag{2}\]

This expression is fully three-dimensional.
To make further progress—and because most laser-cooling setups are arranged along fixed beam axes—we now restrict the motion to one spatial dimension, say the \(x\)-axis. In this 1D geometry:

  • the wave vector becomes \(\vec{k} = k\,\hat{x}\),
  • the velocity becomes \(\vec{v} = v\,\hat{x}\),
  • and therefore \(\vec{k}\cdot\vec{v} = kv\).

This simplifies the Doppler-shifted scattering force and allows us to study how it depends on the sign and magnitude of the atomic velocity.

Two Counterpropagating Beams

We now consider the standard cooling configuration of two counterpropagating laser beams, one with wave vector \(+k\) and one with \(-k\). Neglecting Autler–Townes shifts and interference effects, the total force along the beam axis is the difference of the forces from the two beams:

\[ \begin{aligned} F &= \hbar k\frac{\gamma}{2} \left[ \frac{S_0}{1 + S_0 + \tfrac{4(\delta - kv)^2}{\gamma^2}} - \frac{S_0}{1 + S_0 + \tfrac{4(\delta + kv)^2}{\gamma^2}} \right] \\[6pt] &\approx \hbar k S_0 \frac{kv}{2} \frac{16(\delta/\gamma)} { 1 + \frac{8}{\gamma^2}(\delta^2 + k^2 v^2) + \frac{16}{\gamma^4}(\delta^2 + k^2 v^2)^2 } , \end{aligned} \tag{3}\]

where the second line assumes the low-saturation limit \(S_0 \ll 1\).

Figure 1 shows the resulting force as a function of velocity for red detuning (\(\delta<0\)). Near \(v=0\), the force always opposes the motion:

  • for \(v>0\), the force is negative;
  • for \(v<0\), the force is positive.

Thus the atom is slowed—this is Doppler cooling. Using three orthogonal beam pairs generalizes this to three-dimensional optical molasses.

Code 
import numpy as np
import matplotlib.pyplot as plt

# --- Set up plot styling ---
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.size'] = 14


def scattering_force_single_beam(v, k, gamma, delta, S0, hbar=1):
    """
    Scattering force from a single laser beam including Doppler shift.
    
    Parameters:
    -----------
    v : array-like
        Atomic velocity (in units where recoil velocity = 1)
    k : float
        Wave vector magnitude
    gamma : float
        Natural linewidth
    delta : float
        Laser detuning (negative for red detuning)
    S0 : float
        On-resonance saturation parameter
    hbar : float
        Reduced Planck constant (default=1 for natural units)
    
    Returns:
    --------
    F : array-like
        Scattering force
    """
    detuning_eff = delta - k * v
    F = hbar * k * (gamma/2) * S0 / (1 + S0 + 4*detuning_eff**2/gamma**2)
    return F


def doppler_cooling_force(v, k, gamma, delta, S0, hbar=1):
    """
    Total force from two counterpropagating beams (Doppler cooling).
    
    Parameters:
    -----------
    v : array-like
        Atomic velocity
    k : float
        Wave vector magnitude
    gamma : float
        Natural linewidth
    delta : float
        Laser detuning (negative for red detuning)
    S0 : float
        On-resonance saturation parameter
    hbar : float
        Reduced Planck constant (default=1 for natural units)
    
    Returns:
    --------
    F : array-like
        Total cooling force
    """
    # Force from beam propagating in +x direction (opposes +v)
    F_plus = scattering_force_single_beam(v, k, gamma, delta, S0, hbar)
    
    # Force from beam propagating in -x direction (opposes -v)
    F_minus = scattering_force_single_beam(-v, k, gamma, delta, S0, hbar)
    
    # Total force (note the sign: +k beam gives +F, -k beam gives -F)
    F_total = F_plus - F_minus
    
    return F_total


def doppler_cooling_force_low_sat(v, k, gamma, delta, S0, hbar=1):
    """
    Doppler cooling force in low saturation approximation (S0 << 1).
    This is the analytical approximation from Eq. (eq-coolingforce).
    """
    numerator = hbar * k * S0 * (k * v / 2) * 16 * (delta / gamma)
    denominator = (1 + 8/gamma**2 * (delta**2 + k**2 * v**2) + 
                   16/gamma**4 * (delta**2 + k**2 * v**2)**2)
    return numerator / denominator

def linear_approx(v, k, gamma, delta, S0, hbar=1):
    F = 8 * v * hbar/gamma * (k**2) * S0 * delta / (1 + (2*delta/gamma)**2)**2
    return F


# Set parameters (in natural units where gamma = 1)
gamma = 1.0
k = 1.0  # Wave vector
delta_values = [-0.5]  # Red detuning (negative values)
S0 = 0.1  # Low saturation

# Velocity range (in units of recoil velocity)
v = np.linspace(-10, 10, 500)

# Create figure with subplots
fig, ax = plt.subplots(figsize=(7.5, 5))

# Plot 1: Force vs velocity for different detunings
for delta in delta_values:
    F = doppler_cooling_force(v, k, gamma, delta, S0)
    Fl = scattering_force_single_beam(v, k, gamma, delta, S0)
    Fr = scattering_force_single_beam(v, -k, gamma, delta, S0)
    Flin = linear_approx(v, k, gamma, delta, S0)
    ax.plot(v, F, linewidth=2, label=r'\noindent Total scattering force \\$\delta = -\frac{1}{2}\gamma$')
    ax.plot(v, Fl, 'r:', linewidth=1.5, label=r'Force due to an individual beam')
    ax.plot(v, Fr, 'r:', linewidth=1.5)
    ax.plot(v, Flin, 'g--', linewidth=1, label=r'Low-velocity limit')

ax.axhline(y=0, color='k', linestyle='--', linewidth=0.8, alpha=0.5)
ax.axvline(x=0, color='k', linestyle='--', linewidth=0.8, alpha=0.5)
ax.set_xlabel('Velocity $v$ [in units of $\gamma/k$]', fontsize=13)
ax.set_ylabel('Force $F$ [in units of $\hbar k \gamma$]', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(-5, 5)
ax.set_ylim(-.05, .05)


plt.tight_layout()
plt.show()
Scattering force for two counterpropagating, red-detuned laser beams. The dashed curves show the forces from each beam; the solid curve is the sum.
Figure 1: Scattering force for two counterpropagating, red-detuned laser beams. The dashed curves show the forces from each beam; the solid curve is the sum.

1.3 Low-velocity limit

For \(kv \ll \gamma\) and \(kv \ll |\delta|\), Equation 3 reduces to a linear viscous force,

\[ F = -\alpha v, \qquad \alpha = -\frac{8\,\hbar k^2 S_0(\delta/\gamma)}{\bigl[1+(2\delta/\gamma)^2\bigr]^2}, \]

which shows that, at sufficiently small velocity, the Doppler cooling force behaves like friction.

We start from the low-saturation approximate expression from Equation 3 (one-dimensional motion along the beam axis):

\[ F \approx \hbar k S_0 \frac{k v}{2}\; \frac{16(\delta/\gamma)} { 1 + \frac{8}{\gamma^2}(\delta^2 + k^2 v^2) + \frac{16}{\gamma^4}(\delta^2 + k^2 v^2)^2 } . \]

Now take the low-velocity limit \(k v \ll \gamma\) and \(k v \ll |\delta|\). In this limit we can set \(k^2 v^2 \to 0\) in the denominator, so the denominator reduces to

\[ 1 + \frac{8\delta^2}{\gamma^2} + \frac{16\delta^4}{\gamma^4} = \bigl[\,1 + (2\delta/\gamma)^2 \,\bigr]^2. \]

The numerator simplifies as

\[ \hbar k S_0\cdot\frac{k v}{2}\cdot 16\frac{\delta}{\gamma} = 8\,\hbar k^2 S_0\frac{\delta}{\gamma}\, v . \]

Combining numerator and denominator gives

\[ F \approx \frac{8\,\hbar k^2 S_0(\delta/\gamma)}{\bigl[1+(2\delta/\gamma)^2\bigr]^2}\; v . \]

Writing this in the standard damping form \(F=-\alpha v\) gives

\[ \alpha = -\frac{8\,\hbar k^2 S_0(\delta/\gamma)}{\bigl[1+(2\delta/\gamma)^2\bigr]^2}. \]

Since the factors \(\hbar\), \(k^2\), \(S_0\), and \(\gamma\) are positive, the sign of \(\alpha\) depends entirely on the detuning \(\delta\):

  • For red detuning \(\delta<0\) we have \(\alpha>0\), so \(F=-\alpha v\) is damping (opposes the velocity) and yields cooling.
  • For blue detuning \(\delta>0\) we get \(\alpha<0\), so \(F\) amplifies the motion (heating).

Thus Doppler cooling only occurs for red detuning.

Cooling (energy) rate and time constant

With kinetic energy \(E=\tfrac12 m v^2\) we have

\[ \left(\frac{dE}{dt}\right)_{\text{cool}} = F v = -\alpha v^2 = -\frac{2\alpha}{m}\,E. \tag{4}\]

Hence the energy decays exponentially,

\[ E(t) = E(0)\,e^{- (2\alpha/m)\,t}, \]

with an energy damping rate \(2\alpha/m\) and a characteristic cooling time

\[ \tau_{\rm cool} = \frac{m}{2\alpha}. \]

For red detuning (\(\alpha>0\)), \(\tau_{\rm cool}\) is positive and sets how quickly the atomic energy decreases.

Optimal detuning for maximum damping

To maximize the damping strength \(|\alpha|\), set \(x = \delta/\gamma\) and maximize

\[ |\alpha| \propto \frac{|x|}{\bigl[1+(2x)^2\bigr]^2}. \]

The maximum occurs at

\[ |x| = \frac{1}{2\sqrt{3}} \approx 0.289, \qquad \delta \approx -\frac{\gamma}{2\sqrt{3}}. \]

This detuning gives the fastest damping, but not the lowest achievable steady-state temperature. The minimal Doppler temperature,

\[ T_{\min}=\frac{\hbar \gamma}{2k_B}, \]

is reached at \(\delta = -\gamma/2\). Both optimal detunings are of order \(\gamma\), but correspond to different goals (fast cooling vs. lowest temperature).

1.4 The Doppler Cooling Limit

Each spontaneous emission event imparts a random recoil momentum of magnitude \(\hbar k\) to the atom.
These kicks act like a momentum-space diffusion process. Averaging over many such events gives a heating rate.

To estimate it, note first that the mean kinetic energy of a particle in one dimension is related to temperature via the Maxwell–Boltzmann distribution, \[ E = \frac{1}{2} m\langle v^2\rangle = \frac{1}{2} k_B T . \]

The heating rate is the mean recoil energy per scattering event multiplied by the scattering rate: \[ \left(\frac{dE}{dt}\right)_{\text{heat}} \approx \frac{\hbar^2 k^2}{m}\,\Gamma. \]

The kinetic energy of an atom with momentum \(\vec p\) is \[ E = \frac{p^2}{2m}. \] During a scattering event the atom receives a momentum kick \(\Delta \vec p\) of magnitude \(\hbar \vec k\). Its new momentum is \[ \vec p'=\vec p + \Delta \vec p \] so the new kinetic energy is \[ E = \frac{(\vec p+\Delta \vec p)^2}{2m}. \] The change in kinetic energy is \[ \Delta E = E' - E = \frac{(\vec p+\Delta \vec p)^2-p^2}{2m} = \frac{2\vec p\Delta \vec p + (\Delta p)^2}{2m} \]

This expression has two contributions:

  1. The cross term: \[ \frac{\vec p\Delta \vec p}{m}. \] Because spontaneous emission is isotropic, the direction of the recoil is random, so \[ \langle \Delta \vec p\rangle =0 , \quad \Rightarrow \quad \langle \vec p \Delta \vec p\rangle = 0. \] This term averages to zero.
  2. Diffusion term: \[ \frac{(\Delta p)^2}{2m}. \] This is always positive and therefore always increases the kinetic energy.

Thus the mean energy increase per recoil event is \[ \langle \Delta E\rangle = \frac{\langle (\Delta p)^2 \rangle }{2m}. \] For a photon of wave number \(k\) the recoil momentum has magnitude \(|\Delta p|=\hbar k\), and emission directions are random. Therefore, \[ \langle \Delta E\rangle = \frac{(\hbar k)^2}{2m} \] is the mean recoil energy per spontaneous emission event.

Using the low-velocity form of the scattering rate, \[ \Gamma = \frac{\gamma}{2}\, \frac{S_0}{1 + (2\delta/\gamma)^2}, \] we obtain the approximate heating rate \[ \left(\frac{dE}{dt}\right)_{\text{heat}} \approx \frac{\hbar^2 k^2}{m}\, \frac{\gamma S_0}{1 + (2\delta/\gamma)^2}. \tag{5}\]

This expression implicitly assumes that the atomic velocity distribution is broad enough that recoil events effectively randomize momentum, but not so large that Doppler shifts modify \(\Gamma\) appreciably.

Steady-state temperature

Cooling (from Equation 4) and heating (from Equation 5) balance in steady state:

\[ \left(\frac{dE}{dt}\right)_{\text{cool}} + \left(\frac{dE}{dt}\right)_{\text{heat}} = 0. \]

Solving for the mean-square velocity gives

\[ \langle v^2\rangle = \frac{\hbar k}{4m}\, \frac{1 + (2\delta/\gamma)^2}{2|\delta|/\gamma}. \]

This is minimized for the optimal detuning \(\delta = -\gamma/2\), yielding the Doppler cooling limit: \[ T_{\min} = \frac{m}{k_B}\langle v^2\rangle_{\min} = \frac{\hbar\gamma}{2k_B}. \]

This is the fundamental temperature limit for Doppler cooling based solely on absorption and spontaneous emission.

NoteBeyond the Doppler limit

Colder temperatures can be achieved through mechanisms requiring multilevel structure and polarization gradients, such as Sisyphus or polarization-gradient cooling. These are known as sub-Doppler processes and can reach temperatures near the recoil limit:

\[ T_{\text{recoil}} = \frac{\hbar k^2}{2 m k_B}. \]

2. Magneto-Optical Traps (MOTs)

Laser cooling, as discussed in the previous section, can reduce the velocity spread of an atomic gas. However, optical molasses alone does not confine atoms spatially: the ensemble slowly drifts away. This limitation is a consequence of the Optical Earnshaw Theorem, which states that no purely optical scattering-force configuration can create a stable 3D trap for structureless atoms.

2.1 Optical Earnshaw Theorem

The central assumption is that atoms are point-like dipoles with internal structure independent of position. Under this assumption, the scattering force from a laser field is proportional to the Poynting vector \[ \vec S = \frac{1}{\mu_0}\,\vec E \times \vec B. \]

From Poynting’s theorem, the divergence of \(\vec S\) is
\[ \nabla\!\cdot\vec S = -\, \vec J\cdot\vec E - \frac{1}{2}\frac{d}{dt}\!\left(\varepsilon_0 E^2 + \frac{1}{\mu_0}B^2\right). \]

In the region where atoms are located:

  • there are no free charges or currents, so \(\vec J = 0\),
  • the laser field is (after averaging) time-independent.

Thus, \[ \nabla\!\cdot\vec S = 0. \]

If the force is proportional to \(\vec S\), then \[ \nabla\!\cdot\vec F = 0. \]

A divergence-free vector field cannot point inward on the entire boundary of any finite volume. Therefore:

No configuration of scattering forces alone can create a stable trapping potential.
Any such trap must “leak”—atoms eventually escape.

This motivates adding a mechanism that breaks the assumptions of the Earnshaw theorem. In a magneto-optical trap, this loophole is provided by making the atomic energy levels position dependent using a magnetic field.

2.2 Position-dependent structure: Zeeman effect

A MOT introduces a weak, spatially varying magnetic field produced by anti-Helmholtz coils, creating a near-linear field gradient, \[ B_z(z) \approx B'_z\,z. \]

This magnetic field shifts the energies of magnetic sublevels via the Zeeman effect. For a simple model, consider:

  • a ground state with orbital angular momentum \(\ell=0\) (an s-state),
  • an excited state with \(\ell=1\) (a p-state) with magnetic sublevels \(m_\ell=-1,0,+1\).

Ignoring spin for now, the magnetic sublevels of the excited state shift by \[ \Delta\omega_{\rm Zeeman} = \frac{\mu_B}{\hbar}\,g\,m_\ell\,B_z, \] where \(\mu_B\) is the Bohr magneton and \(g\) is the Landé g-factor.

Thus, at positive \(z\) the \(m_\ell=+1\) state shifts up in energy, while the \(m_\ell=-1\) state shifts down; the opposite happens at negative \(z\). The key consequence: resonance frequencies now depend on position.

This breaks the assumptions of the optical Earnshaw theorem and allows a restoring force.

2.3 Selection rules and polarization

A second ingredient is the ability to drive transitions selectively into specific Zeeman sublevels using circular polarization:

  • A beam propagating in the \(+z\) direction with \(\sigma^+\) polarization drives
    \[ \Delta m_\ell = +1 \quad \Rightarrow\quad m_\ell = +1 \text{ excited state}. \]

  • A counterpropagating beam from \(−z\) with \(\sigma^-\) polarization drives
    \[ \Delta m_\ell = -1 \quad \Rightarrow\quad m_\ell = -1 \text{ excited state}. \]

Because the Zeeman shifts depend on position:

  • At positive \(z\), the \(m_\ell=-1\) state is closer to resonance.
    → The \(\sigma^-\) beam (from the −z direction) scatters more strongly.
    → The resulting force pushes atoms back toward \(z=0\).

  • At negative \(z\), the \(m_\ell=+1\) state is closer to resonance.
    → The \(\sigma^+\) beam (from +z) scatters more strongly.
    → Again the force pushes atoms toward \(z=0\).

Including a slight red detuning makes the Doppler cooling mechanism work simultaneously: atoms feel both a restoring force and a velocity-damping force.

The combination of:

  1. Doppler cooling
  2. Zeeman-shift–induced spatial dependence
  3. Polarization-selective transitions

creates a three-dimensional viscous and restoring force field: the magneto-optical trap.

To visualize how the MOT combines a spatially varying magnetic field with polarization-dependent optical forces, Figure 2 shows (left) the magnetic quadrupole field generated by anti-Helmholtz coils and (right) the resulting Zeeman-shifted excited-state manifold together with the \(\sigma^\pm\) laser couplings.

The linear field gradient, evident in the left panel, provides the position-dependent Zeeman shift that breaks the assumptions of the optical Earnshaw theorem. The right panel illustrates how each circularly polarized beam selectively addresses one of the shifted \(m_\ell=\pm1\) states, thereby producing a restoring force toward the trap center.

Code 
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipk, ellipe


# --- Set up plot styling ---
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.size'] = 14

# Constants
mu0 = 4 * np.pi * 1e-7

# --- FIELD CALCULATION FUNCTIONS (re-used for context) ---
def B_field_single_loop(x_point, z_point, R, current, z_offset):
    # Vectorized Biot-Savart function
    rho = np.abs(x_point)
    z_prime = z_point - z_offset
    
    k_squared_numerator = 4 * R * rho
    k_squared_denominator = (R + rho)**2 + z_prime**2
    
    k_squared = np.divide(k_squared_numerator, k_squared_denominator, 
                          out=np.zeros_like(x_point), 
                          where=k_squared_denominator != 0)
    
    k_squared = np.where(k_squared < 0, 0, k_squared)
    k_squared = np.where(k_squared > 1, 1.0 - 1e-9, k_squared)
    
    K_val = ellipk(k_squared)
    E_val = ellipe(k_squared)
    
    term1 = (mu0 * current) / (2 * np.pi)
    denom_sqrt_K = np.sqrt(k_squared_denominator)
    
    # Bx Component
    m1 = (rho != 0) 
    Bx_numerator = z_prime * (E_val * ( (R**2 + rho**2 + z_prime**2) / ((R - rho)**2 + z_prime**2) ) - K_val)
    Bx_factor = term1 / denom_sqrt_K
    
    Bx = np.zeros_like(x_point)
    Bx[m1] = Bx_factor[m1] * (Bx_numerator[m1] / rho[m1])
    Bx *= np.sign(x_point) 

    # Bz Component
    Bz_factor = term1 / denom_sqrt_K
    Bz_term = K_val + E_val * ( (R**2 - rho**2 - z_prime**2) / ((R - rho)**2 + z_prime**2) )
    Bz = Bz_factor * Bz_term
    
    return Bx, Bz

def B_anti_helmholtz_2D(x, z, R_coil, d_coil, current):
    Bx1, Bz1 = B_field_single_loop(x, z, R_coil, current, -d_coil/2)
    Bx2, Bz2 = B_field_single_loop(x, z, R_coil, -current, d_coil/2)
    return Bx1 + Bx2, Bz1 + Bz2


# --- Plotting Setup ---
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3), width_ratios=[2,3])

# Coil parameters
R_coil = 1.0
d_coil = R_coil
current = 1.0

# Plotting area
xlim = (-2.5 * R_coil, 2.5 * R_coil)
zlim = (-2.5 * R_coil, 2.5 * R_coil)
ax1.set_xlim(xlim)
ax1.set_ylim(zlim)
ax1.set_aspect('equal')

# Create a grid for streamplot
num_points = 100
x_grid = np.linspace(xlim[0], xlim[1], num_points)
z_grid = np.linspace(zlim[0], zlim[1], num_points)
X, Z = np.meshgrid(x_grid, z_grid)

# Calculate the magnetic field at each grid point
Bx_field, Bz_field = B_anti_helmholtz_2D(X, Z, R_coil, d_coil, current)

# --- Define Explicit Symmetric Seed Points (Avoiding Axes) ---

# 1. Define seeds only for the first quadrant (x > 0, z > 0)
num_radial_points = 5 
num_angular_points = 6 

# Radial distances: Start well away from the origin (0.4)
radial_distances = np.linspace(0.4, 2.0, num_radial_points) 

# Angular range: from 15 degrees (pi/12) to 75 degrees (5*pi/12), avoiding the axes (0 and 90)
angles = np.linspace(np.pi / 12, 5 * np.pi / 12, num_angular_points)

# Generate seed points in Quadrant 1
seed_x_q1 = []
seed_z_q1 = []

for r in radial_distances:
    for angle in angles:
        seed_x_q1.append(r * np.cos(angle))
        seed_z_q1.append(r * np.sin(angle))

seed_x_q1 = np.array(seed_x_q1)
seed_z_q1 = np.array(seed_z_q1)

# 2. Reflect Q1 seeds to the other three quadrants to ensure symmetry

# Quadrant 2 (x < 0, z > 0) -> Reflection across z-axis (x -> -x)
seed_x_q2 = -seed_x_q1
seed_z_q2 = seed_z_q1

# Quadrant 3 (x < 0, z < 0) -> Reflection across origin (x -> -x, z -> -z)
seed_x_q3 = -seed_x_q1
seed_z_q3 = -seed_z_q1

# Quadrant 4 (x > 0, z < 0) -> Reflection across x-axis (z -> -z)
seed_x_q4 = seed_x_q1
seed_z_q4 = -seed_z_q1

# 3. Combine all seeds
seed_x_all = np.concatenate([seed_x_q1, seed_x_q2, seed_x_q3, seed_x_q4])
seed_z_all = np.concatenate([seed_z_q1, seed_z_q2, seed_z_q3, seed_z_q4])

seed_points = np.vstack([seed_x_all, seed_z_all])

# --- Streamline Visualization with Fixed Symmetric Seed Points ---
ax1.streamplot(X, Z, Bx_field, Bz_field, 
              color='black', 
              linewidth=1, 
              start_points=seed_points.T, 
              arrowstyle='->', 
              arrowsize=1.5)

# --- Plot Coil Cross-sections ---
circle_radius = 0.1 

# Coil 1 (z = -d_coil/2)
ax1.add_patch(plt.Circle((R_coil, -d_coil/2), circle_radius,color='lightgray', zorder=5)) 
ax1.text(R_coil+0.025, -d_coil/2+0.021, r'$\times$', color='black', fontsize=8, ha='center', va='center', weight='bold', zorder=6)
ax1.add_patch(plt.Circle((-R_coil, -d_coil/2), circle_radius, color='lightgray', zorder=5)) 
ax1.text(-R_coil+0.025, -d_coil/2+0.045, r'$\cdot$', color='black', fontsize=12, ha='center', va='center', weight='bold', zorder=6)

# Coil 2 (z = d_coil/2)
ax1.add_patch(plt.Circle((R_coil, d_coil/2), circle_radius, color='lightgray', zorder=5)) 
ax1.text(R_coil+0.025, d_coil/2+0.045, r'$\cdot$', color='black', fontsize=12, ha='center', va='center', weight='bold', zorder=6)
ax1.add_patch(plt.Circle((-R_coil, d_coil/2), circle_radius, color='lightgray', zorder=5)) 
ax1.text(-R_coil+0.025, d_coil/2+0.021, r'$\times$', color='black', fontsize=8, ha='center', va='center', weight='bold', zorder=6)

ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax1.spines['bottom'].set_visible(False)
ax1.spines['left'].set_visible(False)
ax1.set_xticks([])
ax1.set_yticks([])
ax1.grid(False)


# --- Parameters for the Atomic System and MOT ---
# Unperturbed transition frequency (arbitrary units for plotting, but represents E/hbar)
omega0_level = 1.0 # Energy of the unperturbed excited state (l=1, ml=0)

# Zeeman splitting coefficient (related to Bohr magneton and g-factor)
# This determines the slope of the Zeeman-split levels.
mu_B_prime_g = 0.2 # (mu_B * g_factor * B_gradient_strength) in E/hbar units

# Laser detuning
delta0 = 0.2 # Amount of red-detuning (omega0 - omega)
laser_omega_level = omega0_level - delta0

# --- Plotting Range ---
z_min, z_max = -2.0, 2.0
z_axis = np.linspace(z_min, z_max, 200)

# --- Zeeman Level Calculations ---
# l=0 ground state (ml=0)
E_ground = 0.0 * z_axis

# l=1 excited states
# In an anti-Helmholtz field, B(z) is roughly B'z. So Zeeman splitting is proportional to z.
# Assuming B(z) = B' * z (linear field with gradient B' at origin)
# Energy split: E_ml = E_0 + ml * (mu_B * g) * B(z)
# For simplicity, let mu_B_prime_g * z represent the ml=1 splitting from ml=0 level
# The sign depends on the convention and actual field direction.
# For the plot, ml=+1 increases with z, ml=-1 decreases with z.
E_l1_ml_1 = omega0_level + mu_B_prime_g * z_axis  # m_l = +1
E_l1_ml_0 = omega0_level * np.ones_like(z_axis)    # m_l = 0
E_l1_ml_neg1 = omega0_level - mu_B_prime_g * z_axis # m_l = -1

# --- Plot Energy Levels ---
ax2.plot(z_axis, E_ground, 'k-', linewidth=1.5, label=r'$l=0, m_l=0$')
ax2.plot(z_axis, E_l1_ml_0, 'k-', linewidth=1.5, label=r'$l=1, m_l=0$')
ax2.plot(z_axis, E_l1_ml_1, 'k-', linewidth=1.5, label=r'$l=1, m_l=+1$')
ax2.plot(z_axis, E_l1_ml_neg1, 'k-', linewidth=1.5, label=r'$l=1, m_l=-1$')

ax2.quiver(z_min, -0.2, 4.5, 0, scale=1, scale_units='xy', angles='xy', width=0.003, headwidth=10, headlength=10, headaxislength=9)

# --- Add the fake break here ---
# Coordinates for the rectangle to hide the middle section
break_y_min = 0.5
break_y_max = 0.55
ax2.quiver(0, break_y_max, 0, 2-break_y_max, scale=1, scale_units='xy', angles='xy', width=0.003, headwidth=10, headlength=10, headaxislength=9)
ax2.quiver(0, -0.25, 0, break_y_min+0.22, scale=1, scale_units='xy', angles='xy', width=0.003, headwidth=0, headlength=0, headaxislength=0)
# Add the break symbol "//"
ax2.text(0, (break_y_min + break_y_max) / 2-0.02, r'$//$',
        ha='center', va='center', fontsize=10, color='black', rotation=90)


# --- Plot Laser Frequency (Red-Detuned) ---
ax2.plot([-2,2],[laser_omega_level,laser_omega_level], color='red', linestyle='--', linewidth=1)

# --- Add Transition Arrows and Labels ---
# Arrow for omega0 (unperturbed transition at z=0)
ax2.annotate("", xy=(0.6, omega0_level), xytext=(0.6, E_ground[0]),
            arrowprops=dict(arrowstyle="<->", color="black", lw=1))
ax2.text(0.7, (omega0_level + E_ground[0])/2, r'$\hbar\omega_0$', va='center', ha='left', fontsize=12)
# Arrow for omegaL
ax2.annotate("", xy=(z_min+1, laser_omega_level), xytext=(z_min+1, E_ground[0]),
            arrowprops=dict(arrowstyle="<->", color="black", lw=1))
ax2.text(z_min+1.1, (laser_omega_level + E_ground[0])/2, r'$\hbar\omega_L$', va='center', ha='left', fontsize=12)


# --- Add Labels for Energy Levels and m_l states ---
# Ground state label
ax2.text(z_max + 0.6, -0.25, r'$z, B_z$', va='center', ha='left', fontsize=12)
ax2.text(0.1, 2.05, r'$E$', va='center', ha='left', fontsize=12)


ax2.text(z_max + 0.1, E_l1_ml_1[-1]+0.3, r'$m_\ell$', va='center', ha='left', fontsize=12)
ax2.text(z_max + 0.1, E_ground[0], r'$0$', va='center', ha='left', fontsize=12)
ax2.text(z_max + 0.1, E_l1_ml_0[-1], r'$0$', va='center', ha='left', fontsize=12)
ax2.text(z_max + 0.1, E_l1_ml_1[-1], r'$1$', va='center', ha='left', fontsize=12)
ax2.text(z_max + 0.1, E_l1_ml_neg1[-1], r'$-1$', va='center', ha='left', fontsize=12)
ax2.text(z_min - 0.3, E_ground[0]+0.125, r'$|g\rangle$', va='center', ha='center', fontsize=12)
ax2.text(z_min - 0.3, E_l1_ml_0[-1]+0.125, r'$|e\rangle$', va='center', ha='center', fontsize=12)
ax2.text(z_min - 0.3, E_ground[0]-0.125, r'$\ell=0$', va='center', ha='center', fontsize=12)
ax2.text(z_min - 0.3, E_l1_ml_0[-1]-0.125, r'$\ell=1$', va='center', ha='center', fontsize=12)

# Photon arrows
# --- Parameters for the Sine Wave Arrow ---
start_x = z_min
end_x = z_min+0.6
mid_y = 0.35  # Central y-position for the wave
amplitude = 0.05 # Height of the wave
frequency = 10   # Number of oscillations along the arrow's length
num_points = 100 # Number of points to define the wave (smoother is better)

# 1. Generate the sine wave path
x_path = np.linspace(start_x, end_x, num_points)
y_path = mid_y + amplitude * np.sin(frequency * np.pi * (x_path - start_x) / (end_x - start_x))
y_path[-20:] = y_path[-21]


# 2. Plot the sine wave line
ax2.plot(x_path, y_path, color='red', linewidth=1)
# 3. Add an arrowhead at the end of the sine wave
# We need the last point and a point just before it to calculate the tangent for the arrow direction.
arrow_end_x = x_path[-1]
arrow_end_y = y_path[-1]

# Get a point slightly before the end for tangent calculation
# Using the second to last point
arrow_base_x = x_path[-2]
arrow_base_y = y_path[-2]

# Calculate the direction vector for the arrowhead
dx = arrow_end_x - arrow_base_x
dy = arrow_end_y - arrow_base_y

# Use ax.annotate for the arrowhead, with xytext very close to xy
# The arrow will point from xytext to xy
ax2.annotate("", 
            xy=(arrow_end_x, arrow_end_y), 
            xytext=(arrow_end_x - dx*0.1, arrow_end_y - dy*0.1), # A short segment just before the end
            arrowprops=dict(arrowstyle="simple, head_length=0.2, head_width=0.15",  # <-- INCREASE THESE VALUES
                            color="red", 
                            lw=2, 
                            shrinkA=0, 
                            shrinkB=0))

ax2.text((start_x+end_x)/2+0.3, mid_y+0.05, r'$\sigma^+$', color='red', va='bottom', ha='center', fontsize=12)


# Photon arrows2
# --- Parameters for the Sine Wave Arrow ---
start2_x = z_max
end2_x = z_max-0.6

# 1. Generate the sine wave path
x2_path = np.linspace(start2_x, end2_x, num_points)
y2_path = mid_y + amplitude * np.sin(frequency * np.pi * (x_path - start_x) / (end_x - start_x))
y2_path[-20:] = y_path[-21]


# 2. Plot the sine wave line
ax2.plot(x2_path, y2_path, color='red', linewidth=1)
# 3. Add an arrowhead at the end of the sine wave
# We need the last point and a point just before it to calculate the tangent for the arrow direction.
arrow2_end_x = x2_path[-1]
arrow2_end_y = y2_path[-1]

# Get a point slightly before the end for tangent calculation
# Using the second to last point
arrow2_base_x = x2_path[-2]
arrow2_base_y = y2_path[-2]

# Calculate the direction vector for the arrowhead
dx2 = arrow2_end_x - arrow2_base_x
dy2 = arrow2_end_y - arrow2_base_y

# Use ax.annotate for the arrowhead, with xytext very close to xy
# The arrow will point from xytext to xy
ax2.annotate("", 
            xy=(arrow2_end_x, arrow2_end_y), 
            xytext=(arrow2_end_x - dx2*0.1, arrow2_end_y - dy2*0.1), # A short segment just before the end
            arrowprops=dict(arrowstyle="simple, head_length=0.2, head_width=0.15",  # <-- INCREASE THESE VALUES
                            color="red", 
                            lw=2, 
                            shrinkA=0, 
                            shrinkB=0))
ax2.text((start2_x+end2_x)/2-0.2, mid_y+0.05, r'$\sigma^-$', color='red', va='bottom', ha='center', fontsize=12)


# General settings
ax2.set_xlim([-2.5,2.5])
ax2.set_ylim([-0.5,omega0_level+1])
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)
ax2.spines['bottom'].set_visible(False)
ax2.spines['left'].set_visible(False)
ax2.set_xticks([])
ax2.set_yticks([])


plt.show()
Two-panel figure: left panel shows quadrupole magnetic field lines from anti-Helmholtz coils; right panel shows position-dependent Zeeman-split levels and sigma-plus/minus transitions used for MOT trapping.
Figure 2: Magnetic field of anti-Helmholtz coils (left) and schematic working principle of the MOT (right). The field gradient produces a position-dependent Zeeman shift, while the \(\sigma^\pm\) beams selectively drive the shifted magnetic sublevels to yield a restoring force.

Key Takeaways

  • Optical molasses provides efficient velocity damping through Doppler cooling, but cannot confine atoms spatially. The fundamental limitation arises because, for structureless atoms in static optical fields, the scattering force field is divergence-free, prohibiting a true optical trap.

  • This constraint is formalized in the Optical Earnshaw Theorem, which follows from Poynting’s theorem: with no free charges or currents, the divergence of the Poynting vector vanishes, and thus no purely optical, inward-pointing force field can enclose a stable trapping volume.

  • Magneto-optical traps (MOTs) circumvent Earnshaw’s theorem by giving the atom’s internal structure a position dependence—specifically through a spatially varying Zeeman shift in an inhomogeneous magnetic field.

  • The anti-Helmholtz coil configuration generates a magnetic quadrupole field with a linear gradient. This gradient produces a position-dependent resonance condition, which modifies the scattering force asymmetrically for atoms displaced from the trap center.

  • Circular polarization is essential: a \(\sigma^+\) beam selectively drives transitions to the \(m_\ell = +1\) Zeeman sublevel, while a counterpropagating \(\sigma^-\) beam drives \(m_\ell = -1\). Combined with a red detuning, this ensures that atoms experience a restoring force toward the field zero, yielding both cooling and trapping.

  • The MOT therefore combines:

    1. Doppler cooling (velocity damping)
    2. Zeeman-shift–assisted spatial confinement (restoring force)
      into one of the most robust and widely used methods for preparing cold atomic samples.