dphil-thesis/Chapter2/chapter2.tex

%*******************************************************************************
%*********************************** 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


%********************************** % First Section  ****************************

\section{Ultracold Atoms in Optical Lattices}

%********************************** % Second Section  ***************************

\section{Quantum Optics of Quantum Gases}

Having introduced and described the behaviour of ultracold bosons
trapped and manipulated using classical light, it is time to extend
the discussion to quantized optical fields. We will first 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=1.0\textwidth]{LatticeDiagram}
  \caption[Experimental Setup]{Atoms (green) trapped in an optical
    lattice are illuminated by a coherent probe beam (red). The light
    scatters (blue) in free space or into a cavity and is measured by
    a detector. 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.

\subsection{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})$. We now
proceed in the same way as when deriving the conventional Bose-Hubbard
Hamiltonian in the zero temperature limit. The field operators
$\Psi(\b{r})$ can be expanded using localised Wannier functions of
$V_\mathrm{cl}(\b{r})$ and 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.

\subsection{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$, and thus we neglect any cavity
pumping $\eta_l$. We also limit ourselves to only a single scattered
mode, $a_1$. 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.

\subsection{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 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. 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. Ideally, one could
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. A single-shot success rate of atom detection
will be small. As single-site addressing is challenging, we proceed
with the global scattering.

\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 the 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_1$, where $J_1$ is a constant. 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[-\mathcal{F}[W_0]\left(\frac{2\pi}{d}\right)/\mathcal{F}[W_0](0)].
\end{equation}
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. 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](\pi/d) \sum_m (-1)^m \hat{n}_m
              \cos(\varphi_0) \cos(\varphi_1) \nonumber \\ 
  \hat{B} = & -\mathcal{F}[W_1](\pi/d) \sum_m (-1)^m \hat{B}_m
              \cos(\varphi_0) \sin(\varphi_1).
\end{align}
The corresponding $J_{i,j}$ coefficients are given by
$J_{i,i+1} = -(-1)^i J_2$, where $J_2$ is a constant, and are shown in
Fig. \ref{fig:WannierProducts}d for $\varphi_0=0$. 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.

\subsection{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}.

\begin{table}[!htbp]
  \centering
  \begin{tabular}{|c|c|c|}
    \hline
    Value & Miyake \emph{et al.} & Weitenberg \emph{et al.} \\ \hline
    $\omega_a$ & \multicolumn{2}{|c|}{$2 \pi \cdot 384$ THz}\\ \hline
    $\Gamma$ & \multicolumn{2}{|c|}{$2 \pi \cdot 6.07$ MHz} \\ \hline
    $\Delta_a$ & $2\pi \cdot 30$ MHz & $2 \pi \cdot 85$ MHz \\ \hline
    $I$ & $4250$ Wm$^{-2}$ & N/A \\ \hline
    $\Omega_0$ & 293$\times 10^6$ s$^{-1}$ & 42.5$\times 10^6$ s$^{-1}$ \\ \hline
    $N_K$ & 10$^5$ & 147 \\ \hline \hline
    $n_{\Phi}$ & $6 \times 10^{11}$ s$^{-1}$ & $2 \times 10^6$ s$^{-1}$ \\ \hline
  \end{tabular}
  \caption[Photon Scattering Rates]{ Rubidium data taken from
    Ref. \cite{steck}; Miyake \emph{et al.}: based on
    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 is said to be a few natural
    linewidths. 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.}:
    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. 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$.}
  \label{tab:photons}
\end{table}