qsim.tex

\documentclass[aps,prl,twocolumn,showpacs,amsmath,amssymb,superscriptaddress]{revtex4-1}

\usepackage{graphicx}
\usepackage[justification=centerlast]{caption} % FIXME: 'justified' should be default, but it does not work
\usepackage{subfig}
\usepackage{amsbsy}
\usepackage{amssymb}
\usepackage{esint}
\usepackage{bm}
\usepackage[usenames]{color}
\usepackage[shortcuts]{extdash}

\newcommand{\bogdansremark}[1]{\textcolor{red}{{[}BO: #1{]}}}
\newcommand{\petersremark}[1]{\textcolor{cyan}{{[}PD: #1{]}}}
\newcommand{\andreisremark}[1]{\textcolor{green}{{[}AS: #1{]}}}

% macros for different symbols
\newcommand{\Rb}{$^{87}$Rb}
\newcommand{\figref}[1]{Fig.~\ref{#1}}
\newcommand{\bra}[1]{\mbox{\ensuremath{\langle #1 |}}}
\newcommand{\ket}[1]{\mbox{\ensuremath{| #1 \rangle}}}
\newcommand{\xvec}{\boldsymbol{x}}
\newcommand{\ivec}{\boldsymbol{i}}
\newcommand{\jvec}{\boldsymbol{j}}
\newcommand{\kvec}{\boldsymbol{k}}
\newcommand{\Psivec}{\boldsymbol{\Psi}}

% Magic macro for creating tables in figure environment
\def\imagetop#1{\vtop{\null\hbox{#1}}}

\begin{document}

\title{Quantitative simulations in three-dimensional BEC interferometry}

\newcommand{\swinburneaffiliation}{Centre for Atom Optics and Ultrafast Spectroscopy, Swinburne University of Technology, Hawthorn, VIC 3122, Australia}
\newcommand{\uqaffiliation}{Physics Department, University of Queensland, Queensland, Australia}

\author{B.~Opanchuk}
\affiliation{\swinburneaffiliation}
\author{M.~Egorov}
\affiliation{\swinburneaffiliation}
\author{S.~Hoffmann}
\affiliation{\uqaffiliation}
\author{A.~Sidorov}
\affiliation{\swinburneaffiliation}
\author{P.~D.~Drummond}
\affiliation{\swinburneaffiliation}

\begin{abstract}
We give a first-principles comparison of a calculation of quantum
noise and squeezing in a three-dimensional Bose gas, with nonlinear losses
as compared to recent experimental measurements, using a truncated Wigner simulation.
Without using fitting parameters, we find excellent agreement
between the predicted fringe contrast and experimental data,
in an atom interferometer measurement.
\end{abstract}

\date{\today}

\maketitle

Atom interferometry is an important quantum technology,
at the heart of many suggested future applications of ultra-cold atomic physics.
Bose-Einstein condensates or atom lasers have potential advantages as interferometric detectors and sensors,
provided one can extract atomic phase information.
However, unlike photons, atoms can interact strongly, causing dephasing.
An intimate understanding of quantum many-body dynamics is the key to understanding
the precise nature of interaction-induced dephasing in the measurement process.
This is essential for a quantitative theory of atom interferometry.

In this Letter we give a simple, yet quantitatively accurate
theoretical approach to calculations of atom interferometry at large atom number,
using a truncated Wigner representation~\cite{Drummond1993,Steel1998,Sinatra2002}.
This method extends the conventional Gross-Pitaevski equations
describing a Bose condensate to include quantum noise effects, including noise due to linear and nonlinear losses.
The theory we describe here predicts the results of
interferometric measurements on Bose condensates within experimental errors~\cite{Egorov2010}.
No fitting parameters are required in these comparisons.
As a result, we suggest that the theoretical approach employed here
can be reliably used as a benchmark for predictions.
Importantly, we can clearly demonstrate where phase decay is driven by quantum fluctuations,
and where it is driven by trap inhomogeneity, in order to choose optimal conditions for quantum noise reduction and spin-squeezing.

Quantum phase-diffusion is defined as the phase noise induced by number fluctuations,
which are intrinsically conjugate to phase.
As such, this is a fundamental limitation of any BEC interferometry,
and can only be removed when there are no interactions driving phase evolution.
However, there are other causes of density fluctuations as well as the overall number uncertainty,
which are also important in causing observed phase decays.
The important feature of the approach used here is that it is able to capture
all three significant features of atom interferometry that can result in decoherence.
These are:
\begin{itemize}
\item fluctuations caused by an atomic beam-splitter,
\item linear and nonlinear losses
\item trap inhomogeneity and propagation
\end{itemize}
The results given in this paper require a controllable truncation
of an expansion in the inverse atom number.
They are applicable to simulations where the atom number per lattice point or mode is large.
Such quantitatively predictive theories have a great practical value
in assessing the utility of atom interferometry in new environments.
We expect that multi-configurational Hartree theories should give comparable or greater accuracy,
in principle, as they do not require truncation.
However, MCH methods are currently limited to lossless,
one-dimensional calculations with small atom numbers,
owing to computational complexity issues.

The theory we use is straightforward.
We start by assuming that the trapped BEC has $s$-wave interactions,
together with Markovian losses due to $n$-body collisions.
We employ the Wigner-Moyal quantum phase-space representation~\cite{Gardiner2004}
and the corresponding Lindblad-type master equation,
together with a truncation of third and higher-order derivatives in the equations of motion.
If we regard the commonly used Gross-Pitaevskii equation as a classical,
first approximation to mean-field condensate dynamics,
the truncated Wigner approach is best thought of as the second approximation
in an expansion in inverse particle number.
Such truncations have been shown to be valid in the limit of large particle number~\cite{Drummond1993,Steel1998,Sinatra2002},
provided the effects of unoccupied vacuum modes can be neglected.
They are particularly useful in low-dimensional or trap environments,
where they have been successful in predictions of quantum squeezing
and phase-diffusion effects to high accuracy in photonics~\cite{Corney2008}.
The truncation approximation has been further validated by comparison
with the exact positive-P method,
and with experiment in one-dimensional quantum soliton propagation dynamics.

In the present Letter, we treat an ultra-cold,
interacting multi-component spinor Bose gas in $D$ effective dimensions.
The basic Hamiltonian is easily expressed using quantum fields
$\widehat{\Psi}_{i}^{\dagger}(\xvec),\widehat{\Psi}_{j}(\xvec)$,
where $\widehat{\Psi}_{i}^{\dagger}(\xvec)$ creates a bosonic atom of spin $i$
at location $\xvec$, and $\widehat{\Psi}_{i}(\xvec)$ destroys one;
the commutators are
$[\widehat{\Psi}_{i}(\xvec),\widehat{\Psi}_{j}^{\dagger}(\xvec^\prime)] =
\delta^{(D)}(\xvec-\xvec^\prime)\delta_{ij}\,\,.$
The resulting physics of a dilute, low-temperature Bose gas
is well-described in the $s$-wave scattering limit by an effective Hamiltonian
with contact interactions and external potentials:
\begin{equation}
	\hat{H} / \hbar = \int d^{D}\xvec \left\{
		\widehat{\Psi}_{i}^{\dagger} K_{ij} \widehat{\Psi}_{j} +
		\frac{U_{ij}}{2} \widehat{\Psi}_{i}^{\dagger} \widehat{\Psi}_{j}^{\dagger}
		\widehat{\Psi}_{j} \widehat{\Psi}_{i}
	\right\}.
\end{equation}
Here we omit the field argument $(\xvec)$ after field operators for brevity,
and use the Einstein summation convention of summing over repeated indices.
$K_{ij}$ is the unitary evolution operator:
\begin{equation}
	K_{ij} = \left( -\frac{\hbar}{2m} \nabla^2 + \omega_i + V_i(\xvec) / \hbar \right) \delta_{ij} +
		\Omega_{ij}(t),
\end{equation}
where $V_{i}$ is the external trapping potential for spin $i$,
$\omega_{i}$ is the internal Zeeman energy splitting associated with spin $i$,
and $\Omega_{ij}$ represents a time-dependent coupling
that is used to rotate one spin projection into another.
Thus $\langle \widehat{\Psi}_{i}^{\dagger} \widehat{\Psi}_{i} \rangle$
is the spin $i$ atomic density, $m$ is the atomic mass,
and $U_{ij}$ is the atom-atom interaction term.
For a dilute gas at low enough temperatures,
$U_{ij}$ is related to the $s$-wave scattering length in three dimensions by:
\begin{equation}
	U_{ij}=\frac{4\pi\hbar a_{ij}}{m}.
\end{equation}
Here we implicitly assume a momentum cutoff $k_{c} \ll 1 / a_{ij}$,
otherwise the couplings must be renormalized~\cite{Sinatra2002}.

The truncated Wigner-Moyal quantum phase-space theory is well-suited
to these types of calculation of many-body dynamics.
Although the Wigner function is not a true probability,
its becomes one to a good approximation in the limit of large particle number,
so that a probabilistic (random) sampling of phase-space trajectories
can be used to treat quantum fluctuations.
This means that there are no complexity problems due to exponential growth
of the numbers of quantum states.
Furthermore, it is straightforward to treat a spatially inhomogeneous,
many-mode problem as found in an extended D-dimensional environment.
Initial conditions such as a coherent beam-splitter
as employed in typical atom interferometry experiments can be readily included.

This approach does require an approximation of truncating third-order derivatives
in an otherwise exact functional Fokker-Planck equation.
This approximation is valid when $N \gg M$, where $N$ is the atom number and
$M$ is the number of low-energy modes~\cite{Drummond1993,Sinatra2002,Norrie2006},
where it has been checked in one-dimensional cases.
We note that in free-space calculations in three dimensions,
the ultra-violet (UV) divergence in free-space mode density at large momentum cutoff
makes the truncated Wigner method unreliable,
unless a momentum truncation is used so that $N \gg M$.

In the relevant limits where the technique is applicable,
the equations simply reduce to Gross-Pitaevskii equations with Gaussian fluctuations
of the initial conditions and additional stochastic terms.
Damping and losses can be included via an additional Markovian master equation~\cite{Jack2002}
defined so that,
\begin{equation}
	\frac{d\hat{\rho}}{dt} = -\frac{i}{\hbar} \left[ \hat{H}, \hat{\rho} \right] +
	\sum_{n,\jvec} \kappa_{\jvec}^{(n)}
	\int d^D\xvec \mathcal{L}_{\jvec}^{(n)} \left[ \hat{\rho} \right],
\end{equation}
where $n$ is the number of interacting particles,
and $\jvec$ is a vector indicating the spins that are coupled.

Here we have introduced local Liouville loss terms,
\begin{equation}
	\mathcal{L}_{\jvec}^{(n)} \left[ \hat{\rho}\right] =
	2 \hat{O}_{\jvec}^{(n)} \hat{\rho}\hat{O}_{\jvec}^{(n)\dagger} -
	\hat{O}_{\jvec}^{(n)\dagger}\hat{O}_{\jvec}^{(n)}\hat{\rho} -
	\hat{\rho}\hat{O}_{\jvec}^{(n)\dagger}\hat{O}_{\jvec}^{(n)},
\end{equation}
where the reservoir coupling operators $\hat{O}_{\jvec}^{(n)}$
are all $n$-fold products of local field annihilation operators,
\begin{equation}
	\hat{O}_{\jvec}^{(n)} =
	\hat{O}_{\jvec}^{(n)} \left( \widehat{\Psivec} \right) =
	\widehat{\Psi}_{j_{1}} \left( \xvec \right)
	\widehat{\Psi}_{j_{2}} \left( \xvec \right) \ldots
	\widehat{\Psi}_{j_{n}} \left( \xvec \right),
\end{equation}
describing local $n$-body collision losses.
Temporarily ignoring unitary evolution, the corresponding mean amplitude equations
for an $n$-body loss of this form are
\begin{equation}
	\frac{\partial}{\partial t} \left\langle
		\widehat{\Psi}_{i} \left( \xvec \right)
	\right\rangle =
	-\left\langle \hat{\Gamma}_{i} \left( \xvec \right) \right\rangle ,
\end{equation}
where we have introduced  nonlinear loss operators:
\begin{equation}
	\hat{\Gamma}_{i} \left( \xvec \right) \equiv
	\sum_{n,\jvec} \kappa_{\jvec}^{(n)}
	\frac{\partial\hat{O}_{\jvec}^{(n)\dagger}	\left( \widehat{\Psivec} \right)}
		{\partial \hat{\Psi}_{i}^{\dagger} \left( \xvec \right)}
	\hat{O}_{\jvec}^{(n)} \left( \xvec \right) .
\end{equation}

We proceed by using a stochastic phase-space method
that allows a numerical simulation of the quantum dynamics~\cite{Drummond1993,Steel1998,Hoffmann2008}.
Defining a Wigner function $W\left(\Psi_{1},\Psi_{2}\right)$,
where $\Psi$ is a c-number field corresponding to the quantum field $\hat{\Psi}$,
this has a time-evolution equation:
\begin{equation}
\begin{split}
	\frac{\partial W}{\partial t} & = \int d^D \xvec \,\left\{
		-\frac{\delta}{\delta\Psi_{i}} A_{i} -
		\frac{\delta}{\delta\Psi_{i}^*} A_{i}^* \right. \\
& + \left. \frac{\delta^{2}}{\delta\Psi_{i} \delta\Psi_{j}^{*}}D_{ij} +
		\mbox{O} \left[ \frac{\delta^{3}}{\delta\Psi_{i}^{3}} \right]
	\right\} W,
\end{split}
\end{equation}
% TODO:
%	\petersremark{The other way to proceed is to use a normalized field approach
%	with a normalization so that $\Psi \sim \epsilon\Psi$.
%	Then we can simply use a notation that the terms dropped are $O(\epsilon^3)$.
%	Here $\epsilon = \sqrt(M/N)$, say.
%	Then we must change the later notation as well, for consistency.}
where higher derivative terms of type $\mbox{O}\left[\delta^3/\delta\Psi_{i}^{3}\right]$
will be truncated.
This approximation neglects higher-order terms in an expansion in $1/\sqrt{N}$,
and is therefore valid in the limit of $N \gg M$.
The drift term $A_{i}$ gives Gross-Pitaevskii evolution with nonlinear damping:
\begin{equation}
	A_{i} = -i \left( K_{ij} + U_{ij} \lvert \Psi_{j} \rvert^{2} \right) \Psi_{j} -
	\Gamma_{i} \left( \Psivec \right),
\end{equation}
where $\Gamma_{i}$ is a classical loss term:
\begin{equation}
	\Gamma_{i} = \sum_{n,\jvec} \kappa_{\jvec}^{(n)}
	\frac{
		\partial O_{\jvec}^{(n)*} \left( \Psivec \right)
	}{
		\partial \Psi_{i}^* \left( \xvec \right)
	}
	O_{\jvec}^{(n)} \left( \xvec \right).
\end{equation}

The next terms in the evolution equation give rise to
quantum noise associated with the loss reservoirs,
with:
\begin{equation}
	D_{ij} = \sum_{n,\kvec} \kappa_{\kvec}^{(n)}
	\frac{\partial O_{\kvec}^{(n)*}}{\partial\Psi_{i}^{*}}
	\frac{\partial O_{\kvec}^{(n)}}{\partial\Psi_{j}}.
\end{equation}

This effective Fokker-Planck equation can be transformed
to a coupled stochastic Gross-Pitaevskii equation:
\begin{equation}
	\frac{d\Psi_{j}}{dt} = -i \left( K_{ij} + U_{ij} \lvert \Psi_{i} \rvert ^{2} \right) \Psi_{j} -
	\Gamma_{j} + \sum_{n,\ivec} \gamma_{\ivec,j}^{(n)} \zeta_{\ivec}^{(n)} (\xvec,t).
\end{equation}
Here $\zeta_{\ivec}^{(n)}(\xvec, t)$ is a complex,
stochastic delta-correlated Gaussian noise with
\begin{equation}
	\left\langle
		\zeta_{\ivec}^{(n)} (\xvec,t) \zeta_{\jvec}^{(m)}(\xvec^\prime, t^\prime)
	\right\rangle =
	\delta_{\ivec \jvec} \delta^{nm} \delta^{D} \left(
		\xvec - \xvec^\prime
	\right)
	\delta \left( t - t^\prime \right).
\end{equation}
The multiplicative noise coefficient
\begin{equation}
	\gamma_{\ivec,j}^{(n)} \left( \Psivec \right) =
	\sqrt{\kappa_{\ivec}^{(n)}}
	\frac{\partial O_{\ivec}^{(n)}}{\partial\Psi_{j}}
\end{equation}
is a fluctuation-dissipation term,
so that the Wigner variables remain equivalent to the corresponding operators.

For initial conditions in interferometry it is usually sufficient
to consider a coherent state amplitude $\Psi_{s}^{c}$,
corresponding to a typical initial state with Poissonian number fluctuations.
In this case the initial Wigner amplitude has a Gaussian random distribution,
with $\Psi_{s}(\xvec,t_{0})=\Psi_{s}^{c}(\xvec)+\Delta\Psi_{s}(\xvec)$,
where:
\begin{equation}
	\left\langle
		\Delta\Psi_{s}(\xvec) \Delta\Psi_{u}^{*}(\xvec^\prime)
	\right\rangle =
	\frac{1}{2} \delta_{su} \delta^{D} \left( \xvec - \xvec^\prime\right).
\end{equation}
For greater accuracy, the coherent-state initial condition
can be modified to account for initial ground-state correlations,
or additional number-fluctuations above the Poissonian level.
In numerical simulations, the delta-function correlations mean that
each low-energy mode has a vacuum occupation of a $1/2$ of a particle.
This is because the Wigner representation generates symmetrically ordered correlation functions,
so that when normal ordered correlations are measured\---which is a typical case\---
one has to express them as a sum of symmetrically ordered terms first.
The exact equation depends on the chosen basis and its low-energy subset.
In this work we used the simple basis of plane waves in the volume $V$,
and the atomic density of component $s$ is therefore:
\begin{equation}
\begin{split}
n_{s} (\xvec) & =
\langle \widehat{\Psi}^\dagger_{s} (\xvec) \widehat{\Psi}_{s} (\xvec) \rangle \\
& = \langle \left\{ \widehat{\Psi}^\dagger_{s} (\xvec)
	\widehat{\Psi}_{s} (\xvec) \right\}_{\mathrm{sym}} \rangle - \frac{M}{2V} \\
& = \int \Psi^*_{s} \Psi_{s} W(\Psivec, \Psivec^*) d\Psivec d\Psivec^* - \frac{M}{2V} \\
& = \langle \Psi^*_{s} (\xvec) \Psi_{s} (\xvec) \rangle_{\mathrm{paths}} - \frac{M}{2V}.
\end{split}
\end{equation}
Here we used the fact that the approximate Wigner function is a probability distribution
equivalent to an averaged sum over different simulation paths.

As a specific example of the utility of this method,
we consider recent three-dimensional dipole trap interference experiments,
involving ground-states of \Rb~BEC with quantum numbers ${F=1,\, m_{F}=+1}$
and ${F=2,\, m_{F}=-1}$ respectively~\cite{Egorov2010}.

\begin{figure}
	\begin{tabular}{l l}
	% 0.381in 1.95in for 0.7 size
	\imagetop{\hspace*{0.44in}\includegraphics[width=2.6in]{figures_precreated/ramsey_sequence.pdf}} & \imagetop{(a)} \\
	\imagetop{\includegraphics{figures_generated/long_ramsey_visibility.pdf}} & \imagetop{(b)} \\
%	\vskip 0.2cm
	\imagetop{\hspace*{0.44in}\includegraphics[width=2.6in]{figures_precreated/echo_sequence.pdf}} & \imagetop{(c)} \\
	\imagetop{\includegraphics{figures_generated/long_rephasing_visibility.pdf}} & \imagetop{(d)}
	\end{tabular}

	\caption{(color online).
	Timeline of the experiment for Ramsey (a) and Ramsey with symmetric spin-echo (c),
	(b) and (d) are the corresponding plots of interferometric visibility.
	Classical GPE (red dashed lines) and Wigner calculations (solid blue lines)
	are shown. $N = 5.5 \times 10^4$,
	$\omega_x = \omega_y = 2 \pi \times 97.0\,\mathrm{Hz}$,
	$\omega_z = 2 \pi \times 11.69\,\mathrm{Hz}$,
	$a_{11} = 100.4\,a_0$, $a_{12} = 97.99\,a_0$, $a_{22} = 95.57\,a_0$~\cite{Egorov2010},
	$a_0$ is the Bohr radius.
	Particle losses:
	$\kappa_{111} = 9.0 \times 10^{-31}\,\mathrm{cm^6/s}$~\cite{Mertes2007},
	$\kappa_{12} = 7.6 \times 10^{-15}\,\mathrm{cm^3/s}$,
	$\kappa_{22} = 1.675 \times 10^{-14}\,\mathrm{cm^3/s}$~\cite{Egorov2010}.}

	\label{fig:visibility}
\end{figure}

The experiment consisted of the preparation of a stable BEC with state $\ket{1,+1}$ in a magnetic trap
together with the application of coupling electromagnetic pulses acting as beam-splitters in spin-space.
Two different schemes were used.
The first one was Ramsey interferometry with $\pi/2$-pulse in the beginning,
evolution of resulting non-equilibrium superposition of states $\ket{1,+1}$ and $\ket{2,-1}$,
and the final $\pi/2$-pulse, which transformed the phase difference of the components
into a population difference (see~\figref{fig:visibility},~(a) for the timeline).
The second scheme had an additional $\pi$-pulse (``spin echo'') in the middle of the evolution,
which partially restored interferometric contrast lost due to spatial propagation
(\figref{fig:visibility},~(c)).

% Egorov
By means of applying the first $\pi/2$ pulse a coherent superposition of two different spin-states is created.
Ground state for that consists of state $\ket{2,+1}$ focused in the center of the trap,
while state $\ket{1,-1}$ is pushed towards the edges of it~\cite{Mertes2007}.
However, after the short $\pi/2$ pulse is applied, both components have parabolic density profiles.
Thus, this superposition is not in equilibrium and collective oscillations set on,
both components oscillate about their two-component ground state.
In the process of oscillations local BEC velocities periodically become equal to almost $0$ at the turnover points,
which means that the phase of each is spatially uniform.
Both components are best overlapped at one of the turnover points as well.
As such, interference contrast of the two-component BEC in Ramsey interferometry experiences periodical revivals at that
turnover points with the frequency of collective oscillations~(\figref{fig:visibility},~(b)).
% TODO maybe need a better join of this with the next paragraph, avoid repeating visibility

The fringe visibility $\mathcal{V}$ for Ramsey and spin echo schemes is plotted on~\figref{fig:visibility}.
Simulations of Ramsey scheme experiment (\figref{fig:visibility},~(b)) with Wigner method
using typical experimental parameters
show that quantum noise amplifies the speed of visibility decay noticeably.
In the spin echo scheme (\figref{fig:visibility},~(d)) the central $\pi$-pulse is supposed
to restore visibility by inverting the evolution of the condensate.
In simulations with the Wigner method, visibility is not restored to its initial value even at revival points.
This is the effect of quantum noise, which does not appear in GPE simulations.

\begin{figure}
	\includegraphics{figures_generated/ramsey_squeezing.pdf}

	\caption{(color online).
	Simulation results for the degree of squeezing in a Ramsey sequence.
	Simulation parameters are the same as in~\figref{fig:visibility},
	except for the inter-component scattering length and corresponding loss rate~\cite{Kaufman2009}:
	$a_{12} = 80.0\,a_0$, $\kappa_{12} = 3.65 \times 10^{-13}\,\mathrm{cm^3/s}$ (blue solid line),
	$a_{12} = 85.0\,a_0$, $\kappa_{12} = 1.835 \times 10^{-13}\,\mathrm{cm^3/s}$ (red dashed line),
	$a_{12} = 90.0\,a_0$, $\kappa_{12} = 6.7 \times 10^{-14}\,\mathrm{cm^3/s}$ (green dash-dotted line) and
	$a_{12} = 96.0\,a_0$, $\kappa_{12} = 7.6 \times 10^{-15}\,\mathrm{cm^3/s}$ (black dotted line)}

	\label{fig:squeezing}
\end{figure}

The Wigner method is able to produce not only the average values of different measurables,
but their variances too.
As an example, we can obtain the degree of squeezing $\xi^2$, introduced in~\cite{Wineland1994,Sorensen2001}.
It can be defined as following~\cite{Li2009}:
\begin{equation}
\label{eqn:squeezing}
	\xi^2 = \frac{N \Delta S^2_{\perp, \mathrm{min}}}{\langle S \rangle^2},
\end{equation}
where $N$ is the total number of atoms,
$\Delta S^2_{\perp, \mathrm{min}}$ is the minimal variance of the spin in the plane orthogonal to the mean spin,
and $\langle S \rangle$ is the length of the average spin.
This parameter serves as a measure of atomic entanglement in the condensate,
when the state is entangled if $\xi^2 < 1$~\cite{Sorensen2001}.
It can be represented using moments of field operators $\widehat{\Psi}_1$
and $\widehat{\Psi}_2$~\cite{Li2009}.
The squeezing parameter in a Ramsey sequence for different scattering lengths
close and far from Feshbach resonance is plotted on~\figref{fig:squeezing}.


The important feature of these quantum dynamical simulations
is that they are able to treat large numbers of atoms (55,000 in this case),
while correctly tracking all the quantum noise sources, and also extending the simulations to long time-scales.
Both these features, of large atom numbers and long time-scales,
are essential ingredients to accurate interferometric measurements.
As can be seen from the graphs, there is excellent agreement
with experimental visibility measurements over the entire experimentally accessible time-scale.
The simulations give accurate results despite large, multi-mode dynamical motion in three dimensions
and substantial losses of most of the condensate atoms.
On longer time-scales, the experimental accuracy is limited by technical noise in the RF signal generators,
and we have no data for comparisons.

In conclusion, the Wigner method of simulating BEC dynamics
is quantitatively accurate  for a range of interferometric experiments
involving large numbers of atoms and multiple loss mechanisms.
It allows one to obtain both average values and variances of different measurables,
properly accounting for both quantum effects and nonlinear losses in the condensate.
At the same time, its computational complexity does not directly depend on number of atoms,
which, together with its perfect parallelizability, makes it very fast and scalable.
A disadvantage is that it is quite sensitive to the lattice size, the value of energy cutoff, and
the choice of basis.
Nevertheless,  the simulation results agree very well with experimental data,
allowing the method to serve as a useful predictive framework for experiments.

\bibliography{qsim.bib}

\end{document}