# Momentum distribution and time-of-flight images

## Theory

### Momentum distribution

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}) \,\, .$

### Time of flight images (in experiment)

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.