%*******************************************************************************
%*********************************** Second Chapter ****************************
%*******************************************************************************

%Title of the Second Chapter
\chapter{Quantum Optics of Quantum Gases}  

\ifpdf
    \graphicspath{{Chapter2/Figs/Raster/}{Chapter2/Figs/PDF/}{Chapter2/Figs/}}
\else
    \graphicspath{{Chapter2/Figs/Vector/}{Chapter2/Figs/}}
\fi

\section{Introduction}

In this chapter, we derive a general Hamiltonian that describes the
coupling of atoms with far-detuned optical beams
\cite{mekhov2012}. This will serve as the basis from which we explore
the system in different parameter regimes, such as nondestructive
measurement in free space or quantum measurement backaction in a
cavity. Another interesting direction for this field of research are
quantum optical lattices where the trapping potential is treated
quantum mechanically. However this is beyond the scope of this work.

\mynote{insert our paper citations here}

We consider $N$ two-level atoms in an optical lattice with $M$
sites. For simplicity we will restrict our attention to spinless
bosons, although it is straightforward to generalise to fermions,
which yields its own set of interesting quantum phenomena
\cite{atoms2015, mazzucchi2016, mazzucchi2016af}, and other spin
particles. The theory can be also be generalised to continuous
systems, but the restriction to optical lattices is convenient for a
variety of reasons. Firstly, it allows us to precisely describe a
many-body atomic state over a broad range of parameter values due to
the inherent tunability of such lattices. Furthermore, this model is
capable of describing a range of different experimental setups ranging
from a small number of sites with a large filling factor (e.g.~BECs
trapped in a double-well potential) to a an extended multi-site
lattice with a low filling factor (e.g.~a system with one atom per
site which will exhibit the Mott insulator to superfluid quantum phase
transition). 

\mynote{extra fermion citations, Piazza? Look up Gabi's AF paper.}

As we have seen in the previous section, an optical lattice can be
formed with classical light beams that form standing waves. Depending
on the detuning with respect to the atomic resonance, the nodes or
antinodes form the lattice sites in which atoms accumulate. As shown
in Fig. \ref{fig:LatticeDiagram} the trapped bosons (green) are
illuminated with a coherent probe beam (red) and scatter light into a
different mode (blue) which is then measured with a detector. The most
straightforward measurement is to simply count the number of photons
with a photodetector, but it is also possible to perform a quadrature
measurement by using a homodyne detection scheme. The experiment can
be performed in free space where light can scatter in any
direction. The atoms can also be placed inside a cavity which has the
advantage of being able to enhance light scattering in a particular
direction. Furthermore, cavities allow for the formation of a fully
quantum potential in contrast to the classical lattice trap.

\begin{figure}[htbp!]
  \centering
  \includegraphics[width=\linewidth]{LatticeDiagram}
  \caption[Experimental Setup]{Atoms (green) trapped in an optical
    lattice are illuminated by a coherent probe beam (red), $a_0$,
    with a mode function $u_0(\b{r})$ which is at an angle $\theta_0$
    to the normal to the lattice. The light scatters (blue) into the
    mode $\a_1$ in free space or into a cavity and is measured by a
    detector. Its mode function is given by $u_1(\b{r})$ and it is at
    an angle $\theta_1$ relative to the normal to the lattice. If the
    experiment is in free space light can scatter in any direction. A
    cavity on the other hand enhances scattering in one particular
    direction.}
  \label{fig:LatticeDiagram}
\end{figure}

For simplicity, we will be considering one-dimensional lattices most
of the time. However, the model itself is derived for any number of
dimensions and since none of our arguments will ever rely on
dimensionality our results straightforwardly generalise to 2- and 3-D
systems. This simplification allows us to present a much more
intuitive picture of the physical setup where we only need to concern
ourselves with a single angle for each optical mode. As shown in
Fig. \ref{fig:LatticeDiagram} the angle between the normal to the
lattice and the probe and detected beam are denoted by $\theta_0$ and
$\theta_1$ respectively. We will consider these angles to be tunable
although the same effect can be achieved by varying the wavelength of
the light modes. However, it is much more intuitive to consider
variable angles in our model as this lends itself to a simpler
geometrical representation.

\section{Derivation of the Hamiltonian}
\label{sec:derivation}

A general many-body Hamiltonian coupled to a quantized light field in
second quantized can be separated into three parts,
\begin{equation}
\label{eq:FullH}
  \H = \H_f + \H_a + \H_{fa}.
\end{equation}
The term $\H_f$ represents the optical part of the Hamiltonian,
\begin{equation}
\label{eq:Hf}
  \H_f = \sum_l \hbar \omega_l \ad_l \a_l -
  i \hbar \sum_l \left( \eta_l^* \a_l - \eta_l \ad \right).
\end{equation}
The operators $\a_l$ ($\ad$) are the annihilation (creation) operators
of light modes with frequencies $\omega_l$, wave vectors $\b{k}_l$,
and mode functions $u_l(\b{r})$, which can be pumped by coherent
fields with amplitudes $\eta_l$. The second part of the Hamiltonian,
$\H_a$, is the matter-field component given by
\begin{equation}
\label{eq:Ha}
  \H_a = \int \mathrm{d}^3 \b{r} \Psi^\dagger(\b{r}) \H_{1,a}
  \Psi(\b{r}) + \frac{2 \pi a_s \hbar^2}{m} \int \mathrm{d}^3 \b{r}
  \Psi^\dagger(\b{r}) \Psi^\dagger(\b{r}) \Psi(\b{r}) \Psi(\b{r}).
\end{equation}
Here, $\Psi(\b{r})$ ($\Psi^\dagger(\b{r})$) are the matter-field
operators that annihilate (create) an atom at position $\b{r}$, $a_s$
is the $s$-wave scattering length characterising the interatomic
interaction, and $\H_{1,a}$ is the atomic part of the single-particle
Hamiltonian $\H_1$. The final component of the total Hamiltonian is
the interaction given by 
\begin{equation}
  \label{eq:Hfa}
  \H_{fa} = \int \mathrm{d}^3 \b{r} \Psi^\dagger(\b{r}) \H_{1,fa}
  \Psi(\b{r}),
\end{equation}
where $\H_{1,fa}$ is the interaction part of the single-particle
Hamiltonian, $\H_1$.

The single-particle Hamiltonian in the rotating-wave and dipole
approximation is given by
\begin{equation}
  \H_1 = \H_f + \H_{1,a} + \H_{1,fa},
\end{equation}
\begin{equation}
  \H_{1,a} = \frac{\b{p}^2} {2 m_a} + \frac{\hbar \omega_a}{2} \sigma_z,
\end{equation}
\begin{equation}
  \H_{1,fa} =  - i \hbar \sum_l \left[ \sigma^+ g_l \a_l u_l(\b{r}) - \sigma^- g^*_l
    \ad_l u^*_l(\b{r}) \right].
\end{equation}
In the equations above, $\b{p}$ and $\b{r}$ are the momentum and
position operators of an atom of mass $m_a$ and resonance frequency
$\omega_a$. The operators $\sigma^+ = |g \rangle \langle e|$,
$\sigma^- = |e \rangle \langle g|$, and
$\sigma_z = |e \rangle \langle e| - |g \rangle \langle g|$ are the
atomic raising, lowering and population difference operators, where
$|g \rangle$ and $| e \rangle$ denote the ground and excited states of
the two-level atom respectively. $g_l$ are the atom-light coupling
constants for each mode. It is the inclusion of the interaction of the
boson with quantized light that distinguishes our work from the
typical approach to ultracold atoms where all the optical fields,
including the trapping potentials, are treated classically.

We will now simplify the single-particle Hamiltonian by adiabatically
eliminating the upper excited level of the atom. The equations of
motion for the time evolution of operator $\hat{A}$ in the Heisenberg
picture are given by
\begin{equation}
  \dot{\hat{A}} = \frac{i}{\hbar} \left[\H, \hat{A} \right].
\end{equation}
Therefore, the Heisenberg equation for the lowering operator of a
single particle is
\begin{equation}
  \dot{\sigma}^- = \frac{i}{\hbar} \left[\H_1, \hat{\sigma}^- \right]
  = \hbar \omega_a \sigma^- + i \hbar \sum_l \sigma_z g_l \a_l u_l(\b{r}).
\end{equation}
We will consider nonresonant interactions between light and atoms
where the detunings between the light fields and the atomic resonance,
$\Delta_{la} = \omega_l - \omega_a$, are much larger than the
spontaneous emission rate and Rabi frequencies $g_l \a_l$. Therefore,
the atom will be predominantly found in the ground state and we can
set $\sigma_z = -1$ which is also known as the linear dipole
approximation as the dipoles respond linearly to the light amplitude
when the excited state has negligible population. Moreover, we can
adiabatically eliminate the polarization $\sigma^-$. Firstly we will
re-write its equation of motion in a frame rotating at $\omega_p$, the
external probe frequency, such that
$\sigma^- = \tilde{\sigma}^- \exp(i \omega_p t)$, and similarly for
$\tilde{\a}_l$. The resulting equation is given by
\begin{equation}
  \dot{\tilde{\sigma}}^- = - \hbar \Delta_a \tilde{\sigma}^- - i \hbar
  \sum_l g_l \tilde{\a}_l u_l(\b{r}),
\end{equation}
where $\Delta_a = \omega_p - \omega_a$ is the atom-probe
detuning. Within this rotating frame we will take
$\dot{\tilde{\sigma}}^- \approx 0$ and thus obtain the following
equation for the lowering operator
\begin{equation}
  \label{eq:sigmam}
  \sigma^- = - \frac{i}{\Delta_a} \sum_l g_l \a_l u_l(\b{r}).
\end{equation}
Therefore, by inserting this expression into the Heisenberg equation
for the light mode $m$ given by
\begin{equation}
  \dot{\a}_m = - \sigma^- g_m u^*_m(\b{r})
\end{equation}
we get the following equation of motion
\begin{equation}
  \dot{\a}_m = \frac{i}{\Delta_a} \sum_l g_l g_m u_l(\b{r})
  u^*_m(\b{r}) \a_l.
\end{equation}
An effective Hamiltonian which results in the same optical equations
of motion can be written as
$\H^\mathrm{eff}_1 = \H_f + \H^\mathrm{eff}_{1,a} +
\H^\mathrm{eff}_{1,fa}$. The effective atomic and interaction
Hamiltonians  are
\begin{equation}
\label{eq:aeff}
  \H^\mathrm{eff}_{1,a} = \frac{\b{p}^2}{2 m_a} + V_\mathrm{cl}(\b{r}),
\end{equation}
\begin{equation}
\label{eq:faeff}
  \H^\mathrm{eff}_{1,fa} = \frac{\hbar}{\Delta_a} \sum_{l,m}
  u_l^*(\b{r}) u_m(\b{r}) g_l g_m \ad_l \a_m,
\end{equation}
where we have explicitly extracted
$V_\mathrm{cl}(\b{r}) = \hbar g_\mathrm{cl}^2 |a_\mathrm{cl}
u_\mathrm{cl}(\b{r})|^2 / \Delta_{\mathrm{cl},a}$, the classical
trapping potential, from the interaction terms. However, we consider
the trapping beam to be sufficiently detuned from the other light
modes that we can neglect any scattering between them. A later
inclusion of this scattered light would not be difficult due to the
linearity of the dipoles we assumed.

Normally we will consider scattering of modes $a_l$ much weaker than
the field forming the lattice potential
$V_\mathrm{cl}(\b{r})$. Therefore, we assume that the trapping is
entirely due to the classical potential and the field operators
$\Psi(\b{r})$ can be expanded using localised Wannier functions of
$V_\mathrm{cl}(\b{r})$. By keeping only the lowest vibrational state
we get
\begin{equation}
  \Psi(\b{r}) = \sum_i^M b_i w(\b{r} - \b{r}_i),
\end{equation} 
where $b_i$ ($\bd_i$) is the annihilation (creation) operator of an
atom at site $i$ with coordinate $\b{r}_i$, $w(\b{r})$ is the lowest
band Wannier function, and $M$ is the number of lattice
sites. Substituting this expression in Eq. \eqref{eq:FullH} with
$\H_{1,a} = \H_{1,a}^\mathrm{eff}$ and
$\H_{1,fa} = \H_{1,fa}^\mathrm{eff}$ given by Eq. \eqref{eq:aeff} and
\eqref{eq:faeff} respectively yields the following generalised
Bose-Hubbard Hamiltonian, $\H = \H_f + \H_a + \H_{fa}$,
\begin{equation}
  \H = \H_f + \sum_{i,j}^M J^\mathrm{cl}_{i,j} \bd_i b_j + 
  \sum_{i,j,k,l}^M \frac{U_{ijkl}}{2} \bd_i \bd_j b_k b_l + 
  \frac{\hbar}{\Delta_a} \sum_{l,m} g_l g_m \ad_l \a_m 
  \left( \sum_{i,j}^K J^{l,m}_{i,j} \bd_i b_j \right).
\end{equation}
The optical field part of the Hamiltonian, $\H_f$, has remained
unaffected by all our approximations and is given by
Eq. \eqref{eq:Hf}. The matter-field Hamiltonian is now given by the
well known Bose-Hubbard Hamiltonian
\begin{equation}
  \H_a = \sum_{i,j}^M J^\mathrm{cl}_{i,j} \bd_i b_j + 
  \sum_{i,j,k,l}^M \frac{U_{ijkl}}{2} \bd_i \bd_j b_k b_l,
\end{equation}
where the first term represents atoms tunnelling between sites with a
hopping rate given by
\begin{equation}
  J^\mathrm{cl}_{i,j} = \int \mathrm{d}^3 \b{r} w (\b{r} - \b{r}_i ) 
  \left( -\frac{\b{p}^2}{2 m_a} + V_\mathrm{cl}(\b{r}) \right) w(\b{r}
  - \b{r}_i),
\end{equation}
and $U_{ijkl}$ is the atomic interaction term given by
\begin{equation}
  U_{ijkl} = \frac{4 \pi a_s \hbar^2}{m_a} \int \mathrm{d}^3 \b{r} 
  w(\b{r} - \b{r}_i) w(\b{r} - \b{r}_j) w(\b{r} - \b{r}_k) w(\b{r} - \b{r}_l).
\end{equation}
Both integrals above depend on the overlap between Wannier functions
corresponding to different lattice sites. This overlap decreases
rapidly as the distance between sites increases. Therefore, in both
cases we can simply neglect all terms but the one that corresponds to
the most significant overlap. Thus, for $J_{i,j}^\mathrm{cl}$ we will
only consider $i$ and $j$ that correspond to nearest neighbours.
Furthermore, since we will only be looking at lattices that have the
same separtion between all its nearest neighbours (e.g. cubic or
square lattice) we can define $J_{i,j}^\mathrm{cl} = - J^\mathrm{cl}$
(negative sign, because this way $J^\mathrm{cl} > 0$). For the
inter-atomic interactions this simplifies to simply considering
on-site collisions where $i=j=k=l$ and we define $U_{iiii} =
U$. Finally, we end up with the canonical form for the Bose-Hubbard
Hamiltonian
\begin{equation}
  \H_a = -J^\mathrm{cl} \sum_{\langle i,j \rangle}^M \bd_i b_j + 
  \frac{U}{2} \sum_{i}^M \hat{n}_i (\hat{n}_i - 1),
\end{equation}
where $\langle i,j \rangle$ denotes a summation over nearest
neighbours and $\hat{n}_i = \bd_i b_i$ is the atom number operator at
site $i$.

Finally, we have the light-matter interaction part of the Hamiltonian
given by
\begin{equation}
  \H_{fa} = \frac{\hbar}{\Delta_a} \sum_{l,m} g_l g_m \ad_l \a_m
  \left( \sum_{i,j}^K J^{l,m}_{i,j} \bd_i b_j \right),
\end{equation}
where the coupling between the matter and optical fields is determined
by the coefficients
\begin{equation}
  J^{l,m}_{i,j} = \int \mathrm{d}^3 \b{r} w(\b{r} - \b{r}_i)
  u_l^*(\b{r}) u_m(\b{r}) w(\b{r} - \b{r}_j).
\end{equation}
This contribution can be separated into two parts, one which couples
directly to the on-site atomic density and one that couples to the
tunnelling operators. We will define the operator
\begin{equation}
  \label{eq:F}
  \hat{F}_{l,m} = \hat{D}_{l,m} + \hat{B}_{l,m},
\end{equation}
where $\hat{D}_{l,m}$ is the direct coupling to atomic density
\begin{equation}
  \label{eq:D}
  \hat{D}_{l,m} = \sum_{i}^K J^{l,m}_{i,i} \hat{n}_i,
\end{equation}
and $\hat{B}_{l,m}$ couples to the matter-field via the
nearest-neighbour tunnelling operators
\begin{equation}
  \label{eq:B}
  \hat{B}_{l,m} = \sum_{\langle i, j \rangle}^K J^{l,m}_{i,j} \bd_i b_j,
\end{equation}
where $K$ denotes a sum over the illuminated sites and we neglect
couplings beyond nearest neighbours for the same reason as before when
deriving the matter Hamiltonian.

\mynote{make sure all group papers are cited here} These equations
encapsulate the simplicity and flexibility of the measurement scheme
that we are proposing. The operators given above are entirely
determined by the values of the $J^{l,m}_{i,j}$ coefficients and
despite its simplicity, this is sufficient to give rise to a host of
interesting phenomena via measurement back-action such as the
generation of multipartite entangled spatial modes in an optical
lattice \cite{elliott2015, atoms2015, mekhov2009pra}, the appearance
of long-range correlated tunnelling capable of entangling distant
lattice sites, and in the case of fermions, the break-up and
protection of strongly interacting pairs \cite{mazzucchi2016,
  kozlowski2016zeno}. Additionally, these coefficients are easy to
manipulate experimentally by adjusting the optical geometry via the
light mode functions $u_l(\b{r})$.

It is important to note that we are considering a situation where the
contribution of quantized light is much weaker than that of the
classical trapping potential. If that was not the case, it would be
necessary to determine the Wannier functions in a self-consistent way
which takes into account the depth of the quantum poterntial generated
by the quantized light modes. This significantly complicates the
treatment, but can lead to interesting physics. Amongst other things,
the atomic tunnelling and interaction coefficients will now depend on
the quantum state of light.  \mynote{cite Santiago's papers and
  Maschler/Igor EPJD}

Therefore, combining these final simplifications we finally arrive at
our quantum light-matter Hamiltonian
\begin{equation}
  \label{eq:fullH}
  \H = \H_f -J^\mathrm{cl} \sum_{\langle i,j \rangle}^M \bd_i b_j + 
  \frac{U}{2} \sum_{i}^M \hat{n}_i (\hat{n}_i - 1) + 
  \frac{\hbar}{\Delta_a} \sum_{l,m} g_l g_m \ad_l \a_m \hat{F}_{l,m} -
  i \sum_l \kappa_l \ad_l \a_l,
\end{equation}
where we have phenomologically included the cavity decay rates
$\kappa_l$ of the modes $a_l$. A crucial observation about the
structure of this Hamiltonian is that in the interaction term, the
light modes $a_l$ couple to the matter in a global way. Instead of
considering individual coupling to each site, the optical field
couples to the global state of the atoms within the illuminated region
via the operator $\hat{F}_{l,m}$. This will have important
implications for the system and is one of the leading factors
responsible for many-body behaviour beyond the Bose-Hubbard
Hamiltonian paradigm. 

Furthermore, it is also vital to note that light couples to the matter
via an operator, namely $\hat{F}_{l,m}$, which makes it sensitive to
the quantum state of matter. This is a key element of our treatment of
the ultimate quantum regime of light-matter interaction that goes
beyond previous treatments.

\section{The Bose-Hubbard Hamiltonian}

\subsection{The Model}

The original Hubbard model was developed as a model for the motion of
electrons within transition metals \cite{hubbard1963}. The key
elements were the confinement of the motion to a lattice and the
approximation of the Coulomb screening interaction as a simple on-site
interaction. Despite these enormous simplifications the model was very
succesful \cite{leggett}. The Bose-Hubbard model is an even simpler
variation where instead of fermions we consider spinless bosons. It
was originally devised as a toy model, but in 1998 it was shown by
Jaksch \emph{et.~al.} that it can be realised with ultracold atoms in
an optical lattice \cite{jaksch1998}. Shortly afterwards it was
obtained in a ground-breaking experiment \cite{greiner2002}. The model
has been the subject of intense research since then, because despite
its simplicity it possesses highly nontrivial properties such as the
superfluid to Mott insulator quantum phase transition. Furthermore, it
is one of the most controllable quantum many-body systems thus
providing a solid basis for new experiments and technologies.

The model we have derived is essentially an extension of the
well-known Bose-Hubbard model that also includes interactions with
quantised light. Therefore, it should come as no surprise that if we
eliminate all the quantized fields from the Hamiltonian we obtain
exactly the Bose-Hubbard model
\begin{equation}
  \H_a = -J \sum_{\langle i,j \rangle}^M \bd_i b_j + 
  \frac{U}{2} \sum_{i}^M \hat{n}_i (\hat{n}_i - 1).
\end{equation}
The first term is the kinetic energy term and denotes the hopping of
atoms between neighbouring sites at a hopping rate given by $J$. The
second term describes the on-site repulsion which acts mainly to
suppress multiple occupancies of any site, but at the same time it
lies at the heart of strongly correlated phenomena in the Bose-Hubbard
model. Unfortunately, despite its simplicity, this system is not
solvable for any finite value $U/J$ using known techniques in finite
dimensions \cite{krauth1991, kolovsky2004, calzetta2006}. However,
exact solutions exist for the ground state in the two limits of
$U = 0$ and $J = 0$ which already provide some insight as they
correspond to the two different phases of the quantum phase transition
that is present in the model. In higher dimensions (e.g.~two or three)
the system is reasonably well described by a mean-field approximation
which provides us with a good understanding of the full quantum phase
diagram.

\subsection{The Superfluid and Mott Insulator Ground States}

Before trying to understand the quantum phase transition itself we
will look at the ground states of the two extremal cases first and for
simplicity we will consider only homogenous systems with periodic
boundary conditions here. We will begin with the $U = 0$ limit where
the Hamiltonian's kinetic term can be diagonalised by transforming to
momentum space defined by the annihilation operator
\begin{equation}
  b_\b{k} = \frac{1} {\sqrt{M}} \sum_m b_m e^{i \b{k} \cdot \b{r}_m},
\end{equation}
where $\b{k}$ denotes the wavevector running over the first Brillouin
zone. The Hamiltonian then is given by
\begin{equation}
  \H_a = \sum_\b{k} \epsilon_\b{k} \bd_\b{k} b_\b{k},
\end{equation}
where the free particle dispersion is
\begin{equation}
  \epsilon_\b{k} = -2 J \sum_{d = 1}^D \cos(k_d a),
\end{equation}
and $d$ denotes the dimension out of a total of $D$ dimensions. It
should now be obvious that the ground state is simply the state with
all the atoms in the $\b{k} = 0$ state and thus it is given by
\begin{equation}
  \label{eq:GSSF}
  | \Psi_\mathrm{SF} \rangle = \frac{ ( \bd_{\b{k} = 0} )^N } {\sqrt{N!}} | 0
  \rangle = \frac{1} {\sqrt{N!}} \left( \frac{1} {\sqrt{M}} \sum_m^M
    \bd_m \right)^N | 0 \rangle,
\end{equation}
where $| 0 \rangle$ denotes the vacuum state. This state is called a
superfluid (SF) due to its viscous-free flow properties
\cite{leggett1999}. The macroscopic occupancy of the single-particle
ground state is the signature and defining property of a
Bose-Einstein condensate. Note that the zero momentum states consist of a
superposition of an atom at all the sites of the optical lattice
highlighting the fact that kinetic energy is minimised by delocalising
the particles. This in turn implies that a superfluid state will have
infinite range correlations. It is also important to note that in the
thermodynamic limit the superfluid state is gapless, i.e.~it takes
zero energy to excite the system.

At the other end of the spectrum, i.e.~$J = 0$, the kinetic term
vanishes and the system is dominated by the on-site interactions. This
has the effect of decoupling the Hamiltonian into a sum of identical
single-site Hamiltonians which are already diagonal in the atom number
basis. Therefore, the ground state for $N = gM$, where $g$ is an
integer, is given by
\begin{equation}
  \label{eq:GSMI}
  | \Psi_\mathrm{MI} \rangle = \prod_m^M \frac{1} {\sqrt{g!}}
  (\bd_m)^g | 0 \rangle.
\end{equation}
This state has precisely the same number of bosons, $g$, at each site
and is commonly known as the Mott insulator (MI). Unlike in the
superfluid state, the atoms are all localised to a specific lattice
site and thus not only are there no long-range correlations, there are
actually no two-site correlations in this state at all. Additionally,
this state is gapped, i.e.~it costs a finite amount of energy to
excite the system into its lowest excited state even in the
thermodynamic limit. In fact, the existance of the gap is the defining
property that distinguishes the Mott insulator from the superfluid in
the Bose-Hubbard model. Note that we have only considered a
commensurate case (number of atoms an integer multiple of the number
of sites). If we were to insert another atom on top of the Mott
insulator, the energy of the system would not depend on the lattice
site in which it is inserted. This implies that as soon as $J$ becomes
finite, the ground state would prefer this atom be delocalised and
effectively becomes gapless which means the system is a superfluid and
it will not exhibit a quantum phase transition.

\subsection{Mean-Field Theory}

Now that we have a basic understanding of the two limiting cases we
can now consider the model in between these two extremes. We have
already mentioned that an exact solution is not known, but fotunately
a very good mean-field approximation exists \cite{fisher1989}. In this
approach the interaction term is treated exactly, but the kinetic
energy term is decoupled as
\begin{equation}
  \bd_m b_n = \langle \bd_m \rangle b_n + \bd_m \langle b_n \rangle -
  \langle \bd_m \rangle \langle b_n \rangle = \Phi^* b_n + \Phi \bd_m
  - | \Phi |^2,
\end{equation}
where the expectation values are for the $T = 0$ ground state. Note
that we have assumed that $\Phi = \langle b_m \rangle $ is non-zero
and homogenous. $\Phi$ describes the influence of hopping between
neighbouring sites and is the mean -field order parameter. This
decoupling effect means that we can write the Hamiltonian as $\hat{H}_a
= \sum_m^M \hat{h}_a$, where
\begin{equation}
  \hat{h}_a = -z J ( \Phi \bd \Phi^* b) + z J | \Phi |^2 + \frac{U}{2}
  \n (\n - 1) - \mu \n,
\end{equation}
and $z$ is the coordination number, i.e. the number of nearest
neighbours of a single site. We have also introduced the chemical
potential $\mu$ since we now consider the grand canonical ensemble as
the Hamiltonian no longer conserves the total atom number. The ground
state of the $| \Psi_0 \rangle$ of the overall system will be a
site-wise product of the individual ground states of
$\hat{h}_a$. These can be found very easily using standard
diagonalisation techniques. $\Phi$ is then self-consistently
determined by minimising the energy of the ground state with respect
to $\Phi$.

The main advantage of the mean-field treatment is that it lets us
study the quantum phase transition between the sueprfluid and Mott
insulator phases discussed in the previous sections. The phase
boundaries can be obtained from second-order perturbation theory as 
\begin{equation}
  \left( \frac{\mu} {U} \right)_\pm = \frac{1}{2} \left[ 2g - 1 -
    \left( \frac{zJ}{U} \right) \pm \sqrt{ 1 - 2 (2g + 1) \left( \frac{zJ}{U}
        \right) + \left( \frac{zJ}{U} \right)^2} \right].
\end{equation}
This yields a phase diagram with multiple Mott insulating lobes as
seen in Fig. \ref{fig:BHPhase} where each lobe represents a Mott insualtor with a
different number of atoms per site. The tip of the lobe which
corresponds to the quantum phase transition along a line of fixed
density is obtained by equating the two branches of the above equation
to get
\begin{equation}
  \left( \frac{J} {U} \right) = \frac{1} {z (2g + 1 + \sqrt{4g^2 + 4g} )}.
\end{equation}
It is important to note that the fact that this formalism predicts a
phase transition is in contrast to other mean-field theories such as
the Bogoliubov approximation \cite{PitaevskiiStringari}. This
highlights the fact that the interactions, which are treated exactly
here, are the dominant driver leading to this phase transition and
other strongly correlated effects in the Bose-Hubbard model.

\begin{figure}
  \centering
  \includegraphics[width=0.8\linewidth]{BHPhase}
  \caption[Mean-Field Bose-Hubbard Phase Diagram]{Mean-field phase
    diagram of the Bose-Hubbard model in 1D, i.e.~ $z = 2$, from
    Ref. \cite{StephenThesis}. The shaded regions are the Mott
    insulator lobes and each lobe corresponds to a different on-site
    filling labeled by $n$. The rest of the space corresponds to the
    superfluid phase. The dashed lines are the phase boundaries
    obtained from first-order perturbation theory. The solid lines are
    are lines of constant density in the superfluid
    phase. \label{fig:BHPhase}}
\end{figure}

\subsection{The Bose-Hubbard Model in One Dimension}

The mean-field theory in the previous section is very useful tool for
studying the quantum phase transition in the Bose-Hubbard
model. However, it is effectively an infinite-dimensional theory and
in practice it only works in two dimensions or more. The phase
transition in 1D is poorly described, because it actually belongs to a
different universality class. This is clearly seen from the one
dimensional phase diagram shown in Fig. \ref{fig:1DPhase}.

\begin{figure}
  \centering
  \includegraphics[width=0.8\linewidth]{1DPhase}
  \caption[Exact 1D Bose-Hubbard Phase Diagram]{Exact phase diagram
    of the Bose-Hubbard model in 1D, i.e.~ $z = 2$, from
    Ref. \cite{StephenThesis}. The shaded region corresponds to the
    $n = 1$ Mott insulator lobe. The sharp tip distinguishes it from
    the mean-field phase diagram. The red dotted line denotes the
    reentrance phase transition which occurs for a fixed $\mu$.The
    critical point for the fixed density $n = 1$ transition is denoted
    by an 'x'. \label{fig:1DPhase}}
\end{figure}

Some general conclusions can be obtained by looking at Haldane's
prescription for Luttinger liquids \cite{haldane1981,
  giamarchi}. Without a periodic potential the low-energy physics of
the system is described by the Hamiltonian
\begin{equation}
  \hat{H}_a = \frac{1}{2 \pi} \int \mathrm{d} x \left\{ v K [ \hat{\Pi}(x)
    ]^2 + \frac{v} {K} [\partial_x \hat{\Phi}(x) ]^2 \right\},
\end{equation}
where we have expressed the bosonic field operators in terms of a
density operator $\hat{\rho}(x)$ and a phase operator $\hat{\Phi}(x)$
as $\hat{\Psi}(x) = \sqrt{\hat{\rho}(x)} e^{i \hat{\hat{\Phi}(x)}}$
and $\hat{\Pi}(x)$ is the density fluctuation operator. Provided the
parameters $v$ and $K$ can be correctly determined this Hamiltonian
gives the correct description of the gapless superfluid phase of the
Bose-Hubbard model. Most importantly it gives an expression for the
spatial correlation functions such as
\begin{equation}
  \langle \bd_m b_n \rangle = A \left( \frac{\alpha} {|m - n|} \right)^{K/2},
\end{equation}
where $A$ is some amplitude and $\alpha$ is a necessary cutoff to
regularise the theory at short distances. Unlike the superfluid ground
state in Eq. \eqref{eq:GSSF} this state does not have infinite range
correlations. They decay according to a power law. However, for
non-interacting systems $K = 0$ and long-range order is re-established
as before though it is important to note that in higher dimensions
this long-range order persists in the whole superfluid phase even with
interactions present. In order to describe the phase transition and
the Mott insulating phase it is necessary to introduce a periodic
lattice potential. It can be shown that this system exhibits at
$T = 0$ a Berezinskii-Kosterlitz-Thouless phase transition as the
parameter $K$ is varied with a critical point at $K = \frac{1}{2}$
where $K < \frac{1}{2}$ is a superfluid. Above $K = \frac{1}{2}$ the
value of $K$ jumps discontinuously to $K \rightarrow \infty$ producing
the Mott insulator phase. Unlike the gapless superfluid phase the
spatial correlations decay exponentially as
\begin{equation}
  \langle \bd_m b_n \rangle = B e^{ - |m - n|/\xi},
\end{equation}
where $B$ is some constant  and the correlation length is given by
$\xi = v / \Delta$ where $\Delta$ is the energy gap.

Using advanced numerical methods such as density matrix
renormalisation group (DMRG) calculations it is possible to identify
the critical point by fitting the power-law decay correlations in
order to obtain $K$. The resulting phase transition is shown in
Fig. \ref{fig:1DPhase} and the critical point was shown to be $(U/zJ) = 1.68$. Note
that unusually the phase diagram exhibits a reentrance phase
transition for a fixed $\mu$.

\section{Scattered light behaviour}
\label{sec:a}

Having derived the full quantum light-matter Hamiltonian we will now
look at the behaviour of the scattered light. We begin by looking at
the equations of motion in the Heisenberg picture
\begin{equation}
  \dot{\a_l} = \frac{i}{\hbar} [\hat{H}, \a_l] = 
  -i \left(\omega_l + U_{l,l}\hat{F}_{1,1} \right) \a_l
   - i \sum_{m \ne l} U_{l,m} \hat{F}_{l,m} \a_m
    - \kappa_l \a_l + \eta_l,
\end{equation}
where we have defined $U_{l,m} = g_l g_m / \Delta_a$. By considering
the steady state solution in the frame rotating with the probe
frequency $\omega_p$ we can show that the stationary solution is given
by
\begin{equation}
  \a_l = \frac{\sum_{m \ne l} U_{l,m} \a_m \hat{F}_{l,m} + i \eta_l} 
  {(\Delta_{lp} - U_{l,l} \hat{F}_{l,l}) + i \kappa_l},
\end{equation}
where $\Delta_{lp} = \omega_p - \omega_l$ is the detuning between the
probe beam and the light mode $a_l$. This equation gives us a direct
relationship between the light operators and the atomic operators. In
the stationary limit, the quantum state of the light field
adiabatically follows the quantum state of matter.

The above equation is quite general as it includes an arbitrary number
of light modes which can be pumped directly into the cavity or
produced via scattering from other modes. To simplify the equation
slightly we will neglect the cavity resonancy shift,
$U_{l,l} \hat{F}_{l,l}$ which is possible provided the cavity decay
rate and/or probe detuning are large enough. We will also only
consider probing with an external coherent beam, $a_0$ with mode
function $u_0(\b{r})$, and thus we neglect any cavity pumping
$\eta_l$. We also limit ourselves to only a single scattered mode,
$a_1$ with a mode function $u_1(\b{r})$. This leads to a simple linear
relationship between the light mode and the atomic operator
$\hat{F}_{1,0}$
\begin{equation}
  \label{eq:a}
  \a_1 = \frac{U_{1,0} a_0} {\Delta_{p} + i \kappa} \hat{F} =
  C \hat{F},
\end{equation}
where we have defined $C = U_{1,0} a_0 / (\Delta_{p} + i \kappa)$
which is essentially the Rayleigh scattering coefficient into the
cavity. Furthermore, since there is no longer any ambiguity in the
indices $l$ and $m$, we have dropped indices on $\Delta_{1p} \equiv
\Delta_p$, $\kappa_1 \equiv \kappa$, and $\hat{F}_{1,0} \equiv
\hat{F}$. We also do the same for the operators $\hat{D}_{1,0} \equiv
\hat{D}$, $\hat{B}_{1,0} \equiv \hat{B}$, and the coefficients
$J^{1,0}_{i,j} \equiv J_{i,j}$.

The operator $\a_1$ itself is not an observable. However, it is
possible to combine the outgoing light field with a stronger local
oscillator beam in order to measure the light quadrature
\begin{equation}
  \hat{X}_\phi = \frac{1}{2} \left( \a_1 e^{-i \phi} + \ad_1 e^{i \phi} \right),
\end{equation}
which in turn can expressed via the quadrature of $\hat{F}$,
$\hat{X}^F_\beta$, as
\begin{equation}
  \hat{X}_\phi = |C| \hat{X}_\beta^F = \frac{|C|}{2} \left( \hat{F}
    e^{-i \beta} + \hat{F}^\dagger e^{i \beta} \right),
\end{equation}
where $\beta = \phi - \phi_C$, $C = |C| \exp(i \phi_C)$, and $\phi$ is
the local oscillator phase.

Whilst the light amplitude and the quadratures are only linear in
atomic operators, we can easily have access to higher moments via
related quantities such as the photon number
$\ad_1 \a_1 = |C|^2 \hat{F}^\dagger \hat{F}$ or the quadrature
variance
\begin{equation}
  \label{eq:Xvar}
  ( \Delta X_\phi )^2 = \langle \hat{X}_\phi^2 \rangle - \langle
  \hat{X}_\phi \rangle^2 = \frac{1}{4} + |C|^2 (\Delta X^F_\beta)^2,
\end{equation}
which reflect atomic correlations and fluctuations.

Finally, even though we only consider a single scattered mode, this
model can be applied to free space by simply varying the direction of
the scattered light mode if multiple scattering events can be
neglected. This is likely to be the case since the interactions will
be dominated by photons scattering from the much larger coherent
probe.

\section{Density and Phase Observables}
\label{sec:B}

Light scatters due to its interactions with the dipole moment of the
atoms which for off-resonant light, like the type that we consider,
results in an effective coupling with atomic density, not the
matter-wave amplitude. Therefore, it is challenging to couple light to
the phase of the matter-field as is typical in quantum optics for
optical fields. Most of the exisiting work on measurement couples
directly to atomic density operators \cite{mekhov2012, LP2009,
  rogers2014, ashida2015, ashida2015a}. However, it has been shown
that one can couple to the interference term between two condensates
(e.g.~a BEC in a double-well) by using interference measurements
\cite{cirac1996, castin1997, ruostekoski1997, ruostekoski1998,
  rist2012}. Such measurements establish a relative phase between the
condensates even though the two components have initially well-defined
atom numbers which is phase's conjugate variable. In a lattice
geometry, one would ideally measure between two sites similarly to
single-site addressing \cite{greiner2009, bloch2011}, which would
measure a single term $\langle \bd_i b_{i+1}+b_i
\bd_{i+1}\rangle$. This could be achieved, for example, by superposing
a deeper optical lattice shifted by $d/2$ with respect to the original
one, catching and measuring the atoms in the new lattice
sites. However, a single-shot success rate of atom detection will be
small and as single-site addressing is challenging, we proceed with
our global scattering scheme.

In our model light couples to the operator $\hat{F}$ which consists of
a density component, $\hat{D} = \sum_i J_{i,i} \hat{n}_i$, and a phase
component, $\hat{B} = \sum_{\langle i, j \rangle} J_{i,j} \bd_i
b_j$. In general, the density component dominates,
$\hat{D} \gg \hat{B}$, and thus $\hat{F} \approx \hat{D}$.  However,
it is possible to engineer an optical geometry in which $\hat{D} = 0$
leaving $\hat{B}$ as the dominant term in $\hat{F}$. This approach is
fundamentally different from the aforementioned double-well proposals
as it directly couples to the interference terms caused by atoms
tunnelling rather than combining light scattered from different
sources. Furthermore, it is not limited to a double-wellsetup and
naturally extends to a lattice structure which is a key
advantage. Such a counter-intuitive configuration may affect works on
quantum gases trapped in quantum potentials \cite{mekhov2012,
  mekhov2008, larson2008, chen2009, habibian2013, ivanov2014,
  caballero2015} and quantum measurement-induced preparation of
many-body atomic states \cite{mazzucchi2016, mekhov2009prl,
  pedersen2014, elliott2015}.

\mynote{add citiations above if necessary}

For clarity we will consider a 1D lattice as shown in
Fig. \ref{fig:LatticeDiagram} with lattice spacing $d$ along the
$x$-axis direction, but the results can be applied and generalised to
higher dimensions. Central to engineering the $\hat{F}$ operator are
the coefficients $J_{i,j}$ given by
\begin{equation}
  \label{eq:Jcoeff}
  J_{i,j} = \int \mathrm{d} x \,\,\, w(x - x_i) u_1^*(x) u_0(x) w(x - x_j).
\end{equation}
The operators $\hat{B}$ and $\hat{D}$ depend on the values of
$J_{i,i+1}$ and $J_{i,i}$ respectively. These coefficients are
determined by the convolution of the coupling strength between the
probe and scattered light modes, $u_1^*(x)u_0(x)$, with the relevant
Wannier function overlap shown in Fig. \ref{fig:WannierProducts}a. For
the $\hat{B}$ operator we calculate the convolution with the nearest
neighbour overlap, $W_1(x) \equiv w(x - d/2) w(x + d/2)$, and for the
$\hat{D}$ operator we calculate the convolution with the square of the
Wannier function at a single site, $W_0(x) \equiv w^2(x)$. Therefore,
in order to enhance the $\hat{B}$ term we need to maximise the overlap
between the light modes and the nearest neighbour Wannier overlap,
$W_1(x)$. This can be achieved by concentrating the light between the
sites rather than at the positions of the atoms. 

\mynote{Potentially expand details of the derivation of these
  equations} 

In order to calculate the $J_{i,j}$ coefficients we perform numerical
calculations using realistic Wannier functions
\cite{walters2013}. However, it is possible to gain some analytic
insight into the behaviour of these values by looking at the Fourier
transforms of the Wannier function overlaps,
$\mathcal{F}[W_{0,1}](k)$, shown in Fig.
\ref{fig:WannierProducts}b. This is because for plane and standing
wave light modes the product $u_1^*(x) u_0(x)$ can be in general
decomposed into a sum of oscillating exponentials of the form
$e^{i k x}$ making the integral in Eq. \eqref{eq:Jcoeff} a sum of
Fourier transforms of $W_{0,1}(x)$. We consider both the detected and
probe beam to be standing waves which gives the following expressions
for the $\hat{D}$ and $\hat{B}$ operators
\begin{align}
  \label{eq:FTs}
  \hat{D} = & \frac{1}{2}[\mathcal{F}[W_0](k_-)\sum_m\hat{n}_m\cos(k_-
              x_m +\varphi_-) \nonumber\\ 
            & + \mathcal{F}[W_0](k_+)\sum_m\hat{n}_m\cos(k_+ x_m +\varphi_+)],
              \nonumber\\ 
  \hat{B} = & \frac{1}{2}[\mathcal{F}[W_1](k_-)\sum_m\hat{B}_m\cos(k_- x_m
              +\frac{k_-d}{2}+\varphi_-) \nonumber\\ 
            & +\mathcal{F}[W_1](k_+)\sum_m\hat{B}_m\cos(k_+
              x_m +\frac{k_+d}{2}+\varphi_+)],
\end{align}
where $k_\pm = k_{0x} \pm k_{1x}$,
$k_{(0,1)x} = k_{0,1} \sin(\theta_{0,1}$),
$\hat{B}_m=b^\dag_mb_{m+1}+b_mb^\dag_{m+1}$, and
$\varphi_\pm=\varphi_0 \pm \varphi_1$. The key result is that the
$\hat{B}$ operator is phase shifted by $k_\pm d/2$ with respect to the
$\hat{D}$ operator since it depends on the amplitude of light in
between the lattice sites and not at the positions of the atoms
allowing to decouple them at specific angles.

\begin{figure}[hbtp!]
  \centering
  \includegraphics[width=0.8\linewidth]{BDiagram}
  \caption[Maximising Light-Matter Coupling between Lattice
  Sites]{Light field arrangements which maximise coupling, $u_1^*u_0$,
    between lattice sites. The thin black line indicates the trapping
    potential (not to scale). (a) Arrangement for the uniform pattern
    $J_{i,i+1} = J_1$. (b) Arrangement for spatially varying pattern
    $J_{i,i+1}=(-1)^m J_2$; here $u_0=1$ so it is not shown and $u_1$
    is real thus $u_1^*u_0=u_1$. \label{fig:BDiagram}}
\end{figure}

\begin{figure}[hbtp!]
  \centering
  \includegraphics[width=\linewidth]{WF_S}
  \caption[Wannier Function Products]{The Wannier function products:
    (a) $W_0(x)$ (solid line, right axis), $W_1(x)$ (dashed line, left
    axis) and their (b) Fourier transforms $\mathcal{F}[W_{0,1}]$. The
    Density $J_{i,i}$ and matter-interference $J_{i,i+1}$ coefficients
    in diffraction maximum (c) and minimum (d) as are shown as
    functions of standing wave shifts $\varphi$ or, if one were to
    measure the quadrature variance $(\Delta X^F_\beta)^2$, the local
    oscillator phase $\beta$. The black points indicate the positions,
    where light measures matter interference $\hat{B} \ne 0$, and the
    density-term is suppressed, $\hat{D} = 0$. The trapping potential
    depth is approximately 5 recoil energies.}
  \label{fig:WannierProducts}
\end{figure}

The simplest case is to find a diffraction maximum where
$J_{i,i+1} = J^B_\mathrm{max}$, where $J^B_\mathrm{max}$ is a
constant. This results in a diffraction maximum where each bond
(inter-site term) scatters light in phase and the operator is given by
\begin{equation}
  \hat{B} = J^B_\mathrm{max} \sum_m^K \hat{B}_m .
\end{equation}
This can be achieved by crossing the light modes such that
$\theta_0 = -\theta_1$ and $k_{0x} = k_{1x} = \pi/d$ and choosing the
light mode phases such that $\varphi_+ = 0$. Fig. \ref{fig:BDiagram}a
shows the resulting light mode functions and their product along the
lattice and Fig. \ref{fig:WannierProducts}c shows the value of the
$J_{i,j}$ coefficients under these circumstances. In order to make the
$\hat{B}$ contribution to light scattering dominant we need to set
$\hat{D} = 0$ which from Eq. \eqref{eq:FTs} we see is possible if
\begin{equation}
  \xi \equiv \varphi_0 = -\varphi_1 =
  \frac{1}{2}\arccos\left[\frac{-\mathcal{F}[W_0](2\pi/d)}{\mathcal{F}[W_0](0)}\right].
\end{equation}
Under these conditions, the coefficient $J^B_\mathrm{max}$ is simply
given by $J^B_\mathrm{max} = \mathcal{F}[W_1](2 \pi / d)$. This
arrangement of light modes maximizes the interference signal,
$\hat{B}$, by suppressing the density signal, $\hat{D}$, via
interference compensating for the spreading of the Wannier functions.

Another possibility is to obtain an alternating pattern similar
corresponding to a diffraction minimum where each bond scatters light
in anti-phase with its neighbours giving
\begin{equation}
  \hat{B} = J^B_\mathrm{min} \sum_m^K (-1)^m \hat{B}_m,
\end{equation}
where $J^B_\mathrm{min}$ is a constant. We consider an arrangement
where the beams are arranged such that $k_{0x} = 0$ and
$k_{1x} = \pi/d$ which gives the following expressions for the density
and interference terms
\begin{align}
  \label{eq:DMin}
  \hat{D} = & \mathcal{F}[W_0]\left(\frac{\pi}{d}\right) \sum_m (-1)^m \hat{n}_m
              \cos(\varphi_0) \cos(\varphi_1) \nonumber \\ 
  \hat{B} = & -\mathcal{F}[W_1]\left(\frac{\pi}{d}\right) \sum_m (-1)^m \hat{B}_m
              \cos(\varphi_0) \sin(\varphi_1).
\end{align}
For $\varphi_0 = 0$ the corresponding $J_{i,j}$ coefficients are given
by $J_{i,i+1} = (-1)^i J^B_\mathrm{min}$, where
$J^B_\mathrm{min} = -\mathcal{F}[W_1](\pi / d)$, and are shown in
Fig. \ref{fig:WannierProducts}d. The light mode coupling along the
lattice is shown in Fig. \ref{fig:BDiagram}b. It is clear that for
$\varphi_1 = \pm \pi/2$, $\hat{D} = 0$, which is intuitive as this
places the lattice sites at the nodes of the mode $u_1(x)$. This is a
diffraction minimum as the light amplitude is also zero,
$\langle \hat{B} \rangle = 0$, because contributions from alternating
inter-site regions interfere destructively. However, the intensity
$\langle \ad_1 \a \rangle = |C|^2 \langle \hat{B}^2 \rangle$ is
proportional to the variance of $\hat{B}$ and is non-zero.

Alternatively, one can use the arrangement for a diffraction minimum
described above, but use travelling instead of standing waves for the
probe and detected beams and measure the light quadrature variance.
In this case
$\hat{X}^F_\beta = \hat{D} \cos(\beta) + \hat{B} \sin(\beta)$ and by
varying the local oscillator phase, one can choose which conjugate
operator to measure.

\section{Electric Field Stength}
\label{sec:Efield}

The Electric field operator at position $\b{r}$ and at time $t$ is
usually written in terms of its positive and negative components:
\begin{equation}
  \b{\hat{E}}(\b{r},t) = \b{\hat{E}}^{(+)}(\b{r},t) + \b{\hat{E}}^{(-)}(\b{r},t),
\end{equation}
where
\begin{subequations}
  \begin{equation}
    \b{\hat{E}}^{(+)}(\b{r},t) = \sum_{\b{k}} \hat{\epsilon}_{\b{k}} 
    \mathcal{E}_{\b{k}} \a_{\b{k}} e^{-i \omega_{\b{k}} t + i \b{k} \cdot \b{r}},
  \end{equation}
  \begin{equation}
    \b{\hat{E}}^{(-)}(\b{r},t) = \sum_{\b{k}} \hat{\epsilon}_{\b{k}} 
    \mathcal{E}_{\b{k}} \a^\dagger_{\b{k}} e^{i \omega_{\b{k}} t - i \b{k} \cdot \b{r}},
  \end{equation}
\end{subequations}
$\hat{\epsilon}_{\b{k}}$ is a unit polarization vector, $\b{k}$ is the
wave vector,
$\mathcal{E}_{\b{k}} = \sqrt{\hbar \omega_{\b{k}} / 2 \epsilon_0 V}$,
$\epsilon_0$ is the free space permittivity, $V$ is the quantization
volume and $\a_\b{k}$ and $\a_\b{k}^\dagger$ are the annihilation and
creation operators respectively of a photon in mode $\b{k}$, and
$\omega_\b{k}$ is the angular frequency of mode $\b{k}$.

In Ref. \cite{Scully} it was shown that the operator
$\b{\hat{E}}^{(+)}(\b{r},t)$ due to light scattering from a single
atom located at $\b{r}^\prime$ at the observation point $\b{r}$ is
given by
\begin{equation}
  \label{eq:Ep}
  \b{\hat{E}}^{(+)}(\b{r},\b{r}^\prime,t) = \frac{\omega_a^2 d_A \sin \eta}{4 \pi
    \epsilon_0 c^2 |\b{r} - \b{r}^\prime|} \hat{\epsilon} \sigma^-
  \left( \b{r}^\prime, t - \frac{|\b{r} - \b{r}^\prime|}{c} \right),
\end{equation}
with a similar expression for
$\b{E}^{(-)}(\b{r},\b{r}^\prime,t)$. Eq. \eqref{eq:Ep} is valid only
in the far field, $\eta$ is the angle the dipole makes with
$\b{r} - \b{r}^\prime$, $\hat{\epsilon}$ is the polarization vector
which is perpendicular to $\b{r} - \b{r}^\prime$ and lies in the plane
defined by $\b{r} - \b{r}^\prime$ and the dipole, $\omega_a$ is the
atomic transition frequency, and $d_A$ is the dipole matrix element
between the two levels, and $c$ is the speed of light in vacuum.

We have already derived an expression for the atomic lowering
operator, $\sigma^-$, in Eq. \eqref{eq:sigmam} and it is given by
\begin{equation}
  \sigma^-(\b{r}, t) = - \frac{i} {\Delta_a} \sum_\b{k} g_\b{k}
  a_\b{k}(t) u_\b{k}(\b{r}).
\end{equation}
We will want to substitute this expression into Eq. \eqref{eq:Ep}, but
first we note that in the setup we consider the coherent probe will be
much stronger than the scattered modes. This in turn implies that the
expression for $\sigma^-$ will be dominated by this external
probe. Therefore, we drop the sum and consider only a single wave
vector, $\b{k}_0$, of the dominant contribution from the external
beam corresponding to the mode $\a_0$.

We now use the definititon of the atom-light coupling constant
\cite{Scully} to make the substitution
$g_0 a_0 = -i \Omega_0 e^{-i \omega_0 t} / 2$, where
$\Omega_0$ is the Rabi frequency for the probe-atom system. We now
evaluate the polarisation operator at the retarded time
\begin{equation}
  \sigma^- \left(\b{r}^\prime, t - \frac{|\b{r} - \b{r}^\prime|}{c} \right) =
  -\frac{\Omega_0}{2 \Delta_a} u_0 (\b{r}^\prime) \exp \left[-i \omega_0
    \left( t - \frac{|\b{r} - \b{r}^\prime|}{c} \right)\right].
\end{equation}
We are considering observation in the far field regime so
$|\b{r} - \b{r}^\prime| \approx r - \hat{\b{r}} \cdot
\b{r}^\prime$. We also consider the incoming, $\b{k}_0$, and the
outgoing, $\mathbf{k}_1$, waves to be of the same wavelength,
$|\b{k}_0| = |\b{k}_1| = k = \omega_0 / c$, hence
$(\omega_0 / c) |\b{r} - \b{r}^\prime| \approx k (r - \hat{\b{r}}
\cdot \b{r}^\prime) = \b{k}_1 \cdot (\b{r} - \b{r}^\prime )$ and
\begin{equation}
  \sigma^- \left(\b{r}^\prime, t - \frac{|\b{r} - \b{r}^\prime|}{c} \right) \approx
  -\frac{\Omega_0}{2 \Delta_a} u_0 (\b{r}^\prime) \exp \left[i
    \b{k}_1 \cdot (\b{r} - \b{r}^\prime ) - i \omega_0 t \right].
\end{equation}

The many-body version of the electric field operator is given by
\begin{equation}
  \b{\hat{E}}^{(+)}_N(\b{r},t) = \int \mathrm{d}^3 \b{r}^\prime
  \Psi^\dagger (\b{r}^\prime) \b{\hat{E}}^{(+)}(\b{r},\b{r}^\prime,t) \Psi(\b{r}^\prime),
\end{equation}
where $\Psi(\b{r})$ is the atomic matter-field operator. As before, we
expand this field operator using localized Wannier functions
corresponding to the lattice potential and keeping only the lowest
vibrational state at each site:
$\Psi(\b{r}) = \sum_i b_i w(\b{r} - \b{r}_i)$, where $b_i$ is the
annihilation operator of an atom at site $i$ with the coordinate
$\b{r}_i$. Thus, the relevant many-body operator is
\begin{subequations}
  \begin{equation}
    \b{\hat{E}}^{(+)}_N(\b{r},t) = \hat{\epsilon} C_E
    \sum_{i,j}^K \bd_i b_j \int \mathrm{d}^3 \b{r}^\prime
    w^* (\b{r}^\prime - \b{r}_i)  \frac{u_0
      (\b{r}^\prime)}{|\b{r} - \b{r}^\prime|} e^ {i \b{k}_1
      \cdot (\b{r} - \b{r}^\prime ) - i \omega_0 t }
    w(\b{r}^\prime - \b{r}_j),
  \end{equation}
  \begin{equation}
    C_E = -\frac{\omega_a^2 \Omega_0 d \sin \eta}{8 \pi \epsilon_0 c^2
      \Delta_a}
  \end{equation}
\end{subequations}
where $K$ is the number of illuminated sites. We now assume that the
lattice is deep and that the light scattering occurs on time scales
much shorter than atom tunneling. Therefore, we ignore all terms for
which $i \ne j$ and the remaining integrals become $\int \mathrm{d}^3
\b{r}_0 w^*(\b{r}_0 - \b{r}_i) f (\b{r})
w(\b{r}_0 - \b{r}_i) = f(\b{r}_i)$. The final form of
the many body operator is then
\begin{equation}
  \b{E}^{(+)}_N(\b{r},t) = \hat{\epsilon} C_E
  \sum_{j = 1}^K \hat{n}_j \frac{u_0 (\b{r}_j)}{|\b{r} -
    \b{r}_j|} e^ {i \b{k}_1 \cdot (\b{r} - \b{r}_j
    ) - i \omega_0 t },
\end{equation}
where $\hat{n}_i = b^\dagger_ib_i$ is the atom number operator at site
$i$. We will now consider the incoming probe to be a plane wave,
$u_0(\b{r}) = e^{i \b{k}_0 \cdot \b{r}}$, and we
approximate $|\b{r} - \b{r}_j| \approx r$ in the
denominator. Thus,
\begin{equation}
  \b{E}^{(+)}_N(\b{r},t) = \hat{\epsilon} C_E \frac{e^ {ikr - i\omega_0t}}{r}
  \sum_{j = 1}^K \hat{n}_j e^ {i (\b{k}_0 - \b{k}_1) \cdot \b{r}_j}.
\end{equation}
We note that this is exactly the same form of the optical field as
derived from a generalised Bose-Hubbard model which considers a fully
quantum light-matter interaction Hamiltonian \cite{mekhov2007prl,
  mekhov2007pra, mekhov2012}. We can even express the result in terms of
the $\hat{D} = \sum_j^K u_1^*(\b{r}_j) u_0(\b{r}_j)
\hat{n}_j$ formalism used in those works
\begin{equation}
  \b{E}^{(+)}_N(\b{r},t) = \hat{\epsilon} C_E \frac{e^
    {ikr - i\omega_0t}}{r} \hat{D}.
\end{equation}

The average light intensity at point $\b{r}$ at time $t$ is given
by the formula $I = c \epsilon_0 |E|^2/2$ and yields
\begin{equation}
  \langle I (\b{r}, t) \rangle = c \epsilon_0 \langle
  E^{(-)} (\b{r}, t) E^{(+)} (\b{r}, t) \rangle =
  \frac{c \epsilon_0}{r^2} C_E^2 \langle \hat{D}^* \hat{D} \rangle.
\end{equation}
The scattering is dominated by a uniform background for which $\langle
\hat{D}^* \hat{D} \rangle \approx N_K$ for a superfluid
\cite{mekhov2012}, where $N_K$ is the mean number of atoms in the
illuminated volume. Thus, the approximate number of photons scattered
per second can be obtained by integrating over a sphere
\begin{equation}
  n_{\Phi} = \frac{4 \pi c \epsilon_0}{3 \hbar \omega_a} \left(\frac{\omega_a^2
      \Omega_0 d}{8 \pi \epsilon_0 c^2 \Delta_a}\right)^2 N_K.
\end{equation}
In the Weisskopf-Wigner approximation it can be shown \cite{Scully}
that for a two level system the decay constant, $\Gamma$, is
\begin{equation}
  \Gamma = \frac{\omega_a^3 d^2}{3 \pi \epsilon_0 \hbar c^3}.
\end{equation}
Therefore, we can now express the quantity $n_{\Phi}$ as
\begin{equation}
  n_{\Phi} = \frac{1}{8} \left(\frac{\Omega_0}{\Delta_a}\right)^2 \frac{\Gamma}{2} N_K.
\end{equation}

Estimates of the scattering rate using real experimental parameters
are given in Table \ref{tab:photons}. Rubidium atom data has been
taken from Ref. \cite{steck}. Miyake \emph{et al.} experimental
parameters are from Ref. \cite{miyake2011}. The $5^2S_{1/2}$,
$F=2 \rightarrow 5^2P_{3/2}$, $F^\prime = 3$ transition of $^{87}$Rb
is considered. For this transition the Rabi frequency is actually
larger than the detuning and and effects of saturation should be taken
into account in a more complete analysis. However, it is included for
discussion. The detuning was specified as ``a few linewidths'' and
note that it is much smaller than $\Omega_0$. The Rabi fequency is
$\Omega_0 = (d_\mathrm{eff}/\hbar)\sqrt{2 I / c \epsilon_0}$ which is
obtained from definition of Rabi frequency,
$\Omega = \mathbf{d} \cdot \mathbf{E} / \hbar$, assuming the electric
field is parallel to the dipole, and the relation
$I = c \epsilon_0 |E|^2 /2$. We used
$d_\mathrm{eff} = 1.73 \times 10^{-29}$ Cm \cite{steck}. Weitenberg
\emph{et al.} experimental parameters are based on
Ref. \cite{weitenberg2011, weitenbergThesis}. The
$5^2S_{1/2} \rightarrow 5^2P_{3/2}$ transition of $^{87}$Rb is
used. Ref. \cite{weitenberg2011} gives the free space detuning to be
$\Delta_\mathrm{free} = - 2 \pi \cdot 45$ MHz, but
Ref. \cite{weitenbergThesis} clarifies that the relevant detuning is
$\Delta = \Delta_\mathrm{free} + \Delta_\mathrm{lat}$, where
$\Delta_\mathrm{lat} = - 2 \pi \cdot 40$
MHz. Ref. \cite{weitenbergThesis} uses a saturation parameter
$s_\mathrm{tot}$ to quantify the intensity of the beams which is
related to the Rabi frequency,
$s_\mathrm{tot} = 2 \Omega^2 / \Gamma^2$ \cite{steck,foot}. We can
extract $s_\mathrm{tot}$ for the experiment from the total scattering
rate by
$\Gamma_\mathrm{sc} = (\Gamma/2) (s_\mathrm{tot}) /
(1+s_\mathrm{tot}+(2 \Delta / \Gamma)^2)$. A scattering rate of 60 kHz
per atom \cite{weitenberg2011} gives $s_\mathrm{tot} = 2.5$.

\begin{table}
  \centering
  \begin{tabular}{l c c}
    \toprule
    Value & Miyake \emph{et al.} & Weitenberg \emph{et al.} \\ \midrule
    $\omega_a$ & \multicolumn{2}{ c }{$2 \pi \cdot 384$ THz}\\
    $\Gamma$ & \multicolumn{2}{ c }{$2 \pi \cdot 6.07$ MHz} \\
    $\Delta_a$ & $2\pi \cdot 30$ MHz & $2 \pi \cdot 85$ MHz \\
    $I$ & $4250$ Wm$^{-2}$ & N/A \\ 
    $\Omega_0$ & 293$\times 10^6$ s$^{-1}$ & 42.5$\times 10^6$ s$^{-1}$ \\
    $N_K$ & 10$^5$ & 147 \\ \midrule
    $n_{\Phi}$ & $6 \times 10^{11}$ s$^{-1}$ & $2 \times 10^6$ s$^{-1}$ \\
    \bottomrule
  \end{tabular}
  \caption[Photon Scattering Rates]{Experimental parameters used in
    estimating the photon scattering rates.}
  \label{tab:photons}
\end{table}