Bragg beam splitter

Physics

The Bragg beam splitter splits the wave function into two momentum states (but same internal state). Here, the main physics is sketched. The interested reader is referred to the textbook by Grynberg, Aspect and Fabre.

The formalism describing the Bragg atom-light interaction is based on a semi-classical model. The atoms are described quantum mechanically via a wave function, while the light is described classically (its high intensity ensures a substantial photon presence at all times). The Hamiltonian of the system is:

\[\hat H = \frac{\hat p^2}{2m} - \hat{\vec D} \vec E_L (\vec r, t)\]

where \(\hat p\) is the momentum operator, \(\hat{\vec D}\) is the dipole operator and \(\vec E_L\) is the electric field.

The electric field can be described by two counterpropagating laser beams with frequencies \(\omega_L+\omega_r\) (drives the transition \(| g,0 \rangle \rightarrow | e, \hbar k_L \rangle\)) and \(\omega_r\) (drives the transition \(| e, \hbar k_L \rangle \rightarrow | g, 2\hbar k_L\rangle\)). Note that \(\frac{(2\hbar k_L)^2}{2m} = \hbar \omega_r\). It should be noted that both lasers are detuned by \(\Delta\), such that the transition \(|g,0 \rangle \rightarrow | e, \hbar k_L \rangle\) is unlikely to happen without the stimulated emission directly after it.

Here, a one dimensional wave function \(\Psi (x)\) is considered. Adiabatic elimination of the excited state leads to

\[\hat H = -\frac{\hbar^2}{2m}\nabla^2 + 2 \hbar \Omega \cos ^2 \left( k_L x - \frac{\omega_r}{2} t \right)\]

where \(\Omega\) is the effective Rabi frequency. \(\Omega\) is determined by the laser properties and has typically a Gaussian temporal profile to ensure good velocity selectivity. If the atoms are freely falling, an acceleration term \(\frac{1}{2}a_\text{laser}t^2\) is added to the laser phase to ensure that the laser beams stay resonant to the falling atoms. Also common is an additional constant phase shift \(\Phi_0\).

After the atom-light interaction, the atom is left in a superposition of states \(|g,0\rangle\) and \(|g,2\hbar k_L\rangle\), and typically higher orders like \(|g,-2\hbar k_L\rangle\) and \(|g,4\hbar k_L\rangle\) [Siemß 2020]. Idealized, this sequence applies a momentum transfer of \(2\hbar k_L\) to the atom with a \(50\%\) chance.

Implementation

This example illustrates the implementation of propagating a wave function through a Bragg beam splitter. We propagate a wave function (described by an fa.Array) under the influence of the time-dependent Bragg beam splitter potential.

from typing import Any
import numpy as np
import math
from bokeh.io import output_notebook

import fftarray as fa
from helpers import plt_array, plt_array_values_space_time

# setup bokeh in jupyter notebook
output_notebook(hide_banner=True)

Initialize global constants

from scipy.constants import hbar, pi, atomic_mass
# angular frequency used to initialize the ground state (of quantum harmonic
# oscillator)
omega_x = 2*pi*10 # Hz
# laser pulse parameters
# Rabi frequency. This specific value was found as a binary search to
# optimize a 50/50 split of the two momentum classes for this specific beam
# splitter duration and pulse form.
rabi_frequency = 25144.285917282104 # Hz
phi_0 = 0. # phase jump
bragg_acc = 0. # bragg acceleration
sigma_bs = 25e-6 # temporal pulse width (s)

# Rb87 properties
m = 86.909 * atomic_mass # mass in kg
lambda_L = 780.032 * 1e-9 # wavelength in m
k_L = 2 * pi / lambda_L # wave number in 1/m
hbark: float = hbar * k_L # momentum transfer
hbarkv: float = hbark/m # recoil velocity
w_r = 2 * hbarkv * k_L # recoil frequency in rad/s

Define the time grid

In order to simulate the time evolution of the wave function, we need to discretize time into intervals of length \(\Delta t\) during which the potential is assumed constant. Then, for each time step, we can apply the time evolution operator as \(|\Psi(t+\Delta t)\rangle = e^{-i \hat H \Delta t /\hbar}|\Psi(t)\rangle\) using the split step method.

Here, we choose to sample the Gaussian temporal profile of the Bragg beam splitter potential with Gaussian width \(\sigma_\mathrm{bs}\) with a step size of \(\Delta t = 1\) µs for a duration of \(4\sigma_\mathrm{bs}\). Additionally, we let the wave function freely propagate for \(25\) ms with \(\Delta t = 50\) µs after applying the beam potential to illustrate the separation of both momentum states in position space.

 laser start      intensity peak        laser end            simulation end
------|-----------------|-------------------|----------------------|-----> time
      |------------ 4 sigma_bs -------------|-- free propagation --|

# define how many sigmas of the gauss to sample in each direction before
# reaching zero intensity:
sampling_range_mult = 4. # * sigma_bs
dt_bs = 1e-6 # time step size
# total number of pulse grid steps = gauss_width * scaling_factor / step_size
steps_limit = int(round(sigma_bs * sampling_range_mult / dt_bs))
t_offset = steps_limit*dt_bs
nt_bs = 2*steps_limit # number of time steps for beam splitter
dt_free = 5e-5 # defines time step size for free propagation
nt_free = 50 # number of time steps for free propagation
# time lists
t_list_bs: Any = np.arange(nt_bs)*dt_bs
t_list_free = t_list_bs[-1]+np.arange(1,nt_free+1)*dt_free
t_list = np.concatenate((t_list_bs, t_list_free))

Initialize the wave function

We initialize the wave function as the groundstate of a quantum harmonic oscillator with frequency \(\omega_\mathrm{QHO} = 2\pi\times 10\) Hz.

# coordinate grid
x_dim: fa.Dimension = fa.dim_from_constraints(
    name = "x",
    pos_min = -50e-6,
    pos_max = 50e-6,
    freq_middle = 0.,
    freq_extent = 32*k_L,
    loose_params = ["freq_extent"],
)

# initialize an Array as harmonic oscillator groundstate
x_coords: fa.Array = fa.coords_from_dim(x_dim, "pos")
psi_init: fa.Array = (m * omega_x / (pi*hbar))**(1./4.) * fa.exp(-(m * omega_x * (x_coords**2.)/(2.*hbar)))
psi_init = psi_init / fa.integrate(fa.abs(psi_init)**2)
plt_array(psi_init)

Define the potential

Now, we implement the external potential

\[V = 2 \hbar \Omega(t) \cos ^2 \left( k_L x - \frac{\omega_r}{2} t \right)\]

where the time dependent Rabi frequency is defined by

\[\Omega(t) = \Omega_0 \exp(-t^2/(2\sigma_\mathrm{bs}^2)).\]

from scipy.constants import hbar

def V(ramp: float, t: float) -> fa.Array:
    """Bragg pulse potential.

    Parameters
    ----------
    ramp : float
        The pulse ramp (scaling the rabi frequency).
    t : float
        The global lab time.

    Returns
    -------
    fa.Array
        The potential at time t.
    """
    return (rabi_frequency * ramp * 2. * hbar * fa.cos(
        k_L * (x_coords - 0.5 * bragg_acc * t**2)
        - 0.5 * w_r * t
        + phi_0/2.
    )**2)

Now, we plot the potential in position space at peak intensity (t=t_offset).

from bokeh.plotting import figure, show
# plot the energy trend during the imaginary time evolution
plt = figure(
    width=700, height=400, min_border=50,
    title="The potential versus the wave function's position grid",
    x_axis_label="x [m]",
    y_axis_label="Potential/hbar [Hz]"
)
x_grid = x_dim.values("pos", xp=np)
potential_values = V(ramp=1, t=t_offset).values("pos", xp=np)/hbar
plot = plt.line(x_grid, potential_values)
show(plt)

Propagate the wave function

Now, we implement the iterative procedure of evolving the wave function in time according to the previously defined potential and the subsequent free propagation.

def gauss(t: float, sigma: float):
    """Helper function for the temporal Bragg beam profile."""
    return math.exp(-0.5 * (t / sigma)**2)

# compute the boundary value of the gauss to ensure that the potential is zero
# when we start and stop applying of the potential
gauss_offset = gauss(t = -t_offset, sigma = sigma_bs)

# wavenumber operator squared (used for kinetic operator)
k_sq = ((2*pi*fa.coords_from_arr(psi_init, x_dim.name, "freq"))**2.).into_dtype("complex")

# store probability densities at every time step
psi_data_bs = np.empty((2, len(t_list_bs), x_dim.n))
psi_data_free = np.empty((2, len(t_list_free), x_dim.n))

psi = psi_init

# ---------------------------- Bragg beam splitter --------------------------- #
# propagate the wave function using the Bragg beam potential for times t_list_bs
# (assuming t_list_bs is evenly spaced with dt_bs)
for i, t in enumerate(t_list_bs):
    # evaluate the potential at time t
    ramp = gauss(t=t-t_offset, sigma=sigma_bs) - gauss_offset
    potential = V(ramp, float(t)).into_dtype("complex")

    # apply half kinetic propagator
    psi = psi.into_space("freq") * fa.exp((-1.j * 0.5*dt_bs * hbar / (2*m)) * k_sq)

    # apply potential propagator
    psi = psi.into_space("pos") * fa.exp((-1.j / hbar * dt_bs) * potential)

    # apply half kinetic propagator
    psi = psi.into_space("freq") * fa.exp((-1.j * 0.5*dt_bs * hbar / (2*m)) * k_sq)

    # save probability density in pos and freq space
    psi_data_bs[0,i] = fa.abs(psi.into_space("pos")).values("pos", xp=np)**2
    psi_data_bs[1,i] = fa.abs(psi.into_space("freq")).values("freq", xp=np)**2

# ----------------------------- Free propagation ----------------------------- #
# further freely propagate the wave function for times t_list_free (assuming
# t_list_free is evenly spaced with dt_free)
for i, _t in enumerate(t_list_free):
    # free propagation step
    psi = psi.into_space("freq") * fa.exp((-1.j * dt_free * hbar / (2*m)) * k_sq)

    psi_data_free[0,i] = fa.abs(psi.into_space("pos")).values("pos", xp=np)**2
    psi_data_free[1,i] = fa.abs(psi.into_space("freq")).values("freq", xp=np)**2

psi_final_free = psi

Plot the final wave function

The final wave function is in a superposition of momentum states |0hbark> and |2hbark>.

# plot the final wave function
plt_array(psi_final_free)

Plot the complete time evolution

# every reduce_pts_fac point is shown in time and pos/freq
reduce_pts_fac = 2 # 1 for higher resolution
if nt_bs % nt_free > 0:
    print("Warning: Time grid not equidistant in contour plot.")
psi_init_abs_pos = fa.abs(psi_init.into_space("pos")).values("pos", xp=np)**2
psi_init_abs_freq = fa.abs(psi_init.into_space("freq")).values("freq", xp=np)**2
# list of all position space values
plt_data_psi_final_pos = np.concatenate((
    # initial value
    [psi_init_abs_pos[::reduce_pts_fac]],
    # beam splitter values
    psi_data_bs[0,::reduce_pts_fac*nt_bs//nt_free,::reduce_pts_fac],
    # free propagation values
    psi_data_free[0,::reduce_pts_fac,::reduce_pts_fac]
))
# list of all momentum space values
plt_data_psi_final_freq = np.concatenate((
    # initial value
    [psi_init_abs_freq[::reduce_pts_fac]],
    # beam splitter values
    psi_data_bs[1][::reduce_pts_fac*nt_bs//nt_free,::reduce_pts_fac],
    # free propagation values
    psi_data_free[1,::reduce_pts_fac,::reduce_pts_fac]
))
plt_data_t_list = np.concatenate((
    [0],
    t_list_bs[::reduce_pts_fac] + dt_bs,
    t_list_free[::reduce_pts_fac] + dt_bs
))
plt_array_values_space_time(
    pos_values = np.real(plt_data_psi_final_pos),
    freq_values = np.real(plt_data_psi_final_freq),
    pos_grid = x_dim.values("pos", xp=np),
    freq_grid = x_dim.values("freq", xp=np)*lambda_L,
    time = plt_data_t_list,
    freq_unit = "hbar k_L",
    pos_range = (-1e-5, 4e-5),
    freq_range = (-0.5, 2.5),
)