MOT#
Simulating a MOT with atomsmltr
# - Packages
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib widget < disabled when building
# - atomSmltr
from atomsmltr.atoms import Ytterbium
from atomsmltr.environment.lasers import GaussianLaserBeam
from atomsmltr.environment.fields import MagneticQuadrupoleZ, MagneticOffset
from atomsmltr.environment.zones import Limits, Box
from atomsmltr.simulation import Configuration, ScipyIVP_3D
from atomsmltr.environment.lasers.polarization import CircularLeft, CircularRight
global config#
# - Init global config
config = Configuration(atom=Ytterbium())
main = config.atom.trans["main"]
# - setup magnetic field
B_grad_G_per_cm = 30
B_grad = B_grad_G_per_cm * 1e-4 / 1e-2
mag_quad = MagneticQuadrupoleZ(origin=(0, 0, 0), slope=B_grad, tag="MOT field")
# - setup lasers
# cf. config from Letellier et al. 2023
l399 = {
"wavelength": 399e-9,
"waist": 22e-3,
"power": 100e-3 / 6,
"waist_position": (0, 0, 0),
}
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 = GaussianLaserBeam(**l399)
laser.direction = dir
laser.tag = tag
lasers[tag] = laser
# check that s ~ 0.037
s = lasers["x>"].get_value((0, 0, 0)) / main.Isat
print(f"I = {s:.3f} I_sat")
I = 0.037 I_sat
def plot_force_3D(
sim,
ax1="x",
ax2="vx",
ax1_span=(-0.1, 0.1),
ax2_span=(-10, 10),
ax1_N=100,
ax2_N=101,
):
# - setup grid
# init
axes = {"x": 0, "y": 1, "z": 2, "vx": 3, "vy": 4, "vz": 5}
axi = {v:k for k, v in axes.items()}
lims = {ax: {"min": 0, "max": 0} for ax in axes}
N = {ax: 1 for ax in axes}
# select axis to plot
lims[ax1]["min"] = ax1_span[0]
lims[ax1]["max"] = ax1_span[1]
N[ax1] = ax1_N
lims[ax2]["min"] = ax2_span[0]
lims[ax2]["max"] = ax2_span[1]
N[ax2] = ax2_N
# generate grid
arrays = {ax: np.linspace(lims[ax]["min"], lims[ax]["max"], N[ax]) for ax in axes}
arrays_ord = [arrays[axi[i]] for i in range(6)]
grid = np.meshgrid(*arrays_ord, indexing="ij")
X1 = np.squeeze(grid[axes[ax1]])
X2 = np.squeeze(grid[axes[ax2]])
grid = np.squeeze(np.array([*grid]))
# - compute force
force = sim.get_force(grid.T)
FX, FY, FZ = force.T
# - prepare scales
vfx = np.max(np.abs(FX))
vfy = np.max(np.abs(FY))
vfz = np.max(np.abs(FZ))
# - plot
fig, axes = plt.subplots(1, 3, figsize=(8, 3), tight_layout=True)
axes[0].pcolormesh(X1, X2, FX, cmap="bwr", vmin=-vfx, vmax=vfx)
axes[0].set_title("Fx")
axes[1].pcolormesh(X1, X2, FY, cmap="bwr", vmin=-vfy, vmax=vfy)
axes[1].set_title("Fy")
axes[2].pcolormesh(X1, X2, FZ, cmap="bwr", vmin=-vfz, vmax=vfz)
axes[2].set_title("Fz")
for ax in axes:
ax.set_xlabel(ax1)
ax.set_ylabel(ax2)
ax.grid()
plt.show()
1D test#
laser polarization illustration#
We set the polarization of the laser beams and check their projections on the local magnetic field. Here, we always define the quantization axis with respect to the local magnetic field direction.
Here, the magnetic field gradient is positive along x.
Hence, the laser propagating in the +x direction, labeled “x>” should be “sigma-” for x<0 and “sigma+” for x>0.
# -- Plot illustration
# set laser polarization
lasers["x>"].polarization = CircularRight()
lasers["x<"].polarization = CircularRight()
# compute mag field
x = np.linspace(-0.1, 0.1, 1000)
r = np.array([[x, 0, 0] for x in x])
B = mag_quad.get_value(r)
Bx, By, Bz = B.T
# compute laser projection on quant axis
polar_up = lasers["x>"].get_polarization_quant(B)
polar_up_pi, polar_up_sp, polar_up_sm = polar_up.T
polar_dw = lasers["x<"].get_polarization_quant(B)
polar_dw_pi, polar_dw_sp, polar_dw_sm = polar_dw.T
# plot
fig, axes = plt.subplots(1, 3, figsize=(10,4), tight_layout=True)
axes[0].plot(x * 1e2, Bx * 1e4, label="Bx")
axes[0].plot(x * 1e2, By * 1e4, label="By")
axes[0].plot(x * 1e2, By * 1e4, label="Bz")
axes[0].set_xlabel("x (cm)")
axes[0].set_ylabel("B (G)")
axes[0].legend()
axes[0].grid()
axes[1].plot(x, polar_up_sm, color="C0", label="x>, sigma-")
axes[1].plot(x, polar_up_sp, color="C0", dashes=[2,2], label="x>, sigma+")
axes[1].set_xlabel("x (cm)")
axes[1].set_ylabel("polar proj")
axes[1].legend()
axes[1].grid()
axes[2].plot(x, polar_dw_sm, color="C1", label="x<, sigma-")
axes[2].plot(x, polar_dw_sp, color="C1", dashes=[2,2], label="x<, sigma+")
axes[2].set_xlabel("x (cm)")
axes[2].set_ylabel("polar proj")
axes[2].legend()
axes[2].grid()
plt.show()
good configuration#
# -- Testing with a simulation
# init config
config1D = Configuration(atom=Ytterbium())
# set limit
xlim = Limits(-0.15, 0.15, axis=0, target="position", action="stop", tag="xlim")
# add objects
config1D += xlim, lasers["x<"], lasers["x>"], mag_quad
# setup atomlight
for laser in config1D.list_lasers():
config1D.add_atomlight_coupling(laser=laser, transition="main", detuning=-2*main.Gamma)
# init sim
sim = ScipyIVP_3D(config1D)
# -- Run simulation
# settings
t = np.linspace(0, 0.1, 1000)
v_list = np.linspace(0.5, 30, 20)
u0_list = [(-0.1, 0, 0, vx, 0, 0) for vx in v_list]
# run
sim.u0_list = u0_list
coll = sim.run(t, npools=5, verbose=True)
Show code cell output
100%|██████████| 20/20 [00:00<00:00, 38.72it/s]
# -- compute force
grid = np.mgrid[-0.1:0.1:100j, 0:0:1j, 0:0:1j, -10:30:101j, 0:0:1j, 0:0:1j]
grid = np.squeeze(grid)
pos = grid.T
X, _, _, VX, _ ,_ = grid
force = sim.get_force(pos)
FX, FY, FZ = force.T
# -- plot
plt.figure()
plt.pcolormesh(X * 100, VX, FX, cmap="bwr")
for res in coll:
x = res.y[0]
vx = res.y[3]
if x[-1] > 0.05:
color="C3"
else:
color="C0"
plt.plot(x * 100, vx, color=color, linewidth=0.8)
plt.xlim(-10, 10)
plt.ylim(-10, 30)
plt.xlabel("position (cm)")
plt.ylabel("speed (m/s)")
plt.grid()
plt.show()
bad polarization#
# Invert laser polarization, for fun
lasers["x>"].polarization = CircularLeft()
lasers["x<"].polarization = CircularLeft()
config1D.update_objects([lasers["x>"], lasers["x<"]])
# run
sim.u0_list = u0_list
coll = sim.run(t, npools=3, verbose=True)
Show code cell output
100%|██████████| 20/20 [00:00<00:00, 75.09it/s]
# -- compute force
grid = np.mgrid[-0.1:0.1:100j, 0:0:1j, 0:0:1j, -10:30:101j, 0:0:1j, 0:0:1j]
grid = np.squeeze(grid)
pos = grid.T
X, _, _, VX, _ ,_ = grid
force = sim.get_force(pos)
FX, FY, FZ = force.T
plt.figure()
for res in coll:
x = res.y[0]
vx = res.y[3]
if x[-1] > 0.05:
color="C3"
else:
color="C0"
plt.plot(x * 100, vx, color=color, linewidth=0.8)
plt.pcolormesh(X * 100, VX, FX, cmap="bwr")
plt.xlim(-10, 10)
plt.ylim(-10, 30)
plt.xlabel("position (cm)")
plt.ylabel("speed (m/s)")
plt.grid()
plt.show()
simulation - add offset#
# Invert laser polarization, for fun
lasers["x>"].polarization = CircularRight()
lasers["x<"].polarization = CircularRight()
config1D.update_objects([lasers["x>"], lasers["x<"]])
mag_offset = MagneticOffset((30e-4,0,0))
config1D.rm_all_magnetic_fields()
config1D += mag_offset, mag_quad
# run
sim.u0_list = u0_list
coll = sim.run(t, npools=3, verbose=True)
Show code cell output
0%| | 0/20 [00:00<?, ?it/s]
100%|██████████| 20/20 [00:00<00:00, 25.87it/s]
# -- compute force
grid = np.mgrid[-0.1:0.1:100j, 0:0:1j, 0:0:1j, -10:30:101j, 0:0:1j, 0:0:1j]
grid = np.squeeze(grid)
pos = grid.T
X, _, _, VX, _ ,_ = grid
force = sim.get_force(pos)
FX, FY, FZ = force.T
plt.figure()
plt.pcolormesh(X * 100, VX, FX, cmap="bwr")
for res in coll:
x = res.y[0]
vx = res.y[3]
if x[-1] > 0.05:
color="C3"
else:
color="C0"
plt.plot(x * 100, vx, color=color, linewidth=0.8)
plt.plot(-1, 0, "ok") # new center should be at -1cm, since gradient is 30 G/cm
plt.xlim(-10, 10)
plt.ylim(-10, 30)
plt.xlabel("position (cm)")
plt.ylabel("speed (m/s)")
plt.grid()
plt.show()
3D test#
configuration#
# -- Config
# set laser polarization
# NB : the strong axis of the quadrupole is along z
lasers["x>"].polarization = CircularRight()
lasers["x<"].polarization = CircularRight()
lasers["y>"].polarization = CircularRight()
lasers["y<"].polarization = CircularRight()
lasers["z>"].polarization = CircularLeft()
lasers["z<"].polarization = CircularLeft()
# limits
bsize = 0.15
limits = [i * bsize for j in range(3) for i in (-1,1)]
spatial_limits = Box(*limits, target="position", action="stop", tag="limits")
# config
config3D = Configuration(atom=Ytterbium())
config3D += mag_quad
config3D += [*lasers.values()]
# atom-light
main = config3D.atom.trans["main"]
for laser in config3D.list_lasers():
config3D.add_atomlight_coupling(laser=laser, transition="main", detuning=-2*main.Gamma)
config3D.print_atomlight_info()
────────────────────────
| Atom-light couplings |
────────────────────────
. transition > 'main' :
├── laser 'x>' : detuning=-3.63e+08 (-2.00Γ)
├── laser 'x<' : detuning=-3.63e+08 (-2.00Γ)
├── laser 'y>' : detuning=-3.63e+08 (-2.00Γ)
├── laser 'y<' : detuning=-3.63e+08 (-2.00Γ)
├── laser 'z>' : detuning=-3.63e+08 (-2.00Γ)
└── laser 'z<' : detuning=-3.63e+08 (-2.00Γ)
. transition > 'intercombination' :
└── empty
# init sim
sim = ScipyIVP_3D(config3D)
Illustration : plot forces#
# -- compute force in position x speed space
slim = (-0.05, 0.05)
vlim = (-10, 10)
plot_force_3D(sim, "x", "vx", slim, vlim)
plot_force_3D(sim, "y", "vy", slim, vlim)
plot_force_3D(sim, "z", "vz", slim, vlim)
# -- compute force in position x position
slim = (-0.05, 0.05)
vlim = (-10, 10)
plot_force_3D(sim, "x", "y", slim, slim)
plot_force_3D(sim, "x", "z", slim, slim)
plot_force_3D(sim, "y", "z", slim, slim)
run simulations#
# -- Run simulation
# settings
t = np.linspace(0, 0.1, 1000)
v_list = np.linspace(0.5, 30, 20)
u0_list = [(-0.1, 0, 0, vx, 0, 0) for vx in v_list]
# run
sim.u0_list = u0_list
coll = sim.run(t, npools=5, verbose=True)
Show code cell output
100%|██████████| 20/20 [00:01<00:00, 14.49it/s]
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
for res in coll:
x = res.y[0]
y = res.y[1]
vx = res.y[3]
if x[-1] > 0.05:
color="C3"
else:
color="C0"
axes[0].plot(x * 100, vx, color=color, linewidth=0.8)
axes[1].plot(x * 100, y* 100, color= color, linewidth=0.8)
axes[0].plot(0, 0, "ok")
axes[0].set_xlim(-10, 10)
axes[0].set_ylim(-10, 30)
axes[0].set_xlabel("x (cm)")
axes[0].set_ylabel("vx (m/s)")
axes[0].grid()
axes[1].plot(0, 0, "ok")
axes[1].set_xlim(-10, 10)
axes[1].set_ylim(-10, 10)
axes[1].set_xlabel("x (cm)")
axes[1].set_ylabel("y (cm)")
axes[1].grid()
plt.show()
# -- Run simulation
# settings
t = np.linspace(0, 0.1, 1000)
v = 5
theta_list = np.linspace(0, np.pi/4, 20)
u0_list = [(-0.1, 0, 0, v * np.cos(t), v * np.sin(t), 0) for t in theta_list]
# run
sim.u0_list = u0_list
coll = sim.run(t, npools=5, verbose=True)
Show code cell output
100%|██████████| 20/20 [00:00<00:00, 27.41it/s]
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
for res in coll:
x = res.y[0]
y = res.y[1]
vx = res.y[3]
if x[-1] > 0.05:
color="C3"
else:
color="C0"
axes[0].plot(x * 100, vx, color=color, linewidth=0.8)
axes[1].plot(x * 100, y* 100, color= color, linewidth=0.8)
axes[0].plot(0, 0, "ok")
axes[0].set_xlim(-6, 20)
axes[0].set_ylim(-10, 30)
axes[0].set_xlabel("x (cm)")
axes[0].set_ylabel("vx (m/s)")
axes[0].grid()
axes[1].plot(0, 0, "ok")
axes[1].set_xlim(-10, 10)
axes[1].set_ylim(-1, 5)
axes[1].set_xlabel("x (cm)")
axes[1].set_ylabel("y (cm)")
axes[1].grid()
plt.show()
plot in 3D#
# -- Run simulation
# settings
t = np.linspace(0, 0.05, 1000)
vlist = [0, 1, 5]
u0_list = [(-0.01, -0.01, -0.01, v, 0, 0) for v in vlist]
# run
sim.u0_list = u0_list
coll = sim.run(t, npools=5, verbose=True)
Show code cell output
100%|██████████| 3/3 [00:00<00:00, 8.14it/s]
ax = plt.figure().add_subplot(projection='3d')
for res in coll:
x = res.y[0]
y = res.y[1]
z = res.y[2]
ax.plot(x * 100, y * 100, z * 100)
ax.plot(0, 0, 0, "ok")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
plt.show()
animate#
import matplotlib.animation as animation
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
for res in coll:
x = res.y[0]
y = res.y[1]
z = res.y[2]
ax.plot(x * 100, y * 100, z * 100, color="k", alpha=0.1)
ax.plot(0, 0, 0, "ok")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
lines = []
lines = [ax.plot([], [], [])[0] for _ in coll]
def update_lines(num, lines, coll):
for l, res in zip(lines, coll):
x = res.y[0][:num]
y = res.y[1][:num]
z = res.y[2][:num]
l.set_data_3d(x * 100, y * 100, z*100)
return lines
num_steps = len(res.t)
anim = animation.FuncAnimation(
fig, update_lines, num_steps, fargs=(lines, coll), interval=10, repeat=True,
)
# writer = animation.PillowWriter(fps=30,
# metadata=dict(artist='Me'),
# bitrate=1800)
# anim.save('MOT.gif', writer=writer)
plt.show()
