Physical benchmarks: Doppler temperature#
In this notebook we benchmark atomsmltr simulations on simple physical examples. We focus on the effect of fluctuations in the temperature of a collection of atoms in a three-dimensionnal optical molasses, and check that we do reach the so-called Doppler temperature.
Formulaes and orders of magnitude#
We consider an atom of mass \(m\) illuminated by six pairs of counter-propagating laser beams aligned with the axes of a Cartesian frame. For the sake of simplicity, we will treat the laser beams as plane waves, and focus on the evolution of the atom’s speed. We will not consider any magnetic field, and model our atom as a two-level system. In this setting, there laser polarization has no effect.
The lasers have a wavenumber \(k\) and are detuned by \(\delta\) with respect the atom’s optical transition, which has a natural linewidth \(\Gamma\). We denote \(s_0 = I_0 / I_\mathrm{sat}\) the saturation parameter for one laser beam, where \(I_0\) is the intensity of the laser and \(I_\mathrm{sat}\) the saturation intensity associated to the atom’s transition.
For a red-detuned laser beam (\(\delta < 0\)), the atom will undergo Doppler cooling until it reaches a minimum temperature, known as the Doppler temperature \(T_D\), given by:
Which reaches a minimum of \(T = \hbar\Gamma / 2k_B\) for \(\delta = - \Gamma/2\).
# imports
from atomsmltr.atoms import Ytterbium, Strontium
import scipy.constants as csts
# display examples
for atom in [Ytterbium(), Strontium()]:
for tr in ["main", "intercombination"]:
transition = atom.trans[tr]
print(f">>> {atom.name} ({tr})")
print(f"+ transition Γ : 2π x {transition.Gamma / 2 / csts.pi * 1e-6:.3f} MHz")
print(f"+ transition λ : {transition.wavelength * 1e9:.2f} nm")
print(f"+ transition k : {transition.k * 1e-9:.3f} nm^-1")
print(f"+ Doppler limit : {transition.Doppler_temperature * 1e6:.2f} µK")
print("")
>>> Ytterbium (main)
+ transition Γ : 2π x 28.900 MHz
+ transition λ : 398.91 nm
+ transition k : 0.016 nm^-1
+ Doppler limit : 693.49 µK
>>> Ytterbium (intercombination)
+ transition Γ : 2π x 0.182 MHz
+ transition λ : 555.80 nm
+ transition k : 0.011 nm^-1
+ Doppler limit : 4.37 µK
>>> Strontium (main)
+ transition Γ : 2π x 32.000 MHz
+ transition λ : 460.86 nm
+ transition k : 0.014 nm^-1
+ Doppler limit : 767.88 µK
>>> Strontium (intercombination)
+ transition Γ : 2π x 0.007 MHz
+ transition λ : 689.45 nm
+ transition k : 0.009 nm^-1
+ Doppler limit : 0.18 µK
Test: laser cooling of ytterbium#
We will benchmark Doppler cooling with ytterbium. In this example, we use a set of six laser beams to cool an ytterbium atom. For the sake of simplicity, we will model the laser beams as plane waves with infinite spatial extent, so that we do not have to take care of polarization issues or spatial diffusion.
Setup configuration with atomsmltr#
Let’s define a configuration with a set of six laser beams.
# - IMPORTS
import numpy as np
import matplotlib.pyplot as plt
from scipy import constants as csts
from atomsmltr.simulation import Configuration
from atomsmltr.atoms.collection import Ytterbium
from atomsmltr.environment import PlaneWaveLaserBeam
# - GENERAL CONFIGURATION
# set laser saturation parameter
s0 = 0.05 # I/I_sat
main = Ytterbium().trans["main"] # 399nm transition of Yb
# create laser list
lasers = {}
for axis, direction in zip(["x", "y", "z"], [(1, 0, 0), (0, 1, 0), (0, 0, 1)]):
for head, mult in zip([">", "<"], [1, -1]):
dir = np.array(direction) * mult
tag = axis + head
laser = PlaneWaveLaserBeam(wavelength=main.wavelength)
laser.direction = dir
laser.tag = tag
laser.set_power_from_I(s0 * main.Isat)
lasers[tag] = laser
# generate config
config = Configuration(atom=Ytterbium())
config += [*lasers.values()]
Checkout Doppler temperature at optimal detuning#
We start with a detuning of \(\delta = -\Gamma/2\) to check optimal cooling.
# - PREPARE SIMULATION
# update config parameters
delta = -0.5 * main.Gamma
config.reset_atomlight_coupling()
for laser in config.list_lasers():
config.add_atomlight_coupling(laser=laser, transition="main", detuning=delta)
# initial conditions and time
u0 = np.zeros(shape=(2_000, 6)) # a set of 2e3 atoms at zero speed and position
t = np.linspace(0, 0.5e-3, 500) # integrate over 0.5ms
# - SIMULATE
from atomsmltr.simulation import RK4St
sim = RK4St(config)
res = sim.integrate(u0, t)
# - CHECK RESULT
# extract temperature from speed
v = res.y[:, 3:6, :] # we take only the speed vector
v_norm = np.linalg.norm(v, axis=1)
T = Ytterbium().mass * np.mean(v_norm**2, axis=0) / csts.k / 3
# compute Doppler temperature
T_Doppler = main.get_Doppler_temperature(delta)
# plot
plt.figure(figsize=(8, 5))
plt.plot(res.t * 1e3, T * 1e6, label="atom temperature")
plt.hlines(
T_Doppler * 1e6, 0, res.t.max() * 1e3, color="C1", label="Doppler temperature"
)
plt.xlim(0, res.t.max() * 1e3)
plt.grid()
plt.legend(loc="lower right")
plt.xlabel("time (ms)")
plt.ylabel("temperature (µK)")
plt.tight_layout()
plt.show()
Scanning the detuning#
We now scan the detuning to see the evolution of the atom’s temperature. Since the scattering rate is also impacted, we need to adjust the integration time accordingly.
# - PREPARE SIMULATION
delta_list = -np.linspace(0.2, 2, 10) * main.Gamma
v0 = np.array([0, 0, 0])
N_samples = 2_000
N_photons = 3_000 # number of scattered photons per beam
N_steps = 600
# - RUN
u0 = np.zeros(shape=(N_samples, 6))
T_list = []
for delta in delta_list:
# simulation settings
scatt_rate = main.get_scattering_rate(main.Isat * s0, 0, [1, 0, 0], delta)
Tmax = N_photons / scatt_rate
# print
print(f"+ delta = {delta / main.Gamma:.1f} Gamma | Tmax = {Tmax * 1e3:.2f} ms")
# update config
config.reset_atomlight_coupling()
for laser in config.list_lasers():
config.add_atomlight_coupling(laser=laser, transition="main", detuning=delta)
# simulate
t = np.linspace(0, Tmax, N_steps)
sim = RK4St(config)
res = sim.integrate(u0,t)
# compute temperature
v = res.y[:, 3:6, :] # we take only the speed vector
v_norm = np.linalg.norm(v, axis=1)
T = Ytterbium().mass * np.mean(v_norm**2, axis=0) / csts.k / 3
# store
T_list.append(T)
T_list = np.array(T_list)
+ delta = -0.2 Gamma | Tmax = 0.80 ms
+ delta = -0.4 Gamma | Tmax = 1.12 ms
+ delta = -0.6 Gamma | Tmax = 1.65 ms
+ delta = -0.8 Gamma | Tmax = 2.39 ms
+ delta = -1.0 Gamma | Tmax = 3.34 ms
+ delta = -1.2 Gamma | Tmax = 4.50 ms
+ delta = -1.4 Gamma | Tmax = 5.87 ms
+ delta = -1.6 Gamma | Tmax = 7.46 ms
+ delta = -1.8 Gamma | Tmax = 9.26 ms
+ delta = -2.0 Gamma | Tmax = 11.27 ms
# - PLOT
# theory settings
delta_plot = -np.linspace(0.1, 2.2, 1000) * main.Gamma
T_mean = np.mean(T_list[:, -150:], axis=1)
T_std = np.std(T_list[:, -150:], axis=1)
# figure
fig, ax = plt.subplots(1, 2, figsize=(9,4), tight_layout=True)
# plot temperature evolution
ax[0].plot(T_list.T * 1e3)
ax[0].set_xlabel("integration step")
ax[0].set_xlim(0, N_steps)
# plot final temperature
ax[1].plot(delta_plot / main.Gamma, main.get_Doppler_temperature(delta_plot) * 1e3)
ax[1].errorbar(delta_list / main.Gamma, T_mean * 1e3, yerr= T_std * 1e3, fmt="o")
ax[1].set_xlabel("detuning (units of $\Gamma$)")
# final setup
for cax in ax:
cax.set_ylabel("temperature (mK)")
cax.grid()
plt.show()
Bonus : test vectorization overhead#
Let’s check how adding more atoms impacts the simulation performances
import time
# - PREPARE SIMULATION
# update config parameters - back to optimal conditions
delta = -0.5 * main.Gamma
config.reset_atomlight_coupling()
for laser in config.list_lasers():
config.add_atomlight_coupling(laser=laser, transition="main", detuning=delta)
# initial conditions and time
t = np.linspace(0, 0.6e-3, 600)
sim = RK4St(config)
# scan atom number
N_list = np.logspace(1, 5, 9, dtype=int)
T_mean = []
T_std = []
duration = []
for N_samples in N_list:
print(f"+ N = {N_samples}")
# define sample
u0 = np.zeros(shape=(N_samples, 6))
# integrate
t0 = time.time()
res = sim.integrate(u0, t)
tf = time.time()
# get temperature
v = res.y[:, 3:6, :] # we take only the speed vector
v_norm = np.linalg.norm(v, axis=1)
T = Ytterbium().mass * np.mean(v_norm**2, axis=0) / csts.k / 3
# store
T_mean.append(np.mean(T[-100:]))
T_std.append(np.std(T[-100:]))
duration.append(tf - t0)
print(f" > done in {tf-t0:.1f} seconds")
T_mean = np.array(T_mean)
T_std = np.array(T_std)
duration = np.array(duration)
+ N = 10
> done in 2.2 seconds
+ N = 31
> done in 2.5 seconds
+ N = 100
> done in 2.9 seconds
+ N = 316
> done in 4.0 seconds
+ N = 1000
> done in 6.0 seconds
+ N = 3162
> done in 13.3 seconds
+ N = 10000
> done in 35.7 seconds
+ N = 31622
> done in 114.9 seconds
+ N = 100000
> done in 418.4 seconds
# figure
fig, ax = plt.subplots(1, 2, figsize=(9,4), tight_layout=True)
# time
ax[0].plot(N_list, duration, "o")
ax[0].set_ylabel("integration time (s)")
# result
ax[1].errorbar(N_list, T_mean * 1e3, yerr=T_std * 1e3, fmt="o")
ax[1].axhline(main.Doppler_temperature * 1e3, color="C1")
ax[1].set_ylabel("T (mK)")
for cax in ax:
cax.set_xscale("log")
cax.grid(which="both")
cax.set_xlabel("atom number")
plt.show()
