Documentation:Bosons in an optical lattice

Bandstructure of an homogeneous optical lattice

Theory

At this first moment, we shall look at the simplest case, i.e. a single particle of mass m which experiences a periodic potential $V(\vec{r})$, where

$V(\vec{r}) = \sum_{x_\alpha = x,y,z} V_0^{x_\alpha} \sin^2 (\pi x_\alpha)$

in the units of recoil energy $E_r^\alpha = \frac{\hbar^2}{2m} \left( \frac{2\pi}{\lambda_\alpha} \right)^2$ and lattice spacing $\frac{\lambda_\alpha}{2}$.

The quantum mechanical behaviour of the single particle follows

$\left[\frac{1}{\pi^2} \left( -i \nabla + 2\pi \vec{k} \right)^2 + \sum_{x_\alpha = x,y,z} V_0^{x_\alpha} \sin^2 (\pi x_\alpha) \right] u_k (\vec{r}) = \epsilon_k u_k(\vec{r})$

which is clearly separable to say the x-component:

$\left[\frac{1}{\pi^2} \left( -i \partial_x + 2\pi k_x \right)^2 + V_0^{x} \sin^2 (\pi x) \right] u_{k_x} (x) = \epsilon_{k_x} u_{k_x}(x) \,\,\,\,\, ,\,\,\,\,\, k_x = 0, \frac{1}{L_x} ,\cdots \frac{L_x-1}{L_x}$

In the plane wave basis,

$u_{k_x} (x) = \frac{1}{\sqrt{L_x}} \sum_{m \in \mathbf{Z}} c_m^{(k_x)} e^{i2m\pi x}$

we arrive at a tridiagonal diagonalization problem:

$\left[ 4(m + k_x)^2 + \frac{V_0^x}{2} \right] c_m^{(k_x)} - \frac{V_0^x}{4} c_{m-1}^{(k_x)} - \frac{V_0^x}{4} c_{m+1}^{(k_x)} = \epsilon_{k_x} c_m^{(k_x)}$

The wannier function is defined as:

$w(x) = \frac{1}{\sqrt{L_x}} \sum_{k_x} u_{k_x} (x) e^{i 2\pi k_x x} = \frac{1}{L_x} \sum_{k_x} \sum_{m \in \mathbf{Z}} c_m^{(k_x)} e^{i 2\pi (m+k_x) x}$

and from there, one can calculate the onsite interaction:

$U = g \int | w(x) |^4 dx = \frac{4 \pi a_s \hbar^2}{m} \int | w(x) |^4 dx$

After a little bit of algebra, we arrive at the hopping strength:

$t = -\frac{1}{L_x} \sum_{k_x} \epsilon_{k_x} e^{-i2\pi k_x}$

Finally, the Fourier transform of the wannier function is:

$\tilde{w}(q_x) = \frac{1}{\sqrt{L_x}} \int w(x) e^{-i2\pi q_x x} dx = \frac{1}{\sqrt{L_x}} \sum_{k_x} \sum_{m \in \mathbf{Z}} c_m^{(k_x)} \delta_{q_x, k_x+m}$

Implementation in Python

An example

For instance:

import numpy;
import pyalps.dwa;

V0   = numpy.array([8. , 8. , 8.]);      # in recoil energies
wlen = numpy.array([843., 843., 843.]);  # in nanometer
a    = 114.8;                            # s-wave scattering length in bohr radius
m    = 86.99;                            # mass in atomic mass unit
L    = 200;                              # lattice size (along 1 direction)

band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);

A first glance of the band structure:

>>> band

Optical lattice: 
================
V0    [Er] = 8	8	8	
lamda [nm] = 843	843	843	
Er2nK      = 154.89	154.89	154.89	
L          = 200 
g          = 5.68473

Band structure:
===============
t [nK] : 4.77051	4.77051	4.77051	
U [nK] : 38.7018
U/t    : 8.11272	8.11272	8.11272	

wk2[0 ,0 ,0 ] : 5.81884e-08
wk2[pi,pi,pi] : 1.39558e-08

Well, the values of t (nK), U (nK), and U/t can be obtained via:

>>> numpy.array(band.t())
array([ 4.77050984,  4.77050984,  4.77050984])
>>>
>>> numpy.array(band.U())
array(38.7018197381118)
>>>
>>> numpy.array(band.Ut())
array([ 8.11272192,  8.11272192,  8.11272192])

In momentum ($\vec{q}$) space, the (squared) wannier function $|\tilde{w}(\vec{q})|^2$ can be obtained in the x-direction from:

>>> numpy.array(band.q(0))
array([-5.   , -4.995, -4.99 , ...,  5.985,  5.99 ,  5.995])
>>> 
>>> numpy.array(band.wk2(0))
array([  7.57249518e-15,   7.88189086e-15,   8.20434507e-15, ...,
         1.62988573e-18,   1.56057426e-18,   1.49429285e-18])

and the y- or z- direction by replacing the index 0 to 1 and 2 respectively.

Model hamiltonian

The DWA code simulates the single band boson Hubbard model

$\hat{H} = -t \sum_{\langle i,j \rangle} \hat{b}_i^+ \hat{b}_j + \frac{U}{2} \sum_i \hat{n}_i (\hat{n}_i - 1) - \sum_i ( \mu - V_T ( \vec{r}_i) ) \hat{n}_i$

with hopping strength $t$, onsite interaction strength $U$, and chemical potential $\mu$ at finite temperature $T$ via Quantum Monte Carlo implemented in the directed worm algorithm. Here, $\hat{b}$ ($\hat{b}^+$) is the annihilation (creation) operator, and $\hat{n}_i$ being the number operator at site i. Bosons in an optical lattice are confined, say in a 3D parabolic trapping potential, i.e.

$V_T (\vec{r}_i) = K_x x_i^2 + K_y y_i^2 + K_z z_i^2$.

due to the gaussian beam waists as well as other sources of trapping. At finite temperature T, the physics is essentially captured by the partition function

$Z = \mathrm{Tr} \, \exp \left(-\beta \hat{H} \right)$

and physical quantities such as the local density

$\langle n_i \rangle = \frac{1}{Z} \, \mathrm{Tr} \, \hat{n}_i \exp \left(-\beta \hat{H} \right) = \frac{1}{Z} \sum_{\mathcal{C}} n_i (\mathcal{C}) Z(\mathcal{C})$

for some configuration $\mathcal{C}$ in the complete configuration space, with inverse temperature $\beta = 1/T$. Here, the units will be cleverly normalized later on.

Momentum distribution and time-of-flight images

The momentum distribution $\langle n(\vec{k})\rangle$ of bosons in an optical lattice at equilbrium can be expressed as: $\langle \hat{\psi}^+ (\vec{k}) \, \hat{\psi} (\vec{k}) \rangle = \frac{1}{N} \int d\vec{r} \, d\vec{r}' \langle \hat{\psi}^+ (\vec{r}) \, \hat{\psi} (\vec{r}') \rangle \, e^{i \vec{k} \cdot (\vec{r} - \vec{r}')} = \frac{1}{N} \sum_{i,j} \int d\vec{r} \, d\vec{r}' \, w(\vec{r} - \vec{r}_i) e^{i\vec{k}\cdot(\vec{r} - \vec{r}_i)} \, w(\vec{r}' - \vec{r}_j) e^{-i\vec{k}\cdot(\vec{r}' - \vec{r}_j)} \langle \hat{b}^+_i \hat{b}_j \rangle \, e^{i \vec{k} \cdot (\vec{r}_i - \vec{r}_j)}$

Expressing the fourier transform of the wannier function $w(\vec{r})$ as

$\tilde{w}(\vec{k}) = \frac{1}{\sqrt{N}} \int d\vec{r} \, w(\vec{r}) e^{-i\vec{k}\cdot\vec{r}} \,\, ,$

and the interference term:

$S(\vec{k}) = \sum_{i,j} \langle \hat{b}^+_i \hat{b}_j \rangle \, e^{i \vec{k} \cdot (\vec{r}_i - \vec{r}_j)} \,\, ,$

the momentum distribution reduces to the following form:

$\langle n(\vec{k})\rangle = \left| \tilde{w}(\vec{k})\right|^2 S(\vec{k}) \,\, .$

To measure it in experiments, the optical lattice is momentarily switched off and the bosons expand freely by assumption with momentum gained from the lattice momentum previously. Measurements are performed in our classical world, and therefore semiclassical treatment is already sufficient, i.e.

$\hbar \vec{k} = m \left( \frac{\vec{r} }{t_f} \right)$

where $t_f$ is the time-of-flight taken by the bosons to move from the origin (experiment) to the detector probe at position $\vec{r}$. Here, we assume the simplest picture, ie. 1) no interaction, 2) no collision, and 3) uniform motion.

Of course, at equilibrium, i.e. "after the bosons have travelled for a long time", the time-of-flight image will capture a density of

$\langle n_\mathrm{tof} (\vec{r} ) \rangle = \left( \frac{m}{\hbar t_f} \right)^3 \left| \tilde{w}\left( \frac{m}{\hbar t} \vec{r} \right) \right|^2 S(\vec{k})$

which indirectly probes the original momentum distribution of the lattice bosons. This form is only correct in the far-field limit or in the limit of infinite time-of-flight (tof).

To improve upon the accuracy, we shall correct it with semiclassical dynamics during the time of flight. For the bosons that originate from lattice site at $\vec{r}_j$ to the detector probe at $\vec{r}$, their time-dependent wavefunction must be:

$\Psi(\vec{r},t) = \left( \frac{m}{\hbar t_f} \right)^{3/2} \tilde{w}\left( \frac{m}{\hbar t_f} (\vec{r} - \vec{r}_j) \right) \, \exp \left( \frac{i \epsilon_k t_f}{\hbar} \right)$

where the kinetic energy of the free bosons must be:

$\epsilon_k = \frac{\hbar^2 k^2}{2m} = \frac{\hbar^2}{2m} \left( \frac{m}{\hbar t_f} \right)^2 (\vec{r} - \vec{r}_j )^2$

Therefore, the corrected time-of-flight image must be:

$\langle n_\mathrm{tof} (\vec{r} ) \rangle \approx \left( \frac{m}{\hbar t_f} \right)^3 \left| \tilde{w}\left( \frac{m}{\hbar t} \vec{r} \right) \right|^2 \sum_{i,j} \, \langle \hat{b}^+_i \hat{b}_j \rangle \exp \left( - \frac{i}{\hbar} \frac{\hbar^2}{2m} \left( \frac{m}{\hbar t_f} \right)^2 (\vec{r} - \vec{r}_i )^2 t_f \right) \, \exp \left( \frac{i}{\hbar} \frac{\hbar^2}{2m} \left( \frac{m}{\hbar t_f} \right)^2 (\vec{r} - \vec{r}_j )^2 t_f \right)$

where the dependence of the initial site on the wannier enveloped function is neglected. Finally, we arrive at:

$\langle n_\mathrm{tof} (\vec{k} ) \rangle \approx \left| \tilde{w} ( \vec{k} ) \right|^2 \sum_{i,j} \, \langle \hat{b}^+_i \hat{b}_j \rangle \, e^{i \vec{k} \cdot (\vec{r}_i - \vec{r}_j) - i \gamma_f (\vec{r}_i^2 - \vec{r}_j^2) } \,\, .$

with time-of-flight phase $\gamma_f = \frac{m}{2\hbar t_f}$.

For example, the time-of-flight phase $\gamma$ for Rb-87 bosons (mass 86.99 a.m.u.) with a time-of-flight of 15.5 ms

a. expanded from an isotropic optical lattice of wavelength 765 nm:

>>> pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=765, mass=86.99)
0.006464602556863682

b. expanded from an optical lattice of wavelength 843 nm along x-direction and 765 nm along y- and z- directions:

>>> pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[843,765,765], mass=86.99)
[0.007850080468935228, 0.006464602556863682, 0.006464602556863682]

The following table summarizes the basic option to turn on time-of-flight imaging options.

Option	Default	Remark
tof_phase	0.	time-of-flight phase $\gamma$

Symmetry leads to further simplification

Symmetries:

$\langle n_\mathrm{tof} (\vec{k} ) \rangle = \langle n_\mathrm{tof} (-\vec{k} ) \rangle$

$\langle \hat{b}^+_i \hat{b}_j \rangle = \langle \hat{b}^+_j \hat{b}_i \rangle$

Therefore:

$\langle n_\mathrm{tof} (\vec{k} ) \rangle \approx \left| \tilde{w} ( \vec{k} ) \right|^2 \sum_{i,j} \, \langle \hat{b}^+_i \hat{b}_j \rangle \, \cos( \vec{k} \cdot (\vec{r}_i - \vec{r}_j) ) \, \cos(\gamma_f (r_i^2 - r_j^2)) \,\, .$

Define green function:

$g_f(i,j) = \langle \hat{b}^+_i \hat{b}_j \rangle \, \cos(\gamma_f (r_i^2 - r_j^2))$

In reduced coordinates, $\vec{r}_\alpha = \vec{r}_i - \vec{r}_j$ , we redefine green function:

$g_f(\alpha) = \sum_{i,j : \, \vec{r}_i - \vec{r}_j \, = \, \pm \vec{r}_\alpha} \langle \hat{b}^+_i \hat{b}_j \rangle \, \cos(\gamma_f (r_i^2 - r_j^2))$

and therefore the momentum distribution/ time-of-flight image:

$\langle n_\mathrm{tof} (\vec{k} ) \rangle \approx \left| \tilde{w} ( \vec{k} ) \right|^2 \sum_{\alpha : \, |r_\alpha| \geq 0} \, g_f (\alpha) \cos( \vec{k} \cdot (\vec{r}_\alpha) ) \,\, .$

Here, we distinguish the momentum distribution from the interference

$\langle n_k (\vec{k} ) \rangle = \sum_{\alpha : \, |r_\alpha| \geq 0} \, g_f (\alpha) \cos( \vec{k} \cdot (\vec{r}_\alpha) ) \,\, .$

Please refer to the list of measurement observables provided by the DWA code.