# Variance Control in Weak Value Measurement Pointers

A. D. Parks and J. E. Gray  Electromagnetic and Sensor Systems Department,
Naval Surface Warfare Center, Dahlgren, VA 22448, USA

(April 11, 2011)

###### Abstract

The variance of an arbitrary pointer observable is considered for the general
case that a complex weak value is measured using a complex valued pointer
state. For the typical cases where the pointer observable is either its
position or momentum, the associated expressions for the pointer’s variance
after the measurement contain a term proportional to the product of the weak
value’s imaginary part with the rate of change of the third central moment of
position relative to the initial pointer state just prior to the time of the
measurement interaction when position is the observable - or with the initial
pointer state’s third central moment of momentum when momentum is the
observable. These terms provide a means for controlling pointer position and
momentum variance and identify control conditions which - when satisfied - can
yield variances that are smaller after the measurement than they were before
the measurement. Measurement sensitivities which are useful for estimating
weak value measurement accuracies are also briefly discussed.

###### pacs:

03.65.-w, 03.65.Ca, 03.65.Ta, 06.20.Dk

LABEL:FirstPage1 LABEL:LastPage#1102

##  I Introduction

The weak value $A_{w}$ of a quantum mechanical observable $A$ was introduced
by Aharonov et al A1 ; A2 ; A3  a quarter century ago. This quantity is the
statistical result of a standard measurement procedure performed upon a pre-
and post-selected (PPS) ensemble of quantum systems when the interaction
between the measurement apparatus and each system is sufficiently weak, i.e.
when it is a weak measurement. Unlike a standard strong measurement of $A$
which significantly disturbs the measured system (i.e. it ”collapses” the
wavefunction), a weak measurement of $A$ for a PPS system does not appreciably
disturb the quantum system and yields $A_{w}$ as the observable’s measured
value. The peculiar nature of the virtually undisturbed quantum reality that
exists between the boundaries defined by the PPS states is revealed in the
eccentric characteristics of $A_{w}$, namely that $A_{w}$ can be complex
valued and that the real part $\operatorname{Re}A_{w}$ of $A_{w}$ can lie far
outside the eigenvalue spectral limits of $\widehat{A}$. While the
interpretation of weak values remains somewhat controversial, experiments have
verified several of the interesting unusual properties predicted by weak value
theory R ; P ; RL ; W ; H ; Y ; D .

The pointer of a measurement apparatus is fundamental to the theory of quantum
measurement because the values of measured observables are determined from the
pointer’s properties (e.g. from the pointer’s mean position). Understanding
these properties has become more important in recent years - in large part due
to the increased interest in the theory of weak measurements and weak value
theory. The properties of pointers associated with weak value measurements
have been studied - for example - by Johansen J , Aharonov and Botero AB ,
Jozsa Jo , Di Lorenzo and Egues DE , and Cho et al C .

The purpose of this paper is to extend Jozsa’s work Jo  to obtain the general
expression for the variance associated with an arbitrary pointer observable
when a complex valued pointer state is used to measure a complex weak value
$A_{w}$. For the typical cases where position or momentum are the pointer
observables, the associated expressions each contain a variance control term.
This term is proportional to the product of the imaginary part
$\operatorname{Im}A_{w}$ of $A_{w}$ with the rate of change of the third
central moment of position relative to the initial pointer state just prior to
measurement when the observable is position - or with the initial pointer
state’s third central moment of momentum when momentum is the observable.
Control conditions associated with these terms are identified which - if
satisfied - can yield pointer position and momentum variances after a
measurement that are smaller than they were prior to the measurement. These
results are used to briefly discuss sensitivities associated with weak value
measurements.

##  II Weak Measurements and Weak Values

For the reader’s convenience, this section provides a brief review of weak
measurement and weak value theory. For additional details the reader is
invited to consult references A1 ; A2 ; A3 ; A4 .

Weak measurements arise in the von Neumann description of a quantum
measurement at time $t_{0}$ of a time-independent observable $A$ that
describes a quantum system in an initial fixed pre-selected state
$\left|\psi_{i}\right\rangle={\displaystyle\sum\nolimits_{J}}c_{j}\left|a_{j}\right\rangle$
at $t_{0}$, where the set $J$ indexes the eigenstates
$\left|a_{j}\right\rangle$ of $\widehat{A}$. In this description, the
Hamiltonian for the interaction between the measurement apparatus and the
quantum system is

$\widehat{H}=\gamma(t)\widehat{A}\widehat{p}.$

Here $\gamma(t)=\gamma\delta(t-t_{0})$ defines the strength of the
measurement’s impulsive coupling interaction at $t_{0}$ and $\widehat{p}$ is
the momentum operator for the pointer of the measurement apparatus which is in
the initial normalized state $\left|\phi\right\rangle$. Let $\widehat{q}$ be
the pointer’s position operator that is conjugate to $\widehat{p}$.

Prior to the measurement the pre-selected system and the pointer are in the
tensor product state $\left|\psi_{i}\right\rangle\left|\phi\right\rangle$.
Immediately following the interaction the combined system is in the state

$\left|\Phi\right\rangle=e^{-\frac{i}{\hbar}\int\widehat{H}dt}\left|\psi_{i}\right\rangle\left|\phi\right\rangle={\displaystyle\sum\nolimits_{J}}c_{j}e^{-\frac{i}{\hbar}\gamma
a_{j}\widehat{p}}\left|a_{j}\right\rangle\left|\phi\right\rangle,$

where use has been made of the fact that
$\int\widehat{H}dt=\gamma\widehat{A}\widehat{p}$. The exponential factor in
this equation is the translation operator $\widehat{S}\left(\gamma
a_{j}\right)$ for $\left|\phi\right\rangle$ in its $q$-representation. It is
defined by the action $\left\langle q\right|\widehat{S}\left(\gamma
a_{j}\right)\left|\phi\right\rangle$ which translates the pointer’s
wavefunction over a distance $\gamma a_{j}$ parallel to the $q$-axis. The
$q$-representation of the combined system and pointer state is

$\left\langle
q\right|\left.\Phi\right\rangle={\displaystyle\sum\nolimits_{J}}c_{j}\left\langle
q\right|\widehat{S}\left(\gamma
a_{j}\right)\left|\phi\right\rangle\left|a_{j}\right\rangle.$

When the measurement interaction is strong, the quantum system is appreciably
disturbed and its state ”collapses” with probability $\left|c_{n}\right|^{2}$
to an eigenstate $\left|a_{n}\right\rangle$ leaving the pointer in the state
$\left\langle q\right|\widehat{S}\left(\gamma
a_{n}\right)\left|\phi\right\rangle$. Strong measurements of an ensemble of
identically prepared systems yield $\gamma\left\langle
A\right\rangle\equiv\gamma\left\langle\psi_{i}\right|\widehat{A}\left|\psi_{i}\right\rangle$
as the centroid of the associated pointer probability distribution with
$\left\langle A\right\rangle$ as the measured value of $\widehat{A}$.

A weak measurement of $\widehat{A}$ occurs when the interaction strength
$\gamma$ is sufficiently small so that the system is essentially undisturbed
and the uncertainty $\Delta q$ is much larger than $\widehat{A}$’s eigenvalue
separation. In this case, the pointer distribution is the superposition of
broad overlapping $\left|\left\langle q\right|\widehat{S}\left(\gamma
a_{j}\right)\left|\phi\right\rangle\right|^{2}$ terms. Although a single
measurement provides little information about $\widehat{A}$, many repetitions
allow the centroid of the distribution to be determined to any desired
accuracy.

If a system is post-selected after a weak measurement is performed, then the
resulting pointer state is

$\left|\Psi\right\rangle\equiv\left\langle\psi_{f}\right|\left.\Phi\right\rangle={\displaystyle\sum\nolimits_{J}}c_{j}^{\prime\ast}c_{j}\widehat{S}\left(\gamma
a_{j}\right)\left|\phi\right\rangle,$

where
$\left|\psi_{f}\right\rangle={\displaystyle\sum\nolimits_{J}}c_{j}^{\prime}\left|a_{j}\right\rangle$,
$\left\langle\psi_{f}\right|\left.\psi_{i}\right\rangle\neq 0$, is the post-
selected state at $t_{0}$. Since

$\widehat{S}\left(\gamma
a_{j}\right)={\displaystyle\sum\limits_{m=0}^{\infty}}\frac{\left(-i\gamma
a_{j}\widehat{p}/\hbar\right)^{m}}{m!},$

then

$\left|\Psi\right\rangle={\displaystyle\sum\nolimits_{J}}c_{j}^{\prime\ast}c_{j}\left\\{1-\frac{i}{\hbar}\gamma
A_{w}\widehat{p}+{\displaystyle\sum\limits_{m=2}^{\infty}}\frac{\left(-i\gamma\widehat{p}/\hbar\right)^{m}}{m!}\left(A^{m}\right)_{w}\right\\}\left|\phi\right\rangle\approx\left\\{{\displaystyle\sum\nolimits_{J}}c_{j}^{\prime\ast}c_{j}\right\\}e^{-\frac{i}{\hbar}\gamma
A_{w}\widehat{p}}\left|\phi\right\rangle$

in which case

$\left|\Psi\right\rangle\approx\left\langle\psi_{f}\right|\left.\psi_{i}\right\rangle\widehat{S}\left(\gamma
A_{w}\right)\left|\phi\right\rangle.$

Here

$\left(A^{m}\right)_{w}=\frac{{\displaystyle\sum\nolimits_{J}}c_{j}^{\prime\ast}c_{j}a_{j}^{m}}{{\displaystyle\sum\nolimits_{J}}c_{j}^{\prime\ast}c_{j}}=\frac{\left\langle\psi_{f}\right|\widehat{A}^{m}\left|\psi_{i}\right\rangle}{\left\langle\psi_{f}\right|\left.\psi_{i}\right\rangle},$

with the weak value $A_{w}$ of $\widehat{A}$ defined by

$A_{w}\equiv\left(A^{1}\right)_{w}=\frac{\left\langle\psi_{f}\right|\widehat{A}\left|\psi_{i}\right\rangle}{\left\langle\psi_{f}\right|\left.\psi_{i}\right\rangle}.$
(1)

From this expression it is obvious that $A_{w}$ is - in general - a complex
valued quantity that can be calculated directly from theory and that when the
PPS states are nearly orthogonal $\operatorname{Re}A_{w}$ can lie far outside
$\widehat{A}$’s eigenvalue spectral limits.

For the general case where both $A_{w}$ and $\phi\left(q\right)$ are complex
valued, the mean pointer position and momentum after a measurement are given
by Jo

$\left\langle\Psi\right|\widehat{q}\left|\Psi\right\rangle=\left\langle\phi\right|\widehat{q}\left|\phi\right\rangle+\gamma\operatorname{Re}A_{w}+\left(\frac{\gamma}{\hbar}\right)\operatorname{Im}A_{w}\left(m\frac{d\Delta_{\phi}^{2}q}{dt}\right)$
(2)

and

$\left\langle\Psi\right|\widehat{p}\left|\Psi\right\rangle=\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle+2\left(\frac{\gamma}{\hbar}\right)\operatorname{Im}A_{w}\left(\Delta_{\phi}^{2}p\right),$
(3)

respectively. Here $m$ is the mass of the pointer, $\Delta_{\phi}^{2}p$ is the
pointer’s initial momentum variance, and the time derivative of
$\Delta_{\phi}^{2}q$ is the rate of change of the initial pointer position
variance just prior to $t_{0}$.

##  III Pointer Variance

The mean value of an arbitrary pointer observable $M$ after a measurement of
$A_{w}$ is Jo

$\displaystyle\left\langle\Psi\right|\widehat{M}\left|\Psi\right\rangle$
$\displaystyle=\left\langle\phi\right|\widehat{M}\left|\phi\right\rangle-i\left(\frac{\gamma}{\hbar}\right)\operatorname{Re}A_{w}\left\langle\phi\right|\left[\widehat{M},\widehat{p}\right]\left|\phi\right\rangle+$
(4)
$\displaystyle\left(\frac{\gamma}{\hbar}\right)\operatorname{Im}A_{w}\left(\left\langle\phi\right|\left\\{\widehat{M},\widehat{p}\right\\}\left|\phi\right\rangle-2\left\langle\phi\right|\widehat{M}\left|\phi\right\rangle\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle\right),$

where
$\left\\{\widehat{M},\widehat{p}\right\\}=\widehat{M}\widehat{p}+\widehat{p}\widehat{M}$.
Note that eq.(4) reduces to eq.(3) when $\widehat{M}=\widehat{p}$ and that it
is also in complete agreement with eq.(2) when $\widehat{M}=\widehat{q}$ since
$\left[\widehat{q},\widehat{p}\right]=i\hbar$ and the equations of motion for
$\left\langle\phi\right|\widehat{q}\left|\phi\right\rangle$ and
$\left\langle\phi\right|\widehat{q}^{2}\left|\phi\right\rangle$ yield

$\left\langle\phi\right|\left\\{\widehat{q},\widehat{p}\right\\}\left|\phi\right\rangle=m\frac{d\left\langle\phi\right|\widehat{q}^{2}\left|\phi\right\rangle}{dt}$
(5)

and

$2\left\langle\phi\right|\widehat{q}\left|\phi\right\rangle\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle=m\frac{d\left\langle\phi\right|\widehat{q}\left|\phi\right\rangle^{2}}{dt}.$
(6)

Here the time derivatives are rates of change of the corresponding quantities
just prior to the interaction time $t_{0}$.

The pointer variance for $M$ is easily determined from eq.(4) by subtracting
its square from the expression obtained from eq.(4) when $\widehat{M}$ is
replaced by $\widehat{M}^{2}$. Retaining terms through first order in
$\left(\frac{\gamma}{\hbar}\right)$ yields the following result:

$\Delta_{\Psi}^{2}M=\Delta_{\phi}^{2}M-i\left(\frac{\gamma}{\hbar}\right)\operatorname{Re}A_{w}\mathcal{F}\left(\widehat{M}\right)+\left(\frac{\gamma}{\hbar}\right)\operatorname{Im}A_{w}\mathcal{G}\left(\widehat{M}\right).$
(7)

Here $\Delta_{\phi}^{2}M$ and $\Delta_{\Psi}^{2}M$ are the initial and final
variances, respectively,

$\mathcal{F}\left(\widehat{M}\right)\equiv\left\langle\phi\right|\left[\widehat{M}^{2},\widehat{p}\right]\left|\phi\right\rangle-2\left\langle\phi\right|\widehat{M}\left|\phi\right\rangle\left\langle\phi\right|\left[\widehat{M},\widehat{p}\right]\left|\phi\right\rangle,$

and

$\mathcal{G}\left(\widehat{M}\right)\equiv\left\langle\phi\right|\left\\{\widehat{M}^{2},\widehat{p}\right\\}\left|\phi\right\rangle-2\left\langle\phi\right|\widehat{M}\left|\phi\right\rangle\left\langle\phi\right|\left\\{\widehat{M},\widehat{p}\right\\}\left|\phi\right\rangle-2\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle\left(\Delta_{\phi}^{2}M-\left\langle\phi\right|\widehat{M}\left|\phi\right\rangle^{2}\right).$
(8)

As anticipated from eq.(4), eq.(7) clearly shows that for such a measurement
the pointer variance associated with an arbitrary pointer observable is also
generally effected by both the real and imaginary parts of the weak value.

However, for the typical cases of interest where $\widehat{M}=\widehat{q}$ or
$\widehat{M}=\widehat{p}$, the pointer’s variance is independent of
$\operatorname{Re}A_{w}$ because

$\mathcal{F}\left(\widehat{q}\right)=0=\mathcal{F}\left(\widehat{p}\right).$

Here use has been made of the facts that
$\left[\widehat{q},\widehat{p}\right]=i\hbar$ and
$\left[\widehat{q}^{2},\widehat{p}\right]=2i\hbar\widehat{q}$. Consequently,
for these cases eq.(7) can be written as

$\Delta_{\Psi}^{2}M=\Delta_{\phi}^{2}M+\left(\frac{\gamma}{\hbar}\right)\operatorname{Im}A_{w}\mathcal{G}\left(\widehat{M}\right),\text{
}M=q,p.$ (9)

Now consider $\mathcal{G}\left(\widehat{M}\right)$ in more detail. When
$\widehat{M}=\widehat{q}$, then eq.(8) becomes

$\mathcal{G}\left(\widehat{q}\right)=\left\langle\phi\right|\left\\{\widehat{q}^{2},\widehat{p}\right\\}\left|\phi\right\rangle-2\left\langle\phi\right|\widehat{q}\left|\phi\right\rangle\left\langle\phi\right|\left\\{\widehat{q},\widehat{p}\right\\}\left|\phi\right\rangle-2\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle(\Delta_{\phi}^{2}q-\left\langle\phi\right|\widehat{q}\left|\phi\right\rangle^{2}).$
(10)

From the equation of motion for
$\left\langle\phi\right|\widehat{q}^{3}\left|\phi\right\rangle$ it is found
that

$\frac{d\left\langle\phi\right|\widehat{q}^{3}\left|\phi\right\rangle}{dt}=-\frac{i}{\hbar}\left\langle\phi\right|\left[\widehat{q}^{3},\widehat{H}\right]\left|\phi\right\rangle=-\frac{i}{2m\hbar}\left\langle\phi\right|\left[\widehat{q}^{3},\widehat{p}^{2}\right]\left|\phi\right\rangle=\frac{3}{2m}\left\langle\phi\right|\left\\{\widehat{q}^{2},\widehat{p}\right\\}\left|\phi\right\rangle,$

where $\widehat{H}=\frac{\widehat{p}^{2}}{2m}+V(\widehat{q})$ is the pointer’s
Hamiltonian operator. Applying this result - along with eqs.(5) and (6) - to
eq.(10) yields

$\mathcal{G}\left(\widehat{q}\right)=\frac{2m}{3}\frac{dq_{3}}{dt},$

so that eq.(9) can be compactly written as

$\Delta_{\Psi}^{2}q=\Delta_{\phi}^{2}q+\frac{2\gamma
m}{3\hbar}\operatorname{Im}A_{w}\left(\frac{dq_{3}}{dt}\right).$

Here
$q_{3}\equiv\left\langle\phi\right|\left(\widehat{q}-\left\langle\phi\right|\widehat{q}\left|\phi\right\rangle\right)^{3}\left|\phi\right\rangle$
is the third central moment of $\widehat{q}$ relative to the initial pointer
state and its time derivative is the rate of change of $q_{3}$ just prior to
$t_{0}$.

When $\widehat{M}=\widehat{p}$, then eq.(8) becomes

$\mathcal{G}\left(\widehat{p}\right)=2\left[\left\langle\phi\right|\widehat{p}^{3}\left|\phi\right\rangle-3\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle\left\langle\phi\right|\widehat{p}^{2}\left|\phi\right\rangle+2\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle^{3}\right]=2p_{3},$

where
$p_{3}\equiv\left\langle\phi\right|\left(\widehat{p}-\left\langle\phi\right|\widehat{p}\left|\phi\right\rangle\right)^{3}\left|\phi\right\rangle$
is the third central moment of $\widehat{p}$ relative to the pointer’s initial
state, and eq.(9) assumes the form

$\Delta_{\Psi}^{2}p=\Delta_{\phi}^{2}p+2\left(\frac{\gamma}{\hbar}\right)\operatorname{Im}A_{w}\left(p_{3}\right).$

The quantities $q_{3}$ and $p_{3}$ provide measures of the skewness of the
pointer position and momentum probability distribution profiles. If the
pointer position profile’s skewness is fixed, then $\frac{dq_{3}}{dt}=0$ and
$\Delta_{\Psi}^{2}q=\Delta_{\phi}^{2}q$. Otherwise, $\Delta_{\Psi}^{2}q$ can
be manipulated through the judicious selection of the control term
$\operatorname{Im}A_{w}\left(\frac{dq_{3}}{dt}\right)$. In particular, observe
that $0<\Delta_{\Psi}^{2}q\leq\Delta_{\phi}^{2}q$ when this control term
satisfies the inequality

$\text{ }-\left(\frac{3\hbar}{2\gamma
m}\right)\Delta_{\phi}^{2}q<\operatorname{Im}A_{w}\left(\frac{dq_{3}}{dt}\right)\leq
0.$ (11)

Similarly, the control term $\operatorname{Im}A_{w}\left(p_{3}\right)$ can be
used to manipulate $\Delta_{\Psi}^{2}p$ when it satisfies the inequality

$-\left(\frac{\hbar}{2\gamma}\right)\Delta_{\phi}^{2}p<\operatorname{Im}A_{w}\left(p_{3}\right)\leq
0.$ (12)

Thus, when measuring complex weak values the final pointer position (momentum)
variance can be made smaller than its initial value by choosing
$\operatorname{Im}A_{w}$ or $\frac{dq_{3}}{dt}$ ($p_{3}$) so that condition
(11) ((12)) is satisfied.

##  IV Closing Remarks

Because of the growing interest in the practical application of weak values,
estimating their measurement sensitivities has also become important from both
the experimental and device engineering perspectives. Applying the calculus of
error propagation to the above results defines the measurement sensitivities
$\delta_{q}\operatorname{Re}A_{w}$ and $\delta_{q}\operatorname{Im}A_{w}$ for
determining $\operatorname{Re}A_{w}$ and $\operatorname{Im}A_{w}$ from the
mean pointer position. These sensitivities are the positive square roots of
the following expressions:

$\delta_{q}^{2}\operatorname{Re}A_{w}\equiv\frac{\Delta_{\Psi}^{2}q}{\left|\frac{\partial\left\langle\Psi\right|\widehat{q}\left|\Psi\right\rangle}{\partial\operatorname{Re}A_{w}}\right|^{2}}=\frac{\Delta_{\phi}^{2}q}{\gamma^{2}}+\frac{2m}{3\gamma\hbar}\operatorname{Im}A_{w}\left(\frac{dq_{3}}{dt}\right)$
(13)

(this quantity is obviously undefined when $A_{w}$ is purely imaginary) and

$\delta_{q}^{2}\operatorname{Im}A_{w}\equiv\frac{\Delta_{\Psi}^{2}q}{\left|\frac{\partial\left\langle\Psi\right|\widehat{q}\left|\Psi\right\rangle}{\partial\operatorname{Im}A_{w}}\right|^{2}}=\left(\frac{\hbar}{\gamma
m}\right)^{2}\left(\frac{\Delta_{\phi}^{2}q}{\left|\frac{d\Delta_{\phi}^{2}q}{dt}\right|^{2}}\right)+\frac{2}{3}\left(\frac{\hbar}{\gamma
m}\right)\operatorname{Im}A_{w}\left(\frac{\frac{dq_{3}}{dt}}{\left|\frac{d\Delta_{\phi}^{2}q}{dt}\right|^{2}}\right),\frac{d\Delta_{\phi}^{2}q}{dt}\neq
0$ (14)

(this quantity is obviously undefined when $A_{w}$ is real valued or when
$\frac{d\Delta_{\phi}^{2}q}{dt}=0$ \- in which case the mean position does not
depend upon $\operatorname{Im}A_{w}$). It is clear from eqs.(13) and (14)
that: (i) these measurement sensitivities depend upon the variance control
term $\operatorname{Im}A_{w}\left(\frac{dq_{3}}{dt}\right)$ and that this
dependence vanishes when $q_{3}$ is fixed (or if $A_{w}$ is real valued); (ii)
these measurement accuracies decrease (increase) as the measurement gets
weaker (stronger) - i.e. as $\gamma$ gets smaller (larger); (iii) in principle
- the accuracies associated with measuring both $\operatorname{Re}A_{w}$ and
$\operatorname{Im}A_{w}$ can be arbitrarily increased (for a fixed $\gamma>0$
and $m$) by invoking condition (11) and choosing
$\operatorname{Im}A_{w}\left(\frac{dq_{3}}{dt}\right)=-\left(\frac{3\hbar}{2\gamma
m}\right)\Delta_{\phi}^{2}q+\epsilon$, where $\epsilon$ is a small positive
real number; and (iv) surprisingly, the measurement accuracy for
$\operatorname{Re}A_{w}$ decreases with increasing pointer mass, whereas that
for $\operatorname{Im}A_{w}$ increases.

The sensitivity $\delta_{p}\operatorname{Im}A_{w}$ for determining
$\operatorname{Im}A_{w}$ from the mean pointer momentum is the positive square
root of

$\delta_{p}^{2}\operatorname{Im}A_{w}\equiv\frac{\Delta_{\Psi}^{2}p}{\left|\frac{\partial\left\langle\Psi\right|\widehat{p}\left|\Psi\right\rangle}{\partial\operatorname{Im}A_{w}}\right|^{2}}=\left(\frac{\hbar}{2\gamma}\right)^{2}\left(\frac{1}{\Delta_{\phi}^{2}p}\right)+\left(\frac{\hbar}{2\gamma}\right)\operatorname{Im}A_{3}\left(\frac{p_{3}}{\left(\Delta_{\phi}^{2}p\right)^{2}}\right),\Delta_{\phi}^{2}p\neq
0$ (15)

(this quantity is undefined when $A_{w}$ is real valued). Inspection of
eq.(15) reveals that for such measurements: (i) the sensitivity depends upon
the variance control term $\operatorname{Im}A_{w}\left(p_{3}\right)$; (ii) the
accuracy decreases (increases) as the measurement gets weaker (stronger) -
i.e. as $\gamma$ gets smaller (larger); and (iii) the accuracy can be
arbitrarily increased (for a fixed $\gamma>0$) via eq.(12) by choosing
$\operatorname{Im}A_{w}\left(p_{3}\right)=-\left(\frac{\hbar}{2\gamma}\right)\Delta_{\phi}^{2}p+\epsilon$,
where $\epsilon$ is again a small positive real number.

In closing, it is important to note that the results discussed and developed
above apply when the measurement interaction is instantaneous and the
measurement is read from the pointer immediately after the interaction DE .

##  V Acknowledgement

This work was supported in part by a grant from the Naval Surface Warfare
Center Dahlgren Division’s In-house Laboratory Independent Research Program
sponsored by the Office of Naval Research.

##  VI References

## References

  * (1) Aharonov Y, Albert D, Casher D and Vaidman L 1986 New Techniques and Ideas in Quantum Measurement Theory ed D Greenberger (New York: New York Academy of Sciences) p 417
  * (2) Aharonov Y, Albert D and Vaidman L 1988 Phys. Rev. Lett. 60 1351
  * (3) Aharonov Y and Vaidman L 1990 Phys. Rev. A 41 11
  * (4) Ritchie N, Storey J and Hulet R 1991 Phys. Rev. Lett. 66 1107
  * (5) Parks A, Cullin D and Stoudt D 1998 Proc. R. Soc. A 454 2997
  * (6) Resch K, Lundeen J and Steinberg A 2004 Phys. Lett. A 324 125
  * (7) Wang Q, Sun F, Zhang Y, Li J, Huang Y and Guo G 2006 Phys. Rev. A 73 023814
  * (8) Hosten O and Kwiat P 2008 Science 319 787
  * (9) Yokota K, Yamamoto T, Koashi M and Imoto N 2009 N. J. Phys. 11 033011
  * (10) Dixon P, Starling D, Jordan A and Howell J 2009 Phys. Rev. Lett. 102 173601
  * (11) Johansen L 2004 Phys. Rev. Lett. 93 120402
  * (12) Aharonov Y and Botero A 2005 Phys. Rev. A 72 052111
  * (13) Jozsa R 2007 Phys. Rev. A 76 044103
  * (14) Di Lorenzo A and Egues J 2008 Phys. Rev. A 77 042108
  * (15) Cho Y, Lim H, Ra Y and Kim Y 2010 N. J. Phys. 12 023036
  * (16) Aharonov Y and Rohrlich D 2005 Quantum Paradoxes: Quantum Theory for the Perplexed (Wiley-VCH, Weinheim) p.225


# Dissipative Preparation of Antiferromagnetic Order in the Fermi-Hubbard
Model

J. Kaczmarczyk<EMAIL_ADDRESS>Institute of Science and Technology
Austria, Am Campus 1, 3400 Klosterneuburg, Austria  H. Weimer  Institut für
Theoretische Physik, Leibniz Universität Hannover, Appelstr. 2, 30167
Hannover, Germany  M. Lemeshko<EMAIL_ADDRESS>Institute of
Science and Technology Austria, Am Campus 1, 3400 Klosterneuburg, Austria

###### Abstract

The Fermi-Hubbard model is one of the key models of condensed matter physics,
which holds a potential for explaining the mystery of high-temperature
superconductivity. Recent progress in ultracold atoms in optical lattices has
paved the way to studying the model’s phase diagram using the tools of quantum
simulation, which emerged as a promising alternative to the numerical
calculations plagued by the infamous sign problem. However, the temperatures
achieved using elaborate laser cooling protocols so far have been too high to
show the appearance of antiferromagnetic and superconducting quantum phases
directly. In this work, we demonstrate that using the machinery of dissipative
quantum state engineering, one can observe the emergence of the
antiferromagnetic order in the Fermi-Hubbard model with fermions in optical
lattices. The core of the approach is to add incoherent laser scattering in
such a way that the antiferromagnetic state emerges as the dark state of the
driven-dissipative dynamics. The proposed controlled dissipation channels
described in this work are straightforward to add to already existing
experimental setups.

Hubbard model, ultracold gases, antiferromagnetic phase, lattice fermion
models, dissipative preparation

###### pacs:

67.85.-d, 75.10.Jm, 71.10.Fd

††preprint: New. J. Phys. 18, 093042 (2016); doi:10.1088/1367-2630/18/9/093042

##  I Introduction

Experimental progress with ultracold fermions in optical lattices Lewenstein
_et al._ (2007); Bloch  _et al._ (2008) leads the way to achieving one of the
key goals of quantum simulation Georgescu  _et al._ (2014) – mimicking
realistic condensed matter systems. To date, the experiments covered a broad
range of systems and interaction regimes, from probing the BEC-BCS crossover
in lattices Chin  _et al._ (2006), to the observation of a fermionic Mott
insulator Schneider  _et al._ (2008); Jordens  _et al._ (2008), to studying
short range magnetism Greif  _et al._ (2013) and multiflavor spin dynamics
Krauser  _et al._ (2012), to realizing topological Haldane model Jotzu  _et
al._ (2014) and artificial graphene sheets Uehlinger  _et al._ (2013). These
discoveries pave the way to use ultracold atoms to reveal the properties of
the repulsive Fermi-Hubbard model Esslinger (2010); Le Hur and Rice (2009).
The latter is of particular importance since it represents a playground to get
insight into the physics of high-temperature superconductivity and related
phenomena observed in the cuprates Keimer  _et al._ (2015).

In the case of one particle per site and large on-site interaction, $U$, the
Fermi-Hubbard model exhibits the transition to the Mott-insulating state
Schneider  _et al._ (2008); Jordens  _et al._ (2008) around the temperature
$T\sim U$. If the temperature is decreased further and reaches the so-called
‘Néel temperature,’ $T_{\rm N}\sim 4t^{2}/U$, where $t$ gives the hopping rate
between neighboring sites, the transition to the antiferromagnetic (AF) phase
is expected Hart  _et al._ (2015); Greif  _et al._ (2015). Currently the
temperatures achievable in experiment are slightly above the Néel temperature
where AF correlations can already be observed, for instance, $T/T_{\rm
N}\approx 1.42$ has been reached in Ref. Hart  _et al._ (2015). Ultimately, in
order to study the superconducting phase or other phenomena related to pairing
in high-temperature superconductors, the temperature needs to be substantially
lower. Therefore, due to the experimental limitations inherent to the standard
laser cooling techniques, it is crucial to develop alternative approaches Ho
and Zhou (2009); Bernier  _et al._ (2009); Heidrich-Meisner  _et al._ (2009);
Diehl  _et al._ (2010a); Williams  _et al._ (2010); Medley  _et al._ (2011);
Lee  _et al._ (2011); Colomé-Tatché  _et al._ (2011); Lemeshko  _et al._
(2012); Yi  _et al._ (2012); Carr and Saffman (2013); Murmann  _et al._
(2015); Mazzucchi  _et al._ (2016) to preparation of quantum phases in optical
lattices.

In this work, we propose an efficient scheme for the preparation of
antiferromagnetic order in the Fermi-Hubbard model, based on the ideas of
dissipative state engineering, which have recently emerged in the context of
many-particle systems Diehl  _et al._ (2008); Kraus  _et al._ (2008); Witthaut
_et al._ (2008); Verstraete  _et al._ (2009); Weimer  _et al._ (2010); Diehl
_et al._ (2010a, 2011); Yi  _et al._ (2012); Alharbi and Ficek (2010);
Watanabe and Mäkelä (2012); Kordas  _et al._ (2012); Bermudez  _et al._
(2013); Rao and Mølmer (2013); Carr and Saffman (2013); Lemeshko and Weimer
(2013); Weimer (2013); Lee  _et al._ (2011); Otterbach and Lemeshko (2014);
Lucia  _et al._ (2015); Budich  _et al._ (2015); Weimer (2015a) and have been
implemented experimentally Barreiro  _et al._ (2010); Gommers  _et al._
(2005); Barreiro  _et al._ (2011); Krauter  _et al._ (2011); Schindler  _et
al._ (2013); Lin  _et al._ (2013); Shankar  _et al._ (2013). In such
scenarios, a many-body state of interest (here: states exhibiting AF order) is
prepared as a steady state of the quantum master equation governing the open
system dynamics, as opposed to the ground state of the Hamiltonian. Such
steady state can undergo quantum phase transitions to an ordered state of
matter, which can be classified in close analogy to equilibrium systems Diehl
_et al._ (2010b); Tomadin  _et al._ (2011); Lee  _et al._ (2011); Le Boité
_et al._ (2013); Jin  _et al._ (2013); Sieberer  _et al._ (2013); Lee  _et
al._ (2013); Weimer (2015a); Chan  _et al._ (2015).

We start with a system of fermions in an optical lattice as described by the
Fermi-Hubbard model. The parameters of the Hamiltonian are left intact,
instead we introduce dissipative channels on top of the unitary evolution. As
a result, fermions remain mobile in the optical lattice during the entire
dissipative preparation stage, which should help with retaining coherence
after the dissipation channels are switched off. Furthermore, the dissipation
channels of our scheme are implemented using the level structure of fermionic
40K, currently used in several laboratories Krauser  _et al._ (2012); Greif
_et al._ (2013, 2015); Cheuk  _et al._ (2015); Haller  _et al._ (2015); Edge
_et al._ (2015). Consequently, the presented scheme can be readily implemented
into already existing experimental setups.

Theoretical description of open many-body quantum systems represents a
challenging task and is currently an active field of research Finazzi  _et
al._ (2015); Weimer (2015b, a); Overbeck and Weimer (2016); Bernier  _et al._
(2014); Prosen (2008); Dzhioev and Kosov (2015). In our analysis of the
dissipative Fermi-Hubbard model we use two complementary techniques: the Monte
Carlo wave function (MCWF) Dalibard  _et al._ (1992); Mølmer  _et al._ (1993);
Carmichael (1993) and the variational method Weimer (2015b, a), which is
generalized here to the description of fermionic systems at half-filling. By
using these two methods we demonstrate that a substantial AF magnetization is
present in the system both for an exact solution on a $3\times 3$ lattice, as
well as in the thermodynamic limit.

##  II  The dissipative Fermi-Hubbard model

We start with the Fermi-Hubbard Hamiltonian

$\hat{H}=\sum_{\mathbf{i},\mathbf{j},\sigma}t_{\mathbf{i}\mathbf{j}}\hat{c}^{\dagger}_{\mathbf{i}\sigma}\hat{c}_{\mathbf{j}\sigma}+U\sum_{\mathbf{i}}\hat{n}_{\mathbf{i}\uparrow}\hat{n}_{\mathbf{i}\downarrow},$
(1)

which has been experimentally realized in a range of systems such as 6Li Hart
_et al._ (2015) and 40K Jordens  _et al._ (2008). Our goal is to design
dissipative processes in such a way that the state with an AF order is the
dark state of the dissipative dynamics and the time evolution of the open
system will drive it towards such a dark state.

The dynamics of an open quantum system is governed by the master equation for
the system’s density matrix

$\dot{\rho}=-i\left[\hat{H},\rho\right]+\sum_{\mathbf{i},\mathbf{j},\sigma,\alpha}{}^{{}^{\prime}}\left(\hat{j}^{(\alpha)}_{\mathbf{i}\mathbf{j},\sigma}\,\;\rho\,\;\hat{j}^{(\alpha)\dagger}_{\mathbf{i}\mathbf{j},\sigma}-\frac{1}{2}\left\\{\hat{j}^{(\alpha)\dagger}_{\mathbf{i}\mathbf{j},\sigma}\,\;\hat{j}^{(\alpha)}_{\mathbf{i}\mathbf{j},\sigma},\rho\right\\}\right),$
(2)

where we set $\hbar\equiv 1$ and the primed sum runs over nearest-neighbor
sites. Since we start with a disordered sample, all possible nearest-neighbor
configurations, including $|\uparrow;\uparrow\rangle$,
$|\downarrow;\downarrow\rangle$, and $|\downarrow\uparrow;0\rangle$ will be
present. The jump operators, therefore, need to convert the latter into those
with the local antiferromagnetic order, $|\uparrow;\downarrow\rangle$.

We choose the jump operators to be as follows:

$\displaystyle\hat{j}^{(1)}_{\mathbf{i}\mathbf{j},\uparrow}$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{1}}\,\hat{n}_{\mathbf{i}\uparrow}(1-\hat{n}_{\mathbf{j}\downarrow})\hat{c}^{\dagger}_{\mathbf{i}\downarrow}\hat{c}_{\mathbf{j}\uparrow},\qquad\hat{j}^{(2)}_{\mathbf{i}\mathbf{j},\uparrow}=\left(\hat{j}^{(1)}_{\mathbf{i}\mathbf{j},\uparrow}\right)^{\dagger},$
(3) $\displaystyle\hat{j}^{(1)}_{\mathbf{i}\mathbf{j},\downarrow}$
$\displaystyle=$
$\displaystyle\sqrt{\gamma_{1}}\,\hat{n}_{\mathbf{i}\downarrow}(1-\hat{n}_{\mathbf{j}\uparrow})\hat{c}^{\dagger}_{\mathbf{i}\uparrow}\hat{c}_{\mathbf{j}\downarrow},\qquad\hat{j}^{(2)}_{\mathbf{i}\mathbf{j},\downarrow}=\left(\hat{j}^{(1)}_{\mathbf{i}\mathbf{j},\downarrow}\right)^{\dagger},$
(4) $\displaystyle\hat{j}^{(3)}_{\mathbf{i}\mathbf{j}}$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{2}}\,(1-\hat{n}_{\mathbf{i}\uparrow})\hat{n}_{\mathbf{j}\uparrow}\hat{c}^{\dagger}_{\mathbf{i}\downarrow}\hat{c}_{\mathbf{j}\downarrow},$
(5)

where $\gamma_{1}$, $\gamma_{2}$ are the dissipation rates. The resulting
dissipative dynamics is visualized in figure 1(a): the jump operators
$\hat{j}^{(1)}_{\mathbf{i}\mathbf{j},\sigma}$ (with
$\sigma=\uparrow,\downarrow$) turn the configurations on neighboring sites
from $|\sigma;\sigma\rangle$ to $|\uparrow\downarrow;0\rangle$, whereas
$\hat{j}^{(2)}_{\mathbf{i}\mathbf{j},\sigma}$ act in the opposite direction.
These processes are labelled by the amplitude $\gamma_{1}$ in the figure.
Finally, the jump operators $\hat{j}^{(3)}_{\mathbf{i}\mathbf{j}}$ turn the
configurations $|\uparrow\downarrow;0\rangle$, into those with the AF order,
${|\uparrow;\downarrow\rangle}$ (as labelled by $\gamma_{2}$). As a result,
the dissipative dynamics drives the system towards the AF phase.

We note that the jump operators $\hat{j}^{(3)}_{\mathbf{i}\mathbf{j}}$ break
the SU(2) symmetry, as the down-spin atom becomes more mobile. This, however,
does not constitute a limitation of our scheme. Moreover, it is possible to
reestablish the symmetry by using additional auxiliary states to induce
hopping of the $|\uparrow\rangle$-state atom away from the double-occupancy
configuration.

As we discuss in the following section, such choice of the jump operators is
straighforward to realize in experiment using incoherent laser scattering.

Figure 1: Illustration of the dissipative processes. (a) Action of the jump
operators on the nearest-neighbor sites; and (b)-(c) their implementation
using Raman-assisted hopping. The Raman beams are labelled as
$\Omega_{r1},...,\Omega_{r5}$, the pumping beam by $\Omega$, and the decay
rate is given by $\gamma$; $\Delta_{\uparrow\downarrow}$ gives the hyperfine
splitting between the $|\uparrow\rangle$ and $|\downarrow\rangle$ states,
$\Delta_{X\downarrow}$ gives the Zeeman splitting between the $|X\rangle$ and
$|\downarrow\rangle$ states, while $U$ denotes the on-site interaction between
the $|\uparrow\rangle$ and $|\downarrow\rangle$ states.

##  III Experimental implementation of the jump operators

Simulating the Fermi-Hubbard model requires mapping the two spin states,
$|\uparrow\rangle$ and $|\downarrow\rangle$, onto the fine or hyperfine
components of the ground electronic state manifold of an ultracold atom. The
on-site interaction between the spin components, $U$, can be tuned using a
Feshbach resonance Chin  _et al._ (2010). We exemplify the scheme using the
atomic level structure of fermionic 40K Tiecke (2011), with the levels
$|\uparrow\rangle\equiv|F=\frac{7}{2};m_{F}=-\frac{7}{2}\rangle$ and
$|\downarrow\rangle\equiv|\frac{9}{2};-\frac{7}{2}\rangle$. Furthermore, in
order to realize the dissipative part of the dynamics, we introduce an
auxiliary state, $|X\rangle\equiv|\frac{9}{2};-\frac{9}{2}\rangle$, belonging
to the ${}^{2}S_{1/2}$ manifold, as well as an electronically excited
${}^{2}P_{3/2}$ state, $|e\rangle\equiv|\frac{11}{2};-\frac{9}{2}\rangle$.

The first-stage jump operators, $\hat{j}^{(1)}_{\mathbf{i}\mathbf{j},\sigma}$
and $\hat{j}^{(2)}_{\mathbf{i}\mathbf{j},\sigma}$, can be implemented using
Raman-assisted hopping, as illustrated in figure 1(b). In such a process,
transitions between two quantum states are induced using two laser beams,
which are detuned from some excited state (here the ${}^{2}P_{3/2}$ state).
For example, the Raman beams $\Omega_{r1}$ and $\Omega_{r3}$ realize the
Raman-assisted hopping between the $|\uparrow;\uparrow\rangle$ and
$|0;\uparrow\downarrow\rangle$ states. If the Raman beams,
$\Omega_{r1,...,r3}$, are not phase-locked such hopping processes are
dissipative. Since the Raman-assisted hopping takes place directly between the
initial and final states of the jump operators, the related dissipative
processes are bidirectional. Therefore, we need to avoid populating the
$|\uparrow;\downarrow\rangle$ state at this stage. Otherwise, the jump
operators would also lead from the $|\uparrow;\downarrow\rangle$ state back to
the $|\uparrow;\uparrow\rangle$ and $|\downarrow;\downarrow\rangle$ states.
The on-site interaction energy $U$ (which we assume to be on the order of a
few kHz) can be used for this purpose. We note that this step can be
implemented in a coherent way as well, however, the required phase-locking of
the lasers would introduce an additional complication into the experimental
setup.

In order to implement the second-stage jump operators,
$\hat{j}^{(3)}_{\mathbf{i}\mathbf{j}}$, in a one-directional fashion, we use
Raman-assisted hopping (lasers $\Omega_{r4}$ and $\Omega_{r5}$) from the
$|\uparrow\downarrow;0\rangle$ state to the $|\uparrow;X\rangle$ configuration
with an auxiliary $|X\rangle$ state. This $|X\rangle$ state is then pumped
(laser $\Omega$) to an excited $|e\rangle$ state, which can decay to the
$|\downarrow\rangle$ state completing the process, as illustrated in figure
1(c). The $|e\rangle$ state cannot decay to the $|\uparrow\rangle$ state
because of selection rules on the $F$ quantum number. The Zeeman splitting,
$\Delta_{X\downarrow}$, and the energy, $U+\Delta_{X\downarrow}$,
differentiate among the three states from the lower band in figure 1(c). In
order to resolve between these three states it is sufficient to use selection
rules. To resolve between the $|\uparrow;X\rangle$ and $|\uparrow X;0\rangle$
states as the final states of the Raman-assisted hopping process we need
nonzero on-site interaction between the $|\uparrow\rangle$ and $|X\rangle$
states. The typical values of the background scattering lengths,
$a=105\,a_{0}$, in units of the Bohr radius $a_{0}$ Bohn  _et al._ (1999);
Tiecke (2011), should be sufficient for this purpose. Spontaneous emission
from the $|e\rangle$ state ensures that the resulting dissipative processes
are unidirectional and take place from the $|\uparrow\downarrow;0\rangle$ to
the $|\uparrow;\downarrow\rangle$ state with AF ordering.

In total our dissipative preparation scheme requires six lasers
($\Omega_{r1}$, …, $\Omega_{r5}$, and $\Omega$) in order to induce the
dissipative transitions, in addition to the lasers used to trap the atoms and
prepare them in the $|\uparrow\rangle$ and $|\downarrow\rangle$ states. Due to
the $s$-wave character of the ${}^{2}S_{1/2}$ manifold, the optical lattice
potential for the $|\uparrow\rangle$, $|\downarrow\rangle$, and $|X\rangle$
states will be the same, provided the optical-trapping lasers are sufficiently
far-detuned from the ${}^{2}P_{3/2}$ states (cf. the discussion in Ref. Liu
_et al._ (2004), section II.B).

The values of $\gamma_{1}$ and $\gamma_{2}$ are unrelated to the nearest-
neighbor hopping integral $t$, however, all three of them are proportional to
certain integrals involving two Wannier functions on the nearest-neighboring
sites. As a safe estimate we consider $t$, $\gamma_{1}$, and $\gamma_{2}$ to
be on the same order of magnitude.

Please note that the presented dissipative-preparation scheme is general and
other states can be used as $|\uparrow\rangle$, $|\downarrow\rangle$,
$|X\rangle$, and $|e\rangle$. In particular, it should be possible to use the
states of 40K, for which the Feschbach resonance is already known, i.e.
$|\downarrow\rangle=|\frac{9}{2};-\frac{9}{2}\rangle$,
$|\uparrow\rangle=|\frac{9}{2};-\frac{7}{2}\rangle$ or
$|\downarrow\rangle=|\frac{9}{2};-\frac{9}{2}\rangle$,
$|\uparrow\rangle=|\frac{9}{2};-\frac{5}{2}\rangle$. In the first case one
could use $|X\rangle=|\frac{7}{2},-\frac{7}{2}\rangle$ and
$|e\rangle=|\frac{11}{2},-\frac{11}{2}\rangle$, which requires a two-photon
process from $|X\rangle$ to $|e\rangle$. This could be realized with Raman
beams detuned from the ${}^{2}D_{5/2}$ state or directly as a two-photon
excitation. In the second case one could use
$|X\rangle=|\frac{9}{2},-\frac{7}{2}\rangle$ and
$|e\rangle=|\frac{9}{2},-\frac{9}{2}\rangle$ ($|e\rangle$ is from the
${}^{2}P_{3/2}$ manifold), which would require resolving the
$|X\rangle\leftrightarrow|e\rangle$ transition from
$|\uparrow\rangle\leftrightarrow|\frac{9}{2},-\frac{7}{2}\rangle$.

##  IV Time evolution and steady-state properties

In order to reveal the properties of the system we use two complementary
techniques: the Monte Carlo wave function (MCWF) technique Dalibard  _et al._
(1992); Mølmer  _et al._ (1993); Carmichael (1993) on a $3\times 3$ lattice
and the variational method Weimer (2015b, a) in the thermodynamic limit. While
the latter was originally formulated for bosons, here we extend it to
fermionic systems. In both methods we start from the Jordan-Wigner
transformation in two spatial dimensions Derzhko (2001); Fradkin (1989); Wang
(1991); Azzouz (1993). The related Jordan-Wigner strings restrict the
applicability of our variational scheme to the situation with one particle per
site (half-filling). Experimentally, this is the most interesting regime as it
corresponds to the maximal Néel temperature.

###  IV.1 Monte Carlo wave function

We study the dynamics of the driven-dissipative system governed by equation
(2) using the MCWF method implemented in the QuTiP numerical library Johansson
_et al._ (2012, 2013). We consider only the nearest-neighbor hopping
$t_{\mathbf{i}\mathbf{j}}\equiv-t$ (we use $t$ as the unit of energy
hereafter) and study the time evolution of the half-filled $3\times 3$ lattice
as a function of the parameters $\gamma_{1}$, $\gamma_{2}$, and $U$. Since the
lattice dimensions are odd numbers, we use the antiperiodic boundary
conditions. This is required in the presence of the AF ordering because a
particle hopping to its nearest neighbor across a boundary does not change the
sublattice index, whereas for the same process within the boundaries such a
change occurs. Therefore, when the hopping process takes place across the
boundary, we introduce an additional spin-flip
($\hat{c}^{\dagger}_{\mathbf{i}\sigma}\hat{c}_{\mathbf{j}\overline{\sigma}}$),
which ensures that hopping processes within and across the boundaries are
equivalent. We also introduce analogous corrections to the jump operators in
equations (3)-(5).

The initial states for the MCWF realizations were chosen with randomly-
positioned spin-up or spin-down particles (also allowing for double
occupancies), however, the steady-state properties were found to be
independent on the initial conditions.

Figure 2: Formation of the steady state. (a) Time-evolution of the total
number of spin-up, $N_{\uparrow}$, and spin-down, $N_{\downarrow}$, atoms per
site, the double occupancy probability, $D$, and the spin-structure factor,
$I_{\rm AF}$, for $U=100$, $\gamma_{1}=1$, and $\gamma_{2}=2$. (b) The von
Neumann entropy per particle, $S$, for selected values of
$(U,\gamma_{1},\gamma_{2})$, as labelled in the graph. The results have been
averaged from 256 Monte Carlo realizations. The grey areas around selected
curves correspond to one-sigma confidence interval of the results. The results
($I_{\rm AF}$ and $S$) for the Néel phase are represented by flat dashed
lines.

The time evolution of the system’s properties is shown in figure 2. One can
see that the steady state is reached for $\tau\times t\approx 50-100$, which
corresponds to $\tau=1-2$ seconds for $t=50$ Hz. Even for a large value of the
on-site interaction, $U=100$, a small number of double occupancies is still
present in the system, $D\approx 0.04$, see figure 2(a). These states are
involved in the dissipation processes as an intermediate step towards
preparation of the AF ordered phase, cf. figure 1(a). Non-zero double
occupancies can lead to inelastic losses of atoms (Esslinger, 2010), which
however are more problematic for the attractive Regal  _et al._ (2004), than
for the repulsive Jordens  _et al._ (2008) potassium gas. In the latter case
the inelastic decay time for atoms on doubly occupied sites was reported
Jordens  _et al._ (2008) to exceed $850$ ms. Consequently, such inelastic
losses should not constitute a limitation of our scheme.

To quantify the AF ordering of the system, we evaluate the spin-structure
factor, as defined by

$I(\mathbf{Q})\equiv\frac{4}{N^{2}}\sum_{\mathbf{i},\mathbf{j}}e^{i(\mathbf{R_{i}}-\mathbf{R_{j}})\mathbf{Q}}\langle\sigma^{z}_{\mathbf{i}}\sigma^{z}_{\mathbf{j}}\rangle.$
(6)

Here, $\mathbf{Q}$ is the ordering vector, which for the case of the AF phase
is equal to $\mathbf{Q}_{\rm AF}=[\pi,\pi]$, $\mathbf{R_{i}}$ and
$\mathbf{R_{j}}$ are the lattice site vectors, $N$ is the number of particles
($N=9$ in our system), and
$\sigma_{\mathbf{i}(\mathbf{j})}^{z}=\pm\frac{1}{2}$ is the $z$ component of
the particles’ spin. In the case of the AF ordering, the relevant spin-
structure factor is given by $I_{\rm AF}\equiv I(\mathbf{Q}_{\rm AF})$. For
the steady-state it can be as large as $I_{\rm AF}\approx 0.8$ (cf. also
figure 3), quite close to the fully polarized Néel state, for which $I_{\rm
AF}=1$. Note that the Néel state is an ‘ideal’ AF state, to which the
dissipative processes would drive the system in the limit of
$\gamma_{1},\gamma_{2},U\gg t$. This is because the jump operators suppress
intersite coherence in the system. However, the fact that fermions remain
mobile in the optical lattice during the entire dissipative preparation stage
should help with retaining coherence after the dissipation channels are
switched off.

Figure 2(b) illustrates the decrease of the system’s entropy per particle,
$S\equiv S_{\text{tot}}/N=-1/N\,{\rm Tr}\rho\log\rho$, with time. Due to the
large size of the system’s density matrix $\rho$, in our numerical
calculations we use the equivalent formula, $S_{\text{tot}}=-{\rm Tr}A\log A$.
Here $A_{ij}=\langle\psi_{i}|\psi_{j}\rangle$ is the matrix of overlaps of
wave functions obtained from single realizations of the Monte Carlo algorithm.
The resulting steady-state entropy per particle can be as low as $S\approx
0.1-0.3$ (cf. also figure 3). For the Néel phase the entropy per particle is
equal to $\frac{1}{9}\log 2~{}\approx~{}0.077$ due to the two-fold degeneracy
corresponding to flipping of all spins. The relaxation time (usually below one
second) increases with the increasing on-site interaction, $U$. For larger
systems the relaxation times can be longer, due to possible formation of
domains, as it is the case for coherent preparation strategies.

A slightly different number of spin-up, $N_{\uparrow}$, and spin-down,
$N_{\downarrow}$, atoms in the steady state is related to breaking of the
$SU(2)$ symmetry by the jump operators $\hat{j}^{(3)}_{\mathbf{i}\mathbf{j}}$.
As a result, the spin-down atom becomes more ‘mobile’. Additionally, in the
$3\times 3$ lattice the number of sites is odd. Therefore, inherently in the
steady state there is a spin-direction imbalance with
$N_{\uparrow}>N_{\downarrow}$. For a larger system, as well as for a system
with even number of sites we would have $N_{\uparrow}\approx N_{\downarrow}$.
This, however, does not preclude the formation of the AF order.

Figure 3:  Steady-state properties: spin-structure factor, $I_{\rm AF}$, and
entropy per particle, $S$, as a function of magnitudes of jump operators, (a)
$\gamma_{1}$ and (b) $\gamma_{2}$, and (c) the on-site repulsion energy, $U$.
The error bars of $I_{\rm AF}$ are calculated as the standard error for the
results of single Monte-Carlo realizations, whereas those of $S$ are obtained
using the standard deviation of the Monte-Carlo results averaged over time.

Figure 3 shows the steady-state properties: the entropy per particle and spin-
structure factor as a function of the parameters. In the employed range of
parameters these features turn out to be strongly dependent on $\gamma_{2}$
and $U$, cf. figures 3(b), (c), however, only weakly dependent on
$\gamma_{1}$, cf. figure 3(a). Furthermore, while $I_{\rm AF}$ grows
substantially with increasing $U$ and $\gamma_{2}$, for $\gamma_{1}$ a
saturation effect is observed and increasing the magnitude above
$\gamma_{1}\approx 0.5$ does not improve the efficiency of the scheme
significantly. These observations can be qualitatively understood from figure
1(a). Namely, when the system is close to the AF phase most of the nearest-
neighbor configurations are of the $|\uparrow;\downarrow\rangle$ type. The
processes that drive the system away from the ordered state are related to
coherent hopping from the $|\uparrow;\downarrow\rangle$ state to the
$|\uparrow\downarrow;0\rangle$ state. In the large-$U$ limit, the timescale of
such processes is given by $4t^{2}/U$. Therefore, increasing $U$ reduces the
contribution of the processes that destroy AF ordering. Increase of $I_{\rm
AF}$ with $\gamma_{2}$ is expected, as the related dissipative processes drive
the system directly into the AF ordered state. The saturation effect for
$\gamma_{1}$ can be due to the bidirectional character of the related
dissipative processes. For sufficiently large $\gamma_{1}$, the value of the
spin-structure factor is determined by an interplay between the dissipative
processes related to $\gamma_{2}$ and coherent hopping processes with a time
scale governed by $4t^{2}/U$.

While the values of $U$ required for an efficient preparation of the AF order
are quite large, they are within experimental reach, e.g. $U/t=180$ was
reported in Ref. Jordens  _et al._ (2008). Increasing $U$ even further might
lead to appearance of non-standard terms on top of the Fermi-Hubbard model
Dutta  _et al._ (2015).

###  IV.2 Variational scheme

In order to describe the steady-state properties in the thermodynamic limit,
we use a recently introduced variational principle Weimer (2015b, a). In this
method, one has to minimize a suitable variational norm 111cf. Refs. Weimer
(2015b, a) for details including the discussion of the choice of the
variational norm. of the master equation (2)

$||\dot{\rho}||=||-i[H,\rho]+\mathcal{D}(\rho)||.$ (7)

Here $\mathcal{D}$ is the dissipative part as given by equation (2). The
variational method was originally formulated for bosonic systems where a local
ansatz on the density matrix can be used. For fermionic systems such a
procedure is not possible directly. However, the Jordan-Wigner transformation
Derzhko (2001); Fradkin (1989); Wang (1991); Azzouz (1993) can be used to map
the fermionic creation and annihilation operators to spin operators (see
appendix A for details). The resulting system of spin-$1/2$ particles is
equivalent to a system of hard-core bosons, however, with the fermionic
statistics included. Namely, the equivalence with the starting fermionic
system is ensured by long-range interaction terms – the so-called Jordan-
Wigner strings – manifesting the anticommutation relations of the original
fermionic particles. To apply the variational method we first need to
approximate the ‘non-local’ parts of the Jordan-Wigner strings and, thereby,
transform the master equation (2) to the form with at most two-site
interactions (see appendix B for details of this procedure).

We next consider variational states of the product-state type,
$\rho=\rho_{p}=\prod_{\mathbf{i}}\rho_{\mathbf{i}}$ and minimize the upper
bound of the norm

$||\dot{\rho}||\leq\sum_{\langle\mathbf{i}\mathbf{j}\rangle}{\rm
Tr}|\dot{\rho}_{\mathbf{i}\mathbf{j}}|\to\min,$ (8)

where the reduced two-site operators are defined as
${\dot{\rho}_{\mathbf{i}\mathbf{j}}={\rm
Tr}_{\not{\mathbf{i}}\not{\mathbf{j}}}{\,\dot{\rho}}}$. It is sufficient to
minimize the norm $||\dot{\rho}_{\mathbf{i}\mathbf{j}}||$ of a single bond,
which for the case of an AF order can be expressed as

$\displaystyle||\dot{\rho}_{\mathbf{i}\mathbf{j}}||$ $\displaystyle\equiv$
$\displaystyle||\dot{\rho}_{AB}||={\rm Tr}|\dot{\rho}_{AB}|,$ (9)
$\displaystyle\dot{\rho}_{AB}$ $\displaystyle=$
$\displaystyle-i[H_{AB},\rho_{AB}]+\mathcal{D}_{AB}(\rho_{AB})+$
$\displaystyle\sum_{A^{\prime}}{\rm
Tr}_{A^{\prime}}\left\\{-i[H_{BA^{\prime}},\rho_{ABA^{\prime}}]+\mathcal{D}_{BA^{\prime}}(\rho_{ABA^{\prime}})\right\\}+$
$\displaystyle\sum_{B^{\prime}}{\rm
Tr}_{B^{\prime}}\left\\{-i[H_{B^{\prime}A},\rho_{B^{\prime}AB}]+\mathcal{D}_{B^{\prime}A}(\rho_{B^{\prime}AB})\right\\}.$

Here A and B label the two sublattices and, e.g.,
$\rho_{AB}\equiv\rho_{A}\otimes\rho_{B}$,
$\rho_{ABA^{\prime}}\equiv\rho_{A}\otimes\rho_{B}\otimes\rho_{A^{\prime}}$,
while $\mathcal{D}_{AB}$ gives the dissipative part with the jump operators
acting on the sites A and B. The first two terms of equation (IV.2) correspond
to an exact treatment of a single bond, which already goes beyond the mean-
field description, whereas the next ones describe interaction with the
surrounding sites treated on the mean-field level, as visualized by the dashed
lines in figure 4(a).

Figure 4: Results in the thermodynamic limit. (a) Single link and the
surrounding sites in the presence of AF order. (b) The steady-state magnitude
of the AF spin-structure factor, as a function of the magnitude of the jump
operators. The variational results (red solid line) are compared to the MCWF
method (blue dashed line).

In figure 4(b) we compare the spin-structure factor obtained from our
variational scheme and the MCWF method as a function of magnitude of the jump
operators (which we set equal here without the loss of generality). According
to the results of both methods, the system exhibits substantial ordering (e.g.
with spin-structure factor larger than $0.5$) when
$\gamma_{1}=\gamma_{2}\gtrsim 1$. Although the variational method contains
terms going beyond mean field, its results for AF ordering do not depend on
the value of $U$, as opposed to the exact approach. This happens due to
restriction of the density matrices to the form of product states. Moreover,
the variational method overestimates the AF ordering due to the absence of
fluctuations in our variational manifold. Therefore, the employed approaches
are complementary to each other, and both indicate the formation of an AF
order of substantial magnitude in the steady state.

##  V Conclusions

In this paper, we proposed a scheme for dissipative preparation of
antiferromagnetic order in ultracold fermions trapped in an optical lattice.
We demonstrated that by using a combination of two dissipative processes based
on Raman-assisted hopping it is possible to engineer the dissipative dynamics
in such a way that the AF phase emerges as its dark state. By using a
combination of an exact and variational approaches, we observed the formation
of a strong AF order on the timescales achievable in present-day experiments.

We note that the technique presented here can be readily implemented in the
setups already used to search for the AF order Greif  _et al._ (2015), and
thereby paves the way to an experimental realization of the AF phase in the
Fermi-Hubbard model. While we exemplified the approach using the atomic level
structure of 40K Jordens  _et al._ (2008), the method is general and can be
also applied to other fermionic atoms currently available in laboratory, such
as 6Li Hart  _et al._ (2015), Er Aikawa  _et al._ (2014), Dy Lu  _et al._
(2012), Yb Taie  _et al._ (2010), and Cr Naylor  _et al._ (2015). After
preparation of the AF phase with low entropy it should be possible to explore
the phase diagram of the Hubbard model, including the pseudogap regime, by
coherently removing a fraction of the atoms from the trap thereby introducing
hole carriers into the system. Finally, extending these ideas to single-site
addressable lattices as offered by the fermionic quantum gas microscopes Cheuk
_et al._ (2015); Haller  _et al._ (2015); Parsons  _et al._ (2015); Edge  _et
al._ (2015), opens the door to the preparation of more sophisticated many-
particles states.

###### Acknowledgements.

We acknowledge stimulating discussions with Ken Brown, Tommaso Calarco, Andrew
Daley, Suzanne McEndoo, Tobias Osborne, Cindy Regal, Luis Santos, Michał
Tomza, and Martin Zwierlein. The work was supported by the People Programme
(Marie Curie Actions) of the European Union’s Seventh Framework Programme
(FP7/2007-2013) under REA grant agreement no [291734], by the Volkswagen
Foundation, and by DFG within SFB 1227 (DQ-mat).

##  Appendix A Two-dimensional Jordan-Wigner transformation

The Jordan-Wigner transformation between fermionic operators
$\hat{c}^{\dagger}_{\mathbf{i}\sigma}$, $\hat{c}_{\mathbf{i}\sigma}$ and spin
operators ($\hat{\sigma}^{+},\hat{\sigma}^{-},\hat{\sigma}^{z}$) for spin-1/2
particles is defined as

$\displaystyle\hat{c}^{\dagger}_{\mathbf{i}\sigma}$ $\displaystyle=$
$\displaystyle\hat{\sigma}_{k(\mathbf{i},\sigma)}^{+}e^{+i\pi\sum_{l<k(\mathbf{i},\sigma)}\hat{n}_{l}},$
(11) $\displaystyle\hat{c}_{\mathbf{i}\sigma}$ $\displaystyle=$
$\displaystyle\hat{\sigma}_{k(\mathbf{i},\sigma)}^{-}e^{-i\pi\sum_{l<k(\mathbf{i},\sigma)}\hat{n}_{l}},$
(12)
$\displaystyle\hat{c}^{\dagger}_{\mathbf{i}\sigma}\hat{c}_{\mathbf{i}\sigma}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left(\hat{\sigma}^{z}_{k(\mathbf{i},\sigma)}+1\right)\equiv\hat{n}_{k(\mathbf{i},\sigma)},$
(13)

where the function $k(\mathbf{i},\sigma)$ enumerates the fermionic particles
for a given lattice site, $\mathbf{i}=(i_{x},i_{y})$, and particle’s spin,
$\sigma=\uparrow,\downarrow$. In this manner, fermions are positioned on a
chain, where the position along the chain is given by $k(\mathbf{i},\sigma)$.
The factor $e^{\pm i\pi\sum_{l<k(\mathbf{i},\sigma)}\hat{n}_{l}}$, where the
summation runs over all fermions ‘before’ the one at site $\mathbf{i}$ with
spin $\sigma$, evaluates to $+1(-1)$ for even (odd) number of fermions
‘before’ the given one. Thereby, the anticommutation rules for the fermionic
operators $\hat{c}^{\dagger}_{\mathbf{i}\sigma}$, $\hat{c}_{\mathbf{i}\sigma}$
are fulfilled. We can also use the relation $e^{\pm
i\pi\sum_{l<k(\mathbf{i},\sigma)}\hat{n}_{l}}=\Pi_{l<k(\mathbf{i},\sigma)}(-\hat{\sigma}^{z}_{l})$,
where $-\hat{\sigma}^{z}_{l}$ evaluates to $-1$ when the $l$-th spin is up
(correspondingly, when there is a fermionic particle in the mode $l$) and to 1
in the opposite case.

For a one-dimensional system the function $k(\mathbf{i},\sigma)$ can be chosen
simply as $k(\mathbf{i},\sigma)=2|\mathbf{i}|+\delta_{\sigma,\downarrow}$. In
two dimensions, however, there are a few possibilities to perform the Jordan-
Wigner transformation Derzhko (2001); Fradkin (1989); Wang (1991); Azzouz
(1993) (see Ref. Derzhko (2001) for a review). Here, we consider an
$N_{x}\times N_{y}$ system with lattice site coordinates
$i_{x(y)}\in\\{1,2,...,N_{x(y)}\\}$ and use the following function to
enumerate particles (with $m\in\mathbb{N}$)

$k(\mathbf{i},\sigma)=\begin{cases}2\left[N_{x}(i_{y}-1)+i_{x}-1\right]+\delta_{\sigma,\downarrow},&i_{y}=2m+1,\\\
2\left[N_{x}(i_{y}-1)+N_{x}-i_{x}\right]+\delta_{\sigma,\downarrow},&i_{y}=2m.\end{cases}$
(14)

The chain of fermions and the related Jordan-Wigner string goes from left to
right in the first row of the system, then in the second row goes to the left
and forms a zig-zag (cf. figure 5(a)).

Figure 5: Illustration of the Jordan-Wigner transformation for a $N_{x}\times
N_{y}$ system. (a) Jordan-Wigner string (in green) from site $(1,1)$ to site
$\mathbf{i}$. (b) and (c) illustration of the hopping terms and the related
Jordan-Wigner strings along (b) and across (c) the ’rows’ formed by the
strings. (d) The hopping term as in (c) with approximation of the ’non-local’
(dashed) part of the string. The remaining components of the string act only
on the sites A and B (as visualized by the solid green arrow) and the hopping
term has the same form as in (b).

##  Appendix B Variational method for fermions

By performing the Jordan-Wigner transformation as described in appendix A, we
can transform all terms of the master equation (2), which we exemplify here by
considering the hopping,
$\hat{c}^{\dagger}_{\mathbf{i}\sigma}\hat{c}_{\mathbf{j}\sigma}$. Using the
transformation (11)-(13) the hopping is expressed as

$\hat{c}^{\dagger}_{\mathbf{i}\sigma}\hat{c}_{\mathbf{j}\sigma}=\hat{\sigma}_{k(\mathbf{i},\sigma)}^{+}\left[\Pi_{k(\mathbf{i},\sigma)<l<k(\mathbf{j},\sigma)}(-\hat{\sigma}^{z}_{l})\right]\hat{\sigma}_{k(\mathbf{j},\sigma)}^{-}.$
(15)

If the hopping takes place between two sites in the same row, the Jordan-
Wigner string is local (as illustrated in figure 5(b)). For example, if
$\mathbf{j}=\mathbf{i}+(1,0)$ and $k(\mathbf{i},\sigma)<k(\mathbf{j},\sigma)$,
we have

$\displaystyle\hat{c}^{\dagger}_{\mathbf{i}\uparrow}\hat{c}_{\mathbf{j}\uparrow}$
$\displaystyle=$
$\displaystyle\hat{\sigma}_{k(\mathbf{i},\uparrow)}^{+}\left[-\hat{\sigma}^{z}_{k(\mathbf{i},\downarrow)}\right]\hat{\sigma}_{k(\mathbf{j},\uparrow)}^{-},$
(16)
$\displaystyle\hat{c}^{\dagger}_{\mathbf{i}\downarrow}\hat{c}_{\mathbf{j}\downarrow}$
$\displaystyle=$
$\displaystyle\hat{\sigma}_{k(\mathbf{i},\downarrow)}^{+}\left[-\hat{\sigma}^{z}_{k(\mathbf{j},\uparrow)}\right]\hat{\sigma}_{k(\mathbf{j},\downarrow)}^{-}.$
(17)

If the hopping takes place between two sites in different rows, the Jordan-
Wigner string includes also sites between the $\mathbf{i},\mathbf{j}$ pair and
the (left or right) edge of the system (as illustrated in figure 5(c)). Note
that the number of these ‘non-local’ sites is always even due to our choice of
the enumerating function $k(\mathbf{i},\sigma)$. Consequently, at half-
filling, i.e. with one particle per site on average, there is an even number
of particles, $2N_{1}$, along the ‘non-local’ part of the Jordan-Wigner
string. In such case, the string evaluates to the factor $(-1)^{2N_{1}}=1$. In
our variational method we use this value as an approximation for the ‘non-
local’ part of the string, whereas the local terms are retained. Such a
procedure is similar to the ‘mean-field’ treatment of the strings applied,
e.g. in Refs. Derzhko (2001); Derzhko  _et al._ (2003); Wang (1991). As a
result, the hopping between rows is expressed in the same way as in equations
(16)-(17), which is visualized in figure 5(d).

For the jump operators we use the same approximation. Thereby, the master
equation (2) acquires the form with at most two-site interaction terms (site
here refers to the original fermionic site), for which the variational method
can be readily applied. Note that in the MCWF calculations the ‘non-local’
part of the strings needs to be preserved.

Explicitly, the contribution from the hopping term to the norm of a single
bond, $||{\dot{\rho}_{AB}||=||{\rm Tr}_{\not{A}\not{B}}{\,\dot{\rho}}}||$, on
the example of the first term in equation (IV.2), is given by

$\displaystyle-i[H_{AB},\rho_{AB}]$ $\displaystyle=$ $\displaystyle-
it[\hat{\sigma}_{k(A,\uparrow)}^{+}\hat{\sigma}^{z}_{k(A,\downarrow)}\hat{\sigma}_{k(B,\uparrow)}^{-}$
$\displaystyle+\hat{\sigma}_{k(A,\downarrow)}^{+}\hat{\sigma}^{z}_{k(B,\uparrow)}\hat{\sigma}_{k(B,\downarrow)}^{-}+{\rm
h.c.},\rho_{AB}],$

where we assumed $k(A,\sigma)<k(B,\sigma)$. The treatment of other terms in
the norm is analogous and, similarly, results in the appearance of parity
factors ($-\hat{\sigma}^{z}_{k(A,\sigma)}$ and
$-\hat{\sigma}^{z}_{k(B,\sigma)}$), which originate from the fermionic
anticommutation rules.

Finally, let us note that the results of our variational method do not depend
on the system size, $N_{x}\times N_{y}$ (however, they make sense only for
$N_{x},N_{y}\geq 4$) and hence correspond to the thermodynamic limit.

## References

  * Lewenstein  _et al._ (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Adv. Phys. 56, 243 (2007).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * Chin  _et al._ (2006) J. K. Chin, D. E. Miller, Y. Liu, C. Stan, W. Setiawan, C. Sanner, K. Xu, and W. Ketterle, Nature 443, 961 (2006).
  * Schneider  _et al._ (2008) U. Schneider, L. Hackermüller, S. Will, T. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science 322, 1520 (2008).
  * Jordens  _et al._ (2008) R. Jordens, N. Strohmaier, K. Gunter, H. Moritz, and T. Esslinger, Nature (London) 455, 204 (2008).
  * Greif  _et al._ (2013) D. Greif, T. Uehlinger, G. Jotzu, L. Tarruell, and T. Esslinger, Science 340, 1307 (2013).
  * Krauser  _et al._ (2012) J. S. Krauser, J. Heinze, N. Flaschner, S. Gotze, O. Jurgensen, D.-S. Luhmann, C. Becker, and K. Sengstock, Nat. Phys. 8, 813 (2012).
  * Jotzu  _et al._ (2014) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature 515, 237 (2014).
  * Uehlinger  _et al._ (2013) T. Uehlinger, G. Jotzu, M. Messer, D. Greif, W. Hofstetter, U. Bissbort, and T. Esslinger, Phys. Rev. Lett. 111, 185307 (2013).
  * Esslinger (2010) T. Esslinger, Annu. Rev. Condens. Matter. Phys. 1, 129 (2010).
  * Le Hur and Rice (2009) K. Le Hur and T. M. Rice, Ann. Phys. 324, 1452 (2009).
  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Hart  _et al._ (2015) R. A. Hart, P. M. Duarte, T.-L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, D. A. Huse, and R. G. Hulet, Nature (London) 519, 211 (2015).
  * Greif  _et al._ (2015) D. Greif, G. Jotzu, M. Messer, R. Desbuquois, and T. Esslinger, Phys. Rev. Lett. 115, 260401 (2015).
  * Ho and Zhou (2009) T.-L. Ho and Q. Zhou, Proc. Nat. Acad. Sci. 106, 6916 (2009).
  * Bernier  _et al._ (2009) J.-S. Bernier, C. Kollath, A. Georges, L. De Leo, F. Gerbier, C. Salomon, and M. Köhl, Phys. Rev. A 79, 061601 (2009).
  * Heidrich-Meisner  _et al._ (2009) F. Heidrich-Meisner, S. R. Manmana, M. Rigol, A. Muramatsu, A. E. Feiguin, and E. Dagotto, Phys. Rev. A 80, 041603 (2009).
  * Diehl  _et al._ (2010a) S. Diehl, W. Yi, A. J. Daley, and P. Zoller, Phys. Rev. Lett. 105, 227001 (2010a).
  * Williams  _et al._ (2010) J. R. Williams, J. H. Huckans, R. W. Stites, E. L. Hazlett, and K. M. O’Hara, Phys. Rev. A 82, 011610 (2010).
  * Medley  _et al._ (2011) P. Medley, D. M. Weld, H. Miyake, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 106, 195301 (2011).
  * Lee  _et al._ (2011) T. E. Lee, H. Häffner, and M. C. Cross, Phys. Rev. A 84, 031402 (2011).
  * Colomé-Tatché  _et al._ (2011) M. Colomé-Tatché, C. Klempt, L. Santos, and T. Vekua, New J. Phys. 13, 113021 (2011).
  * Lemeshko  _et al._ (2012) M. Lemeshko, R. V. Krems, and H. Weimer, Phys. Rev. Lett. 109, 035301 (2012).
  * Yi  _et al._ (2012) W. Yi, S. Diehl, A. J. Daley, and P. Zoller, New J. Phys. 14, 055002 (2012).
  * Carr and Saffman (2013) A. W. Carr and M. Saffman, Phys. Rev. Lett. 111, 033607 (2013).
  * Murmann  _et al._ (2015) S. Murmann, F. Deuretzbacher, G. Zürn, J. Bjerlin, S. M. Reimann, L. Santos, T. Lompe, and S. Jochim, Phys. Rev. Lett. 115, 215301 (2015).
  * Mazzucchi  _et al._ (2016) G. Mazzucchi, S. F. Caballero-Benitez, and I. B. Mekhov, Sci. Rep. 6, 31196 (2016).
  * Diehl  _et al._ (2008) S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, Nat. Phys. 4, 878 (2008).
  * Kraus  _et al._ (2008) B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307 (2008).
  * Witthaut  _et al._ (2008) D. Witthaut, F. Trimborn, and S. Wimberger, Phys. Rev. Lett. 101, 200402 (2008).
  * Verstraete  _et al._ (2009) F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Nat. Phys. 5, 633 (2009).
  * Weimer  _et al._ (2010) H. Weimer, M. Muller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nat. Phys. 6, 382 (2010).
  * Diehl  _et al._ (2011) S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Nat. Phys. 7, 971 (2011).
  * Alharbi and Ficek (2010) A. F. Alharbi and Z. Ficek, Phys. Rev. A 82, 054103 (2010).
  * Watanabe and Mäkelä (2012) G. Watanabe and H. Mäkelä, Phys. Rev. A 85, 023604 (2012).
  * Kordas  _et al._ (2012) G. Kordas, S. Wimberger, and D. Witthaut, Europhys. Lett. 100, 30007 (2012).
  * Bermudez  _et al._ (2013) A. Bermudez, T. Schaetz, and M. B. Plenio, Phys. Rev. Lett. 110, 110502 (2013).
  * Rao and Mølmer (2013) D. D. B. Rao and K. Mølmer, Phys. Rev. Lett. 111, 033606 (2013).
  * Lemeshko and Weimer (2013) M. Lemeshko and H. Weimer, Nat. Commun. 4, 2230 (2013).
  * Weimer (2013) H. Weimer, Mol. Phys. 111, 1753 (2013).
  * Otterbach and Lemeshko (2014) J. Otterbach and M. Lemeshko, Phys. Rev. Lett. 113, 070401 (2014).
  * Lucia  _et al._ (2015) A. Lucia, T. S. Cubitt, S. Michalakis, and D. Pérez-García, Phys. Rev. A 91, 040302 (2015).
  * Budich  _et al._ (2015) J. C. Budich, P. Zoller, and S. Diehl, Phys. Rev. A 91, 042117 (2015).
  * Weimer (2015a) H. Weimer, Phys. Rev. A 91, 063401 (2015a).
  * Barreiro  _et al._ (2010) J. T. Barreiro, P. Schindler, O. Guhne, T. Monz, M. Chwalla, C. F. Roos, M. Hennrich, and R. Blatt, Nat. Phys. 6, 943 (2010).
  * Gommers  _et al._ (2005) R. Gommers, S. Bergamini, and F. Renzoni, Phys. Rev. Lett. 95, 073003 (2005).
  * Barreiro  _et al._ (2011) J. T. Barreiro, M. Muller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, Nature (London) 470, 486 (2011).
  * Krauter  _et al._ (2011) H. Krauter, C. A. Muschik, K. Jensen, W. Wasilewski, J. M. Petersen, J. I. Cirac, and E. S. Polzik, Phys. Rev. Lett. 107, 080503 (2011).
  * Schindler  _et al._ (2013) P. Schindler, M. Muller, D. Nigg, J. T. Barreiro, E. A. Martinez, M. Hennrich, T. Monz, S. Diehl, P. Zoller, and R. Blatt, Nat. Phys. 9, 361 (2013).
  * Lin  _et al._ (2013) Y. Lin, J. P. Gaebler, F. Reiter, T. R. Tan, R. Bowler, A. S. Sorensen, D. Leibfried, and D. J. Wineland, Nature 504, 415 (2013).
  * Shankar  _et al._ (2013) S. Shankar, M. Hatridge, Z. Leghtas, K. M. Sliwa, A. Narla, U. Vool, S. M. Girvin, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Nature 504, 419 (2013).
  * Diehl  _et al._ (2010b) S. Diehl, A. Tomadin, A. Micheli, R. Fazio, and P. Zoller, Phys. Rev. Lett. 105, 015702 (2010b).
  * Tomadin  _et al._ (2011) A. Tomadin, S. Diehl, and P. Zoller, Phys. Rev. A 83, 013611 (2011).
  * Le Boité  _et al._ (2013) A. Le Boité, G. Orso, and C. Ciuti, Phys. Rev. Lett. 110, 233601 (2013).
  * Jin  _et al._ (2013) J. Jin, D. Rossini, R. Fazio, M. Leib, and M. J. Hartmann, Phys. Rev. Lett. 110, 163605 (2013).
  * Sieberer  _et al._ (2013) L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl, Phys. Rev. Lett. 110, 195301 (2013).
  * Lee  _et al._ (2013) T. E. Lee, S. Gopalakrishnan, and M. D. Lukin, Phys. Rev. Lett. 110, 257204 (2013).
  * Chan  _et al._ (2015) C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Phys. Rev. A 91, 051601 (2015).
  * Cheuk  _et al._ (2015) L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V. Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein, Phys. Rev. Lett. 114, 193001 (2015).
  * Haller  _et al._ (2015) E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Nat Phys 11, 738 (2015).
  * Edge  _et al._ (2015) G. J. A. Edge, R. Anderson, D. Jervis, D. C. McKay, R. Day, S. Trotzky, and J. H. Thywissen, Phys. Rev. A 92, 063406 (2015).
  * Finazzi  _et al._ (2015) S. Finazzi, A. Le Boité, F. Storme, A. Baksic, and C. Ciuti, Phys. Rev. Lett. 115, 080604 (2015).
  * Weimer (2015b) H. Weimer, Phys. Rev. Lett. 114, 040402 (2015b).
  * Overbeck and Weimer (2016) V. R. Overbeck and H. Weimer, Phys. Rev. A 93, 012106 (2016).
  * Bernier  _et al._ (2014) J.-S. Bernier, D. Poletti, and C. Kollath, Phys. Rev. B 90, 205125 (2014).
  * Prosen (2008) T. Prosen, New J. Phys. 10, 043026 (2008).
  * Dzhioev and Kosov (2015) A. A. Dzhioev and D. S. Kosov, J. Phys. A: Math. Theor. 48, 015004 (2015).
  * Dalibard  _et al._ (1992) J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
  * Mølmer  _et al._ (1993) K. Mølmer, Y. Castin, and J. Dalibard, J. Opt. Soc. Am. B 10, 524 (1993).
  * Carmichael (1993) H. J. Carmichael, Phys. Rev. Lett. 70, 2273 (1993).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Tiecke (2011) T. G. Tiecke,  (2011), www.tobiastiecke.nl/archive/PotassiumProperties.pdf.
  * Bohn  _et al._ (1999) J. L. Bohn, J. P. Burke, C. H. Greene, H. Wang, P. L. Gould, and W. C. Stwalley, Phys. Rev. A 59, 3660 (1999).
  * Liu  _et al._ (2004) W. V. Liu, F. Wilczek, and P. Zoller, Phys. Rev. A 70, 033603 (2004).
  * Derzhko (2001) O. Derzhko, J. Phys. Stud. 5, 49 (2001), arXiv:cond-mat/0101188.
  * Fradkin (1989) E. Fradkin, Phys. Rev. Lett. 63, 322 (1989).
  * Wang (1991) Y. R. Wang, Phys. Rev. B 43, 3786 (1991).
  * Azzouz (1993) M. Azzouz, Phys. Rev. B 48, 6136 (1993).
  * Johansson  _et al._ (2012) J. Johansson, P. Nation, and F. Nori, Comp. Phys. Comm. 183, 1760 (2012).
  * Johansson  _et al._ (2013) J. Johansson, P. Nation, and F. Nori, Comp. Phys. Comm. 184, 1234 (2013).
  * Regal  _et al._ (2004) C. A. Regal, M. Greiner, and D. S. Jin, Phys. Rev. Lett. 92, 040403 (2004).
  * Dutta  _et al._ (2015) O. Dutta, M. Gajda, P. Hauke, M. Lewenstein, D.-S. Lühmann, B. A. Malomed, T. Sowiński, and J. Zakrzewski, Rep. Prog. Phys. 78, 066001 (2015).
  * Note (1) Cf. Refs. Weimer (2015b, a) for details including the discussion of the choice of the variational norm.
  * Aikawa  _et al._ (2014) K. Aikawa, A. Frisch, M. Mark, S. Baier, R. Grimm, J. L. Bohn, D. S. Jin, G. M. Bruun, and F. Ferlaino, Phys. Rev. Lett. 113, 263201 (2014).
  * Lu  _et al._ (2012) M. Lu, N. Q. Burdick, and B. L. Lev, Phys. Rev. Lett. 108, 215301 (2012).
  * Taie  _et al._ (2010) S. Taie, Y. Takasu, S. Sugawa, R. Yamazaki, T. Tsujimoto, R. Murakami, and Y. Takahashi, Phys. Rev. Lett. 105, 190401 (2010).
  * Naylor  _et al._ (2015) B. Naylor, A. Reigue, E. Maréchal, O. Gorceix, B. Laburthe-Tolra, and L. Vernac, Phys. Rev. A 91, 011603 (2015).
  * Parsons  _et al._ (2015) M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Phys. Rev. Lett. 114, 213002 (2015).
  * Derzhko  _et al._ (2003) O. Derzhko, T. Verkholyak, R. Schmidt, and J. Richter, Physica A 320, 407 (2003).