a

r
X
i
v
:
0
8
1
2
.
3
7
6
7
v
1


[
c
o
n
d
-
m
a
t
.
m
e
s
-
h
a
l
l
]


1
9

D
e
c

2
0
0
8
Cold atoms in double-well optical lattices
V.I. Yukalov
1
and E.P. Yukalova
2
1
Bogolubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, Dubna 141980, Russia
2
Department of Computational Physics, Laboratory of Information Technologies,
Joint Institute for Nuclear Research, Dubna 141980, Russia
Abstract
Cold atoms, loaded into an optical lattice with double-well sites, are considered.
Pseudospin representation for an eective Hamiltonian is derived. The system in
equilibrium displays two phases, ordered and disordered. The second-order phase
transition between the phases can be driven either by temperature or by changing
the system parameters. Collective pseudospin excitations have a gap disappearing
at the phase-transition point. Dynamics of atoms is studied, when they are loaded
into the lattice in an initially nonequilibrium state. It is shown that the temporal
evolution of atoms, contrary to their equilibrium thermodynamics, cannot be de-
scribed in the mean-eld approximation, since it results in a structurally unstable
dynamical system, but a more accurate description is necessary taking account of
attenuation eects.
PACS: 03.75.Lm, 03.75.Be, 03.67.Lx
1
I Introduction
Degenerate cold gases, Bose [18] as well as Fermi [9,10], possess many interesting prop-
erties. Loading cold atoms into optical lattices yields highly controllable systems that
could be employed for a variety of applications [1115]. Recently, a novel type of optical
lattices has been experimentally realized, each site of which is formed by a double-well
potential [1621].
The phase diagram of cold bosons in double-well optical lattices has been studied
in Refs. [22,23], where the Hubbard model is used and the main attention is payed to
the peculiarities of the superuid-insulator phase transition in such lattices, as compared
to this transition in the standard single-well optical lattices [1114]. In Ref. [24], the
superuid phase is studied, aiming at nding a type of superuid with the broken time-
reversal symmetry.
The aim of the present paper is to study a dierent regime that can be realized in
the double-well optical lattices. We consider the lattices that are in the insulating state,
far from the boundary of the insulator-superuid transition. In that state, the jumps
of atoms between dierent lattice sites are suppressed, while the tunneling between the
wells of a single-site double well can be of principal importance. We show that in such
insulating double-well optical lattices there exists another phase transition, the order-
disorder phase transition. The investigation of this dierent physical regime in insulating
double-well optical lattices, far from the superuid phase-transition, but close to another,
order-disorder phase transition, is the main motivation for the present work.
The order-disorder phase transition, and the related ordered and disordered states,
can be the most clearly realized for the case of one atom per a double well. Therefore,
we consider exactly this case. An additional argument for considering the lattices with
one atom per site is that the optical lattices with small ling factors, such as one or two,
seem to be good candidates for quantum information processing [11,19].
The advantage of dealing with the case of one atom per a double well is twofold.
The most important, as is emphasized above, is that this is the setup allowing for a
clear realization of the order-disorder phase transition. Another, technical, convenience
is that in this case the pseudospin representation can be invoked. Some types of spin
Hamiltonians using cold atoms in optical lattices can be met in literature. It is straight-
forward to obtain a spin Hamiltonian for spinor condensates [25]. The Hubbard model
for a system of two-component bosons can be reduced to a pseudospin representation by
a second-order perturbation theory in the tunneling parameter [2628]. In our case of
an insulating double-well optical lattice, the pseudospin representation can be introduced
without perturbation theory, by means of an exact canonical transformation, similar to the
canonical transformations used for deriving the spin representations for superconductors
and ferromagnets [29,30].
The key point of realizing the ordered and disordered states in a double-well optical
lattice is the existence of suciently long-range atomic interactions. Such interactions of
dypolar type arise between polar molecules [31], Rydberg atoms [32], and between atoms
with large magnetic moments [33]. The existence of long-range atomic interactions is
another point making our consideration principally dierent from the earlier theoretical
works [2224] on the double-well lattices.
2
Characterizing collective atomic states in the double-well lattices, we, rst, describe
their equilibrium properties. However, we keep in mind that in experiments atoms need
to be loaded into a lattice, and their initial state after the loading may be nonequilibrium.
If so, how then atoms would relax to their equilibrium state? Another way, when atoms
can happen to be in a nonequilibrium states, is if the lattice parameters are varied after
the atoms have been loaded. If this variation is suciently fast, atoms again occur to
be in a nonequilibrium state, from which they should relax to an equilibrium one. In
order that our consideration of the insulating double-well optical lattices would be more
complete, and keeping in mind the possibility of realizing nonequilibrium states, we study
not only the equilibrium properties of such lattices, but also the relaxational dynamics of
atoms from initially nonequilibrium states to their equilibrium.
The paper is organized as follows. First, we derive an eective Hamiltonian for the
considered insulating system and show that it allows for a convenient pseudospin repre-
sentation (Sec. II). Then we study the equilibrium thermodynamics of the model, which
exhibits the existence of two phases, ordered and disordered (Sec. III). Collective excita-
tions, corresponding to pseudospin waves, are described in Sec. IV. Dynamics of atoms,
loaded in the lattice in an initially nonequilibrium state, is considered in Sec. V. Lyapunov
stability and structural stability of solutions to the equations of motion are investigated
in Sec. VI, where the principal importance of taking into account attenuation eects is
demonstrated. The main results are summarized in Sec. VII.
II Model Hamiltonian
We start with the general form of the energy Hamiltonian

H =
_

(r)H
L
(r)(r) dr +
1
2
_

(r)

(r

)(r r

)(r

)(r) drdr

, (1)
where (r) = (r, t) is a eld operator, in which the time dependence, for brevity, is
omitted. Keeping in mind the case of an insulating lattice, with the unity lling factor,
the statistics of atoms is not as important. So, atoms can be either bosons or fermions.
In the lattice Hamiltonian
H
L
(r)

2
2m
+ V (r) , (2)
the lattice potential V (r + a
i
) = V (r) is periodic over the lattice {a
i
: i = 1, 2, . . . , N
L
}
and enjoys the double-well structure around each of the lattice sites a
i
. The interaction
potential (r) = (r), in general, is a sum of the short-range interactions, whose
strength can be regulated by the Feschbach resonance techniques [2,10,15,34,35], and of a
long-range interaction of the dipolar type, such that exists between polar molecules [31],
Rydberg atoms [32], and atoms with large magnetic moments [33]. If, instead of atoms,
we consider ions, then these exists the long-range Coulomb interaction. The strength of
the interaction potential is assumed to be such that the intersite interactions, at least
between the nearest neighbors, cannot be neglected.
The eld operator can be expanded over Wannier functions
(r) =

nj
c
nj
w
n
(r a
j
) , (3)
3
where n is a band index and j enumerates the lattice sites. Then Hamiltonian (1) trans-
forms into

H =

ij

mn
E
mn
ij
c

mi
c
nj
+
1
2

{j}

{n}

n
1
n
2
n
3
n
4
j
1
j
2
j
3
j
4
c

n
1
j
1
c

n
2
j
2
c
n
3
j
3
c
n
4
j
4
, (4)
where
E
mn
ij

_
w

m
(r a
i
) H
L
(r) w
n
(r a
j
) dr (5)
and
n
1
n
2
n
3
n
4
j
1
j
2
j
3
j
4
is the corresponding matrix element of the interaction potential.
It can be shown (see Appendix A), that
E
mn
ij
=
mn
E
n
ij
(6)
is diagonal with respect to the band indices. Here,
E
n
ij
=
ij
E
n
+ (1
ij
)J
n
ij
, (7)
with
E
n
=
_
w

n
(r) H
L
(r) w
n
(r) dr , J
n
ij
=
_
w

n
(r a
ij
) H
L
(r) w
n
(r) dr .
Keeping in mind an insulating lattice implies that the intersite hopping is small, such that

J
ij
E
n

1 (i = j) . (8)
Assuming that each lattice site contains just one atom, we impose the unipolarity
conditions

n
c

nj
c
nj
= 1 , c
nj
c
nj
= 0 . (9)
Taking into account the above consideration, for Hamiltonian (4) we have

H =

nj
E
n
c

nj
c
nj
+
1
2

i=j

mnm

V
mnm

ij
c

mi
c

nj
c
m

j
c
n

i
, (10)
where
V
mnm

ij

mnm

ijji

mnn

ijij
, (11)
with the upper or lower sign for bosons or fermions, respectively.
When a single-well potential is transformed into a double-well potential, then each
spectrum level of a particle, it would possess in the former potential, splits into two lines
for the particle in a double-well potential. The spectrum splitting is connected with the
particle tunneling between the wells of a double-well potential. The splitting magnitude
depends on the characteristics of the double-well potential and can be regulated in a
wide range [36]. Thus, for a double-well lattice, one cannot limit oneself by considering
solely the lowest energy band, but at least two energy levels must be taken into account,
describing the level splitting and interwell tunneling. In what follows, we take into account
two lowest energy levels, so that the index n takes two values n = 1, 2.
4
The symmetry properties of the ground-state wave function and of that for the rst
excited state are known [36] to be dierent. Enumerating the ground state with n = 1
and the excited state with n = 2, one has
w
1
(r) = w
1
(r) , w
2
(r) = w
2
(r) . (12)
Also, both functions w
n
(r), for n = 1, 2, can be taken to be real. Because of the symmetry
property (12), the matrix elements of the type as V
1112
ij
and V
2221
ij
become zero.
For what follows, it is convenient to introduce the notation
E
0

1
2
(E
1
+ E
2
) . (13)
And let us dene the interaction matrix elements
A
ij

1
4
_
V
1111
ij
+ V
2222
ij
+ 2V
1221
ij
_
, B
ij

1
2
_
V
1111
ij
+ V
2222
ij
2V
1221
ij
_
,
C
ij

1
2
_
V
2222
ij
V
1111
ij
_
, I
ij
2V
1122
ij
. (14)
The quantity
E
2
E
1
+

j(=i)
C
ij
(15)
is the tunneling frequency characterizing the tunneling between the wells of a double-well
potential.
The convenience of dealing with the two-level case is that it allows for the introduction
of the pseudospin representation. The pseudospin operators can be dened as
S
x
j
=
1
2
_
c

1j
c
1j
c

2j
c
2j
_
, S
y
j
=
i
2
_
c

1j
c
2j
c

2j
c
1j
_
,
S
z
j
=
1
2
_
c

1j
c
2j
+ c

2j
c
1j
_
, (16)
which gives
c

1j
c
1j
=
1
2
+ S
x
j
, c

2j
c
2j
=
1
2
S
x
j
,
c

1j
c
2j
= S
z
j
iS
y
j
, c

2j
c
1j
= S
z
j
+ iS
y
j
. (17)
The physical meaning of the pseudospin operators (16) can be claried by introducing
the left, c
jL
, and the right, c
jR
, location operators
c
jL

1

2
(c
1j
+ c
2j
) , c
jR

1

2
(c
1j
c
2j
) , (18)
characterizing the left or right location of an atom in the left or right well of a double-well
potential. The pseudospin operators (16), expressed through the location operators (18),
become
S
x
j
=
1
2
_
c

jL
c
jR
+ c

jR
c
jL
_
, S
y
j
=
i
2
_
c

jL
c
jR
c

jR
c
jL
_
,
5
S
z
j
=
1
2
_
c

jL
c
jL
c

jR
c
jR
_
. (19)
This representation demonstrates that S
x
j
characterizes the tunneling intensity between
the left and right wells of a double-well potential centered at the j-site; the operator S
y
j
corresponds to the Josephson current between the wells; while S
z
j
is the displacement
operator describing imbalance between the wells.
Finally, for Hamiltonian(10), we obtain the pseudospin form

H = E
0
N +
1
2

i=j
A
ij

j
S
x
j
+

i=j
B
ij
S
x
i
S
x
j

i=j
I
ij
S
z
i
S
z
j
. (20)
The rst two terms here are not of the operator type, hence, can be omitted. The third
term describes the tunneling between the wells of a double-well potential at the j-site, with
the tunneling frequency . The fourth and fth terms characterize particle interactions,
with the transverse strength B
ij
and longitudinal strength I
ij
. The values of the system
parameters depend on the properties of particle interactions and on the features of the
lattice potential. Generally, these parameters can be varied in rather wide ranges.
In order to illustrate how the tunneling frequency (15) can be varied, we may take the
double-well potential, in the vicinity of the lattice site a
j
= 0, in the form
V (r) V
0
_
r
x
r
0
_
2
_
_
r
x
r
0
_
2
2
_
+ V
H
(r
y
, r
z
) ,
where V
H
(r
y
, r
z
) is a harmonic potential in the y- and z-directions. The tunneling fre-
quency essentially depends on the parameters of the double-well potential V (r), its
depth V
0
and the interwell distance r
0
. These parameters enter the variable

1
_
2mr
2
0
V
0
,
which varies in the interval 0 < < . For intermediate values of , the tunneling fre-
quency can be calculated numerically, while for small and large , asymptotic expressions
are available [36] yielding
6V
0
exp
_

_
( 1) ,
V
0

2/3
( 1) .
These expressions show that the tunneling frequency can be regulated, by changing r
0
and V
0
, in a wide range 0 < < . This means that varying the shape of the double-
well potential, which is achievable in experiments [1621], one can regulate the tunneling
properties of the system.
III Equilibrium Phases
Thermodynamics of the system with Hamiltonian (20), for suciently long-range inter-
actions, can be accurately described in the mean-eld approximation [37] corresponding
6
to the equality
S

i
S

j
= < S

i
> S

j
+ S

i
< S

j
> < S

i
>< S

j
> (i = j) , (21)
where the angle brackets imply statistical averaging. We introduce the notation for the
average interaction strengths
A
1
N
L

i=j
A
ij
, B
1
N
L

i=j
B
ij
, I
1
N
L

i=j
I
ij
. (22)
Also, we dene the eective tunneling frequency

eff
2B < S
x
j
> (23)
and the eective mean eld
h
eff

_

2
eff
+ 4I
2
< S
z
j
>
2
. (24)
Using the standard methods of dealing with pseudospin Hamiltonians [38], we nd
< S
x
j
> =

eff
2h
eff
tanh
_
h
eff
2T
_
, < S
y
j
> = 0 ,
< S
z
j
> = < S
z
j
>
I
h
eff
tanh
_
h
eff
2T
_
, (25)
where T is temperature. The system free energy reads as
F = NE
0
+
N
2
_
A2B < S
x
j
>
2
+2I < S
z
j
>
2
_
NT ln
_
2cosh
h
eff
2T
_
, (26)
in which N = N
L
because of the unity lling factor. From the latter equation, one can
calculate all thermodynamics characteristics.
We shall be mainly interested in the properties of the average tunneling intensity
x
2
N
L

j
< S
x
j
> , (27)
average Josephson current
y
2
N
L

j
< S
y
j
> , (28)
and the average well imbalance
z
2
N
L

j
< S
z
j
> . (29)
The properties of the system depend on the parameters dened in Eq. (22), which
characterize the eective interaction strength. The parameter A enters the free energy
(26) additively, hence, it does not play here an important role. The parameter B is
7
composed of the eective interactions B
ij
given in Eq. (14). The symmetric terms V
1111
ij
,
V
2222
ij
, and V
1221
ij
are close to each other, because of which the transverse interactions B
ij
are small, so that the parameter B is substantially smaller than the parameter I. The
main role among interactions is played by the longitudinal interaction strength I. Here
we assume that I is positive, which is the often situation for the exchange interactions
[30,37,38]. Then, as follows from Hamiltonian (20), the nearest spins S
z
i
and S
z
j
tend to
align parallel in order to reduce the system energy. Therefore only the ferromagnetic-type
order is possible in the system [30,37,38].
Generally, the exchange interactions I
ij
could be negative. Then the longitudinal term
in Hamiltonian (20) would enter with the sign plus, which would imply that the nearest
spins prefer to align antiparallel to each other in order to lower the system energy. In
that case, the antiferromagnetic-type order could be the sole possibility. In the present
paper, we limit the consideration by a positive parameter I > 0, which means that only
a ferromagnetic-type order can arise.
To simplify the formulas, we use the dimensionless quantities, such as the dimensionless
tunneling frequency


I + B
(30)
and the dimensionless transverse interaction
b
B
I + B
. (31)
Also, we dene the dimensionless eld
h
h
eff
I + B
, (32)
which, with the use of the above notations, becomes
h =
_
( bx)
2
+ (1 b)
2
z
2
. (33)
From denitions (27) to (29) and Eqs. (25), we have the tunneling intensity
x =
bx
h
tanh
_
h
2T
_
, (34)
Josephson current in equilibrium
y = 0 , (35)
and the well imbalance
z = z
1 b
h
tanh
_
h
2T
_
, (36)
where temperature T is measured in units of I + B.
We may notice that Eqs. (34) to (36) are invariant under the replacement
x x , , z z . (37)
8
Therefore, without the loss of generality, we can consider only the case with x 0, 0,
and z 0.
Equation (36) shows that there can exist two thermodynamic phases, ordered and
disordered, when, respectively,
z = 0 (ordered) ,
z = 0 (disordered) .
In the ordered phase, the mean well imbalance is nonzero, and one has
x = (z = 0) , (38)
with z given by the equation
1 b
h
tanh
_
h
2T
_
= 1 . (39)
From Eqs. (33) and (38), it follows that
h = (1 b)

2
+ z
2
.
The ordered phase exists, when two conditions are valid, the tunneling frequency is not
too high,
0 < 1 , (40)
and the temperature is lower than the critical temperature
T
c
=
(1 b)
2 artanh
, (41)
so that
0 T < T
c
. (42)
When at least one of conditions (40) or (42) is not valid, that is, either 1 or T T
c
,
the ordered phase cannot exist and transfers to the disordered phase.
In the disordered phase, the mean well imbalance is zero, and one has
x = tanh
_
bx
2T
_
(z = 0) , (43)
which denes the tunneling intensity x.
At zero temperature T = 0, the ordered phase is described by the equations
x = , z =

1
2
( < 1) . (44)
While the disordered phase is characterized by the expressions
x = 1 , z = 0 ( > 1) . (45)
The quantum phase transition between the ordered and disordered phases occurs at the
tunneling frequency = 1.
9
IV Collective Excitations
Collective excitations in the system represented by the pseudospin Hamiltonian (20),
correspond to pseudospin waves. The Heisenberg equations of motion for the pseudospin
operators yield
dS
x
i
dt
= 2S
y
i

j(=i)
I
ij
S
z
j
,
dS
y
i
dt
= S
z
i
2S
x
i

j(=i)
I
ij
S
z
j
2S
z
i

j(=i)
B
ij
S
x
j
,
dS
z
i
dt
= S
y
i
+ 2S
y
i

j(=i)
B
ij
S
x
j
. (46)
For describing collective excitations, we employ the Fourier transformation for the
pseudospin operators
S

j
=
1
N
L

k
e
ika
j
,

k
=

j
S

j
e
ika
j
, (47)
where k pertains to the rst Brillouin zone. Similarly, we expand the interaction functions
I
ij
=
1
N
L

k
I
k
e
ika
ij
, B
ij
=
1
N
L

k
B
k
e
ika
ij
, (48)
in which a
ij
= a
i
a
j
. Then, Eqs. (46) acquire the form
d
x
k
dt
=
2
N
L

p
I
kp

y
p

z
kp
,
d
y
k
dt
=
z
k

2
N
L

p
I
kp

x
p

z
kp

2
N
L

p
B
kp

z
p

x
kp
,
d
z
k
dt
=
y
k
+
2
N
L

p
B
kp

y
p

x
kp
, (49)
where

k
=

k
(t).
To nd the spectrum of collective excitations, we resort to the random-phase approx-
imation. For this purpose, we look for the solutions to Eqs. (49) in the form

k
(t) = <

k
> +

k
e
it
, (50)
describing the deviations from the equilibrium values
<

k
> =
k0
N
L
< S

j
> .
Linearizing Eqs. (49) with respect to

k
yields for the spectrum of pseudospin waves

2
k
= ( Bx)[ (B + I
k
)x] + I(I + B
k
)z
2
. (51)
10
In deriving Eq. (51), we have taken into account that
lim
k0
I
k
= I , lim
k0
B
k
= B ,
with I and B dened in Eq. (22).
To consider the spectrum property in the long-wave limit, we use the asymptotic
equalities
I
k
I
1
2

i
I
ij
(k a
ij
)
2
, B
k
B
1
2

i
B
ij
(k a
ij
)
2
, (52)
where k 0. We dene the gap by the equation
_

I + B
_
2
= ( bx)( x) + (1 b)z
2
. (53)
And let us introduce a matrix G
ij
, given by the equality
G
ij
I + B
= ( bx)xI
ij
(1 b)z
2
B
ij
. (54)
Then the long-wave limit of spectrum in Eq. (51) gives

2
k

2
+
1
2

i
G
ij
(k a
ij
)
2
. (55)
The spectrum is quadratic, which is typical of spin waves.
Specifying the gap and the matrix G
ij
for the ordered and disordered phases, we
have for the ordered phase, when z = 0,

I + B
=

1 b z ,
G
ij
I + B
= (1 b)
_

2
I
ij
z
2
B
ij
_
,
and for the disordered phase, with z = 0,

I + B
=
_
( bx)( x) ,
G
ij
I + B
= ( bx)xI
ij
.
The gap disappears at the critical point of the order-disorder phase transition. Approach-
ing this point from the side of the ordered phase, one has z 0, while approaching the
point from the side of the disordered phase, one has x . In both these cases, 0
at the critical point.
V Nonequilibrium Loading
Let us suppose that atoms are loaded into a double-well lattice in an initially nonequlib-
rium state. For describing their equilibration process, one needs to study the temporal
behavior of the variables dened in Eqs. (27) to (29) and characterizing the tunneling
11
intensity x = x(t), Josephson current y = y(t), and the well imbalance z = z(t). The evo-
lution equations for these quantities can be obtained by averaging the operator equations
(46). For this purpose, one often employs the mean-eld approximation. This approx-
imation does not take into account the relaxation processes due to atomic interactions.
However, such processes could be important in order that equilibrium would be possible.
Therefore, we shall resort to a more accurate approximation, called the local-eld approx-
imation. This approximation, suggested by Wangness [39], treats atomic interactions as
occurring in a local eld formed by other particles, which results in the appearance of
relaxation characterized by the attenuation parameters,
1
and
2
whose values can be
calculated through the given atomic interactions. The local-eld approximation assumes
the existence of local equilibrium [40], so that at each moment of time the variables tend to
relax to the locally-equilibrium state. Averaging Eqs. (46) in the frame of the local-eld
approximation [39] yields the system of equations
dx
dt
= (1 b)yz
2
(x x
t
) ,
dy
dt
= ( x)z
2
(y y
t
) ,
dz
dt
= (bx )y
1
(z z
t
) , (56)
where the time variable is measured in units of 1/(I+B), while the attenuation parameters

1
and
2
, in units of I + B, and the local elds are
x
t
=
bx
h
tanh
_
h
2T
_
, y
t
= 0 ,
z
t
=
1 b
h
z tanh
_
h
2T
_
. (57)
Here T is the temperature, in units of I+B, corresponding to a would be equilibrium state,
the type of the latter being dened by the system parameters. The explicit expression
for h is given in Eq. (33). The evolution equations (56) are supplemented with the initial
conditions
x(0) = x
0
, y(0) = y
0
, z(0) = z
0
, (58)
which, generally, represent a nonequilibrium state.
Keeping in mind applications to ultracold atoms, we may set temperature to zero.
Then, Eqs. (56) possess two types of stationary solutions. One type is dened by the
xed point
x

1
= , y

1
= 0 , z

1
=

1
2
, (59)
which corresponds to the ordered phase. While another type is given by the xed point
x

2
= 1 , y

2
= 0 , z

2
= 0 , (60)
corresponding to the disordered phase. The question remains whether an atomic system
loaded into the double-well lattice, in an initially nonequilibrium state, would relax to
one of the stationary solutions.
12
VI Stability Analysis
The stability of stationary solutions can be studied by the Lyapunov stability analysis.
To this end, we calculate, in the standard way, the Jacobian matrix associated with the
dynamical system (56). Then we nd the eigenvalues of the Jacobian matrix, denoted
by , which give the characteristic exponents. The real parts Re dene the Lyapunov
exponents describing the stability properties. For simplicity, we set
1
=
2
= .
Accomplishing this procedure for the xed point (59), we nd that it is stable when
< 1. One of the characteristic exponents is exactly

1
= . (61)
The expressions for two other characteristic exponents are rather cumbersome, because
of which we write them down to the rst order in , resulting in

2,3

2 b
2
2(1 b)
i
eff
, (62)
with the eective frequency

eff
=
_
(1 b) (1
2
) ( < 1) . (63)
This shows that the xed point (59) is a stable focus, when < 1 and is small. Recall
that by its denition (31), one always has b < 1.
The xed point (60) is stable for > 1. Then, one of the characteristic exponents is
exactly the same as Eq. (61). Two other characteristic exponents, to the rst order in ,
are

2,3

2 1 b
2( b)
i

eff
, (64)
with the eective frequency

eff
=
_
( 1)( b) ( > 1) . (65)
Thus, the xed point (60) is also a stable focus for > 1 and small .
The value of the tunneling frequency = 1 is a bifurcation point related to the
dynamical phase transition. At this point, the dynamical system is neutral, with the
characteristic exponents

1
= ,
2
= 0 ,
3
= ( = 1)
for any .
An important observation, following from the above analysis, is that the dynamical
system is structurally stable only in the presence of the damping parameter > 0. Setting
the latter to zero leads to the characteristic exponents for the ordered xed point (59)

1
= 0 ,
2,3
= i
eff
( = 0, < 1)
13
and for the disordered xed point (60)

1
= 0 ,
2,3
= i

eff
( = 0, > 1)
where
eff
and

eff
are given by Eqs. (63) and (65). These characteristic exponents
demonstrate that the dynamical system is structurally unstable.
But the case of no attenuation corresponds to the usual mean-eld approximation.
Hence the above consideration shows that the mean-eld approximation results in a struc-
turally unstable dynamical system, thus, this approximation cannot correctly describe the
evolutional processes of atoms in a double-well lattice.
In order to illustrate in an explicit way the dierence between the temporal behavior of
solutions to the evolution equations (56) and the related mean-eld approximations, with
no attenuation, we solve Eqs. (56) numerically. The initial conditions (58) are chosen so
that to correspond to an initially nonequilibrium state.
In Fig. 1, the relaxation to the ordered stationary solution (59) is shown for = 1
and < 1. Taking the attenuation parameter equal to one in dimensionless units
implies that it is of order I + B in dimensional units. For these parameters, the xed
point (59) becomes a stable node. Contrary to this, the mean-eld approximation, with
= 0, exhibits permanently oscillating solutions. The curves for x(t) and z(t) always
oscillate around the given initial conditions, which is not a correct behavior. In turn, the
latter is caused by the structural instability of the dynamical system in the mean-eld
approximation. The relaxation to the disordered stationary solution (60) is demonstrated
in Fig. 2 for > 1 and = 1, when Eq. (60) corresponds to a stable node. But the mean-
eld approximation, with = 0, again exhibits an incorrect oscillatory behavior depending
on the choice of the initial conditions (58). Again, it is the structural instability that is
responsible for the incorrect dynamics in the mean-eld approximation.
Full equations (56) can also exhibit the oscillatory behavior of their solutions. This
happens when the xed points are the stable focuses, which occurs for , as follows
from Eqs. (62) to (65). In particular, if 1, than, according to Eq. (65), the oscillation
period approximately is 2/. This situation is illustrated in Fig. 3.
VII Discussion
We have considered a system of atoms in a double-well optical lattice. The case of
an insulating lattice with the unity lling factor is studied. This type of double-well
optical lattices is of special interest, being a convenient setup for realizing the ordered
and disordered states of atoms in a double-well lattice.
The principal dierence of the present paper, as compared to the earlier theoreti-
cal works on the double-well lattices [2224,41], where the superuidity-insulator phase
transition is studied for atoms with local interactions, is that we consider an insulating
double-well lattice, with atoms possessing long-range interactions. An eective Hamilto-
nian in the pseudospin representation is derived. The system exhibits two thermodynamic
phases, ordered and disordered. The phase transition between the phases can be driven
either by temperature or by varying the system parameters. For instance, one can vary
the strength of atomic interactions or the shape of the double-well potential.
14
The spectrum of collective excitations has the form typical of that for pseudospin
waves. In the long-wave limit, the spectrum is quadratic, with a gap. The latter disappears
at the phase transition point.
Nonequilibrium properties of the system are investigated. Physically, the situation
can correspond to atoms loaded into a lattice in an initially nonequilibrium state. Such a
nonequilibrium state can also be prepared by disturbing the lattice by external elds.
The relaxation of the system, from an initially nonequilibrium state to equilibrium,
cannot be described by the mean-eld approximation. This is rather clear from the
physical point of view, since the mean-eld approximation does not take into account
attenuation eects, hence, is not able to describe the equilibration process in principle.
From the mathematical point of view, as we show, the mean-eld approximation results in
a structurally unstable dynamical system. To correctly describe the process of relaxation,
the damping eects, caused by atomic interactions, must be included. This can be done,
e.g., by employing the local-eld approximation.
The long-range interactions that would be sucient for realizing the eects described
above could be of dipolar type, such that occur for polar molecules [31,42], Rydberg atoms
[32], and quantum gases with dipolar atomic interactions (see review articles [33,43]). A
good candidate would be the gas of cold atoms
52
Cr, possessing large magnetic moments
of 6
B
and, as a result, suciently strong dipolar interactions [44,45]. Because of the
existence of several species exhibiting long-range dipolar intreations, the latter can be of
dierent strength. In addition, these interactions can be eectively modulated by external
magnetic or electric elds, thus, tuning the interaction strength in a wide range [4649].
In the dynamics of atoms, relaxing from a nonequilibrium state to their equilibrium,
an important role is played by the attenuation parameter . The calculation of the latter,
for systems with dipolar interactions, is a well known procedure, described in detail in
literature on magnetic resonance [5055]. In the case of particles with magnetic moment
, interacting through magnetic dipolar forces, an accurate value of the attenuation pa-
rameter is given by the expression = z
0

2
/ h, where is the mean particle density
and z
0
is the number of nearest neighbors. This value of is of the order of the interac-
tion strength. That is why setting, in dimensionless units, to one, as we have done in
numerical calculations, is absolutely natural.
Taking, for illustration, the magnetic moment of
52
Cr, equal to = 6
B
, with
B
=
0.927 10
20
erg/G being the Bohr magneton, the typical atomic densities in a trap
(10
12
10
15
) cm
3
, and accepting for the number of nearest neighbors z
0
10,
we have (1 10
3
) s
1
. This means that the relaxation time T
rel
1/ is of order
T
rel
(10
3
1) s. The lifetime of trapped atomic systems can vary in the range of
t
exp
(1 100) s. Therefore, the relaxation time T
rel
in many cases is much smaller than
t
exp
, which implies that taking into account the attenuation is of critical importance. The
simple mean-led approximation without taking account of the relaxation eects would
be qualitatively incorrect.
In conclusion, it is worth noting that the shape of the double-well potential can be
easily regulated. Hence, it is possible to organize any required sequence of transitions
between the ordered and disordered states in the double-well lattice during the system
lifetime.
15
Acknowledgement
We are grateful for nancial support to the Russian Foundation for Basic Research
(Grant 08-02-00118).
16
Appendix A
The eigenproblem for the lattice Hamiltonian (2),
H
L
(r)
nk
(r) = E
nk

nk
(r) ,
denes the Bloch functions
nk
(r) and the Bloch spectrum E
nk
. This eigenproblem,
invoking the relation between the Bloch and Wannier functions

nk
(r) =
1

N
L

j
w
n
(r a
j
) e
ika
j
,
can be rewritten as
H
L
(r) w
n
(r a
j
) =
1
N
L

ik
E
nk
e
ika
ij
w
n
(r a
i
) ,
where a
ij
a
i
a
j
. Using this in Eq. (5) gives Eq. (6), with
E
n
ij
=
1
N
L

k
E
nk
e
ika
ij
.
Representing the latter in form (7) yields
E
n
=
1
N
L

k
E
nk
and
J
n
ij
=
1
N
L

k
E
nk
e
ika
ij
(i = j) .
These properties are used in deriving Hamiltonian (10).
17
References
[1] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Clarendon, Oxford,
2003).
[2] V.I. Yukalov, Laser Phys. Lett. 1, 435 (2004).
[3] J.O. Andersen, Rev. Mod. Phys. 76, 599 (2004).
[4] K. Bongs and K. Sengstock, Rep. Prog. Phys. 67, 907 (2004).
[5] V.I. Yukalov and M.D. Girardeau, Laser Phys. Lett. 2, 375 (2005).
[6] A. Posazhennikova, Rev. Mod. Phys. 78, 1111 (2006).
[7] V.I. Yukalov, Laser Phys. Lett. 4, 632 (2007).
[8] N.P. Proukakis and B. Jackson, J. Phys. B 41, 203002 (2008).
[9] S. Giorgini, L.P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
[10] W. Ketterle and M.W. Zwierlein, Riv. Nuovo Cimento 31, 247 (2008).
[11] D. Jaksch and P. Zoller, Ann. Phys. (N.Y.) 315, 52 (2005).
[12] O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006).
[13] C. Moseley, O. Fialko, and K. Ziegler, arXiv:0707.1979 (2007).
[14] I. Block, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2007).
[15] V.A. Yurovsky, M. Olshani, and D.S. Weiss, Adv. At. Mol. Opt. Phys. 55, 61 (2008).
[16] J. Sebby-Strabley, M. Anderlini, P.S. Jessen, and J.V. Porto, Phys. Rev. A 73, 033605
(2006).
[17] J. Sebby-Strabley, B.L. Brown, M. Anderlini, P.J. Lee, W.D. Phillips, and J.V. Porto,
Phys. Rev. Lett. 98, 200405 (2007).
[18] P.J. Lee, M. Anderlini, B.L. Brown, J. Sebby-Strabley, W.D. Phillips, and J.V. Porto,
Phys. Rev. Lett. 99, 020402 (2007).
[19] M. Anderlini, P.J. Lee, B.L. Brown, J. Sebby-Strabley, W.D. Phillips, and J.V. Porto,
Nature. 448, 452 (2007).
[20] S. Folling, S. Trotzky, P. Cheinet, M. Feld, R. Saers, A. Widera, T. M uller, and I.
Bloch, Nature 448, 1029 (2007).
[21] P. Cheinet, S. Trotzky, M. Feld, U. Schnorrberger, M. Moreno-Cardoner, S. Folling,
and I. Bloch, arXiv:0804.3372 (2008).
18
[22] I. Danshita, J.E. Williams, C.A.R. Sa de Melo, and C.W. Clark, Phys. Rev. A 76,
043606 (2007).
[23] I. Danshita, C.A.R. Sa de Melo, and C.W. Clark, Phys. Rev. A 77, 063609 (2008).
[24] V.M. Stojanovic, C. Wu, W.V. Liu, and S. Das Sarma, Phys. Rev. Lett. 101, 125301
(2008).
[25] H. Pu, W. Zhang, and P. Meystre, Phys. Rev. Lett. 87, 140405 (2001).
[26] A.B. Kuklov and B.V. Svistunov, Phys. Rev. Lett. 90, 100401 (2003).
[27] L.M. Duan, E. Demler, and M.D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
[28] E. Altman, W. Hofstetter, E. Demler, and M.D. Lukin, New J. Phys. 5, 113 (2003).
[29] N.N. Bogolubov, Lectures on Quantum Statistics (New York, Gordon and Breach,
1967), Vol. 1.
[30] N.N. Bogolubov, Lectures on Quantum Statistics (New York, Gordon and Breach,
1970), Vol. 2.
[31] J. Doyle, B. Friedrich, R.V. Krems, and F. Masnou-Seeuws, editors, Ultracold Polar
Molecules: Formation and Collisions, Eur. Phys. J. D 31, N2 (2004).
[32] T.F. Gallagher, Rydberg Atoms (Cambridge University, Cambridge, 1994).
[33] A. Griesmaier, J. Phys. B 40, 91 (2007).
[34] R.A. Duine and H.T.C. Stoof, Phys. Rep. 396 115 (2004).
[35] T. Kohler, K. Goral, and P.S. Julienne, Rev. Mod. Phys. 78, 1311 (2006).
[36] V.I. Yukalov and E.P. Yukalova, J. Phys. A 29, 6429 (1996).
[37] N.N. Bogolubov, Method for Studying Model Hamiltonians (Pergamon, Oxford,
1972).
[38] V.I. Yukalov and A.S. Shumovsky, Lectures on Phase Transitions (World Scientic,
Singapore, 1990).
[39] R.K. Wangness, Phys. Rev. 98, 927 (1955).
[40] S. De Groot and P. Mazur, Non-Equilibrium Thermodynamics (North-Holland, Am-
sterdam, 1962).
[41] I. Danshita, J.E. Williams, C.A.R. Sa de Melo, and C.W. Clark, Laser Phys. 18, 318
(2008).
[42] K.K. Ni, S. Ospelkaus, M. H. de Miranda, A Peer, B. Neyenhius, J.J. Zirbel, S.
Kotochigova, P.S. Julienne, D.S. Jin, and J. Ye, Science 322, 231 (2008).
19
[43] M.A. Baranov, Phys. Rep. 464, 71 (2008).
[44] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Phys. Rev. Lett. 94,
160404 (2005).
[45] Q. Beauls, R. Chicireanu, T. Zanon, B. Laburthe-Tolra, E. Marechal, L. Vernac,
J.C. Keller, and O. Gorceix, Phys. Rev. A 77, 061601 (2008).
[46] S. Giovanazzi, A. Gorlitz, and T. Pfau, Phys. Rev. Lett. 89, 130401 (2002).
[47] K. Glaum and A. Pelster, arXiv:cond-mat/0609374 (2007).
[48] T. Lahaye, J. Metz, B. Frohlich, T. Koch, M. Meister, A. Griesmaier, T. Pfau, H.
Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. Lett. 101, 080401 (2008).
[49] N.G. Parker, C. Ticknor, A.M. Martin, and D.H. ODell, arXiv:0810.2028 (2008).
[50] N. Bloembergen, Nuclear Magnetic Relaxation (Benjamin, New York, 1961).
[51] C.H. Pole and H.A. Farach, Relaxation in Magnetic Resonance (Academic, New York,
1971).
[52] C.P. Slitcher, Principles of Magnetic Resonance (Springer, Berlin, 1980).
[53] A. Abragam and M. Goldman, Nuclear Magnetism: Order and Disorder (Clarendon,
Oxford, 1982).
[54] V.I. Yukalov and E.P. Yukalova, Phys. Part. Nucl. 35, 348 (2004).
[55] V.I. Yukalov, Phys. Rev. B 71, 184432 (2005).
20
Figure Captions
Fig. 1. Dimensionless variables describing the tunneling intensity x(t), Josephson
current y(t), and well population imbalance z(t) as functions of dimensionless time for
= 0.1 and b = 0.5. Initial conditions are x
0
= 0.66, y
0
= 0.75, and z
0
= 0. The case of
no attenuation ( = 0) is shown by the dashed curve. The case with attenuation ( = 1)
is represented by the solid line.
Fig. 2. Dimensionless variables x(t), y(t), and z(t) as functions of dimensionless time
for = 1.5 and b = 0.5. Initial conditions are x
0
= 0.33, y
0
= 0.5, and z
0
= 0.8. The
attenuation parameters are: = 0 (dashed line) and = 1 (solid line).
Fig. 3. Population imbalance for large tunneling = 100, with = 1 and b = 0.5 as
a function of dimensionless time. Initial conditions are x
0
= 0, y
0
= 0, and z
0
= 1.
21
Figure 1: Dimensionless variables describing the tunneling intensity x(t), Josephson cur-
rent y(t), and well population imbalance z(t) as functions of dimensionless time for = 0.1
and b = 0.5. Initial conditions are x
0
= 0.66, y
0
= 0.75, and z
0
= 0. The case of no at-
tenuation ( = 0) is shown by the dashed curve. The case with attenuation ( = 1) is
represented by the solid line.
22
Figure 2: Dimensionless variables x(t), y(t), and z(t) as functions of dimensionless time
for = 1.5 and b = 0.5. Initial conditions are x
0
= 0.33, y
0
= 0.5, and z
0
= 0.8. The
attenuation parameters are: = 0 (dashed line) and = 1 (solid line).
23
Figure 3: Population imbalance for large tunneling = 100, with = 1 and b = 0.5 as a
function of dimensionless time. Initial conditions are x
0
= 0, y
0
= 0, and z
0
= 1.
24