\chapter{Optical Lattice}

\noindent This chapter gives an introduction of the quantum effects of identical particles in an optical lattice.


\section{Atomic gases in an optical lattice}


\subsection{Optical lattice}

\begin{figure}[htp]
\centering
\includegraphics[width=7cm]{./figures/optical-lattice-m-yamashita.jpg}
\caption{Atomic gases in a 2D optical lattice. Counter-propagating laser light form a standing wave, or an optical lattice, due to wave interference effect. The lattice intensity is $V_0$, the wavelength $\lambda$ and thus spacing $d=\lambda/2$. {\it{This figure is modified from its original version, courtesy of M. Yamashita.}}\cite{yamashita1}.}
\label{bosonsol_optical_lattice_experiment}
\end{figure}

The physics of wave interference, generated by counter-propagating laser light, dictates the presence of a standing wave. \cite{quantum_optics}
The environment of standing waves is called an optical lattice with intensity $V_0$ and spacing $d=\frac{\lambda}{2}$, which is conceptually illustrated in figure \ref{bosonsol_optical_lattice_experiment}.
For instance, the potential of an isotropic 3D optical lattice reads
%
\begin{equation}
V({\bf{x}}) = \sum_{x^i=x,y,z} V_0 \, \sin^2 (k x^i)
\end{equation}
%
with wavevector $k=\frac{2\pi}{\lambda}=\frac{\pi}{d}$. \cite{greiner1}


\subsection{Single particle in an optical lattice}

To better conceptualise the physics, let us hypothetically consider a single particle of mass {\it{m}} in an optical lattice.
Within non-relativisitic regime, its quantum mechanical nature is captured by the Schr\"dingier equation in reciprocal space $\hat{H}_{\bf{k}} u_{\bf{k}} = E_{\bf{k}} u_{\bf{k}}$. \cite{solid_state}
Due to the periodicity of the optical lattice, the Bloch Hamiltonian is given by
%
\begin{equation}
\hat{H} = \frac{\hbar^2}{2m} \left( -i \nabla + {\bf{k}} \right)^2 + V({\bf{x}})   \,\, ,
\end{equation}
%
where the {\bf{k}}-points in reciprocal space take values depending individually on different lattice geometries.
Figure \ref{bosonsol_single_particle_energy_band_structure} illustrates the energy band structure $E_{\bf{k}}$ of an optical lattice, where one could clearly observe the opening of a energy band gap $\Delta$ as the shallow lattice (small $V_0$) becomes deep (large $V_0$). 
%
\begin{figure}[htp]
\centering
\includegraphics[width=14cm]{./figures/single-particle-energy-band-structure.jpg}
\caption{Single particle energy band structure $E_{\bf{k}}$ of a simple cubic lattice. Energy band gap $\Delta$ opens up with increasing lattice intensity $V_0$, i.e. from $1.0 E_R$ (left), to $2.2 E_R$ (centre), and finally $3.4 E_R$ (right). 
Here, the high-symmetry {\bf{k}}-points of a simple cubic lattice are namely $\Gamma(0,0,0)$, $X(\frac{\pi}{d},0,0)$, $M(\frac{\pi}{d},\frac{\pi}{d},0)$ and $R(\frac{\pi}{d},\frac{\pi}{d},\frac{\pi}{d})$. In addition, energy is expressed in units of recoil energy $E_R=\frac{h^2}{2m\lambda^2}$. Typically, for a Rb-87 particle in an optical lattice with $\lambda=800nm$, the recoil energy is $E_R = 172 nK$.}
\label{bosonsol_single_particle_energy_band_structure}
\end{figure}
%

In the limit of larger $V_0$ such that $\Delta \gg k_BT$ where T is the thermodynamical temperature,  only the ground state band becomes relevant.
By taking the Fourier transform of ground state Bloch functions $u^{(0)}_{\bf{k}}({\bf{x}})$, i.e.
%
\begin{equation}
w({\bf{x}} - {\bf{x}}_i) = \frac{1}{\sqrt{\Omega}} \sum_{\bf{k}} u^{(0)}_{\bf{k}} ({\bf{x}}) e^{i{\bf{k}} \cdot ({\bf{x}} - {\bf{x}}_i)} \,\, ,
\end{equation} 
%
where $\Omega$ is the volume of the primitive unit cell, we obtain an orthonormal basis of Wannier functions $\{w({\bf{x}} - {\bf{x}}_i)\}$ that are maximally localized at every individual lattice site. \cite{jaksch1}
Next, we define the hopping strength of a particle from site {\it{i}} to {\it{j}} as
%
\begin{equation}
t_{ij} = - \int d{\bf{x}} \, w({\bf{x}} - {\bf{x}}_i) \left( -\frac{\hbar^2}{2m} \nabla^2 + V({\bf{x}}) \right) w({\bf{x}} - {\bf{x}}_j)
\end{equation}
%
and the onsite energy $\epsilon_i = t_{ii}$.
Often, only the hopping strength between nearest neighbouring sites is relevant to the physics, and is simply denoted by $t$ for the case of homogeneous lattice.


\subsection{Many particles in an optical lattice}

An atomic gas is a {\it{statistical}} collection of many identical particles.
By {\it{statistics}}, nature has only 2 different genres, namely {\it{symmetric statistics}} for one group of particles called bosons, whereas anti-symmetric statistics for another group called fermions.
By symmetric or antisymmetric statistics, we mean that the many-body wavefunction remains symmetrical or antisymmetrical under particle permutation, i.e. $\psi(\cdots, {\bf{x}}_i , \cdots , {\bf{x}}_j, \cdots) = \pm \psi(\cdots, {\bf{x}}_j , \cdots , {\bf{x}}_i, \cdots)$ respectively. \cite{quantum_mechanics}
To properly handle particle statistics in real space, we work with field operator $\hat{\psi}({\bf{x}})$ in the second quantization language.
Here, the following description in the single band picture is only valid for the case of large $V_0$, otherwise generalizable to include multiband effects. 
\\ \\
\noindent For single component (spinless) bosons, \cite{quantum_mechanics}
%
\begin{equation}
\hat{\psi}({\bf{x}}) = \sum_i w({\bf{x}} - {\bf{x}}_i) \, \hat{b}_i \,\, ,
\end{equation}
%
where $w({\bf{x}} - {\bf{x}}_i)$ is the wannier function at lattice site {\it{i}} and the summation extends over to all sites.
In the bosonic Fock basis $\{|n_i\rangle \, | \,  n_i = 0, 1, 2, \cdots \}$, the annihilation (creation) operators $\hat{b}_i$ ($\hat{b}^\dag_i$) at site {\it{i}} are defined as
%
\begin{eqnarray}
\hat{b}_i \,|n_i\rangle &=& \sqrt{n_i} \,|n_i-1\rangle \\
\hat{b}^\dag_i \,|n_i\rangle &=& \sqrt{n_i+1} \,|n_i+1\rangle \\
\hat{b}^\dag_i \hat{b}_i \,|n_i\rangle &=&  n_i \,|n_i\rangle
\end{eqnarray}
%
where they satisfy the following commutation relations
%
\begin{equation}
 [\hat{b}_i,\hat{b}_j]=0 \,\,\,\,\, , \,\,\,\,\, [\hat{b}^\dag_i,\hat{b}^\dag_j]=0  \,\,\,\,\, , \,\,\,\,\, [\hat{b}_i,\hat{b}^\dag_j]=\delta_{ij}\,\,.
\end{equation}
\\ \\
\noindent For 2-component (spin $\sigma=\uparrow,\downarrow$) fermions, \cite{quantum_mechanics}
%
\begin{equation}
\hat{\psi}_\sigma({\bf{x}}) = \sum_i w_\sigma ({\bf{x}} - {\bf{x}}_i) \hat{c}_{i\sigma}
\end{equation}
%
where the annihilation (creation) operator $\hat{c}_{i\sigma}$ ($\hat{c}^\dag_{i\sigma}$) of a spin-$\sigma$ fermion at site $i$ is defined by the following anticommutation relations
%
\begin{equation}
\{\hat{c}_{i\sigma} , \hat{c}_{j\sigma'} \} = 0  \,\,\,\,\,  , \,\,\,\,\,
\{\hat{c}^\dag_{i\sigma} , \hat{c}^\dag_{j\sigma'} \} = 0 \,\,\,\,\, ,\,\,\,\,\,
\{\hat{c}_{i\sigma} , \hat{c}^\dag_{j\sigma'} \} = \delta_{ij} \delta_{\sigma\sigma'}   \,\, .
\end{equation}
%
Here, these operators act on the fermionic Fock basis $\{|n_{i\sigma}\rangle \, | \,  n_{i\sigma} = 0, 1 \}$, where the ordering becomes very important.
For example, the following state of a 4-site lattice is defined as
%
\begin{equation}
| \uparrow \downarrow, 0 , \downarrow, \uparrow \rangle = \hat{c}^\dag_{1\uparrow}  \hat{c}^\dag_{1\downarrow}  \hat{c}^\dag_{3\downarrow}  \hat{c}^\dag_{4\uparrow}  | 0 \rangle   \,\,\,\,\,\,\,\, (\,\, = - \hat{c}^\dag_{1\uparrow}  \hat{c}^\dag_{1\downarrow}   \hat{c}^\dag_{4\uparrow} \hat{c}^\dag_{3\downarrow}  | 0 \rangle \,\,)
\end{equation}
%
by convention.
The first 2 anticommutation relations ensure that
%
\begin{eqnarray}
\hat{c}_{i\sigma} \hat{c}_{i\sigma} &=& 0  \\
\hat{c}^\dag_{i\sigma} \hat{c}^\dag_{i\sigma} &=& 0 
\end{eqnarray}
%
or no 2 identical fermions with the same state can be annihilated or created respectively.
This is in fact the Pauli exclusion principle.
\\ \\
\noindent Two identical particles in an optical lattice interact through short-range potential, hereby modelled by contact interaction in real space, i.e.
%
\begin{equation}
U({\bf{x}}, {\bf{x}}') = \frac{4\pi a_s \hbar^2}{m} \delta({\bf{x}}, {\bf{x}}')  \,\, ,
\end{equation}
%
where the scattering mechanism is effectively described by the s-wave scattering length $a_s$, which can be easily tuneable in experiments via Feshbach resonance technique.  \cite{greiner1,chin1} 
Next, we define interaction strength among site {\it{i}}, {\it{j}}, {\it{k}}, {\it{l}} as
%
\begin{equation}
U_{ijkl} = \frac{4\pi a_s \hbar^2}{m}  \int d{\bf{x}} \, w({\bf{x}} - {\bf{x}}_i) \, w({\bf{x}} - {\bf{x}}_j)  \, w({\bf{x}} - {\bf{x}}_l) \, w({\bf{x}} - {\bf{x}}_k) \,\, .
\end{equation}
%
Often, only the onsite interaction strength is relevant to the physics, and is simply denoted by U for the case of homogeneous lattice.
\\ \\
\noindent Finally, the many-particle hamiltonian of the optical lattice system is
%
\begin{equation}
\label{many_body_hamiltonian_1}
\hat{H} 
 = \int d{\bf{x}} \, \hat{\psi}^\dag ({\bf{x}}) \left( -\frac{\hbar^2}{2m} \nabla^2 + V({\bf{x}}) \right) \hat{\psi} ({\bf{x}}) \,\,
 +\,\, \frac{1}{2} \int d{\bf{x}} \, d{\bf{x}}' \hat{\psi}^\dag ({\bf{x}}) \hat{\psi}^\dag ({\bf{x}}') \, U({\bf{x}}, {\bf{x}}') \, \hat{\psi} ({\bf{x}}') \hat{\psi} ({\bf{x}})
\end{equation}
%

\subsection{Bosons in an optical lattice}

For single component bosons in an optical lattice, Hamiltonian \ref{many_body_hamiltonian_1} can be rewritten as
%
\begin{equation}
\label{many_body_hamiltonian_2}
\hat{H} = - \sum_{ij} t_{ij} \, \hat{b}^\dag_i \, \hat{b}_j \,\, + \,\, \frac{1}{2} \sum_{ijkl} U_{ijkl} \,  \hat{b}^\dag_i \, \hat{b}^\dag_j \, \hat{b}_l \, \hat{b}_k
\end{equation}
%
which can be simplified to the {\it{boson Hubbard Hamiltonian}} \cite{fisher1} for homogeneous lattice in the limit of large $V_0$, i.e.
%
\begin{equation}
\label{boson_hubbard_model}
\hat{H} - \mu \hat{N} = -t \sum_{\langle i,j \rangle} \hat{b}^\dag_i  \hat{b}_j + \frac{U}{2} \sum_i n_i(n_i-1) - \mu \sum_i n_i  
\end{equation}
% 
where the summation $\langle i,j \rangle$  extends over all nearest neighbouring lattice sites, and the chemical potential $\mu$ determines the total number of bosons in the system.
%
\begin{figure}[h]
\centering
\includegraphics[width=7cm]{./figures/boson-hubbard-model-phase-diagram.jpg}
\caption{Phases of boson Hubbard model for 3D lattices. \cite{trotsky} At lower temperatures $T/t$, an homogeneous bosonic optical lattice exhibits superfluid (Mott insulator) phase in the limit of small (large) interaction $U/t$. At higher temperatures $T/t$, the system becomes normal fluid.}
\label{boson_hubbard_model_phase_diagram}
\end{figure}
%
\\
\noindent The boson Hubbard model \ref{boson_hubbard_model} has been completely solved at least numerically \cite{trotsky}, with its phase diagram illustrated in figure \ref{boson_hubbard_model_phase_diagram} for different interactions $U/t$ and temperatures $T/t$. 
The finite temperature properties of the system are summarized in the partition function
%
\begin{equation}
Z = \mathrm{Tr} \exp \left[-\beta(\hat{H} - \mu \hat{N}) \right]
\end{equation} 
%
where $\beta = 1/k_B T$. \cite{statistical_mechanics}
The density
%
\begin{equation}
\langle n \rangle = \frac{1}{\beta} \frac{\partial}{\partial \mu} \log Z = \frac{1}{Z} \mathrm{Tr} \,  \hat{n} \exp( -\beta \hat{H} )
\end{equation}
%
of the bosonic system is a physical quantity measurable in experiments, up to single-site resolution. \cite{chin2, bloch1}
%
\begin{figure}[h]
\centering
\includegraphics[width=7cm]{./figures/honeycomb-lattice-towards-integer-filling.jpg}
\caption{Average density $\langle n \rangle$ vs temperature $T/t$ for a $20^2$ honeycomb lattice at interaction $U/t=8$. At constant chemical potential $\mu/t=2.5$ (red) and $\mu/t=4.5$ (yellow), the average density $\langle n \rangle$ converges to an integer value of 1 in the limit of lower temperatures $T/t$. When this happens, the system is a Mott insulator. These values are computed from numerical exact Quantum Monte Carlo methods.}
\label{mott_insulator_illustration}
\end{figure}
%
Figure \ref{mott_insulator_illustration} studies how $\langle n \rangle$ varies with $T/t$, thereby illustrating a Mott-insulating phase, as evident in the integer-valued $\langle n \rangle$ \cite{fisher1}, in the limit towards lower temperatures.
\\ \\
\noindent By looking at the lattice momentum density 
%
\begin{equation}
\langle n_{\bf{k}} \rangle = \langle \hat{b}^\dag_{\bf{k}} \hat{b}_{\bf{k}} \rangle = \frac{1}{N} \sum_{ij} \langle \hat{b}_i^\dag \hat{b}_j \rangle e^{i{\bf{k}}\cdot ({\bf{x}}_i - {\bf{x}}_j )} \,\, ,
\end{equation}
where $N$ is the total number of bosons (see figure \ref{momentum_density}), we get a hint of the existence of a superfluid phase characterised by a large condensate fraction $\langle n_{{\bf{k}} = 0} \rangle/N$ for 3D optical lattices.
%
\begin{figure}[h]
\centering
\includegraphics[width=7cm]{./figures/square-lattice-momentum-density.jpg}
\caption{Cross-sectional lattice momentum density profile $\langle n_{\bf{k}} \rangle $ of a $20^2$ square lattice at interaction $U/t=8$, average density $\langle n \rangle = 1$ and temperatures $T/t=0.5$ (blue), $T/t=2.5$ (red). These values are computed from numerical exact Quantum Monte Carlo methods.}
\label{momentum_density}
\end{figure}
%
\\ \\
\noindent In continuum, the momentum density defined by
%
\begin{equation}
\langle n({\bf{k}}) \rangle = \langle \hat{\psi}^\dag ({\bf{k}})  \hat{\psi} ({\bf{k}})  \rangle = \left| w({\bf{k}}) \right|^2 \langle n_{\bf{k}} \rangle 
\end{equation}
%
with $w({\bf{k}}) = \int d{\bf{x}} \, w({\bf{x}}) e^{-i{\bf{k}}\cdot{\bf{x}}}$, can be experimentally measured by time-of-flight (TOF) images
%
\begin{equation}
n_{\mathrm{TOF}}(k_x,k_y) = \int dk_z \left| w({\bf{k}}) \right|^2 \frac{1}{N} \sum_{ij} \langle \hat{b}_i^\dag \hat{b}_j \rangle e^{i{\bf{k}}\cdot ({\bf{x}}_i - {\bf{x}}_j )}  e^{-\frac{im}{2\hbar t_{\mathrm{TOF}}}(x_i^2 - x_j^2)} \,\, ,
\end{equation}
%
along the direction-of-sight $k_z$.
\\ \\
\noindent In fact, the first piece of indirect evidence of a superfluid-insulator phase transition was derived from the TOF images obtained from the experiment by M. Greiner {\it{et al}} in 2002. \cite{greiner2} 
(See figure \ref{greiner_nature}.) 
%
\begin{figure}[h]
\centering
\includegraphics[width=10cm]{./figures/greiner-nature.jpg}
\caption{Superfluid-insulator phase transition in a gas of ultracold bosonic atoms. \cite{greiner2} Time-of-flight (TOF) images represented by optical density (OD) of up to 200,000 Rb-87 bosons trapped in a 3D optical lattice. As $V_0$ increases from (a) 0 $E_r$, (b) 3 $E_r$, (c) 7 $E_r$, (d) 10 $E_r$, (e) 13 $E_r$, (f) 14 $E_r$, (g) 16 $E_r$, to (h) 20 $E_r$, a superfluid-insulator phase transition has been indirectly observed.} 
\label{greiner_nature}
\end{figure}
%
7 years later, the superfluid-insulator phase transition phenomenon had been numerically confirmed by the boson Hubbard model. \cite{trotsky} 
(See figure \ref{experiment-comparison-time-of-flight}.)
%
\begin{figure}
\centering
\includegraphics[width=10cm]{./figures/experiment-comparison-time-of-flight.jpg}
\caption{Comparison between theory and experiment. Time-of-flight (TOF) images represented by optical density (OD) of 280,000 Rb-87 bosons trapped in a 3D optical lattice, as obtained from (a) experiments with intensity $V_0 = 8E_R$ at various temperatures, as well as from (b) ab-initio Quantum Monte Carlo simulations with interaction $U/t = 8.11$ (from band structure calculations) at best-fitted temperatures. As the temperature increases, the system transits from the superfluid phase (indicated with sharp Bragg peaks) to the normal-fluid phase (indicated with blurred TOF imaging. An excellent agreement has been observed.\cite{trotsky}}
\label{experiment-comparison-time-of-flight}
\end{figure}
%


\subsection{Fermions in an optical lattice}

For 2 component fermions in an optical lattice, Hamiltonian \ref{many_body_hamiltonian_1} can be rewritten as
%
\begin{equation}
\hat{H} = -\sum_{ij\sigma}  t_{ij} \, \hat{c}^\dag_{i\sigma} \hat{c}_{j\sigma} + \frac{1}{2} \sum_{ijkl\sigma\sigma'} U_{ijkl} \, \hat{c}^\dag_{i\sigma} \hat{c}^\dag_{j\sigma'} \hat{c}_{l\sigma'} \hat{c}_{k\sigma}
\end{equation}
%
which can be simplified to the {\it{Fermi Hubbard Hamiltonian}} \cite{hubbard1} for homogeneous lattice in the limit of large $V_0$, i.e.
%
\begin{equation}
\label{fermi_hubbard_model}
\hat{H} = -t \sum_{\langle i, j \rangle , \sigma} \hat{c}^\dag_{i\sigma} \hat{c}_{j\sigma} + U \sum_i n_{i\uparrow} n_{i\downarrow}
\end{equation}
%
where the summation $\langle i,j \rangle$  extends over all nearest neighbouring lattice sites, and $\sigma=\uparrow,\downarrow$.
Unlike the case for bosons, the Fermi Hubbard model \ref{fermi_hubbard_model}, being the simplest correlation model for fermions, remains unsolved in general even till today.
Embarassingly, only with a particle-filling of half, i.e. $\langle n \rangle = 1$, is the model numerically solved, with the finite temperature phase diagram illustrated in figure \ref{fermi_hubbard_model_phase_diagram}. \cite{muramatsu1}
%
\begin{figure}[h]
\centering
\includegraphics[width=7cm]{./figures/phase-diagram-half-filled-fhm.jpg}
\caption{Phases of Fermi Hubbard model with repulsive interactions $U/t$ for 3D lattices at half-filling, i.e. $\langle n \rangle = 1$. \cite{muramatsu1} At lower temperatures $T/t$, an homogeneous half-filled Fermi Hubbard model exhibits antiferromagnetic phase. At higher temperatures $T/t$, the system becomes a Fermi liquid (Mott insulator) at small (large) interactions $U/t$.}
\label{fermi_hubbard_model_phase_diagram}
\end{figure}
%
The phase boundary is determined theoretically by the Neel temperature $T_N$ which onsets antiferromagnetism in the model system.
The local magnetization $\langle m_i \rangle$ at site $i$ is defined as
%
\begin{equation}
\langle m_i \rangle = \frac{\langle n_{i\uparrow} \rangle - \langle n_{i\downarrow} \rangle }{2}
\end{equation}
%
where 
%
\begin{equation}
\langle n_{i\sigma} \rangle = \frac{\mathrm{Tr} \, \hat{n}_{i\sigma}  \exp (-\beta \hat{H})}{\mathrm{Tr}  \exp (-\beta \hat{H})} \,\, .
\end{equation}
%
At low temperatures $T/t$, the model system is antiferromagnetic such that
%
\begin{equation}
\langle m_i \rangle = -\langle m_j \rangle = m
\end{equation}
%
for any 2 adjacent sites $i$ and $j$, and $m$ being the sublattice magnetization. 
In practice, $T_N$ is determined via finite-size scaling analysis of the specific heat $C$ in the magnetic phase transition.  
(See figure \ref{specific_heat_finite_size_scaling}.)
%
\begin{figure}[h]
\centering
\includegraphics[width=7cm]{./figures/specific-heat-finite-size-scaling.jpg}
\caption{Specific heat $C$ as a function of temperature $T$ (in units of $t$) at interaction $U=6$ of a 3D Fermi Hubbard model system of various lengths L. Finite size scaling analysis determines the value of $T_N$. \cite{muramatsu1}}
\label{specific_heat_finite_size_scaling}
\end{figure}
%
\\ \\
\noindent State-of-the-art experiments for fermions in an optical lattice are still pretty hot in the Mott insulating phase, far from the antiferromagnetism phase. 
Tough challenges still remain technically in cooling the fermionic system down beyond Neel transition, due to Pauli exclusion principle.
Besides, potential 3-body recombination losses to form molecules destablizes antiferromagnetism, therefore adding more challenges to existing experiments aiming at studying magnetism.
Hopefully one day after the challenges are overcome, the Fermi Hubbard model can then be finally solved by a fermionic optical lattice emulator.  



\newpage

\section{Bosons in an optical lattice trap}


\subsection{Parabolic trapping}
\label{gaussian_beam_waists}

\begin{figure}[h]
\centering
\includegraphics[width=6cm]{./figures/gaussian-waist-optical-potential.jpg}
\caption{Schematic illustration of the Gaussian confinement of bosons.\cite{greiner1} Top: A pair of counter-propagating laser beams interfere in the $x$-direction to form an enveloped standing wave of intensity $V_0$, with trapping envelope of waist $w_0$. Bottom: Near the centre, the optical potential is almost uniformly periodic, i.e. $V_0 \sin^2(kx)$,  with an additional parabolic trapping term $\frac{V_0}{w_0^2}x^2$.}
\label{gaussian_waist}
\end{figure}

\noindent Bosons in an optical lattice are confined in a trapping potential $V_T({\bf{x}})$, otherwise they would fly apart.
In the Greiner experiment, the confinement has been realised with a tight Gaussian laser focus\cite{greiner1}, as schematically illustrated by figure \ref{gaussian_waist}.
Here, 3 pairs of counter-propagating gaussian laser beams set up a 3D optical lattice which is radially confined by the trapping envelope
%
\begin{equation}
V_0 \exp \left( -\frac{2r^2}{w_0^2} \right)
\end{equation}
%
of waist (or $1/e^2$-radius) $w_0$. \cite{greiner1} 
In current experiments, bosons are trapped in the vicinity close to the center of the gaussian trap, therefore effectively confined by the first order parabolic trapping term.
Taking into account also all other parabolic trapping, i.e. $V_T({\bf{x}}) = V_T {\bf{x}}^2$, where $V_T$ is the strength of the parabolic trapping, the  trapped boson Hubbard hamiltonian reads
%
\begin{equation}
\label{trapped_boson_hubbard_hamiltonian}
\hat{H} = -t \sum_{\langle i,j \rangle} \hat{b}^\dag_i  \hat{b}_j + \frac{U}{2} \sum_i n_i(n_i-1) - \sum_i (\mu - V_T {\bf{x}}_i^2)n_i  \,\, .
\end{equation}
%


\subsection{Local density approximation}

For a weak parabolic trap (small $V_T$), the local density approximation (LDA) readily estimates the density profile, i.e.
%
\begin{equation}
\langle n({\bf{x}}) \rangle \approx \langle n^{LDA}(\mu_i ; {\bf{x}}) \rangle
\end{equation}
%
where $\mu_i = \mu - V_T {\bf{x}}_i^2$ uniquely determines the density $\langle n^{LDA} \rangle$ as if the system is locally homogeneous.
%
\begin{figure}[h]
\centering
\includegraphics[width=15cm]{./figures/LDA-density-profiles.jpg}
\caption{LDA density profiles of a 2D square lattice with different interactions $U/t=8$ (left), $U/t=16$ (centre), $U/t=24$ (right) at temperature $T/t=1.0$ and parabolic trapping strength $V_T/t = 0.05$. The density is taken to be unity at the trap centre. With increasing interactions $U/t$, a Mott plateau emerges. These values are computed from numerical exact Quantum Monte Carlo methods.}
\label{LDA_density_profiles}
\end{figure}
%
In the validity regime of LDA, we get a first insight on how the Mott plateau emerges with increasing interactions $U/t$, as illustrated in figure \ref{LDA_density_profiles}.
\begin{figure}
\centering
\includegraphics[width=15cm]{./figures/coexistence-of-phases.jpg}
\caption{Coexistence of phases for bosons trapped in an 2D optical lattice trap with trapping strength $V_T/t=0.05$, at interaction $U/t=24$ and temperature $T/t=1$. Left: At different chemical potential $\mu/U$,  $\langle n_{{\bf{k}} = 0} \rangle /N$ takes different values, indicating the different phases of the bosonic system. With few particles, the system exhibits superfluid (SF) phase. At integer filling, the system possesses Mott character. When $\langle n \rangle=1$, the bosonic system is in the Mott insulator (MI) phase. Right: LDA density profile at its cross-section $y=0$. Around the centre, a Mott plateau emerges. At the wings, the bosons have superfluid character. See figure. These values are computed from numerical exact Quantum Monte Carlo methods.}
\label{coexistence_of_phases}
\end{figure}
%
Within the LDA approximation, optical lattice bosons exhibit a coexistence of different phases in a parabolic trap.
At low temperatures $T/t$ and large interactions $U/t$, the Mott plateau is surrounded by the superfluid or normal wings, depending on the local densities. (See figure \ref{coexistence_of_phases}.)
\\ \\

\subsection{Experiment validation}

In general, the numerical exact density profile $\langle n({\bf{x}})\rangle$ can be theoretically obtained by full-scale Quantum Monte Carlo (QMC) simulations \cite{prokofev1,pollet1}, where the QMC method will be discussed in detail in chapter \ref{chapter_QMC}.
%
\begin{figure}
\label{experiment-comparison-density-profile}
\centering
\includegraphics[width=8cm]{./figures/experiment-comparison-density-profile.jpg}
\caption{Comparison between theory and experiment. Average density profile $\langle n \rangle$ with respect to radial distance $r$ in units of lattice spacing $d$. Blue dots - optical lattice experiment\cite{chin1}: 9400 Cs-133 bosons are confined in a 2D square optical lattice with intensity $V_0=5E_R$ and trapping strength $V_T=0.01nK$. Band structure calculations give the hopping strength $t=4.18nK$ and onsite repulsion strength $U=10.79nK$.  Red line - numerical-exact Quantum Monte Carlo simulation\cite{shiang1}:  Single-band boson Hubbard model with interaction $U/t = 2.58$, temperature $T/t=5.98$ and trapping strength $V_T/t=0.00239$, thereby exhibiting normal-fluid phase. An excellent agreement has been observed.}
\end{figure}
%
Recent experimental advancement has enabled density probing with single-site resolution, \cite{chin2, bloch1} thus making direct comparison between theory and experiment possible.
Figure  \ref{experiment-comparison-density-profile} illustrates such a direct comparison over radial density profiles of 2D bosonic optical lattices, \cite{shiang1} which shows excellent agreement.
Quantum engineering of bosons in an optical lattice becomes therefore feasible at least in the best possible scenario.


\section{Optical lattice as a quantum emulator}

This section gives the motivation of fermionic optical lattice being a quantum emulator for single-band strongly correlated materials.

\subsection{Electronic properties of materials}

The static (non-relativistic) electronic properties of any material in nature is completely described by the Schr\"odinger's equation in atomic units, i.e.
%
\begin{equation}
\label{electronic_hamiltonian}
\hat{H} = - \sum_i \frac{ \nabla_i^2 }{2}  - \sum_{i,\alpha} \frac{Z_\alpha}{r_{i\alpha}} + \sum_{i<j} \frac{1}{r_{ij}}  
\end{equation} 
%
correct up to the Born-Oppenheimer's approximation. \cite{solid_state}
For the first 2 terms, the summation $i$ ($\alpha$) is over all electrons (nuclei), whereas for the last term, the summation $i< j$ is over all independent electron pairs.
Here, $Z_\alpha$ is the atomic number of nucleus $\alpha$, and $r$ is the inter-particle distance. 
In fact, the complete set of eigen-solutions for equation \ref{electronic_hamiltonian} decodes the entire information regarding the quantum electronic system.
In practice, however, it proves to be a formidable (or impossible) task, and one is forced to make reasonable approximations.
\\ \\
\noindent Separating into the 1-electron and 2-electron component respectively, i.e.
%
\begin{eqnarray}
\hat{H}_0 &=& - \sum_i \frac{ \nabla_i^2 }{2}  - \sum_{i,\alpha} \frac{Z_\alpha}{r_{i\alpha}} \\
\hat{H}_1 &=&  \sum_{i<j} \frac{1}{r_{ij}}
\end{eqnarray}
%
one can apply perturbative methods to approximate the eigen-solutions of the electronic system.  


\subsection{Effective valence electron in a crystal lattice.}

%
\begin{figure}[h]
\centering
\includegraphics[width=6cm]{./figures/pseudopotential.png}
\caption{Illustrating pseudopotential. \cite{pseudopotential1}}
\label{fig:pseudopotential}
\end{figure}
%
\noindent Neglecting electron-electron correlations $\hat{H}_1$, the free electron wavefunction $\phi({\bf{r}})$ satisfies the following equation, i.e.
%
\begin{equation}
\label{free_electron}
\left( -\frac{1}{2} \nabla^2 + V({\bf{x}})  \right) \phi({\bf{x}}) = \epsilon \, \phi({\bf{x}})
\end{equation}
%
whereby the crystal symmetry of the material ensures the lattice periodicity of
%
\begin{equation}
V({\bf{x}}) = -\sum_\alpha \frac{Z_\alpha}{| {\bf{x}} - {\bf{x}}_\alpha |}
\end{equation}
%
i.e. $V({\bf{x}}) = V({\bf{x}} + {\bf{x}}_\alpha)$ for all nuclei $\alpha$. \cite{solid_state}
\\ \\
\noindent In practice, only the valence electrons are relevant to the physics of the material, and therefore the wavefunction behaviour of the valence electron becomes unimportant in the core region ($r < r_c$).
Figure \ref{fig:pseudopotential} illustrates how the wavefunction of a valence electron gets smoothed out in the core region with the replacement of a {\it{pseudopotential}}. \cite{fulde1}
For the sake of simplicity, only alkali metals with 1 valence electron per nucleus will be considered throughout this chapter.
In the presence of lattice symmetry, equation \ref{free_electron} can be transformed into
%
\begin{equation}
\label{free_electron_bloch}
\left[ \frac{1}{2} (- i\nabla + {\bf{k}} )^2 + V({\bf{x}}) \right] u_{\bf{k}} ({\bf{x}}) = \epsilon_{\bf{k}} \, u_{\bf{k}} ({\bf{x}})
\end{equation}
%
with $V({\bf{x}})$ being the periodic pseudopotential, and ${\bf{k}}$ the reciprocal lattice vectors.
The eigensolution of equation \ref{free_electron_bloch} defines all the properties of the effective valence electron in a crystal lattice up to the pseudopotential approximation.


\subsection{Correlation of single band valence electrons. \\Fermi Hubbard model}

In the limit of large pseudopotential $V({\bf{x}})$, the ground state Bloch band $\epsilon^{(0)}_{\bf{k}}$ is well separated from other excited bands with a large band gap $\Delta$.
Therefore, it is already sufficient to consider only the ground state Bloch functions $u^{(0)}_{\bf{k}} ({\bf{x}})$, which correspond to the valence band functions up to the pseudopotential approximation.
Taking the Fourier transform
%
\begin{equation}
w({\bf{x}} - {\bf{x}}_\alpha) = \frac{1}{\sqrt{\Omega}} \sum_{\bf{k}} u^{(0)}_{\bf{k}} ({\bf{x}}) \, e^{i{\bf{k}}\cdot({\bf{x}} - {\bf{x}}_\alpha)}
\end{equation}
%
with $\Omega$ being the volume of the primitive unit cell, we construct an orthonormal Wannier basis $\{ w({\bf{x}} - {\bf{x}}_\alpha)\}$, in which the hopping matrix $t_{ij}^\sigma$ and the interaction matrix $U_{ijkl}^{\sigma\sigma'}$ are defined as
%
\begin{eqnarray}
t^\sigma_{ij} &=& - \int d{\bf{r}}  \, w_\sigma({\bf{x}} - {\bf{x}}_i) \left( -\frac{1}{2} \nabla^2 + V({\bf{x}})  \right) w_\sigma({\bf{x}} - {\bf{x}}_j)  \nonumber \\
U_{ijkl}^{\sigma\sigma'} &=&  \int d{\bf{x}} \, d{\bf{x}}'  w_\sigma({\bf{x}} - {\bf{x}}_i) \, w_{\sigma'}({\bf{x}}' - {\bf{x}}_i) \,  \frac{1}{| {\bf{x}} - {\bf{x}}' |} \,  w_{\sigma'}({\bf{x}}' - {\bf{x}}_l) \, w_{\sigma}({\bf{x}} - {\bf{x}}_k) \nonumber
\end{eqnarray}
%
for spin indices $\sigma,\sigma' = \, \uparrow, \downarrow$. 
These characterize the strength of which valence electrons hop across sites, as well as their interaction among different sites in the crystal lattice.
In this limit, the wannier functions are sharply localized, and only the nearest-neighbor hopping strength $t$ and the onsite interaction U are relevant to the physics. 
\\ \\
\noindent The hamiltonian for single band valence electrons reads
%
\begin{equation}
\label{correlation}
\hat{H} = -\sum_{ij\sigma}  t_{ij} \, \hat{c}^\dag_{i\sigma} \hat{c}_{j\sigma} + \frac{1}{2} \sum_{ijkl\sigma\sigma'} U_{ijkl} \, \hat{c}^\dag_{i\sigma} \hat{c}^\dag_{j\sigma'} \hat{c}_{l\sigma'} \hat{c}_{k\sigma}
\end{equation}
%
which can be simplified to the {\it{Fermi Hubbard Hamiltonian}}, i.e.
%
\begin{equation}
\hat{H} = -t \sum_{\langle i, j \rangle , \sigma} \hat{c}^\dag_{i\sigma} \hat{c}_{j\sigma} + U \sum_i n_{i\uparrow} n_{i\downarrow}  \,\, .
\end{equation}
%
This exactly coincides with the hamiltonian for fermions in a deep optical lattice.
Therefore, the fermionic optical lattice is a natural candidate for being a quantum emulator for single band strongly correlated materials.
This has been the goal of the Optical Lattice Emulator (OLE) project for more than a decade, \cite{ole} aiming at solving the impossibly hard Fermi Hubbard hamiltonian completely, thus providing scientists more insights regarding strongly correlated materials.


