\chapter{Implementing quantitative simulations}
\label{appendix_dwa}

\section{Optical lattice bandstructure}
\label{chapter_bandstructure}

A three-dimensional optical lattice is setup by 3 pairs of counter-propagating laser beams in orthogonal directions with lattice intensity $\vec{V}_0 = (V_{0x}, V_{0y}, V_{0z})$ and wavelength $\vec{\lambda} = (\lambda_x, \lambda_y, \lambda_z)$.
In units of lattice spacing $\vec{d} = \vec{\lambda}/2$, and recoil energy $\vec{E}_R = \frac{h^2}{2m} \left( \frac{1}{\lambda_x^2} , \frac{1}{\lambda_y^2} , \frac{1}{\lambda_z^2} \right)$, the single particle Bloch hamiltonian \cite{jaksch1} reads
%
\begin{equation}
\hat{H}_k = \frac{1}{\pi^2} \left( -i\nabla + 2\pi \vec{k} \right)^2 + \sum_{\alpha=x,y,z} V_{0\alpha} \sin^2 (\pi x_\alpha) 
\end{equation}
%
which is clearly separable.
Along a direction, say $x$, the eigenvalue problem is reduced to
%
\begin{equation}
\label{evalue1}
\left[ \frac{1}{\pi^2} \left( -i\partial_x + 2\pi k_x \right)^2 + V_{0x} \sin^2 (\pi x) \right] u_{k_x} (x) = \epsilon_{k_x} u_{k_x} (x)
\end{equation}
%
with $k_x = 0, \frac{1}{L} , \cdots, \frac{L-1}{L}$ for a periodic lattice of size $L^3$.
Expanding in the plane wave basis, 
%
\begin{equation}
u_{k_x} (x) = \frac{1}{\sqrt{L}}\sum_{G = 0, \pm 1 , \pm 2 \cdots} c_G^{(k_x)} e^{i2\pi G x}  \,\, ,
\end{equation}
%
equation \ref{evalue1} becomes the following tridiagonal eigenvalue problem
%
\begin{equation}
\label{evalue1}
\left( 4 \left( G + k_x \right)^2 + \frac{V_{0x}}{2} \right) c_G^{(k_x)} \, - \, \frac{V_{0x}}{4} c_{G-1}^{(k_x)} \,- \, \frac{V_{0x}}{4} c_{G+1}^{(k_x)} = \epsilon_{k_x} c_G^{(k_x)}  \,\, .
\end{equation}
%
The wannier function is defined as
%
\begin{equation}
w_x(x) = \frac{1}{\sqrt{L}} \sum_{k_x} u_{k_x} (x) e^{i2\pi k_x x} = \frac{1}{L} \sum_{k_x} \sum_G c_G^{(k_x)} e^{i2\pi (G+k_x)x}  \,\, .
\end{equation}
%
Important bandstructure parameters are the hopping strength
%
\begin{equation}
t_x = -\frac{1}{L} \sum_{k_x} \epsilon_{k_x} e^{-i2\pi k_x}
\end{equation}
%
(likewise for $t_y$ and $t_z$), and the onsite interaction strength
%
\begin{equation}
U = g \int |w(\vec{r})|^4 \, d\vec{r}  = \frac{4\pi a_s \hbar^2}{m} \int |w(\vec{r})|^4 \, d\vec{r}
\end{equation}
%
where $a_s$ is the s-wave scattering length easily tunable via Feshbach resonance technique.
Finally, the Fourier transform of the wannier function is
%
\begin{equation}
\tilde{w}_x(q_x) = \frac{1}{\sqrt{L}} \int w_x (x) e^{-i2\pi q_x x} \, dx = \frac{1}{\sqrt{L}} \sum_{k_x} \sum_G c_G^{(k_x)} \delta_{q_x,G+k_x} \,\, .  
\end{equation}
%
\\ \\
\noindent Easy implementation is available in the Python interface, for example:
{\color{blue}{
\begin{lstlisting}
>>> import numpy;
>>> import pyalps.dwa;
>>>
>>> V0   = numpy.array([8.805, 8. ,  8. ]); #lattice strength [Er]
>>> wlen = numpy.array([765., 843., 843.]); #laser wavelength [nm]
>>> a    = 101;    #s-wave scattering length [bohr radius]
>>> m    = 86.99;  #mass [a.m.u.]
>>> L    = 160;    #lattice of size L^3
>>>
>>> band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
>>>
>>> band

Optical lattice:
================
V0    [Er] = 8.805	8	8	
lamda [nm] = 765	843	843	
Er2nK      = 188.086	154.89	154.89	
L          = 160
g          = 5.51132

Band structure:
===============
t [nK] : 4.77257	4.77051	4.77051	
U [nK] : 38.7027
U/t    : 8.1094	8.1129	8.1129	
\end{lstlisting}
}}


\section{Expression for finite time-of-flight}
\label{chapter_tof}

To probe for momentum distribution, the optical lattice is momentarily turned off, and the bosons expand freely with momentum gained from the previous lattice momentum.
Measurements are performed in our classical world, and therefore semiclassical treatment is already sufficient, i.e.
%
\begin{equation}
\hbar \vec{k} = m \left( \frac{\vec{r}}{t_f} \right)
\end{equation}
%
where $t_f$ is the actual time of flight taken by the bosons to move from the origin (experiment) to the detector probe at position $\vec{r}$.
Here, we assume that
%
\begin{enumerate}
\item there is no interaction and collision among the bosons as they fly apart;
\item the bosons move at constant speed during their entire flight.
\end{enumerate}
%
In the limit of infinite $t_f$, the time-of-flight image captures
%
\begin{equation}
\langle n_f (\vec{r}) \rangle = \left( \frac{m}{\hbar t_f} \right)^3 \left| \tilde{w} \left( \frac{m}{\hbar t_f} \vec{r} \right) \right|^2 S(\vec{k})  \,\, .
\end{equation}
%
To start with the correction due to dynamics, we consider the time evolution of the time-dependent wannier wavepacket that orginates from site j
%
\begin{equation}
\mathcal{W}_j (\vec{r},t_f) =  \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)  \,\, ,
\end{equation}
%
where 
%
\begin{equation}
\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  \,\, .
\end{equation}
%
Ignoring the initial-site dependence on the wannier enveloped function, the time-of-flight image captures
%
\begin{eqnarray}
\langle n_f (\vec{r}) \rangle 
 &\approx&  \left( \frac{m}{\hbar t_f} \right)^3 \left| \tilde{w} \left( \frac{m}{\hbar t_f} \vec{r} \right) \right|^2 \sum_{i,j} \langle \hat{b}_i^\dag \hat{b}_j \rangle \nonumber \\
 & & \times  \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)  \nonumber \\
 & & \times  \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) 
\end{eqnarray}
%
at finite $t_f$, therefore arriving at the final expression
%
\begin{equation}
\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)} \,\, .  
\tag{\ref{tof_eq}}
\end{equation}
%
The time-of-flight phase $\gamma_f$ can be easily evaluated in Python, for example:
{\color{blue}{
\begin{lstlisting}
>>> import pyalps.dwa
>>> pyalps.dwa.tofPhase(time_of_flight=15.5, wavelength=[843,765,765], mass=86.99)
[0.007850080468935228, 0.006464602556863682, 0.006464602556863682]
\end{lstlisting}
}}

\section{Obtaining time-of-flight image}
\label{chapter_tof_image}

This section illustrates how one can easily obtain the time-of-flight images, say figure \ref{fig_tof_image} for example, using the Python interface.
%
{\color{blue}{
\begin{lstlisting}
import numpy;
import numpy.fft;
import pyalps;
import pyalps.dwa;
import matplotlib.pyplot;
from mpl_toolkits.mplot3d.axes3d import Axes3D;

resultFile       = pyalps.getResultFiles()[0];
L                = int(pyalps.getParameters(resultFile)[0]['L']);
green_tof        = pyalps.getMeasurements(resultFile, observable='Green Function:TOF')[0]['mean']['value'].reshape([L,L,L]);
momentum_density = numpy.fft.fftn(green_tof).real; 

V0   = numpy.array([8.8, 8.,  8.]);
wlen = numpy.array([765., 843., 843.]);
a    = 101;
m    = 86.99;
band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);

q_z   = numpy.array(band.q(2));
wk2_z = numpy.array(band.wk2(2));
wk2_z = numpy.array([q_z, wk2_z]).transpose();
wk2_z0= wk2_z[wk2_z[:,0] == 0.][0][1];
wk2_z = numpy.transpose(wk2_z[wk2_z[:,1]/wk2_z0 > 1e-4]);
q_z   = wk2_z[0]
wk2_z = wk2_z[1]

momentum_density = numpy.tile(momentum_density, reps=(1,1,2*((q_z[q_z[:] >= 0].size / L)+1)));

dummy = numpy.zeros(momentum_density.shape[2]);
dummy[dummy.size/2:dummy.size/2 + q_z[q_z[:] >= 0].size] = wk2_z[q_z[:] >= 0]
dummy[dummy.size/2-q_z[q_z[:] < 0].size:dummy.size/2]  = wk2_z[q_z[:] < 0]    

momentum_density = numpy.tensordot(momentum_density, dummy, axes=([2],[0]))
momentum_density = numpy.tile(momentum_density, reps=(4,4));

q_x   = numpy.array(band.q(0));
wk2_x = numpy.array(band.wk2(0));
wk2_x = wk2_x[(q_x[:] >= -2.)*(q_x[:] < 2.)];
q_x   = q_x[(q_x[:] >= -2.)*(q_x[:] < 2.)]

q_y   = numpy.array(band.q(1));
wk2_y = numpy.array(band.wk2(1));
wk2_y = wk2_y[(q_y[:] >= -2.)*(q_y[:] < 2.)];
q_y   = q_y[(q_y[:] >= -2.)*(q_y[:] < 2.)]
 
wk2 = numpy.outer(wk2_x, wk2_y);

q_x = numpy.array([q_x] * wk2.shape[1], float);
q_y = numpy.array([q_y] * wk2.shape[0], float).transpose();
tof_image = wk2 * momentum_density;

mag = L/10;
q_x = q_x[0:q_x.shape[0]:mag, 0:q_x.shape[1]:mag];
q_y = q_y[0:q_y.shape[0]:mag, 0:q_y.shape[1]:mag];
tof_image = tof_image[0:tof_image.shape[0]:mag, 0:tof_image.shape[1]:mag] * mag * mag

fig = matplotlib.pyplot.figure();
ax = fig.add_subplot(1, 1, 1, projection='3d') 
surf = ax.plot_surface( q_y, q_x, tof_image, 
                        rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=10);
fig.show();
\end{lstlisting}
}}

\noindent Further detailed information can be found in reference \cite{alps_webpage}.

\newpage

\section{Obtaining density profile}
\label{chapter_density_profile}

This section illustrates how one can easily obtain the density profiles, say figure \ref{fig_density_profile} for example, using the Python interface.
%
{\color{blue}{
\begin{lstlisting}
import numpy;
import pyalps;
import pyalps.dwa;
import matplotlib.pyplot;
from mpl_toolkits.mplot3d.axes3d import Axes3D;

resultFile = pyalps.getResultFiles()[0];
L          = int(pyalps.getParameters(resultFile)[0]['L']);
density    = pyalps.getMeasurements(resultFile, observable='Local Density')[0]['mean']['value'].reshape([L,L,L]);
cross_section_density = density[:,:,L/2];
column_integrated_density = numpy.sum(density, axis=2);
x = numpy.array([range(L)] * L, float) - (L-1)/2.;
y = numpy.transpose(x);

# Plotting cross sectional density profile
fig1 = matplotlib.pyplot.figure();
ax = fig1.add_subplot(1, 1, 1, projection='3d') 
surf = ax.plot_surface( x, y, cross_section_density, 
                        rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=False)
fig1.colorbar(surf, shrink=0.5, aspect=10);
fig1.show();

# Plotting column integrated density profile
fig2 = matplotlib.pyplot.figure();
ax = fig2.add_subplot(1, 1, 1, projection='3d') 
surf = ax.plot_surface( x, y, column_integrated_density, 
                        rstride=1, cstride=1, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=False)
fig2.colorbar(surf, shrink=0.5, aspect=10);
fig2.show();
\end{lstlisting}
}}

\noindent Further detailed information can be found in reference \cite{alps_webpage}.

\newpage

\section{Fixing $\langle N \rangle$ by tuning chemical potential}
\label{chapter_tuning_chemical_potential}

In this section, we shall work through part of the example of quantifying the error budget due to a $\pm 5 \%$ fluctuation in lattice strength $V_{0x}$ to illustrate how easily the chemical potential can be conviently tuned in the Python interface.
(See table \ref{table:errorbudgets}.) 

\subsection{Initial simulations}

The first step is to scan a certain range of chemical potentials $\mu$ in the initial simulations.
\\ \\
\noindent The following prepares the parameter files for the initial simulations in the Python interface:
%
{\color{blue}{
\begin{lstlisting}
import numpy;
import pyalps;
import pyalps.dwa;

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

ps = [];
for ratio in [0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04]:
  V0   = numpy.array([ratio*8.805, 8. ,  8. ]);
  wlen = numpy.array([765., 843., 843.]);
  a    = 101;
  m    = 86.99;
  L    = 160;
 
  band = pyalps.dwa.bandstructure(V0, wlen, a, m, L);
  p = {};
  p.update({'ratio_V0x' : ratio , 'U' : band.Ut()[2] , 'tx_t' : band.t()[0]/band.t()[2]})
  ps.append(p);

params=[]
for p in ps:
  for ratio in [0.8, 0.82, 0.84, 0.88, 0.9, 0.92, 0.94, 0.96, 0.98, 1.00, 1.02, 1.04, 1.06, 1.08, 1.10, 1.12, 1.14, 1.16, 1.18, 1.2]:
    params.append(
      { 'LATTICE'         : 'inhomogeneous simple cubic lattice'   
      , 'MODEL'           : 'boson Hubbard'                        
      , 'L'               : 160                             
      , 'Nmax'            : 20
      , 'ratio_V0x'       : p['ratio_V0x']                 
      , 't'               : 1.
      , 'tx_t'            : p['tx_t']                 
      , 'U'               : p['U']               
      , 'T'               : 1.                 
      , 'mu_homogeneous'  : ratio*4.0265                
      , 'mu'              : 'mu_homogeneous - (0.0073752*(x-(L-1)/2.)*(x-(L-1)/2.) + 0.0036849*(y-(L-1)/2.)*(y-(L-1)/2.) + 0.0039068155*(z-(L-1)/2.)*(z-(L-1)/2.))'
      , 'tof_phase'       : str(tof_phase) 
      , 'SWEEPS'          : 100000             
      , 'SKIP'            : 100                 
      }
    )

h5_infiles = pyalps.writeInputH5Files("errorbudget.V0",params);
\end{lstlisting}
}}
\noindent The next step is thermalization, i.e.
{\color{blue}{
\begin{lstlisting}
for taskfile in h5_infiles:
  pyalps.dwa.recursiveRun( 
    "pyalps.runApplication('dwa', taskfile)" ,  
    cmd_lang = 'python' ,
    follow_up_script = "pyalps.dwa.extract_worldlines(infile=pyalps.input2output(taskfile), outfile=taskfile)" ,
    break_if = "pyalps.dwa.thermalized(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Thermalizing')" ,
    loc = locals()
  );
\end{lstlisting}
}}
\noindent followed by convergence of $\langle N \rangle$, i.e.
{\color{blue}{
\begin{lstlisting}
for taskfile in h5_infiles:  
  pyalps.dwa.recursiveRun( 
    "pyalps.runApplication('dwa', pyalps.input2output(taskfile))" ,  
    cmd_lang = 'python' ,
    break_if = "pyalps.dwa.converged(pyalps.input2output(taskfile), 'Total Particle Number', simplified=True)" ,
    write_status = "pyalps.dwa.write_status(pyalps.input2output(taskfile), 'Converging')" ,
    loc = locals() 
  );
\end{lstlisting}
}}

\newpage

\subsection{Tuning chemical potential}

Next, we prepare the new parameter files for a new interpolated value of chemical potentials:
{\color{blue}{
\begin{lstlisting}
resultFiles = pyalps.getResultFiles(prefix='errorbudget.V0');
data = pyalps.loadMeasurements(resultFiles, 'Total Particle Number');
params = pyalps.paramsAtFixedY(data, x='mu_homogeneous', y='Total Particle Number', foreach=['ratio_V0x'], fixedY=280000);

h5_infiles = pyalps.writeInputH5Files("errorbudget.V0.1",params);
\end{lstlisting}
}}
\noindent and we perform the new simulations similar to before.
This process is repeated until one is satisfied with the accuray of the interpolated chemical potential.

\subsection{Turning on measurements}

After one is satisfied with the correct chemical potential, the measurements for density profile and green function can be turned on
{\color{blue}{
\begin{lstlisting}
for h5_infile in h5_infiles:
  pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Local Density]', 1);
  pyalps.dwa.switchParameter(h5_infile, 'MEASURE[Green Function]', 1);
\end{lstlisting}
}}
\noindent for the simulations to run and collect statistics till satisfaction.
\\ \\
\noindent Further detailed information can be found in reference \cite{alps_webpage}.
