\chapter{Quantitative simulations}

This chapter is devoted to quantitative simulations of bosons trapped in an actual-size 3D optical lattice with simulation parameters mimicked as realistically to the experiments as possible, which are quantitatively modelled by the boson Hubbard model \cite{fisher}, i.e.
%
\begin{equation}
\hat{H} - \mu N = - \sum_{\langle i,j \rangle} t_{ij} \, \hat{b}_i^\dag \hat{b}_j + \sum_i \frac{U_i}{2} n_i (n_i -1) - \sum_i \mu_i n_i  \,\, ,
\end{equation}
%
and efficiently solved by the Quantum Monte Carlo (directed) worm algorithm \cite{prokofev1} (\cite{pollet1}).
The spirit behind this algorithm has been discussed in detail previously in Chapter \ref{chapter_QMC}, and its implementation fully outlined in Appendix \ref{appendix_dwa}.
\\ \\
\noindent This chapter is based on one of my publications \cite{implementing_dwa}, which highlights this state-of-art algorithmic implementation to be one of the fastest, if not the fastest, in the world till date.
In particular, current simulations are ultra robust, which are capable of handling up to 10 million lattice sites.
Of which, some manage to converge results in no longer than a couple of hours.


\section{Mimicking the experiments}

Mimicking the experiments by Bloch \cite{trotsky}, the optical lattice anisotropy gives rise to direction-dependent band structure parameters (refer to appendix \ref{chapter_bandstructure}), i.e. 
%
\begin{eqnarray}
t_{ij} &=& t_x \,\, , \,\,  t_y \,\, , \mathrm{or} \,\, t_z \nonumber \\
U_i &=& U \nonumber \\
\mu_i &=& \mu - (V_{Tx} x_i^2 +V_{Ty} y_i^2 + V_{Tz} z_i^2) \nonumber
\end{eqnarray}
%
which become increasingly valid in the close vicinity towards the centre of the optical lattice trap. 

\subsection{Momentum distribution and time-of-flight image}

The momentum distribution of bosons in an optical lattice at equilibrium is 
%
\begin{eqnarray}
\langle n(\vec{k}) \rangle 
 &=& \langle \hat{\psi}^\dag (\vec{k}) \hat{\psi} (\vec{k}) \rangle   \\
 &=& \frac{1}{L^3}  \int d\vec{r} \, d\vec{r}'   \langle \hat{\psi}^\dag (\vec{r}) \hat{\psi} (\vec{r}') \rangle \, e^{i\vec{k}\cdot (\vec{r} - \vec{r}')}  \nonumber \\
 &=& \frac{1}{L^3} \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^\dag \hat{b}_j \rangle \, e^{i\vec{k}\cdot(\vec{r}_i - \vec{r}_j)}  \nonumber   \,\, .
\end{eqnarray}
%
Defining the Fourier transform of the wannier function $w(\vec{r})$ as
%
\begin{equation}
\tilde{w}(\vec{k}) = \frac{1}{\sqrt{L^3}} \int d\vec{r} \, w(\vec{r}) \, e^{-i\vec{k}\cdot \vec{r}}  \,\, ,
\end{equation}
%
and the interference term 
%
\begin{equation}
S(\vec{k}) = \sum_{i,j} \langle \hat{b}_i^\dag \hat{b}_j \rangle \, e^{i\vec{k}\cdot(\vec{r}_i - \vec{r}_j)} \,\,
\end{equation}
%
the momentum distribution becomes
%
\begin{equation}
\langle n(\vec{k}) \rangle = |\tilde{w}(\vec{k})|^2 S(\vec{k})  \,\, .
\end{equation}
%
which is indirectly probed in experiements by a time-of-flight image \cite{trotsky}, i.e.
%
\begin{equation}
\label{tof_eq}
\langle n_f (\vec{k}) \rangle \approx |\tilde{w}(\vec{k})|^2 \sum_{i,j} \langle \hat{b}_i^\dag \hat{b}_j \rangle \, e^{i\vec{k}\cdot (\vec{r}_i - \vec{r}_j) - i \gamma_f (r_i^2 - r_j^2)}  \,\, ,
\end{equation}
%
having time-of-flight phase $\gamma_f = \frac{m}{2\hbar t_f}$ with $t_f$ being the actual time of flight. 
(See Appendix \ref{chapter_tof} for its derivation.)
\\ \\
\noindent Exploiting the symmetries
%
\begin{eqnarray}
\langle n_f (\vec{k}) \rangle &=& \langle n_f (-\vec{k}) \rangle \\
\langle \hat{b}_i^\dag \hat{b}_j \rangle &=& \langle \hat{b}_j^\dag \hat{b}_i \rangle 
\end{eqnarray}
%
equation \ref{tof_eq} reduces to
%
\begin{equation}
\langle n_f (\vec{k}) \rangle \approx |\tilde{w}(\vec{k})|^2 \sum_{|r_\alpha| \geq 0} g_f(\vec{r}_\alpha) \cos (\vec{k} \cdot \vec{r}_\alpha)  \,\, ,
\end{equation}
%
with time-of-flight Green function defined as
%
\begin{equation}
g_f(\vec{r}_\alpha) = \sum_{i,j \, : \, \vec{r}_i - \vec{r}_j = \pm \vec{r}_\alpha}   \langle \hat{b}_i^\dag \hat{b}_j \rangle  \cos ( \gamma_f (r_i^2 - r_j^2) )  \,\, .
\end{equation}
%
Last but not least, the time-of-flight image observed in experiments are column-integrated along the line-of-sight, say the z-direction:
%
\begin{equation}
\label{tof_eq2}
\langle n_f (k_x, k_y) \rangle = \int dk_z \langle n_f (\vec{k}) \rangle \,\, ,
\end{equation}
%
as illustrated in figure \ref{fig_tof_image}. 
%
\begin{figure}[h]
\includegraphics[width=8cm]{./figures/parm-validation-task1-tof-image.eps}
\includegraphics[width=8cm]{./figures/parm-validation-task2-tof-image.eps}
\caption{Time of flight images (equation \ref{tof_eq2}) obtained from QMC-DWA simulations mimicking optical lattice experiments as realistically as possible, with interaction strength $U/t=8.11$ (left) , $27.5$ (right) at temperature $T/t=1$.
The anisotropic optical lattice, with lattice strength $\vec{V}_0 = (8.8 E_{rx}, 8 E_{ry} , 8 E_{rz})$ (left), $(12.64 E_{rx}, 11.75 E_{ry} , 11.75 E_{rz})$ (right), and laser wavelength $\vec{\lambda} = (765,843,843)\mathrm{nm}$, confines $2.8 \times 10^5$ (left), $9.4 \times 10^4$ (right) bosons in a parabolic trap $\vec{V}_T = (17.1, 10.9, 11.3)\mathrm{Hz}$ (left), $(19.9, 13.0, 13.4)\mathrm{Hz}$ (right).
The horizontal axes are $k_x$ and $k_y$ in units of $2\pi$, and the vertical axis is the time-of-flight distribution $\langle n_f (k_x, k_y)\rangle$ in unit of inverse momentum area resolution $(\Delta k_x \Delta k_y)^{-1}$, taking experimental value $(\Delta k_x \Delta k_y) \approx 0.1^2$ in units of $(2\pi)^2$.
See appendix \ref{chapter_tof_image} for further details.
}
\label{fig_tof_image}
\end{figure}
%
\\ \\
\noindent Important quantities that derive from the time-of-flight distribution include the condensate fraction
%
\begin{equation}
f_c = \langle n_f (0,0) \rangle  \,\, ,
\end{equation}
%
and the visibility
%
\begin{equation}
\mathcal{V} = \frac{\langle n_f (0,0) \rangle - \langle n_f (\pi,\pi) \rangle}{\langle n_f (0,0) \rangle + \langle n_f (\pi,\pi) \rangle}  \,\, .
\end{equation}
%

\subsection{Density profile}

The density profile of bosons in an optical lattice at equilibrium
%
\begin{equation}
\langle n({\bf{r}}) \rangle = \frac{1}{Z} \, \mathrm{Tr} \, n({\bf{r}}) \exp [ -\beta (\hat{H} - \mu N)] \,\, 
\tag{\ref{eq_local_density}}
\end{equation}
%
\begin{figure}[h]
\includegraphics[width=8cm]{./figures/parm-validation-task1-cross-section-density.eps}
\includegraphics[width=8cm]{./figures/parm-validation-task2-cross-section-density.eps}
\includegraphics[width=8cm]{./figures/parm-validation-task1-column-integrated-density.eps}
\includegraphics[width=8cm]{./figures/parm-validation-task2-column-integrated-density.eps}
\caption{Cross-sectional (top) and column-integrated (bottom) density profiles obtained from QMC-DWA simulations mimicking optical lattice experiments as realistically as possible, with interaction strength $U/t=8.11$ (left) , $27.5$ (right) at temperature $T/t=1$.
The anisotropic optical lattice experiments take the same setup parameters as in figure \ref{fig_tof_image}.
See appendix \ref{chapter_density_profile} for further details.
}
\label{fig_density_profile}
\end{figure}
%
is usually measured cross-sectional $\langle n(x,y,z=0) \rangle$, or column-integrated 
%
\begin{equation}
\label{eq_column_integrated_local_density}
\langle n(x,y) \rangle = \int \langle n(\vec{r}) \rangle \, dz 
\end{equation}
%
in experiments, as illustrated in figure \ref{fig_density_profile}.
%

\section{Beyond parabolic trapping}

In this section, we attempt to quantify the non-parabolic effect of gaussian traps due to the $1/e^2$-waist.  

\subsection{Quantifying error budgets}

The first step is to quantify the error budgets for a certain 3D optical lattice experiemnt due to the various imperfections, as illustrated in table \ref{table:errorbudgets}.
%
\begin{table} [!htp]
\centering
\begin{tabular}{| c | c | c | c | c | c |}
\hline
$\pm 5\%$ fluctuations in  & \multicolumn{5}{|c|}{estimated fluctuation in}     \\
 \cline{2-6} 
& $\langle E \rangle$  & $\langle n(0,0,0) \rangle$  & $\langle n(0,0) \rangle$  &  $f_c$ & $\mathcal{V}$  \\
\hline\hline
$\langle N \rangle$  & < 6\% & < 3\% & < 3\% & < 4\% & < 1\%      \\
$^*a$                       & < 2\% & < 5\% & < 2\% & < 5\% & < 1\%      \\
$^*V_{0x}$              & < 1\% & < 3\% & < 1\% & < 4\% & < 1\%      \\
$^*\lambda_x$        & < 1\% & < 2\% & < 1\% & < 6\% & < 1\%      \\
\hline
\end{tabular}
\caption{Quantifying error budgets in energy $\langle E \rangle$, density at trap centre $\langle n(0,0,0) \rangle$, column-integrated density at trap centre $\langle n(0,0)\rangle$, condensate fraction $f_c$, and visibility $\mathcal{V}$, for a $\pm 5\%$ fluctuation in total particle number $\langle N \rangle$, s-wave scattering length $a$, lattice strength $V_{0x}$, and laser wavelength $\lambda_x$.
The anisotropic optical lattice experiments take the same setup parameters as in figure \ref{fig_tof_image}.
The estimated fluctuation is computed via QMC-DWA simulations, where the simulations indicated by (*) are performed such that the total particle number is fixed at $\langle N \rangle = 2.8 \times 10^5$. 
See appendix \ref{chapter_tuning_chemical_potential} for further details.
}
\label{table:errorbudgets}
\end{table}

\subsection{A first remark of waist correction}

The three-dimensional optical lattice is radially confined by the trapping envelope
%
\begin{equation}
V_0 \exp \left( -\frac{2r^2}{w_0^2} \right)
\tag{\ref{eq_trapping_envelope}}
\end{equation}
%
of waist (or $1/e^2$-radius) $w_0$ \cite{greiner1}. 
Consider the same experiment as figure \ref{fig_tof_image}, where a bosonic cloud spans over a radius $L \sim 50$ in an optical lattice with spacing $d \sim 400 \mathrm{nm}$.
A Gaussian trap with a typical waist of $w_0 \sim 150 \mathrm{\mu m}$ will induce an approximate $4.5\%$ change over the lattice strength $V_0$ from the trap centre.
Table \ref{fig_tof_image} gives us a first idea of the most extreme variations over the physical observables we may measure in the experiments.
Special attention should be paid to the spatial dependence of the lattice strength $V_0$, and thus also the interaction strength $U/t$, in particularly to the next generation of experiments which involve a bigger bosonic cloud size whereby the spatial dependence cannot be neglected anymore.
\\ \\
\noindent In this chapter, however, we confine ourselves in the regime of current experiments, where we neglect the spatial dependence of $V_0$ due to Gaussian beam waist.

\subsection{Corrections to parabolic trapping}

For isotropic optical lattices, the correction to the trapping potential due to waist effects is 
%
\begin{equation}
\label{trap_potential_correction}
\Delta V_T (\vec{r}) = V_0 - \frac{2V_0}{w_0^2} r^2 - V_0 e^{-\frac{2r^2}{w_0^2}}  \,\, ,
\end{equation}
%
which has minimal effect on the measurements for current experiments, as illustrated in figure \ref{fig_waist} for instance.
However, the waist effects will become more prominent in future experiments that involve bigger bosonic clouds. 
%
\begin{figure}[h]
\centering
\includegraphics[width=8cm]{./figures/waist/parm-waist-cross-section-1d-density.eps}
\caption{Cross-sectional density profiles obtained from QMC-DWA simulations for $2.8 \times 10^5$ bosons trapped in an isotropic optical lattice with interaction strength $U/t=8.1$ at temperature $T/t=1$. 
The lattice strength is taken to be 8.35 $E_r$, the laser wavelength 843 nm, and the s-wave scattering length 101 $a_B$.
Blue: The trapping potential is assumed to be parabolic, with trapping frequency 10.5 Hz. 
Red: The trapping potential is corrected according to equation \ref{trap_potential_correction} with $w_0 = 150 \mu m$.
Waist correction reduces the density at the trap center $\langle n(0,0,0) \rangle$ by approximately 1\%, therefore only spreading out slightly in the wings of the bosonic cloud which is however statistically irrelevant.
}
\label{fig_waist}
\end{figure}
%



